[gmx-users] Virtual sites causing npt simulation instabilities

[Dallas Warren]

Always good practice to use Berendsen when turning on constant
pressure first, it is much more robust than P-R.  Let that run for
awhile, let it equilibrate, then turn on P-R.
Catch ya,

Dr. Dallas Warren
Drug Delivery, Disposition and Dynamics
Monash Institute of Pharmaceutical Sciences, Monash University
381 Royal Parade, Parkville VIC 3052
dallas.warren at monash.edu
---------------------------------
When the only tool you own is a hammer, every problem begins to resemble a nail.


On 31 May 2018 at 06:40, Eisenhart, Andrew (eisenhaw)
<eisenhaw at mail.uc.edu> 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
>
> -----------------------------------
>
> [ 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
>
>
>
> --
> 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.