[gmx-users] TIP4P water model

[abdus sabuj]

Hi gmx-users,

I am trying to equilibrate my system with water , but my simulation is
breaking down, showing LINCS warnings. I have used the following .top file
for TIP4P.

[ defaults ]
; nbfunc    comb-rule   gen-pairs   fudgeLJ fudgeQQ
1           2        no      0.0 0.0


[ atomtypes ]
;TIP4P
;name bond_type mass    charge   ptype sigma         epsilon
TIPW     8  15.999400     0.000     A    0.315365   0.648520
TIPH1    1  1.0080000     0.520     A    0.00000    0.000000
TIPH2    1  1.0080000     0.520     A    0.000000   0.000000
TIPM     0  0.0000000    -1.040     D    0.000000   0.000000


[ moleculetype ]
; molname   nrexcl
SOL     2

[ atoms ]
; id    at type res nr  residu name at name cg nr   charge
1       TIPW        1       SOL      OW     1       0.0
2       TIPH1       1       SOL     HW1     1       0.52
3       TIPH2       1       SOL     HW2     1       0.52
4       TIPM        1       SOL      MW     1      -1.04

[ bonds ]
; i j   funct   length  force.c.
1   2   1   0.09572      502416.0
1   3   1   0.09572      502416.0

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

[ exclusions ]
1   2   3   4
2   1   3   4
3   1   2   4
4   1   2   3


[ virtual_sites3 ]
; Vsite from    funct   a             b
4   1   2   3   1       0.128012065  0.128012065


[ system ]
; Name
alkane1/water

[ molecules ]
; Compound        #mols
SOL               800

And the options I used  for NPT equlibration in my .mdp file ....


;
integrator               = md
nsteps                   = 2500000
dt                                   = 0.002
nstenergy                = 10000
nstxout                  = 10000
nstvout                  = 10000

vdw-type                 = cut-off


continuation    = yes                ; first 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
lincs-warnangle = 90                    ; warning after rotating the bond


ns_type         = grid          ; search neighboring grid cels
nstlist         = 10                 ; 10 fs
rlist           = 1.4           ; short-range neighborlist cutoff (in nm)
rcoulomb        = 1.4           ; short-range electrostatic cutoff (in nm)
rvdw            = 1.4           ; short-range van der Waals cutoff (in nm)
coulombtype     = PME           ; Particle Mesh Ewald for long-range
electrostatics
pme_order       = 4                 ; cubic interpolation
fourierspacing  = 0.16          ; grid spacing for FFT

; set temperature to 300K
tcoupl                   = v-rescale
tc-grps                  = system
tau-t                    = 0.4
ref-t                    = 303
; and pressure to 1 bar
pcoupl                   = berendsen
ref-p                    = 1
compressibility  = 4.5e-5
tau-p                    = 5
; generate initial velocities
gen-vel                  = no
;gen-temp                 = 300
;gen_seed                = -1
; Periodic boundary conditions
pbc                 = xyz               ; 3-D PBC
; Dispersion correction
DispCorr        = EnerPres      ; account for cut-off vdW scheme

; COM motion removal
; These options remove motion of the protein/bilayer relative to the
solvent/ions
nstcomm         = 1
comm-mode       = Linear
comm-grps       = system


I know, I am doing something wrong, but could not figure it out.

Thanks,
Sabuj