[gmx-users] Virtual sites causing npt simulation instabilities

[Justin Lemkul]

On 5/30/18 4:40 PM, Eisenhart, Andrew (eisenhaw) wrote:
> 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

The particle type for these should be V. This could be the problem - if 
you've encoded an incorrect particle type then the update routines will 
not propagate them correctly.

-Justin

> -----------------------------------
>
> [ 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
>
>
>

-- 
==================================================

Justin A. Lemkul, Ph.D.
Assistant Professor
Virginia Tech Department of Biochemistry

303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com

==================================================