[gmx-users] Amber force fields giving LINCS errors
Justin A. Lemkul
jalemkul at vt.edu
Thu Sep 9 23:04:22 CEST 2010
TJ Mustard wrote:
<snip>
> *Positional restraint mdp file:*
>
>
>
> define = -DPOSRES
> constraints = all-bonds
> integrator = md
> dt = 0.004 ; ps
Unless you're using virtual sites (which, per your commands below, you are not)
this time step is inappropriately large and is a likely source of instability.
With plain constraints, 2 fs is more appropriate.
> nsteps = 2500 ; total ps = dt*nsteps
> nstcomm = 1
> nstxout = 200 ; output coordinates every ps = dt*nstxout
> nstvout = 1000 ; output velocities every ps = dt*nstvout
> nstfout = 0
> nstenergy = 10
> nstlog = 10
> nstlist = 10
> ns_type = grid
> rlist = 0.9
> coulombtype = PME
> rcoulomb = 0.9
> rvdw = 1.0
> fourierspacing = 0.12
> fourier_nx = 0
> fourier_ny = 0
> fourier_nz = 0
> pme_order = 6
> ewald_rtol = 1e-5
> optimize_fft = yes
> ; Berendsen temperature coupling is on
> Tcoupl = v-rescale
> tau_t = 0.1 0.1 0.1 0.1 0.1
> tc_grps = protein RNA SOL NA CL
Also a major concern - never couple solvent and ions separately! See here:
http://www.gromacs.org/Documentation/Terminology/Thermostats
> ref_t = 310 310 310 310 310
> ; Pressure coupling is on
> pcoupl = berendsen
> pcoupltype = isotropic
> tau_p = 0.5
> compressibility = 4.5e-5
> ref_p = 1.0
> ; Generate velocities is on at 310K (core body temp)
> gen_vel = yes
> gen_temp = 310.0
> gen_seed = 173529
>
>
>
>
>
> *main molecular dynamics mdp file:*
>
>
>
<snip>
> define = -DPOSRES
Restraints during data collection, or did you paste the wrong file?
<snip>
> ; Generate velocities is on at 310K (core body temp)
> gen_vel = yes
Re-generating velocities after equilibration defeats the whole purpose of
equilibrating. Per your commands below, you are not preserving the velocities
(grompp -t with .trr or .cpt file) obtained during your PR equilibration.
> gen_temp = 310.0
> gen_seed = 173529
>
>
>
>
>
> *And my script for running the whole process:*
>
>
>
> pdb2gmx -f receptor.pdb -o receptor.gro -p receptor.top
>
> editconf -bt cubic -f receptor.gro -o receptor.gro -c -d 1.5
>
> genbox -cp receptor.gro -cs spc216.gro -o receptor_b4ion.gro -p receptor.top
>
> grompp -f em.mdp -c receptor_b4ion.gro -p receptor.top -o receptor_b4ion.tpr
>
> genion -s receptor_b4ion.tpr -o receptor_b4em.gro -neutral -conc 0.0001
> -pname NA -nname CL -g receptor_ion.log -p receptor.top
>
>
>
> *Here I select the SOL for ions...*
>
>
>
> grompp -f em.mdp -c receptor_b4em.gro -p receptor.top -o receptor_em.tpr
>
> mdrun -v -s receptor_em.tpr -c "$base"_after_em.gro -g emlog.log
>
> grompp -f pr.mdp -c receptor_after_em.gro -p receptor.top -o receptor_pr.tpr
>
> mdrun -v -s receptor_pr.tpr -o receptor_pr.trr -e pr.edr -c
> receptor_after_pr.gro -g prlog.log -cpi state_pr.cpt -cpo state_pr.cpt
>
> grompp -f md.mdp -c receptor_after_pr.gro -p receptor.top -o receptor_md.tpr
>
> mdrun -s receptor_md.tpr -o receptor_md.trr -c receptor_after_pr.gro -g
> md.log -e md.edr -cpi state_md.cpt -cpo state_md.cpt
>
> *Where receptor is of course my protein/RNA pdb name*
>
>
>
> I always make sure that the pdb doesn't give me notes or warnings or
> errors of course in the pdb2gmx step. Most of the time I minimize to my
> computers "machine precision."
>
> **
>
> *This is the return I get from the "EM" step:*
>
>
>
> Stepsize too small, or no change in energy.
> Converged to machine precision,
> but not to the requested precision Fmax < 200
>
> Double precision normally gives you higher accuracy.
> You might need to increase your constraint accuracy, or turn
> off constraints alltogether (set constraints = none in mdp file)
>
> writing lowest energy coordinates.
>
> Steepest Descents converged to machine precision in 324 steps,
> but did not reach the requested Fmax < 200.
> Potential Energy = -1.9868888e+07
> Maximum force = 8.7276396e+03 on atom 19080
This is a very high Fmax, indicating that you have bad clashes in your system
that, when combined with the parameters above, surely leads to an unstable
system. The fact that Gromos96 works is a bit curious, but you have too many
potential pitfalls listed here to specifically diagnose the difference based on
this information alone.
-Justin
> Norm of force = 3.8592850e+01
>
>
>
>
>
> *I will get this error sometime in the "PR" step of my script:*
>
>
>
> step 0
> Step 11 Warning: pressure scaling more than 1%, mu: 1.03087 1.03087 1.03087
>
> Step 11 Warning: pressure scaling more than 1%, mu: 1.03087 1.03087 1.03087
>
> Step 21 Warning: pressure scaling more than 1%, mu: 0.849053 0.849053
> 0.849053
>
> Step 21 Warning: pressure scaling more than 1%, mu: 0.849053 0.849053
> 0.849053
>
> Step 22, time 0.088 (ps) LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.005723, max 0.026407 (between atoms 6545 and 6542)
> bonds that rotated more than 30 degrees:
> atom 1 atom 2 angle previous, current, constraint length
>
>
>
>
>
>
>
> I hope this is enough information. Any help would be much appreciated.
>
>
>
> TJ Mustard
> Email: mustardt at onid.orst.edu
> Cell: 509-879-4173
>
> On September 9, 2010 at 9:11 PM "Justin A. Lemkul" <jalemkul at vt.edu> wrote:
>
> >
> >
> > TJ Mustard wrote:
> > >
> > > First off I am using gromacs 4.5. I will also post all of my files and
> > > errors if they help.
> > >
> > >
> > >
> > > If I run a protein in GROMOS96 all my md runs complete succesfully. But
> > > if I change to any of the AMBER force fields I get LINCS errors in my
> > > positional restraint md run. I have tried using shake, 1 fs step sizes,
> > > -heavyh, and many more. Does anyone know what is going on here?
> > >
> >
> > A complete (but not overly lengthy) post will save everyone a lot of
> time.
> > Based on the information you've provided here, I see now way to
> diagnose the
> > problem. The most important information to post would be your .mdp file.
> > Certain settings can influence stability. A description of the hardware,
> > compilers used, etc. can also be useful.
> >
> > -Justin
> >
> > >
> > >
> > > The reason I want to use AMBER is the fact that I want to run md on the
> > > 30s rybosome and amber converts RNA much easier than GROMOS force
> fields.
> > >
> > >
> > >
> > >
> > >
> > > Thank you in advance,
> > >
> > >
> > >
> > > TJ Mustard Email: mustardt at onid.orst.edu
> > > Cell: 509-879-4173
> > >
> > >
> > >
> > >
> >
> > --
> > ========================================
> >
> > Justin A. Lemkul
> > Ph.D. Candidate
> > ICTAS Doctoral Scholar
> > MILES-IGERT Trainee
> > Department of Biochemistry
> > Virginia Tech
> > Blacksburg, VA
> > jalemkul[at]vt.edu | (540) 231-9080
> > http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
> >
> > ========================================
> > --
> > gmx-users mailing list gmx-users at gromacs.org
> > http://lists.gromacs.org/mailman/listinfo/gmx-users
> > Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> > Please don't post (un)subscribe requests to the list. Use the
> > www interface or send it to gmx-users-request at gromacs.org.
> > Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >
>
>
>
>
>
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
More information about the gromacs.org_gmx-users
mailing list