[gmx-users] Virtual sites causing npt simulation instabilities

[Mark Abraham]

Hi,

If that fixes the issue, then please let us know - we should consider
fixing grompp to check that the inputs have appropriate consistency.

Mark

On Thu, May 31, 2018 at 1:48 PM Justin Lemkul <jalemkul at vt.edu> wrote:

>
>
> 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
>
> ==================================================
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>