[gmx-users] SWM4-NDP water model in GROMACS

[Mikhail Stukan]

Dear gromacs users,

Does anyone have experience in implementation of polarizable SWM4-NDP water model by Lamoureux et al. (Chem.Phys.Lett. 418,245 (2006)) in gromacs?
I have performed several simulations with systems of different size (from 1024 up to 8192 water molecules)  with dispersion correction term and the cut-off 15A for both Ewald and LJ  (according to the original paper), but on the regular basis I have higher density than the one reported in the original paper. Namely, average molecular volume (calculated by g_energy) is 29.67 (density 1008.25) instead of 29.91 (density 1000.18).  Results are  same for 4.5.1 and 4.0.7 version of gromacs. Other parameters such as dielectric constant and surface tension looks fine, though they have much higher error bars.

Below are  the .itp and .mdp files  I am using.

Thank you very much in advance.

Mikhail

.itp
----------------------------------------------------------------------------------------------------------
[ defaults ]
; LJ non-bonded  Lorenz-Berthelot  Generate 1-4 pairs   FudgeLJ    FudgeQQ
       1             2                    yes           1.0        1.0

[ atomtypes ]
;name        mass      charge     ptype   c6     c12
   WO    15.99940    1.71636e+00    A     0.0    0.0
   WH     1.00800    0.55733e+00    A     0.0    0.0
   WV    0.0       -1.11466e+00    V     0.0    0.0
   WS     0.0       -1.71636e+00    S     0.0    0.0

[ nonbond_params ]
;                             sigma             epsilon
WO      WO      1          0.318395e+00      882.57296e-03

 [ moleculetype ]
; molname     nrexcl
SW4           2

[ atoms ]
; id   at type       res nr        residu name   at name              cg nr  charge
1      WO     1      SW4           OW1           1       1.71636e+00
2      WH     1      SW4           HW2           1      0.55733e+00
3      WH     1      SW4           HW3           1      0.55733e+00
4      WV     1      SW4           VW            1      -1.11466e+00
5      WS     1      SW4           SW            1      -1.71636e+00

[ polarization ]
; See notes above.   alpha (nm^3)
1      5      1      0.97825e-03

[ settles ]
; dHH = 0.15139 gives HOH agle equal to 104.52 degree
; i    funct  dOH    dHH
 1     1      0.09572       0.15139


;[ constraints ]
;1   2    1   0.09572
;1   3    1   0.09572
;3   2    1   0.15139


[ virtual_sites3 ]
; The position of the virtual site is computed as follows:
;
;             O
;
;             V
;
;      H             H
;
;
; Virtual site from         funct  a             b
4      1      2      3      2      0.5    0.024034

[ exclusions ]
; iatom excluded from interaction with i
1      2      3      4      5
2      1      3      4      5
3      1      2      4      5
4      1      2      3      5
5      1      2      3      4

.mdp
===================================================================
integrator          =  md
dt                  =  0.001              ; time step
nsteps              =  1000000                   ; number of steps
comm_mode           =  Linear;Angular            ; Remove center of mass translation and rotation
nstcomm             =  1003               ; reset c.o.m. motion
nstxout             =  10000              ; write coords
nstvout             =  10000                  ; write velocities
nstxtcout           =  2000                  ; write coords to xtc-trajectory file
nstlog              =  2000               ; print to logfile
nstlist             =  20                 ; update pairlist
ns_type             =  grid ;simple              ; pairlist method
;================== Polarizable model parameters =======================
emtol             =  0.1              ;the convergency criterion for maximum force
niter               =  20                ;maximum number of iterations for the shell particle optimization
;=================  Cutt off specification =============================
pbc               = xyz            ; periodic boundary conditions
optimize_fft        =  yes                ; perform FFT optimization at start
coulombtype         =  PME
rcoulomb            =  1.5                ; cut-off for coulomb
rlist               =  1.5                ; cut-off for ns
vdw-type            =  cut-off
rvdw                =  1.5                ; cut-off for vdw
dispcorr            =  enerpres
;=======================================================================
Tcoupl              =  v-rescale ; berendsen;nose-hoover;      ;        ; temperature coupling
tc-grps             =  SW4
ref_t             =  298.15
tau_t             =  0.1
Pcoupl              =  berendsen          ; pressure bath
Pcoupltype          =  isotropic          ; pressure geometry
tau_p               =  0.5                ; p-coupoling time
compressibility     =  4.5e-5                 ;
ref_p               =  1.01325
;=======================================================================
gen_vel             =  no                 ; generate initial velocities
gen_temp            =  300                ; initial temperature
gen_seed            =  1903               ; random seed





=============================================
Dr Mikhail Stukan
Schlumberger Dhahran Carbonate Research Center,
Dhahran Techno Valley  - KFUPM,
P.O. Box 39011, Dammam / Doha Camp  31942,
Kingdom of Saudi Arabia
Tel: +966 3 331 0300 ext 6182
Fax:+966 3 330 0845
mstukan at slb.com<mailto:mstukan at slb.com>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100914/434e81a2/attachment.html>