[gmx-users] Virtual sites causing npt simulation instabilities

Justin Lemkul jalemkul at vt.edu
Wed Jun 6 17:28:46 CEST 2018



On 6/6/18 11:13 AM, Eisenhart, Andrew (eisenhaw) wrote:
> new link with sole mdp file
>
> https://files.fm/u/cmyg863c

This system doesn't even survive energy minimization for me, suggesting 
that the topology is not stable. The machine I'm on uses GROMACS 2016.3. 
Please upgrade to at least that version, if not 2018.1, because 5.1.2 is 
not even the terminal version in that series.

Your topology specifies a SETTLE geometry for an SPC-like water, which 
immediately leads to each of your waters (which have O-H distances of 
0.181 nm and H-H of 0.286 nm) to shrink inward, also shifting your 
virtual sites. That minimization simply crashes. I changed the SETTLE 
geometry to agree with your initial geometry, and while the atomic sites 
are OK, the qh1 and qh2 virtual sites expand outward, and though EM 
finishes, the max force is on the order of 10^8.

Check your virtual site construction, as this may be the source of the 
buildup in force. Use "nstxout = 1" in your em.mdp to see this happening 
more clearly.

Your cutoffs are bizarrely short but don't seem to be a problem. Once 
you've sorted out your topology issues, make sure your physical model is 
sane with such tiny cutoffs.

-Justin

> ________________________________
> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Justin Lemkul <jalemkul at vt.edu>
> Sent: Tuesday, June 5, 2018 9:30 PM
> To: gmx-users at gromacs.org
> Subject: Re: [gmx-users] Virtual sites causing npt simulation instabilities
>
>
>
> On 6/5/18 9:28 PM, Eisenhart, Andrew (eisenhaw) wrote:
>> Sorry about that, he’s a another link to the same files with the fixed em.mdp. I hope the 7zip format is ok.
> Please use tgz or just upload a plain-text em.mdp.
>
> -Justin
>
>> https://files.fm/u/r65e3kzq
> [https://files.fm/images/files.fm-fb.png]<https://files.fm/u/r65e3kzq>
>
> 2018-06-06_r65e3kzq<https://files.fm/u/r65e3kzq>
> files.fm
> Fast, Secure and Easy cloud file hosting, storage and safe sharing. FTP alternative. Free signup. Unlimited download traffic via torrents.
>
>
>> Andrew E. Eisenhart
>> Graduate Research Assistant | University of Cincinnati
>> 330.383.5061
>> ________________________________
>> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Justin Lemkul <jalemkul at vt.edu>
>> Sent: Tuesday, June 5, 2018 7:05:51 PM
>> To: gmx-users at gromacs.org
>> Subject: Re: [gmx-users] Virtual sites causing npt simulation instabilities
>>
>>
>>
>> On 6/5/18 7:02 PM, Eisenhart, Andrew (eisenhaw) wrote:
>>> Here's the link to the files, if there's anything else I can send let me know. I've been using 5.1.1 for the initial set up and energy minimization, and 5.1.2 for the nvt and npt simulations.
>>>
>> em.mdp is some corrupted binary file. Please upload a corrected version.
>>
>> -Justin
>>
>>> <https://files.fm/u/c9uef9qh>
>>>
>>> https://files.fm/u/c9uef9qh
>>>
>>>
>>>
>>> Thanks
>>>
>>> Andrew
>>>
>>> ________________________________
>>> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Justin Lemkul <jalemkul at vt.edu>
>>> Sent: Tuesday, June 5, 2018 12:54:50 PM
>>> To: gmx-users at gromacs.org
>>> Subject: Re: [gmx-users] Virtual sites causing npt simulation instabilities
>>>
>>>
>>>
>>> On 6/5/18 9:10 AM, Eisenhart, Andrew (eisenhaw) wrote:
>>>> Hello again,
>>>>
>>>> I have done a few things since I last updated this list.
>>>>
>>>>
>>>> I generated a new system of only 1000 water instead of the original 2000, I used a more careful approach in my energy minimization, minimizing first with all bonds and angles constrained then again without those hash constraints. This succeeded in getting the maximum force in my system to below 100kcal without lincs warnings.
>>>>
>>>>
>>>> Then I proceeded to equilibrate the temperature as before with no issues, and finally I attempted to equilibrate the box size  using the Berendsen barostat running for 2ns. The box size increases as before and still does not converge. The system stays at around 300K, but after a while (~1ns) it seems to fill the box like its a gas visually at least.
>>>>
>>>> So any ideas that I could try would be appreciated.
>>> If you can assemble a tarball that contains starting coordinates (not
>>> minimized), topology, all .mdp files, and a run script that gives the
>>> exact sequence of what you're doing, and post it somewhere to share, I
>>> will try to take a look at it and see if I can reproduce the behavior.
>>>
>>> What version of GROMACS are you using?
>>>
>>> -Justin
>>>
>>>> Thanks,
>>>>
>>>> Andrew
>>>>
>>>> ________________________________
>>>> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Eisenhart, Andrew (eisenhaw) <eisenhaw at mail.uc.edu>
>>>> Sent: Thursday, May 31, 2018 1:51:51 PM
>>>> To: gmx-users at gromacs.org
>>>> Subject: Re: [gmx-users] Virtual sites causing npt simulation instabilities
>>>>
>>>> Hey,
>>>>
>>>> Thanks for everyone's response.
>>>>
>>>>
>>>> I have tried using the berendsen barostat, with it I see similar behavior with the unit cell jumping up in size and then steadily increasing in size. I have some longer runs using it running now (2-3ns), so I'll be able to see it it eventually converges.
>>>>
>>>>
>>>> Nice catch with the particle type being incorrect, but with it corrected I'm still having the same problems (ie the oscillation of the box size) it has yet to throw lincs warnings though (even at ts = 2fs) so I'm going to simulate for a longer amount of time to see if it will converge.
>>>>
>>>>
>>>> Thanks again for everyone's reply, below is my .itp in its entirety maybe there is an error in the exceptions(which im still not 100% sure about) or something else that will stand out.
>>>>
>>>>
>>>> Andrew
>>>>
>>>>
>>>> ;
>>>>
>>>> [ atomtypes ]
>>>> ; name  bond_type    mass    charge   ptype          sigma      epsilon
>>>>       qm1   qm1         0.000   -2.516    V              0.000      0.000
>>>>       qh1   qh1         0.000    1.258    V              0.000      0.000
>>>>
>>>> [ moleculetype ]
>>>> ; molname   nrexcl
>>>> M3      3
>>>>
>>>> [ 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
>>>> #ifndef FLEXIBLE
>>>>
>>>> [ settles ]
>>>> ; OW    funct   doh        dhh
>>>> 1       1       0.09572    0.15139
>>>> #else
>>>> [ bonds ]
>>>> ; i j   funct   length  force.c.
>>>> 1   2   1   0.09572 502416.0 0.09572    502416.0
>>>> 1   3   1   0.09572 502416.0 0.09572    502416.0
>>>>
>>>> [ angles ]
>>>> ; i j   k   funct   angle   force.c.
>>>> 2   1   3   1   104.52  628.02  104.52  628.02
>>>> #endif
>>>>
>>>> [ exclusions ]
>>>> 1   2   3   4   5   6
>>>> 2   1   3   4   5   6
>>>> 3   1   2   4   5   6
>>>> 4   1   2   3   5   6
>>>> 5   1   2   3   4   6
>>>> 6   1   2   3   4   5
>>>>
>>>> [ 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
>>>>
>>>>
>>>> ________________________________
>>>> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Mark Abraham <mark.j.abraham at gmail.com>
>>>> Sent: Thursday, May 31, 2018 7:54:09 AM
>>>> To: gmx-users at gromacs.org
>>>> Subject: Re: [gmx-users] Virtual sites causing npt simulation instabilities
>>>>
>>>> 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.
>>>>>
>>>> --
>>>> 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.
>>>> --
>>>> 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.
>>> --
>>> ==================================================
>>>
>>> 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.
>> --
>> ==================================================
>>
>> 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.
> --
> ==================================================
>
> 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.

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

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

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



More information about the gromacs.org_gmx-users mailing list