RE: AMBER: Questions/observations about Nudged Elastic Band (NEB)

From: Pietro Amodeo <>
Date: Fri, 29 Sep 2006 11:21:43 +0200 (CEST)

Dear Ross,

thank you for the fast and exhaustive (as much as possible, due the
prototype nature of the NEB implementation) answer.

>> 1) after a careful reading of both AMBER9 manuals and Ross Walker's
>> Tutorial 7 from "Amber Workshop Tutorials" I understood that
>> NEB approach
>> can be performed using multiple frames from an approximate path as
>> starting guess for min. energy paths. In the case of our simulated
>> molecule, we obtained the desired conformational transition during a
>> standard MD simulation using implicit GB/SA solvation. Now,
>> we'd like to
>> obtain a NEB path using frames from this MD simulation. What
>> MD protocol
>> (in particular, temperature vs. time profile and frequency of frame
>> sampling from original MD run) is suggested in this case?
> This is the million dollar question. The problem is that the simulated
> annealing version of NEB in Amber 9 is new and because of this should be
> considered still a prototype method. That is I don't think we have a full
> understanding yet of the sampling issues involved or what is the best
> protocol to use. I don't think you can do much damage by over sampling in
> your simulated annealing but under sampling (running the simulated
> annealing
> for too short a time) will cause problems. Where the balance between this
> lies is very much dependent on your system size, degrees of freedom,
> number
> of images etc etc. Suffice to say that the larger the system is and the
> more
> images you have the longer you will need to run for. How long though would
> only be a guess.
> I would suggest you start your system at 0K and then do something along
> the
> follow lines. Note the large values of gamma_ln (this is required for NEB)
> to give the system viscosity. For a regular MD simulation these values are
> orders of magnitude too large:
> Step 1 -> Heat linearly from 0 to 300K dt=0.0005, skmin=10, skmax=10,
> gamma_ln=100,
> nstlim = 40,000 (20ps)
> Step 2 -> Heat linearly to 1000K, dt=0.002, skmin=20, skmax=20,
> gamma_ln=1000,
> nstlim = 25,000 (50ps)
> Step 3 -> Equilibrate at 1000K, dt=0.002, skmin=25, skmax=25,
> gamma_ln=1000,
> nstlim = 50,000 (100ps)
> Step 4 -> Cool 1000K to 300K linearly, dt=0.002, skmin=25, skmax=25,
> gamma_ln=1000,
> nstlim = 50,000 (100ps)
> Step 5 -> Equilibrate at 300K, dt=0.002, skmin=25, skmax=25,
> gamma_ln=1000,
> nstlim = 50,000 (100ps)
> Step 6 -> Cool 300K -> 0K linearly, dt=0.002, skmin=25, skmax=25,
> gamma_ln=1000,
> nstlim = 50,000 (100ps)
> Step 7 -> 0K, dt=0.002, skmin=25, skmax=25, gamma_ln=5000, nstlim = 50,000
> (100ps),
> vv=1, vfac=0.1, ntwx=25000
> I have had success with this protocol. Although the time lengths are
> really
> just a guess. You may want to increase the length of sampling in the
> middle
> if you can afford it. You should also probably run 10 or so independent
> runs
> with different random seeds to test how well your sampling is converged.
> You
> can then see how many different reactions paths you get and compare them.
I'll try your protocol, possibly with systematic exploration of some of
the involved parameters, and, if you are interested, I can send you the
main results/issues found in the analysis.

> In terms of building your initial guess at the pathway there ar eno hard
> and
> fast rules. You just have to provide the code with sufficient input
> structures. And indeed this is very much a "guess" at the pathway so the
> structures really just need to be reasonable representations of points
> along
> the pathway. So I would just pick what you thing are representative of the
> transition from your MD simulation and string them together.

I just wondered if some "quasi-minimization-only" protocol could be
meaningful for opportunely-sampled structures from a plain MD simulation
involving the transition under examination.

> You may also want to just try it with half the images piled up on one end
> and half on the other end (I.e. no guess at the pathway). We have had some
> success with this. You obviously need to equilibrate for longer but if you
> are lucky it won't get trapped in a local minimum due to the high
> temperature.

Calculations with this approach are already running, as we were interested
in a comparison between the two approaches and for some of the systems
under examination we haven't the MD trajectory including the transition.

>> 2) looking at the source code of AMBER9 implementation of the quenched
>> velocity Verlet algorithm used for NEB minimization (see below), we
>> observed that, differently from what stated in the manual,
>> when the dot
>> product v dot f<0 , velocities are simply zeroed and not
>> scaled according
>> to eq. (6.34) v=x(v dot f)f and, consequently, the input
>> scaling factor
>> VFACT is never used (or, equivalently, the program behaves as if VFACT
>> value were hardwired to 0). I've found no explanation for
>> this discrepance
>> either in docs, or in man. upgrades, or in mailing list, or
>> in bugfixes.
>> >>>> File: runmd.f
>> >>>> Subroutine: quench
>> > if (dotproduct>0.0d0) then
>> > v(1:3*natom) = dotproduct*f(1:3*natom)*force
>> > else
>> > v(1:3*natom) = 0.0d0
>> > end if
> This is a typo in the manual. It should state:
> if(v.f>=0) then v = (v.f)f
> if(v.f<0) then v = 0
> This is what the code does and is what should happen. I will update the
> manual errata page.


> Note you may also want to read the following paper describing the NEB
> implementation in Amber 9:
> Mathews, D.H., Case, D.A. J. Mol. Biol. (2006) 357, 1683-1693

We already read the paper and most of the information was desumed from it.

Thank you again,


Dr. Pietro Amodeo
Istituto di Chimica Biomolecolare del CNR
Comprensorio "A. Olivetti", Edificio 70
Via Campi Flegrei 34
I-80078 Pozzuoli (Napoli) - Italy
Phone      +39-0818675072
Fax        +39-0818041770
The AMBER Mail Reflector
To post, send mail to
To unsubscribe, send "unsubscribe amber" to
Received on Sun Oct 01 2006 - 06:07:16 PDT
Custom Search