Hi Vickie,
I actually recall something similar happening to me several years ago when I
was doing some QM/MM simulations. Printing every step (or every 5), my
minimization RMS would spike every 25 steps like clockwork. I believe the
way I got around this was to make "nsnb" larger than the length of the
simulation. Does this work for you?
I don't see any REAL reason why this should fix your problem, but nsnb is
the only variable that has a default value of 25, and who knows what's
ACTUALLY happening when the pairlist gets rebuilt.
Another comment is that the Amber minimizer has never gotten a terribly
large amount of praise that I'm aware of (being nice about it). Perhaps you
should try xmin (ntmin=3, described in the manual). I think this is a
better minimization method.
My final comment is in regards to the specific application you're doing...
I'm not sure how meaningful the comparisons that you're doing are. In
principle what you're comparing is the 0 degree binding enthalpy
comparisons. These are only really useful insofar as they correlate well to
binding free energies, since that's what's actually measured in the lab.
Furthermore, you have no guarantee that you're getting the global minimum,
and since the minimizations are not completing and you're not at a minimum,
the comparisons aren't really the 0K binding enthalpies either. A more
appropriate method along these lines is MM/PBSA, but the combination of
MM/PBSA and the MD required to make the snapshots is certainly more
expensive than a minimization. In any case, my 2cents.
Hope this helps,
Jason
(P.S. I was most of the way through this email when Professor Case's
response popped up, so I just finished it up and sent it out with the
duplications intact :) )
On Fri, Aug 20, 2010 at 3:55 AM, Vickie Roberts <vickie.scripps.edu> wrote:
>
> Referring back to a problem described by:
>
> Balvinder Singh (Wed Sep 08 2004 - 15:49:04 PDT)
>
> I have rediscovered this problem in minimizations of 30 complexes
> of 2 proteins that I minimized with the following protocol: initial
> minimiation with strong restraints in vacuo, then with restraints
> using GB, then unrestrained with GB. The problem shows up in the GB
> minimizations, mostly in the final, long, unconstrained one. Unlike the
> industrious Balvinder Singh, I only printed
> out the energies every 10 steps. My goal was to have all 30
> complexes minimized to the same point (as measured by drms), to see if I
> could distinguish 'good' from 'bad' complexes by comparing their energies
> (all complexes are of the same 2 molecules, 2 proteins made up of standard
> amino acids with no other cofactors).
>
> I used Amber 10, Forcefield was leaprc.ff94
>
> For the final unrestrained minimization, I used
>
> GB Energy Minimization
> &cntrl
> imin = 1, maxcyc = 10000,
> ntmin = 1, ncyc = 20,
> ntpr = 10,
> ntb = 0,
> igb = 1,
> cut = 999.9
> drms = 0.00001
>
>
> The drms is rather stringent, so none of the minimizations got there, but,
> looking at my minimization, a criterion of
> drms = 0.01 or even 0.001 is not good enuf,
> the energy could still drop a few kcal/mol further.
>
> My problem is that the bond (and to a lesser extent valence) energy spikes
> at steps xx30 and xx80. Its likely that, as Balvinder Singh found,
> they also occur at steps xx05 and xx55, i.e. every 25 steps, but I did
> print out the energies at those steps. The spikes do not occur all the
> time,
> in fact they don't appear in the first 2 thousand steps or so. But then, as
> the
> minimization goes on, both the size and frequency of the spikes increases.
> I
> have a plot of the energies for all 30 complexes (attached) that shows the
> problem gets steadily worse as the minimization proceeds. My bond
> energy after about step 1000 was generally about 168, except for
> steps xx30 and xx80, where it would go up a small amount (say 2-4 kcal)
> or a very large amount (50 to more than 600 kcal). For example:
>
>
> NSTEP ENERGY RMS GMAX NAME NUMBER
> 9420 -1.1843E+04 9.2814E-03 2.8792E-01 CD 1625
>
> BOND = 168.3836 ANGLE = 739.8685 DIHED =
> 1302.8428
> VDWAALS = -2835.5918 EEL = -23010.7186 EGB =
> -3978.4321
> 1-4 VDW = 940.2447 1-4 EEL = 14830.0193 RESTRAINT =
> 0.0000
>
>
> NSTEP ENERGY RMS GMAX NAME NUMBER
> 9430 -1.1056E+04 1.8257E+01 6.1400E+02 CD 1625
>
> BOND = 797.9902 ANGLE = 889.1781 DIHED =
> 1306.3816
> VDWAALS = -2836.4958 EEL = -23002.0935 EGB =
> -3975.0551
> 1-4 VDW = 940.5651 1-4 EEL = 14823.1138 RESTRAINT =
> 0.0000
>
>
> NSTEP ENERGY RMS GMAX NAME NUMBER
> 9440 -1.1843E+04 1.6190E-02 3.0229E-01 CD2 2609
>
> BOND = 168.4287 ANGLE = 739.9444 DIHED =
> 1302.8351
> VDWAALS = -2835.6302 EEL = -23010.5702 EGB =
> -3978.4084
> 1-4 VDW = 940.1704 1-4 EEL = 14829.8270 RESTRAINT =
> 0.0000
>
> So you get this huge hit to the RMS (9.2814E-03 at step 9420 going to
> 1.8257E+01 at step 9430) that only gets partially repaired at step
> 9440. Seems like this must be making the minimization much less efficient.
> After 10000 steps of a 4800 atom system (all standard protein residues) I
> am
> at
>
>
> NSTEP ENERGY RMS GMAX NAME NUMBER
> 10000 -1.1844E+04 2.2508E-02 6.7566E-01 CD 4525
>
> BOND = 168.4205 ANGLE = 739.7665 DIHED =
> 1302.4882
> VDWAALS = -2835.9819 EEL = -23010.4237 EGB =
> -3978.8126
> 1-4 VDW = 940.0621 1-4 EEL = 14830.2729 RESTRAINT =
> 0.0000
>
>
>
> (and yes, step 9980 had a big spike in bond and valence energy)
>
> The responses that Balvinder Singh got were not very satisfying.
> Ross Walker replied that something was probably wrong with the system
> or the installation of Amber. I agree that something is wrong, but
> I dont think there is a problem with either the installation or
> my particular system.
> Instead i am wondering if its a problem that is only seen if you
> 1) have ntpr set so that you look at energies of steps xx05, xx30,
> xx55, or xx80 AND
> 2) do a sufficiently long (conjugent gradient) minimization.
>
> I have not tried the AMBER test cases, but I have been running AMBER
> a long time, previously versions in the Case lab. Unfortunately,
> all of my old outputs had ntpr=100 or greater, so if the problem
> occurred i would not see it.
>
> Could it be the system? Well, i know that my protein-protein
> complexes had some close contacts, though they relaxed well during
> the initial restrained minimizations. Just in case, i took a look at the
> atoms with the GMAX for the bad steps, thinking they might be
> in the interface.
>
> I looked at steps ending in 30 and 80 that showed no jumps after
> the bond energy had hit 167. The GMAX atom for these steps
> varied a lot in type - H, N, C - residue, etc.
>
> I looked at about 50 steps numbered xx30 and xx80 that showed funny
> spikes,
> from little (2-4) to big jumps, the corresponding GMAX atoms ALL were
> either
>
> ASP CG or GLU CD
>
> So the same atom type. Most were NOT in the interface
> and they could be from either molecule.
>
> I ran my vacuum minimization for just 200 steps, and did not see the
> spike problem there, but Balvinder Singh had the problem with both
> igb = 1 and using a water box, I believe, so it may not be GB
> dependent.
>
> The other response was from Andreas Svrcek-Seiler
>
> ''Sometimes the minimization routine does a trial step, which
> might be "too large", resulting in energy spikes. Whether (and how often)
> that happens depends on the algorithm as well as on implementation
> details, however it is normal and harmless.
> regards
> Andreas''
>
> Since this problem happened for both me and Balvinder Singh, both using
> a very typical and straightforward protocol, with conjugent gradient
> minimization, i am not convinced. Especially since the problem becomes
> worse as the system becomes more minimized. The advantage of GB is that
> it takes solvation into account, so comparing relative energies of
> my complexes might be useful. However, knowing that my complexes had
> large and varying jumps in energy at step 9980 makes the whole thing
> seem less convincing. Its too bad, because the energy comparison does
> appear to favor correct (RMSD < 4A) over incorrect (RMSD > 20A) complexes.
>
> I am happy to send coords and parameter files, but it would be nice if
> someone tried a familiar ~5000 atom system, 10000 steps of standard
> minimization, with GB on, printing out energies every 5 steps (to catch the
> xx05 and xx55 spikes too). If the spikes show up there, then start to
> worry about if its more general problem. If the problem does not show up,
> then maybe I can get some help figuring out why my system has the
> problem.
>
> Attached is a plot of the total energy for all 30 minimizations saved
> every 10 steps. It makes it very clear that the problem gets worse
> as the minimizations proceed. All complexes have final minimized energies
> between -11,700 and -11,900. ALL of the energies that are not close to
> the minimized energy occur at steps xx30 and xx80.
>
> Thanks,
> vickie
>
>
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
>
--
Jason M. Swails
Quantum Theory Project,
University of Florida
Ph.D. Graduate Student
352-392-4032
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Aug 20 2010 - 05:00:07 PDT
