[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