AMBER Archive (2006)

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

From: Ross Walker (ross_at_rosswalker.co.uk)
Date: Thu Sep 28 2006 - 10:29:29 CDT

Dear Pietro

> 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.

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.

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.
 
> 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

All the best
Ross

/\
\/
|\oss Walker

| HPC Consultant and Staff Scientist |
| San Diego Supercomputer Center |
| Tel: +1 858 822 0854 | EMail:- ross_at_rosswalker.co.uk |
| http://www.rosswalker.co.uk | PGP Key available on request |

Note: Electronic Mail is not secure, has no guarantee of delivery, may not
be read every day, and should not be used for urgent or sensitive issues.

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu