[gmx-users] Re: [gmx-developers] flexible water model
Nilesh Dhumal
ndhumal at andrew.cmu.edu
Wed May 11 22:13:13 CEST 2011
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.
I checked the temp when I use half value of Kbond and it goes max. ~320.
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
>
>
>
More information about the gromacs.org_gmx-users
mailing list