[gmx-users] Re: [gmx-developers] flexible water model
Nilesh Dhumal
ndhumal at andrew.cmu.edu
Fri May 13 20:29:29 CEST 2011
Hello Justin,
I used 0.1fs timestep and I could get dielectric constant around 56 which
looks ok. To get the dielectric constant around ~75, should I reduce the
timestep around 0.05fs.
I used temperature 298k for my simulation and average tmep. I got is
297.83 K. I am doing NPT and used 1.0 bar as reference pressure in my
.mdp. I am getting average pressure around 2.58 bar.
I run the equilibration for 1ns.
Energy Average RMSD Fluct. Drift
Tot-Drift
-------------------------------------------------------------------------------
Temperature 297.83 10.9687 10.9643 -0.00165807
-1.07774
Pressure (bar) 2.57551 4246.51 4246.51 -0.0360116
-23.4075
Why there is change in pressure? Is there any wrong with simulation?
Nilesh
On Wed, May 11, 2011 4:18 pm, Justin A. Lemkul wrote:
>
>
> 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
>
>
> ========================================
> --
> 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