[gmx-users] Re: [gmx-developers] flexible water model
Nilesh Dhumal
ndhumal at andrew.cmu.edu
Wed May 11 21:55:41 CEST 2011
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.
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
coulombtype = PME
vdwtype = cut-off
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
>
>
>
More information about the gromacs.org_gmx-users
mailing list