[gmx-users] CHARMM TIP3-Water with GMX

[Stephane Abel]

Hi all


I am doing some tests with the CHARMM port in GROMACS. Before to start 
more extensive simulations with this ff
I have performed a short run of 1000 TIP3-CHARMM water molecules in a 
cubic during at T=300 K and P=1.015 bar with v-rescale and PR as 
thermostat and barostat. I have used the VERSION 
4.0.99_development_20090927 of GMX.


############### npt.mdp #####################


title        = Bulk TIP3 equili
define        = -DCHARMM_TIP3P ; CHARMM_TIP3_WATER with LJ interactions 
in H atoms
; Run parameters
integrator    = md        ; leap-frog integrator
nsteps        = 500000    ; 2 * 500000 = 1000 ps
dt        = 0.002        ; 2 fs
nstcalcenergy   = 5
nstcomm         = 5 
; Output control
nstxout        = 1000        ; save coordinates every 2 ps
nstvout        = 1000        ; save velocities every 2 ps
nstenergy    = 1000        ; save energies every 2 ps
nstlog        = 1000        ; update log file every 2 ps
energygrps      = system

; Bond parameters
continuation    = yes        ; continuation dynamics run
constraint_algorithm = lincs    ; holonomic constraints
constraints    = all-bonds    ; all bonds (even heavy atom-H bonds) 
constrained
lincs_iter    = 1        ; accuracy of LINCS
lincs_order    = 4        ; also related to accuracy
; Neighborsearching
ns_type        = grid        ; search neighboring grid cels
nstlist        = 5        ; 5*2fs 10 fs
rlist        = 1.0        ; short-range neighborlist cutoff (in nm)
rcoulomb    = 1.0        ; short-range electrostatic cutoff (in nm)
;vdW
rvdw        = 1.2        ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype    = PME        ; Particle Mesh Ewald for long-range 
electrostatics
pme_order    = 4        ; cubic interpolation
fourierspacing    = 0.12        ; grid spacing for FFT
; Temperature coupling is on
tcoupl        = V-rescale    ; modified Berendsen thermostat
tc-grps        = system        ; one coupling groups - more accurate
tau_t        = 0.1        ; time constant, in ps
ref_t        = 300             ; reference temperature, one for each 
group, in K
; Pressure coupling is off
pcoupl        = Parrinello-Rahman     ; pressure coupling in NPT
pcoupltype      = isotropic     ; uniform scaling of box vectors
tau_p           = 1.0           ; time constant, in ps
ref_p           = 1.0135        ; reference pressure, in bar
compressibility = 4.5e-5        ; isothermal compressibility of water, 
bar^-1
; Periodic boundary conditions
pbc        = xyz        ; 3-D PBC
; Dispersion correction
DispCorr    = EnerPres    ; account for cut-off vdW scheme
; Velocity generation
gen_vel        = no        ; assign velocities from Maxwell distribution


###############


After 1ns of MD, I can reproduce the density (1.014 g/cm3) and RDF of 
oxygen-oxygen water given in JCTC paper of Bjelkmar.


and with g_energy_mpi -f -nmol 1000 npt.edr command  I obtain


Statistics over 500001 steps [ 0.0000 thru 1000.0000 ps ], 12 data sets
All averages are over 100001 frames

Energy                      Average       RMSD     Fluct.  Tot-Drift
-------------------------------------------------------------------------------
LJ (SR)                     4.46356   0.197047   0.196872 -0.0287239  
(kJ/mol)
LJ (LR)                  -0.0808426 0.00072743 0.000725833 0.000166929  
(kJ/mol)
Disper. corr.             -0.109557 0.000959736 0.000957692 0.000216851  
(kJ/mol)
Coulomb (SR)               -41.7696   0.309073   0.308661  0.0552643  
(kJ/mol)
Coul. recip.               -3.69787  0.0144802  0.0144802 3.31364e-05  
(kJ/mol)
Potential                  -41.1943    0.19617   0.196016  0.0269571  
(kJ/mol)
Kinetic En.                 7.47811   0.135378   0.135378 0.000498793  
(kJ/mol)
Total Energy               -33.7162   0.237225   0.237093   0.027456  
(kJ/mol)
Temperature                 299.952     5.4301    5.43009  0.0200087  (K)
Pressure                    1.45811    429.519    429.013   -72.2231  (bar)
Volume                      29.4835   0.258304   0.257752  0.0584779  (nm^3)
Density                     1014.72    8.88909    8.87016   -2.00846  
(kg/m^3)


The computed RMSD of the water oxygen atoms with the command g_msd_mpi 
-s npt.tpr -f npt.trr -n index.ndx -b 300 -e 1000



Diffusion : D[        OW] 5.3529 (+/- 0.2862) 1e-5 cm^2/s


The translational diffusion are nearly corrects compared with the 
literature  (5.86 1e-5 cm^2/s at 298 K, by Mark et Nilsson , J. Phys. 
Chem. A, 2001, 105 (43), pp 9954)



But i noticed that my average pressure seem to high (1.5 bar) for 1ns of 
MD and by consequence I am not very confident with the parameters I used 
in my mdp. Probably this run was also to short.

Moreover, since I am not a power user of  GROMACS, I don't know how to 
translate in GROMACS directives the sentence given in JCTC paper " using 
a short-range cutoff of 1.2 nm, and van der Waals interactions were 
switched off between 1.0 to 1.2 nm". And by consequence what are the 
corrects coulomb and vdw parameters for simulations of CHARMM ff ?


Thank you in advance for your advices and comments.


Stephane