Project

General

Profile

Bug #222

Wrong force constant values in flexible water models

Added by Bert empty about 12 years ago. Updated about 12 years ago.

Status:
Closed
Priority:
Low
Assignee:
Erik Lindahl
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

1> From the original literature "A molecular dynamics simulation of a water model with intramolecular degrees of freedom" (http://dx.doi.org/10.1080/00268978700100141)

force constant k_OH was given as 4637 kJ/mol/angstrom^2, but spc.itp or spce.itp gives a value as 345000kJ/mol/nm^2.

-------------------------------from spc.itp-------------------------------
#ifdef FLEXIBLE
[ bonds ]
; i j funct length force.c.
1 2 1 0.1 345000 0.1 345000
1 3 1 0.1 345000 0.1 345000

[ angles ]
; i j k funct angle force.c.
2 1 3 1 109.47 383 109.47 383
#else
--------------------------------------------------------------------------

2> There was a flexible version of tip3p in tip3p.itp, and I think it was for tip3p/fs model. The original literature (THE JOURNAL OF CHEMICAL PHYSICS 124, 024503 2006) for tip3p/fs also gives different values of force constants compared with tip3p.itp.
------------------------------from tip3p.itp---------------------------------
#ifdef FLEXIBLE
[ bonds ]
; i j funct length force.c.
1 2 1 0.09572 502416.0 0.09572 502416.0
1 3 1 0.09572 502416.0 0.09572 502416.0

[ angles ]
; i j k funct angle force.c.
2 1 3 1 104.52 628.02 104.52 628.02
------------------------------------------------------------------------------

the correct one would be
#ifdef FLEXIBLE
[ bonds ]
; i j funct length force.c.
1 2 1 0.09572 443153 0.09572 443153
1 3 1 0.09572 443153 0.09572 443153

[ angles ]
; i j k funct angle force.c.
2 1 3 1 104.52 285 104.52 285


I think it needs some revisions in time.

History

#1 Updated by Bert empty about 12 years ago

1> From the original literature "A molecular dynamics simulation of a water model with intramolecular degrees of freedom" (http://dx.doi.org/10.1080/00268978700100141)

force constant k_OH was given as 4637 kJ/mol/angstrom^2, but spc.itp or spce.itp gives a value as 345000kJ/mol/nm^2.

-------------------------------from spc.itp-------------------------------
#ifdef FLEXIBLE
[ bonds ]
; i j funct length force.c.
1 2 1 0.1 345000 0.1 345000
1 3 1 0.1 345000 0.1 345000

[ angles ]
; i j k funct angle force.c.
2 1 3 1 109.47 383 109.47 383
#else
--------------------------------------------------------------------------

2> There was a flexible version of tip3p in tip3p.itp, and I think it was for tip3p/fs model. The original literature (THE JOURNAL OF CHEMICAL PHYSICS 124, 024503 2006) for tip3p/fs also gives different values of force constants compared with tip3p.itp.
------------------------------from tip3p.itp---------------------------------
#ifdef FLEXIBLE
[ bonds ]
; i j funct length force.c.
1 2 1 0.09572 502416.0 0.09572 502416.0
1 3 1 0.09572 502416.0 0.09572 502416.0

[ angles ]
; i j k funct angle force.c.
2 1 3 1 104.52 628.02 104.52 628.02
------------------------------------------------------------------------------

the correct one would be
#ifdef FLEXIBLE
[ bonds ]
; i j funct length force.c.
1 2 1 0.09572 443153 0.09572 443153
1 3 1 0.09572 443153 0.09572 443153

[ angles ]
; i j k funct angle force.c.
2 1 3 1 104.52 285 104.52 285


I think it needs some revisions in time.

#2 Updated by Erik Lindahl about 12 years ago

Hi Bert,

Interesting. However, I wouldn't classify these as bugs since the parameter files in question predate the paper (TIP3p) or came from very early versions of Gromos that predated the other paper (SPC) - the equilibrium bond and angles of TIP3p/fs (at least the Day/Voth 2002 model) are for instance not exactly the same as for vanilla TIP3p.

The "FLEXIBLE" option is meant to enable more efficient energy minimization by temporarily turning constraints off (and using the default bond/angle values from the FF), rather than being alternative water models - one shouldn't run normal simulations with it (although it's technically possible, of course).

Now, if it is only a matter of adjusting e.g. the bond/angle potential it could be neat to support both models in a single file, so we'll consider doing it for the head branch after 4.0!

Cheers,

Erik

Also available in: Atom PDF