[gmx-users] Lincs all-bonds causing instability in otherwise stable system
Mark Abraham
mark.j.abraham at gmail.com
Fri May 22 14:00:14 CEST 2015
Hi,
After further thought, the base problem is that there must be at least one
rigid triangle, and LINCS can't do this because it couples constraints too
much. The water models use [settles] to solve the particular case of a
single triangle, and TIP4P and TIP5P water models are implemented with
virtual sites on such a triangle. So with a bit of pencil and paper, I'm
sure you could e.g. have a Me-Me-O rigid triangle and use a 3out virtual
site (see manual) for the sulfur. But you'd only do this if the need to use
SHAKE was hurting your simulation performance.
Mark
On Thu, May 21, 2015 at 10:16 AM Cara Kreck <
cara.kreck at student.curtin.edu.au> wrote:
> Thanks Mark,
>
> A grompp warning about problematic constraints sounds like a good idea.
>
> DMSO is actually tetrahedral, not planar. Is that possible (and hopefully
> straightforward) to do with virtual sites? If it is, is there an example I
> could look at?
>
> Thanks,
>
> Cara
>
>
> > From: mark.j.abraham at gmail.com
> > Date: Wed, 20 May 2015 11:47:43 +0000
> > To: gmx-users at gromacs.org; gromacs.org_gmx-users at maillist.sys.kth.se
> > Subject: Re: [gmx-users] Lincs all-bonds causing instability in
> otherwise stable system
> >
> > Hi,
> >
> > Sorry that was so painful for you. The LINCS documentation in the manual
> > suggests that that combination of bonds becoming constraints won't work
> > well. Perhaps we should add a 3-edge cycle detector to grompp?
> >
> > If you actually want a rigid planar DMSO, then three real atoms and a
> > virtual site would be the best implementation.
> >
> > Mark
> >
> > On Wed, May 20, 2015 at 6:41 AM Cara Kreck <
> cara.kreck at student.curtin.edu.au>
> > wrote:
> >
> > > Hi everyone,
> > >
> > > I never got a reply to my message but I figured out the problem by
> myself.
> > > Just in case anyone else runs into a similar problem I thought I should
> > > explain what was wrong and share the solution.
> > >
> > > I was using a DMSO topology from the ATB that uses extra bonds to fix
> the
> > > geometry instead of angle or dihedral terms:
> > >
> > > [ moleculetype ]
> > > ; Name nrexcl
> > > DMSO 3
> > > [ atoms ]
> > > ; nr type resnr resid atom cgnr charge mass total_charge
> > > 1 SDmso 1 DMSO S 1 0.128 32.0600
> > > 2 ODmso 1 DMSO O 1 -0.448 15.9994
> > > 3 CDmso 1 DMSO C 1 0.160 15.0350
> > > 4 CDmso 1 DMSO C 1 0.160 15.0350 ; 0.000
> > > ; total charge of the molecule: 0.000
> > > [ bonds ]
> > > ; ai aj funct c0 c1
> > > 1 2 2 0.1530 8.0400e+06
> > > 1 3 2 0.1938 4.9500e+06
> > > 1 4 2 0.1938 4.9500e+06
> > > 2 3 2 0.2794 2.3900e+06
> > > 2 4 2 0.2794 2.3900e+06
> > > 3 4 2 0.2912 2.1900e+06
> > >
> > > This topology is fine to use with SHAKE but LINCS doesn't seem to be
> able
> > > to handle it. When I removed the extra bonds my simulations were able
> to
> > > run with all bonds constrained by LINCS. Then I found appropriate
> angle and
> > > dihedral parameters in the GROMOS ffbonded.itp file to control the
> geometry
> > > again. My topology file now looks like this:
> > >
> > > [ moleculetype ]
> > > ; Name nrexcl
> > > DMSO 3
> > > [ atoms ]
> > > ; nr type resnr resid atom cgnr charge mass total_charge
> > > 1 SDmso 1 DMSO S 1 0.128 32.0600
> > > 2 ODmso 1 DMSO O 1 -0.448 15.9994
> > > 3 CDmso 1 DMSO C 1 0.160 15.0350
> > > 4 CDmso 1 DMSO C 1 0.160 15.0350 ; 0.000
> > > ; total charge of the molecule: 0.000
> > > [ bonds ]
> > > ; ai aj funct c0 c1
> > > 1 2 2 0.1530 8.0400e+06
> > > 1 3 2 0.1938 4.9500e+06
> > > 1 4 2 0.1938 4.9500e+06
> > > [ angles ]
> > > ; ai aj ak funct angle fc
> > > 3 1 4 2 97.40 469.00
> > > 3 1 2 2 106.75 503.00
> > > 4 1 2 2 106.75 503.00
> > > [ dihedrals ]
> > > ; GROMOS improper dihedrals
> > > ; ai aj ak al funct angle fc
> > > 1 3 4 2 2 35.26439 334.84617 ; tetrahedral centre
> > >
> > > I ran an energy minimisation of a single molecule with the new topology
> > > and its geometry overlapped the old one perfectly. So the problem was
> > > difficult to diagnose but easy to fix. Especially since I was able to
> > > energy minimise the individual molecules with all bonds constrained
> but the
> > > constraints went haywire when everything was combined in the full
> system.
> > > Hopefully I can at least save someone else from wasting 3 weeks trying
> to
> > > get a similar topology to work with LINCS.
> > >
> > > Cara
> > >
> > >
> > > > From: cara.kreck at student.curtin.edu.au
> > > > To: gromacs.org_gmx-users at maillist.sys.kth.se
> > > > Date: Fri, 15 May 2015 14:31:16 +0800
> > > > Subject: [gmx-users] Lincs all-bonds causing instability in otherwise
> > > stable system
> > > >
> > > > I forgot to include an example of my mdp files. I've tried varying
> the
> > > timestep, running with and without pressure coupling, and obviously
> > > changing the constraints as outlined in the previous message:
> > > >
> > > > integrator = md
> > > > dt = 0.001 ; 1fs
> > > > nsteps = 100000 ; 100ps
> > > > comm_grps = DOPC !DOPC
> > > > nstxout = 1000
> > > > nstvout = 0
> > > > nstlog = 1000
> > > > nstenergy = 1000
> > > > energygrps = DOPC !DOPC
> > > > nstcalcenergy = 5
> > > > nstlist = 5
> > > > ns_type = grid
> > > > pbc = xyz
> > > > rlist = 0.8
> > > > coulombtype = Reaction-Field
> > > > rcoulomb = 1.4
> > > > epsilon_rf = 62
> > > > vdwtype = Cut-off
> > > > rvdw = 1.4
> > > > tcoupl = berendsen
> > > > tc-grps = DOPC !DOPC
> > > > tau_t = 0.1 0.1
> > > > ref_t = 303 303
> > > > ;Pcoupl = berendsen
> > > > ;pcoupltype = semiisotropic
> > > > ;tau_p = 1.0 1.0
> > > > ;compressibility = 4.6e-5 4.6e-5
> > > > ;ref_p = 1.0 1.0
> > > > gen_vel = no
> > > > constraints = all-bonds
> > > > constraint_algorithm = shake
> > > > continuation = yes
> > > >
> > > >
> > > > > From: cara.kreck at student.curtin.edu.au
> > > > > To: gromacs.org_gmx-users at maillist.sys.kth.se
> > > > > Date: Fri, 15 May 2015 14:17:28 +0800
> > > > > Subject: [gmx-users] FW: Lincs all-bonds causing instability in
> > > otherwise stable system
> > > > >
> > > > > Hi,
> > > > >
> > > > > I am having trouble constraining all-bonds with lincs in a system
> that
> > > is stable using shake all-bonds or lincs h-bonds in gromacs 4.6.7. The
> > > previous step using lincs h-bonds gave these energies:
> > > > >
> > > > > Step Time Lambda
> > > > > 200000 200.00000 0.00000
> > > > >
> > > > > Energies (kJ/mol)
> > > > > G96Bond G96Angle Proper Dih. Improper Dih.
> > > LJ-14
> > > > > 1.76317e+04 2.43132e+04 1.39308e+04 1.55159e+03
> > > -2.91002e+03
> > > > > Coulomb-14 LJ (SR) LJ (LR) Coulomb (SR)
> Coulomb
> > > (LR)
> > > > > 4.07560e+05 -4.09017e+03 -1.26700e+04 -5.62217e+05
> > > -3.06259e+03
> > > > > RF excl. Potential Kinetic En. Total Energy
> > > Temperature
> > > > > -1.24406e+05 -2.44368e+05 9.29185e+04 -1.51450e+05
> > > 3.02522e+02
> > > > > Pressure (bar) Constr. rmsd
> > > > > -1.98735e+02 1.95240e-06
> > > > >
> > > > > However, when I switch to lincs all-bonds the 0 step kinetic
> energy,
> > > temperature, and pressure are dramatically increased, even though
> > > continuation=yes and gen_vel=no. The system then quickly crashes. I
> tried
> > > using GMX_MAXCONSTRWARN=-1 and particle decomposition to get past the
> > > errors but then it simply stalls without crashing or outputting
> anything.
> > > > >
> > > > > Step Time Lambda
> > > > > 0 0.00000 0.00000
> > > > >
> > > > > Grid: 8 x 8 x 14 cells
> > > > > Energies (kJ/mol)
> > > > > G96Angle Proper Dih. Improper Dih. LJ-14
> > > Coulomb-14
> > > > > 2.43132e+04 1.39308e+04 1.55159e+03 -2.91002e+03
> > > 4.07560e+05
> > > > > LJ (SR) LJ (LR) Coulomb (SR) Coulomb (LR)
> RF
> > > excl.
> > > > > -4.09017e+03 -1.26700e+04 -5.62217e+05 -3.06259e+03
> > > -1.24406e+05
> > > > > Potential Kinetic En. Total Energy Temperature
> Pressure
> > > (bar)
> > > > > -2.62000e+05 4.92862e+05 2.30862e+05 1.93783e+03
> > > 1.04858e+04
> > > > > Constr. rmsd
> > > > > 1.44601e-02
> > > > >
> > > > > When I instead switch to shake all-bonds (with particle
> decomposition)
> > > I get a similar spike in the kinetic energy etc. but the system is
> able to
> > > recover without crashing. Unfortunately gromacs won't let me use
> pressure
> > > coupling with shake due to my twin-range cut-offs though, so I need to
> find
> > > a way to get lincs to work. I tried switching to lincs all-bonds after
> > > 100ps with shake all-bonds, and there was no longer the spike in
> energies
> > > except for the pressure:
> > > > >
> > > > > Step Time Lambda
> > > > > 100000 100.00000 0.00000
> > > > >
> > > > > Energies (kJ/mol)
> > > > > G96Angle Proper Dih. Improper Dih. LJ-14
> > > Coulomb-14
> > > > > 2.46619e+04 1.38286e+04 1.50388e+03 -2.92380e+03
> > > 4.11368e+05
> > > > > LJ (SR) LJ (LR) Coulomb (SR) Coulomb (LR)
> RF
> > > excl.
> > > > > -3.50689e+03 -1.26695e+04 -5.64100e+05 -2.81505e+03
> > > -1.24427e+05
> > > > > Potential Kinetic En. Total Energy Temperature
> Pressure
> > > (bar)
> > > > > -2.59079e+05 7.73473e+04 -1.81732e+05 3.04114e+02
> > > -4.16208e+02
> > > > >
> > > > > Step Time Lambda
> > > > > 0 0.00000 0.00000
> > > > >
> > > > > Energies (kJ/mol)
> > > > > G96Angle Proper Dih. Improper Dih. LJ-14
> > > Coulomb-14
> > > > > 2.46619e+04 1.38286e+04 1.50388e+03 -2.92380e+03
> > > 4.11368e+05
> > > > > LJ (SR) LJ (LR) Coulomb (SR) Coulomb (LR)
> RF
> > > excl.
> > > > > -3.50690e+03 -1.26695e+04 -5.64100e+05 -2.81533e+03
> > > -1.24427e+05
> > > > > Potential Kinetic En. Total Energy Temperature
> Pressure
> > > (bar)
> > > > > -2.59079e+05 7.73979e+04 -1.81681e+05 3.04313e+02
> > > 5.61153e+02
> > > > > Constr. rmsd
> > > > > 2.45638e-04
> > > > >
> > > > > However, it still crashes with a domain decomposition error. When I
> > > switch back to particle decomposition it stalls again. I then tried
> > > decreasing the time-step to 0.1fs and I got lots of these errors before
> > > stalling:
> > > > >
> > > > > WARNING: Listed nonbonded interaction between particles 2485 and
> 2490
> > > > > at distance 3f which is larger than the table limit 3f nm.
> > > > >
> > > > > This is likely either a 1,4 interaction, or a listed interaction
> inside
> > > > > a smaller molecule you are decoupling during a free energy
> calculation.
> > > > > Since interactions at distances beyond the table cannot be
> computed,
> > > > > they are skipped until they are inside the table limit again. You
> will
> > > > > only see this message once, even if it occurs for several
> interactions.
> > > > >
> > > > > IMPORTANT: This should not happen in a stable simulation, so there
> is
> > > > > probably something wrong with your system. Only change the
> > > table-extension
> > > > > distance in the mdp file if you are really sure that is the reason.
> > > > >
> > > > > There were 156 inconsistent shifts. Check your topology
> > > > > There were 190 inconsistent shifts. Check your topology
> > > > >
> > > > >
> > > > > Whereas with domain decomposition and a 0.1 fs time-step it still
> > > crashes with a domain decomposition error. I don't understand why a
> system
> > > that is perfectly fine with lincs h-bonds and shake all-bonds can be so
> > > problematic with lincs all-bonds. Any suggestions?
> > > > >
> > > > > Thanks,
> > > > >
> > > > > Cara
> > > >
> > > >
> > > > --
> > > > 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.
>
> --
> 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.
>
More information about the gromacs.org_gmx-users
mailing list