[gmx-users] Virtual sites causing npt simulation instabilities
Justin Lemkul
jalemkul at vt.edu
Wed Jun 6 01:06:05 CEST 2018
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
==================================================
More information about the gromacs.org_gmx-users
mailing list