[gmx-users] Re: [gmx-developers] flexible water model

[Justin A. Lemkul]

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

========================================