[gmx-users] Re: [gmx-developers] flexible water model
Justin A. Lemkul
jalemkul at vt.edu
Wed May 11 22:03:13 CEST 2011
Nilesh Dhumal wrote:
> Hello Justin,
>
> I used 0.5 fs time step and still I got dielectric constant ~2.
>
> This is the md.mdp file I used. I checked the temperature it didn't go
> more than ~ 320 K.
>
Still, that's unacceptable if your ref_t is 295. The fact that the temperature
is still rocketing up suggests the same instability I found.
> title = cpeptide MD
> cpp = /usr/bin/cpp
> integrator = md
> dt = 0.001 ; ps !
> nsteps = 6500000 ; total 5 ps.
> nstcomm = 1000
> nstxout = 1000
> nstvout = 1000
> nstfout = 1
> nstlist = 1
> ns_type = grid
> rlist = 0.6
> rcoulomb = 0.6
> rvdw = 0.7
These cutoffs make no sense. The paper you've cited used 0.9 nm. If you're
trying to reproduce previous work, use the same settings.
> coulombtype = PME
> vdwtype = cut-off
Similarly here, you haven't applied dispersion correction, which you should, per
the methods used in the JCP article.
-Justin
> pbc = xyz
> fourierspacing = 0.12
> fourier_nx = 0
> fourier_ny = 0
> fourier_nz = 0
> pme_order = 4
> ewald_rtol = 1e-5
> optimize_fft = yes
> ; Berendsen temperature coupling is on
> Tcoupl = v-rescale
> tau_t = 0.1
> tc-grps =system
> ref_t = 295
> ; Pressure coupling is on
> Pcoupl = parrinello-rahman
> pcoupltype = isotropic
> tau_p = 0.5
> compressibility = 4.5e-5
> ref_p = 1.0
> ; Generate velocites is on at 300 K.
> gen_vel = yes
> gen_temp = 295.0
> gen_seed = 173529
>
>
> Nilesh
>
>
> On Tue, May 10, 2011 6:14 pm, Justin A. Lemkul wrote:
>
>> The gmx-developers list is not the forum for these types of questions. I
>> am replying via gmx-users, which is where the discussion should stay.
>>
>> I took a few minutes to dig into this. My conclusion is that your system
>> is not stable. I would encourage you to analyze the temperature and
>> pressure of your systems that are giving odd results. When I used the
>> topology you provided, the temperature spiked to over 600 K and the box
>> began to oscillate extensively as if the system were about to explode.
>> The resulting epsilon value was about 1.6.
>>
>>
>> If I reduce the timestep to 0.5 fs, the results are much more reasonable,
>> and using the Ka and Kb values in the paper (not halved!) I get a much
>> more sensible result, even after just 100 ps. It looks to me as if the
>> combination of these parameters and a 1-fs timestep is not entirely
>> stable. I know the original authors used dt = 1 fs, but the point
>> remains.
>>
>> -Justin
>>
>>
>> Nilesh Dhumal wrote:
>>
>>> Hello,
>>>
>>>
>>> I am using flexible water model for my system. I am referring a paper
>>> J.
>>> Chem. Phys. 124, 024503 (2006). Author have used Amber type force field.
>>>
>>>
>>> i.e. 1/2 Kbond (r-req)**2 + 1/2 Kangle(theta-thetaeq)**2.
>>>
>>> Kbond = 443153.3808 kJ/mol nm**2
>>>
>>>
>>> Kangle= 317.5656 kJ/mol rad**2.
>>>
>>>
>>> I am using olss-aa force field parameters in Gromacs VERSION 4.0.7.
>>>
>>>
>>> I checked some papers in which author have used oplsaa force field in
>>> Gromacs. 1/2 factor is not in opls force field if I compare opls and
>>> amber.
>>>
>>>
>>> I didn’t get the proper dielectric constant for water when I used the
>>> parameters reported in paper
>>>
>>> (Kbond = 443153.3808 kJ/mol nm**2,Kangle= 317.5656 kJ/mol rad**2)
>>>
>>>
>>> I half the value of Kb and get the proper dielectric constant (~80) for
>>> water reported in paper. If I half Kangle then I don’t get proper
>>> value. Below are the results for the dielectric constant of water.
>>> Bond length is nm.
>>>
>>>
>>> Here I have done some analysis. The original value reported in J.
>>> Chem.
>>> Phys. 124, 024503 2006, paper are
>>>
>>>
>>> Kbond = 443153.3808 kJ/mol nm**2
>>> Kangle = 317.5656 kJ/mol rad**2.
>>>
>>>
>>>
>>> bond length Kbond angle Kangle dielectric constant
>>> 0.1012
>>> 443153.3808 113.24 317.5656 ~1.9 : orginal value
>>>
>>>
>>>
>>> 0.1012 221576.6904 113.24 317.5656 ~80 : 1/2 (Kbond)
>>>
>>>
>>>
>>> 0.1012 443153.3808 113.24 158.7828 ~1.58 : 1/2 (kangle)
>>>
>>>
>>>
>>> 0.1012 221576.6904 113.24 317.5656 ~1.9 : 1/2
>>> (Kbond)&(Kangle)
>>>
>>>
>>>
>>> Here I pasted spc_fw.itp file used for the simulation.
>>>
>>>
>>> ; Derived from parsing of runfiles/alat.top.orig
>>> [ defaults ]
>>> ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
>>> ;1 3 yes 0.5 0.5
>>> ; comb-rule 3 is square-root sigma, the OPLSAA version
>>>
>>>
>>> [ atomtypes ]
>>> ; full atom descriptions are available in ffoplsaa.atp
>>> ; name bond_type mass charge ptype sigma epsilon
>>> opls_1001 OW 8 15.99940 -0.820 A 3.1655e-01
>>> 6.503e-01
>>> opls_1002 HW 1 1.00800 0.410 A 0.00e+00 0.00e+00
>>>
>>>
>>> [ bondtypes ]
>>> ; i j func b0 kb
>>> OW HW 1 0.1012 443153.3808 ; J. Chem. Phys.
>>> (2006),124,024503
>>> [ angletypes ]
>>> ; i j k func th0 cth
>>> HW OW HW 1 113.24 317.5656 ; J. Chem. Phys.
>>> (2006),124,024503
>>>
>>>
>>> [ moleculetype ]
>>> ; Name nrexcl
>>> WAT 3
>>>
>>>
>>> [ atoms ]
>>> ; nr type resnr residue atom cgnr charge mass typeB
>>> chargeB massB 1 opls_1001 1 WAT OW 1 -0.82
>>> 15.99940 ;
>>> 2 opls_1002 1 WAT HW1 1 0.41 1.008 ;
>>> 3 opls_1002 1 WAT HW2 1 0.41 1.008 ;
>>>
>>>
>>> [ bonds ]
>>> ; i j funct
>>> 1 2 1
>>> 1 3 1
>>>
>>>
>>> [ angles ]
>>> ; i j k funct
>>> 2 1 3 1
>>>
>>>
>>> In water.top file, I included spc_fw.itp file.
>>>
>>>
>>> ; Include water topology
>>> #include "spc_fw.itp"
>>>
>>>
>>>
>>> I run the simulation 6.5 ns for collecting data and I have total 256
>>> water molecules.
>>>
>>>
>>> Is there anything wrong with my .itp file?
>>>
>>>
>>> I am using Gromacs VERSION 4.0.7.
>>>
>>>
>>>
>>>
>>> Thanks
>>>
>>>
>>> Nilesh
>>>
>>>
>>>
>>>
>> --
>> ========================================
>>
>>
>> Justin A. Lemkul
>> Ph.D. Candidate
>> ICTAS Doctoral Scholar
>> MILES-IGERT Trainee
>> Department of Biochemistry
>> Virginia Tech
>> Blacksburg, VA
>> jalemkul[at]vt.edu | (540) 231-9080
>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>>
>>
>> ========================================
>> --
>> gmx-users mailing list gmx-users at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-request at gromacs.org. Can't post? Read
>> http://www.gromacs.org/Support/Mailing_Lists
>>
>>
>>
>
>
>
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
More information about the gromacs.org_gmx-users
mailing list