[gmx-users] Re: [gmx-developers] flexible water model
Justin A. Lemkul
jalemkul at vt.edu
Wed May 11 22:18:37 CEST 2011
Nilesh Dhumal wrote:
> I have general question. If this parameters are not giving the proper
> results then How come I get the proper results with same parameter if I
> use 1/2 value of Kbond with the same parameters.
>
Random chance, I imagine. Your .mdp settings are not correct and the cutoffs
are unreasonably short for any modern force field. If you're seeking to
reproduce others' work, then you have to use the same protocol. You're not
doing that.
> I checked the temp when I use half value of Kbond and it goes max. ~320.
>
This would also indicate instability. You managed to get a "correct" result for
one observable from an incorrect simulation wherein the thermodynamic
observables do not correspond to what you have set.
The way I found all of this out was first and foremost by watching the
trajectory. After about 30 ps, the water molecules freeze in place and vibrate
as if they are about to shear apart. At this exact time, the value of epsilon
(from epsilon.xvg, output by g_dipoles) plummets down, after having steadily
increased for the first 30 ps, as expected. The temperature spikes and the
pressure begins to oscillate by +/- 10^5, far higher than what it should for a
system of a few hundred water molecules.
Watch your trajectories, make sure you're getting the right average temperature
and pressure, and that all your other energetic terms stabilize, then analyze
the result of the simulation. Otherwise, any inferred output is meaningless.
-Justin
> Nilesh
>
> On Wed, May 11, 2011 4:03 pm, Justin A. Lemkul wrote:
>
>> 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
>>
>>
>> ========================================
>> --
>> 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