[gmx-users] implicit solvent LINCS errors
Joshua L. Phillips
jphillips7 at ucmerced.edu
Wed Jun 1 20:21:01 CEST 2011
This isn't surprising. Pressure coupling without explicit solvent
doesn't work the way that you are anticipating. You might want to look
at what happened to the system using ngmx or VMD to see what went wrong.
Unless you need to simulate in a box for some reason, you should turn
off PBC as Justin suggested:
ns_type = simple
>>>>>>> rlist = 0
> >>>>>>> nstlist = 0
> >>>>>>> rvdw = 0
> >>>>>>> rcoulomb = 0
> >>>>>>> rgbradii = 0
> >>>>>>> pbc = no
> >>>>>>> comm-mode = angular
>
and use the -pd option for mdrun.
-- Josh
On Wed, 2011-06-01 at 12:28 -0500, Michael D. Daily wrote:
> Scratch what I said last night about pressure coupling - it crashes
> after ~200K timesteps due to fluctuations in the cell volume. I'm
> sticking with NVT for now.
>
> On 5/31/2011 6:14 PM, Michael Daily wrote:
> > The following mdp file produces a successful dynamics run out to 100K
> > steps / 200 ps. What I discovered is this:
> >
> > Using the md integrator, it is necessary to turn off pressure
> > coupling. However, pressure coupling works with sd (Langevin)
> > integrator.
> >
> > Mike
> >
> > ---
> >
> > ; title and include files
> > title = 1EX6-S35P_md1
> > ;cpp = cpp
> > ;include =
> > -I/home/yoo2/myusr/gromacs-3.3.3/share/gromacs/top/ -I./
> > define =
> > ; integrator and input/output setting up
> > ;integrator = md
> > integrator = sd
> >
> > nsteps = 1000000 ; 2 ns
> > dt = 0.002
> > nstxout = 5000
> > nstvout = 5000
> > nstenergy = 500
> > nstxtcout = 500
> > nstlog = 500
> > xtc_grps = System
> > energygrps = System
> > comm_mode = Linear
> >
> > ;implicit solvent
> > implicit_solvent = GBSA
> > gb_algorithm = OBC
> > gb_saltconc = 0.15 ; note - this feature is not functional yet
> > rgbradii = 2.0 ; need larger cut-off than atomistic models
> >
> > ; neighbor searching and vdw/pme setting up
> > nstlist = 10
> > ns_type = grid
> > pbc = xyz
> > ;rlist = 1.4
> > rlist = 2.0
> >
> > ;coulombtype = pme
> > coulombtype = Cut-Off
> > fourierspacing = 0.1
> > pme_order = 6
> > ;rcoulomb = 1.4
> > rcoulomb = 2.0
> >
> > ;vdwtype = switch
> > vdwtype = Cut-Off
> > rvdw_switch = 1.0
> > ;rvdw = 1.2
> > rvdw = 2.0
> >
> > ; cpt control
> > tcoupl = v-rescale
> > tc-grps = System
> > tau_t = 0.4
> > ref_t = 300.0
> > Pcoupl = parrinello-rahman
> > pcoupltype = isotropic
> > tau_p = 1.0
> > compressibility = 4.5e-5
> > ref_p = 1.0
> >
> > ; velocity & temperature control
> > gen_vel = yes
> > gen_temp = 300.0
> > annealing = no
> > constraints = hbonds
> > constraint_algorithm = lincs
> > morse = no
> >
> > ---
> >
> >
> > On 5/30/11 8:45 PM, Michael D. Daily wrote:
> >> Thanks for the recommendations everyone. I tried all of the mdp
> >> changes recommended by Justin (increase rlist, rvdw, etc to 2.0;
> >> change T-coupling to v-scale, and eliminate P-coupling). When I
> >> increased the distance cutoffs, it ran about 30 ps then crashed
> >> instead of crashing immediately. Only when I also turned off
> >> P-coupling did it keep running long-term (now it's up to 500K+ steps
> >> / 1 ns), so apparently the problem was a lack of control of system
> >> volume. Reducing the timestep from 2 fs to 1fs did not work by
> >> itself and was not necessary ultimately to get it to run.
> >>
> >> Josh, thanks for the physical explanation of why the protein might
> >> explode. My starting structure was a minimized xtal structure so I
> >> wouldn't call it "extended." I'll try again with 4.5.4 before I
> >> scale this up any further.
> >>
> >> Mike
> >>
> >>
> >> On 5/30/2011 8:02 PM, Mark Abraham wrote:
> >>> On 31/05/2011 10:54 AM, Justin A. Lemkul wrote:
> >>>>
> >>>>
> >>>> Joshua L. Phillips wrote:
> >>>>> I've found that I often get LINCS warnings like this when starting
> >>>>> from
> >>>>> highly extended conformations when using implicit solvent. The GBSA
> >>>>> surface tension combined with the lack of viscosity (due to the
> >>>>> absence
> >>>>> of explicit water) allows the protein to change conformation much
> >>>>> faster
> >>>>> than LINCS likes by default.
> >>>>>
> >>>>> Normally, I just run a short vacuum simulation (keep the same
> >>>>> settings
> >>>>> as Justin suggested but set GBSA = No) to let the system relax a
> >>>>> little
> >>>>> more before starting a GBSA run (EM is often just not enough).
> >>>>> Usually
> >>>>> only about 50ps. Also, doing this with position restraints can
> >>>>> help slow
> >>>>> down the collapse, and usually results in just enough collapse
> >>>>> that the
> >>>>> following GBSA run will satisfy the default LINCS settings.
> >>>>>
> >>>>> Another option might be to run a short simulation with GBSA turned
> >>>>> on,
> >>>>> but using Langevin dynamics to include some additional friction that
> >>>>> should slow down the initial collapse.
> >>>>>
> >>>>
> >>>> That's generally a good idea. We often use the sd integrator when
> >>>> doing GBSA calculations.
> >>>>
> >>>>> Also, you could fudge with the environment variables associated with
> >>>>> LINCS, but this seems a little dangerous compared to the above two
> >>>>> suggestions.
> >>>>>
> >>>>
> >>>> I wouldn't say "a little dangerous," I'd say "very dangerous" :)
> >>>> If the constraints are failing, it's not necessarily (or usually)
> >>>> their fault. The system is unstable, and trying to just override
> >>>> this will give you a trajectory, but you should be concerned that
> >>>> this trajectory is being ruled by spurious forces.
> >>>>
> >>>>> I would be interested in peoples' opinions and suggestions on how to
> >>>>> handle this issue as it is quite common when starting from highly
> >>>>> extended structures using GBSA. Some proteins are more susceptible
> >>>>> than
> >>>>> others...
> >>>>>
> >>>>
> >>>> Reducing the timestep is a good option. If the dynamics are
> >>>> occurring so fast that the constraints can't keep up, reducing dt
> >>>> can help.
> >>>
> >>> Agreed. Even using all-vs-all kernels, I've observed similar
> >>> symptoms as the OP lately with excessive rotation of v-site terminal
> >>> methyl groups. I can generate stable trajectories by starting my
> >>> equilibration with 0.5 fs timesteps, and switching later.
> >>>
> >>> Mark
> >>>
> >>>>
> >>>> -Justin
> >>>>
> >>>>> I also tend to find that version 4.5.4 is much more stable than even
> >>>>> 4.5.3 for implicit solvent simulations.
> >>>>>
> >>>>> -- Josh
> >>>>>
> >>>>> On Mon, 2011-05-30 at 18:06 -0500, Michael D. Daily wrote:
> >>>>>> Thanks Justin, this is very helpful. I'll attempt these fixes
> >>>>>> tomorrow.
> >>>>>>
> >>>>>> Mike
> >>>>>>
> >>>>>> On 5/30/2011 5:50 PM, Justin A. Lemkul wrote:
> >>>>>>>
> >>>>>>> Michael D. Daily wrote:
> >>>>>>>> Hi all,
> >>>>>>>>
> >>>>>>>> I'm trying to run implicit solvent calculations in gromacs 4.5
> >>>>>>>> with the charmm forcefield. I am able to minimize successfully
> >>>>>>>> and compile for
> >>>>>>> When troubleshooting, it is always advisable to try the latest
> >>>>>>> version (4.5.4) to see if the problem is reproducible. If a
> >>>>>>> pertinent bug has been fixed, there's no use troubleshooting the
> >>>>>>> broken version.
> >>>>>>>
> >>>>>>>> mdrun, but soon after starting, mdrun complains about excessive
> >>>>>>>> rotation in LINCS (see the error printed below that). I also
> >>>>>>>> include my mdp file at the bottom. Can anyone advise me as to
> >>>>>>>> the possible cause of such errors, as it is difficult to
> >>>>>>>> diagnose given that grompp worked fine.
> >>>>>>>>
> >>>>>>> For reference:
> >>>>>>>
> >>>>>>> http://www.gromacs.org/Documentation/Terminology/Blowing_Up#Diagnosing_an_Unstable_System
> >>>>>>>
> >>>>>>>
> >>>>>>> If grompp worked, that just means your coordinate and topology
> >>>>>>> matched and there were no internal conflicts within the .mdp
> >>>>>>> file. It is no guarantee that the resulting simulation will
> >>>>>>> actually work, unfortunately.
> >>>>>>>
> >>>>>>>> --- lincs error ---
> >>>>>>>>
> >>>>>>>> Step 1, time 0.002 (ps) LINCS WARNING
> >>>>>>> Typically an instant LINCS failure indicates insufficient
> >>>>>>> minimization. You said you minimized successfully, but what
> >>>>>>> does this mean? What values did you achieve for Fmax and Epot?
> >>>>>>>
> >>>>>>>> relative constraint deviation after LINCS:
> >>>>>>>> rms 0.000780, max 0.020692 (between atoms 880 and 881)
> >>>>>>>> bonds that rotated more than 30 degrees:
> >>>>>>>> atom 1 atom 2 angle previous, current, constraint length
> >>>>>>>> 606 607 36.7 1.0527 0.1080 0.1080
> >>>>>>>> 614 615 35.6 0.8972 0.1125 0.1111
> >>>>>>>> 614 616 75.3 0.1054 0.1121 0.1111
> >>>>>>>> 880 881 58.0 0.1068 0.1134 0.1111
> >>>>>>>> 880 882 50.4 0.9168 0.1121 0.1111
> >>>>>>>> 889 890 55.3 0.1066 0.1122 0.1111
> >>>>>>>> 889 891 35.3 0.8588 0.1118 0.1111
> >>>>>>>>
> >>>>>>>> ---- mdp file ------------
> >>>>>>>>
> >>>>>>>> ; title and include files
> >>>>>>>> title = 1EX6-S35P_md1
> >>>>>>>> cpp = cpp
> >>>>>>>> include =
> >>>>>>>> -I/home/yoo2/myusr/gromacs-3.3.3/share/gromacs/top/ -I./
> >>>>>>> Any reason you're including files from an ancient Gromacs version?
> >>>>>>>
> >>>>>>>> define =
> >>>>>>>> ; integrator and input/output setting up
> >>>>>>>> integrator = md
> >>>>>>>> nsteps = 1000000 ; 2 ns
> >>>>>>>> ;nsteps = 5000 ; 2 ns
> >>>>>>>> dt = 0.002
> >>>>>>>> nstxout = 5000
> >>>>>>>> nstvout = 5000
> >>>>>>>> nstenergy = 500
> >>>>>>>> nstxtcout = 500
> >>>>>>>> nstlog = 500
> >>>>>>>> xtc_grps = System
> >>>>>>>> energygrps = System
> >>>>>>>> comm_mode = Linear
> >>>>>>>>
> >>>>>>>> ;implicit solvent
> >>>>>>>> implicit_solvent = GBSA
> >>>>>>>> gb_algorithm = Still
> >>>>>>>> gb_saltconc = 0.15
> >>>>>>> FYI, gb_saltconc is nonfunctional. Don't expect it to do
> >>>>>>> anything :)
> >>>>>>>
> >>>>>>>> rgbradii = 1.0
> >>>>>>>>
> >>>>>>>> ; neighbor searching and vdw/pme setting up
> >>>>>>>> nstlist = 10
> >>>>>>>> ns_type = grid
> >>>>>>>> pbc = xyz
> >>>>>>>> ;rlist = 1.4
> >>>>>>>> rlist = 1.0
> >>>>>>>>
> >>>>>>>> ;coulombtype = pme
> >>>>>>>> coulombtype = Cut-Off
> >>>>>>>> fourierspacing = 0.1
> >>>>>>>> pme_order = 6
> >>>>>>>> ;rcoulomb = 1.4
> >>>>>>>> rcoulomb = 1.1
> >>>>>>>>
> >>>>>>>> ;vdwtype = switch
> >>>>>>>> vdwtype = Cut-Off
> >>>>>>>> rvdw_switch = 1.0
> >>>>>>>> ;rvdw = 1.2
> >>>>>>>> rvdw = 1.1
> >>>>>>>>
> >>>>>>> All of these are potentially problematic. Running implicit
> >>>>>>> simulations typically requires longer cutoffs than would
> >>>>>>> normally be needed for explicit solvent simulations. Try
> >>>>>>> rlist=rvdw=rcoulomb=2.0 nm.
> >>>>>>>
> >>>>>>>> ; cpt control
> >>>>>>>> tcoupl = nose-hoover
> >>>>>>> A better choice for initial equilibration would be either
> >>>>>>> V-rescale or Berendsen. I know this can be an issue in explicit
> >>>>>>> solvent, when velocities can oscillate a lot at the outset of a
> >>>>>>> simulation using Nose-Hoover and the simulation box can explode;
> >>>>>>> I don't know if this is such a big deal with implicit, but it
> >>>>>>> can't hurt to try.
> >>>>>>>
> >>>>>>>> tc-grps = System
> >>>>>>>> tau_t = 0.4
> >>>>>>>> ref_t = 300.0
> >>>>>>>> Pcoupl = parrinello-rahman
> >>>>>>> I don't know how an implicit box will respond to pressure
> >>>>>>> coupling, but it would be better to try NVT first and see if
> >>>>>>> it's stable, then try NPT and see if things break down.
> >>>>>>>
> >>>>>>> One option that might be advantageous is to use the all-vs-all
> >>>>>>> kernels for a speed upgrade. You can accomplish this with:
> >>>>>>>
> >>>>>>> rlist = 0
> >>>>>>> nstlist = 0
> >>>>>>> rvdw = 0
> >>>>>>> rcoulomb = 0
> >>>>>>> rgbradii = 0
> >>>>>>> pbc = no
> >>>>>>> comm-mode = angular
> >>>>>>>
> >>>>>>> You'd have to run with mdrun -pd (particle decomposition), but
> >>>>>>> the end result can be quite fast and you avoid potential
> >>>>>>> periodicity effects.
> >>>>>>>
> >>>>>>> -Justin
> >>>>>>>
> >>>>>>
> >>>>>> --
> >>>>>> Michael D. Daily, Ph.D.
> >>>>>> Postdoctoral Fellow
> >>>>>> Qiang Cui group
> >>>>>> Department of Chemistry
> >>>>>> University of Wisconsin-Madison
> >>>>>>
> >>>>>
> >>>>
> >>>
> >>
> >>
> >
> >
>
>
> --
> Michael D. Daily, Ph.D.
> Postdoctoral Fellow
> Qiang Cui group
> Department of Chemistry
> University of Wisconsin-Madison
>
--
Joshua L. Phillips
Ph.D. Candidate - School of Engineering
University of California, Merced
jphillips7 at ucmerced.edu
More information about the gromacs.org_gmx-users
mailing list