## Bug #222

### Wrong force constant values in flexible water models

**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