[gmx-users] Virtual sites causing npt simulation instabilities

[Eisenhart, Andrew (eisenhaw)]

Hello all,

I am having an issue with the current systems I am working on. The system consists of 2000 6point water molecules (3 chargeless masses for Ow, hw1, and hw2; and 3 massless point charges for the electrostatics). The three point charges I have implemented as virtual interaction sites using the section below in my topology file.


[ virtual_sites3 ]
; Vsite from            funct   a       b
4   1   2   3   1   0.29869481802   0.29869481802
; Vsite (3fad) funct        theta       d
5   3   2   1   3     18.06    0.025
6   2   3   1   3     18.06    0.025


this is loosely based on the tip4p model. Since the two partial charges off of the hydrogens differ only in location I have them defined as a single type with sigma and epsilon defined as 0. I have also done this to the partial charge off of the oxygen. This is all seen below in my atomtypes and atoms sections.


[ atomtypes ]
; name  bond_type    mass    charge   ptype          sigma      epsilon
  qm1   qm1         0.000   -2.516    A              0.000      0.000
  qh1   qh1         0.000    1.258    A              0.000      0.000

-----------------------------------

[ atoms ]
; id    at type res nr  residu name at name cg nr   charge
1       opls_113        1       M3      OW     1       0.0
2       opls_114        1       M3     HW1     1       0.0
3       opls_114        1       M3     HW2     1       0.0
4       qm1             1       M3      qm     1      -2.516
5       qh1             1       M3     qh1     1       1.258
6       qh1             1       M3     qh2     1       1.258

The problem I am seeing right now depends on the parameters I use in the mdp files. During thermostat equlibration nothing is amiss, but once the volume of the box is allowed to relax the problems show. Using the parrinello-rahman barostat and 0.002fs timestep I crash with lincs warnings. Reducing the timestep to 0.001fs stops the lincs warnings, but the box size never converges It  instead expands to around 10x its original size. So I think that I must have a problem with the implementation of my virtual sites. Any guidance would be greatly appreciated. Thanks below are some other parameters I have set in my input files.

Andrew


integrator          =  md
nsteps              =  500000
nstcomm             =  100
nstxout             =  0
nstvout             =  0
nstfout             =  0
nstlog              =  1000
nstenergy           =  1000
nstxtcout           =  1000
nstlist             =  5
ns_type             = grid
pbc                 =  xyz
coulombtype     = pme
rcoulomb        = 0.5
rlist           = 0.5
vdw-type        = cut-off
rvdw            =0.5
constraint_algorithm = lincs
constraints     = all-bonds
lincs_iter      = 1
lincs_order     = 8
fourierspacing       = 0.10
pme_order        = 6
ewald_rtol       = 1e-06
ewald_geometry       = 3d
optimize_fft         = yes
tcoupl          = Nose-Hoover
tc-grps         = System
tau_t           = 0.4
ref_t           = 300
; Pressure coupling is on
pcoupl          = parrinello-rahman
pcoupltype      = isotropic
tau_p           = 2.0
ref_p           = 1.0
compressibility = 4.5e-5
refcoord-scaling = com
DispCorr        = EnerPres
comm-mode       = linear