[gmx-users] Lincs all-bonds causing instability in otherwise stable system
Cara Kreck
cara.kreck at student.curtin.edu.au
Thu May 21 10:15:56 CEST 2015
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.
More information about the gromacs.org_gmx-users
mailing list