> On Feb 10, 2015, at 10:49 AM, Marek Maly <marek.maly.ujep.cz> wrote:
>
> Dear David and Jason,
>
> first of all thank you for the really quick response !
>
> Since in this case it seems that there is really some memory
> bug inside "mmpbsa_py_nabnmode" code, which might be not so "trivial"
> to identify (or even to correct/solve) especially for the ordinary Amber
> user as I am,
> I am sending you relevant files off this mailing list to allow you i)
> reproduce the issue which I reported ii) analyze the problem (for sure
> 1000% more effectively than I could).
>
> Of course that on both mentioned machines we run 64bit Linux and
> compilation was done
> with these compilers:
>
> gcc version 4.4.7 20120313 (Red Hat 4.4.7-4) (GCC) (my machine)
> gcc version 4.8.3 (Gentoo 4.8.3 p1.1, pie-0.5.9) (cluster)
>
> using 64bit libs (if not there would be inter alia impossible that the
> nmode first phase
> (already after the initial minimization) is working for several hours with
> 24GB RAM requirements).
>
> Regarding drms, that 0.05 was set just to not waste much time with
> minimization and analyze the problem with nmode phase. I know that
> recommended values you mentioned (i.e. at least 1e-6), but in case of
> bigger systems it might be "killing" requirement ... so in such a cases I
> am using 0.01.
>
> The problem is that in the actual version of mmpbsa_py , all the work
> (minimization + nmode) on the
> given frame is done using just 1 CPU core (or even less if the
> multithreading is activated).
>
> So only way in the actual version how to overcome this "weakness" is
> probably to minimize
> selected snaps in parallel manner (before mmpbsa_py analysis) using sander
> to reach drms 1e-6 or less in some acceptable time (especially for the
> bigger molecular systems). Am I right ?
Yes, but then you can’t use those frames to compute implicit solvent energies, since they will no longer have a Boltzmann distribution (you have taken them *off* the free energy surface and dropped them to the potential energy surface).
> Anyway why exactly so low drms values (at least 1e-6) are recommended
> before nmode analysis ?
This has to do with the very nature of the normal mode approximation itself. It is based on diagonalizing the Hessian matrix and treating each of the eigenvectors as an orthogonal “vibrational” mode that is approximately harmonic at the minimum, and using the frequency (eigenvalue) -- which is related to the force constant of the harmonic well -- to calculate the vibrational entropy for that mode.
Now there are 6 modes (5 for linear molecules) that represent purely translational or purely rotational motions -- these should have eigenvalues of 0, ideally, at the true minimum (since there is no potential energy well to rigid body translation or rotation in the absence of external fields). Other than those 6 modes, there will be modes with low frequencies (small force constants) that contribute greatly to the entropy, up to those with very high force constants that contribute very little.
So the modes that contribute *most* to the entropy are those that are closest to 0. These modes are the *most* sensitive to deviations from the minimum, and very small changes to these eigenvalues can result in much *larger* changes to their contributions to the entropy. That’s one overarching problem with normal modes -- the most imprecise eigenvalues contribute the most to the final answer. (Another, of course, being that biomolecules are “soft” and globular and therefore do not rotate rigidly.)
If you are far enough away from the minimum, then more than just those 6 omitted modes will have negative eigenvalues (imaginary frequencies) -- in those cases, nab omits those modes from the entropy calculation altogether... but those modes *would* be among the largest contributors if the structure was more fully minimized, so you can see where the big errors begin to creep in.
> I just for the curiosity did entropy analysis of much smaller complex
> (1174 atoms):
For every atom you add, you add *three* additional degrees of freedom, and three additional possible modes. Therefore, I don’t think that comparisons like this are that scalable.
> a) with drms = 0.01
>
> b) with drms = 1e-6
>
> The mmpbsa_py results (TdS) for these cases are:
>
> a)
> DELTA S total= -49.9888 9.2120
> 3.2569
>
> Total calculation time: 31.836 min.
> Total normal mode calculation time: 16.327 min.
>
>
> b)
> DELTA S total= -48.2373 8.3468
> 2.9510
>
> Total calculation time: 1.146 hr.
> Total normal mode calculation time: 0.532 hr.
>
> As we can see the results are quite similar (difference cca 2%) but the
I encourage you to look at the actual nmode output files. Look at the summary of the eigenvectors and eigenvalues -- how many are omitted? Is it just the 6? How much do those “smaller” modes differ in their entropy contributions between the drms=0.01 and drms=1e-6 cases? Certainly if nmode needed to omit more than the 6 initial modes, you have not minimized enough...
> time requirements
> are here 2.2x higher in case b). The interesting point here is that for
> better minimized
> system (b) not only the initial minimization part lasted longer (which is
> easy understable)
> but also the "Total normal mode calculation time" is significanly (cca 2x)
> longer here why ?
I think what you are seeing may be an artifact. The minimization time is actually included in the normal mode timing, since they are done in the same program at the same time. So all of the time difference should have been attributed to the normal mode calculation time. If you use MPI, though, then the normal mode calculation time may be misleading. The MPI timings are reported individually for *just* the master node. This isn’t a big deal, since usually one frame takes the same amount of time as any other. For normal modes, this isn’t true, since some frames may take *much* longer to minimize. If this frame(s) is being done by a processor *besides* an arbitrarily defined “master”, then the master thread will spend most if its time waiting for the other CPUs to finish. So the amount of time reported by the master on normal modes will not reflect the amount of time spent by the other CPUs. But I never really thought it was that important to get the timers working in parallel...
>
> Just for the curiosity I started mmpbsa_py minimization with drms=1e-6
> condition also in case of
> my problematic big systems (cca 12k atoms), but I am afraid that it will
> take week/(maybe weeks) to reach drms=1e-6 considering already known
> times for reaching drms=0.01.
Your estimates are probably good -- weeks is probably not too far off. You may want to consider the quasi-harmonic approximation (which is PCA-based), as there is no minimization (you build a mass-weighted covariance matrix as a pseudo-Hessian).
HTH,
Jason
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Feb 12 2015 - 06:00:02 PST
