[gmx-users] Regular vs. CHARMM TIP3P water model

[Nicolas SAPAY]

Hello everybody,

I have many difficulties to reproduce TIP3P simulation results with CHARMM
TIP3P. Regular TIP3P gives systematically a lower density than its CHARMM
counterpart, independantly from the cutoff for non-bonded interactions,
the version of GROMACS (4.5.1 or 4.0.7) or the double precision.
Regular 961,067 +/-0,756 g/L
CHARMM  980,860 +/-0,492 g/L

The Enthalpy of vaporization follows a similar scheme:
regular -39,992 +/-0,021 kJ/mol
CHARMM  -40,665	+/-0,009 kJ/mol

In fact, CHARMM gives results closer to what I should obtain at 300K and 1
bar for TIP3P. I suspect an issue with the bond constraints, but I can't
locate precisely where it is. Settle parameters are exactly the same (as
well as constraints for flexible water). Did someone ever face a similar
problem?

Rgards,
Nicolas

P.S. I My system is just a box of 1728 water molecules pre-equilibrated
for 500 ps at 300 K and 1 bar.

P.S. II I'm using the following simulation parameters:
constraints              = hbonds
constraint_algorithm     = LINCS
continuation             = no
lincs-order              = 4
lincs-iter               = 1
lincs-warnangle          = 30
morse                    = no
(LINCS or SHAKE should not make any difference since TIP3P is normally
treated by SETTLE).

Other parameters are:
; COUPLING
; Temperature coupling
tcoupl                   = Berendsen
nsttcouple               = -1
nh-chain-length          = 10
tc_grps                  = System
tau_t                    = 0.1
ref_t                    = 300
; Pressure coupling
pcoupl                   = Berendsen
pcoupltype               = isotropic
nstpcouple               = -1
tau_p                    = 1
compressibility          = 4.5e-5
ref_p                    = 1

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist                  = 10
ns_type                  = grid
pbc                      = xyz
rlist                    = 1.15

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype              = PME
rcoulomb-switch          = 0
rcoulomb                 = 1.15
vdwtype                  = Cutoff
rvdw                     = 1.0
fourierspacing           = 0.10
pme_order                = 6
ewald_rtol               = 1.0e-6
ewald_geometry           = 3d
epsilon_surface          = 0
optimize_fft             = yes

; RUN CONTROL PARAMETERS
integrator               = md
dt                       = 0.002
nsteps                   = 200000
comm_mode                = Linear
nstcomm                  = 1
comm_grps                = System


-- 
[ Nicolas Sapay - Post-Doctoral Fellow ]
CERMAV - www.cermav.cnrs.fr
BP53, 38041 Grenoble cedex 9, France
Phone: +33 (0)4 76 03 76 44/53