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.kr...@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.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. > -- 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.
