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.kr...@student.curtin.edu.au
> To: gromacs.org_gmx-users@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.kr...@student.curtin.edu.au
> > To: gromacs.org_gmx-users@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-requ...@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-requ...@gromacs.org.
