(Note: These tutorials are meant to provide
illustrative examples of how to use the AMBER software suite to carry out
simulations that can be run on a simple workstation in a reasonable period of
time. They do not necessarily provide the optimal choice of parameters or
methods for the particular application area.)
Copyright Ross Walker 2006
AMBER ADVANCED WORKSHOP
TUTORIAL A5 - SECTION 1
The Nudged Elastic Band Approach to Finding the Lowest
Energy
Pathway Between two States
By Ross Walker
1. Background
In the nudged elastic band method1,2 the path for a conformational change is approximated with a series of images of the molecule describing the path. Minimisation of the entire system, but with the end point structures fixed, provides a minimum energy path. Each image between the end points is connected to the previous and next image by 'springs' along the path that keep each image from sliding down the energy landscape onto adjacent images. NEB is derived from the plain elastic band method3 which added the spring forces to the potential of energy surface and minimised the energy of the system. The plain elastic band method found low energy paths, but tended to cut corners in the energy landscape. NEB prevents corner cutting by truncating the spring forces in directions perpendicular to the tangent of the path. Furthermore, the forces from the molecular potential are truncated along the path, so that images remain evenly spaced along the path.This leads to:
F = F┴ + F║
F┴ = -ÑV(P) + ((ÑV(P))×t)t
F║ = [(ki+1|Pi+1 – Pi| - ki|Pi – Pi-1|)×t] t
where, if N is the number of atoms per image, F is the force on image i, Pi is the 3N dimensional position vector of image i, ki is the spring constant between image i-1 and image i, V is the potential described by the force field, and t is the 3N dimensional tangent unit vector that describes the path. The simplest definition of t is:
t = NORM(Pi – Pi-1)
This definition leads to instability in the path caused by kinks that occur where the magnitude of F║ is much larger than the magnitude of F┴. A more stable tangent definition was derived to prevent kinks in the path4 that depends upon the energies, E, of adjacent images:
If (Ei+1 > Ei > Ei-1) then t = NORM(Pi+1 – Pi)
If (Ei+1 < Ei < Ei-1) then t = NORM(Pi – Pi-1)
Or, if Ei is a local minima or maxima, then:
If (Ei+1 > Ei-1) then t = NORM[(Pi+1 – Pi)E+ + (Pi – Pi-1)E-]
If (Ei+1 < Ei-1) then t = NORM[(Pi+1 – Pi)E- + (Pi – Pi-1)E+]
where:
E+ = MAX[|Ei+1 – Ei|,|Ei-1 – Ei|]
E- = MIN[|Ei+1 – Ei|,|Ei-1 – Ei|]
The spring constants can be constant or they can be scaled to move the images closer together in the regions of transition states5:
If (Ei > Eref) then ki = kmax – Dk(Emax – Ei)/(Emax – Eref)
If (Ei ≤ Eref) then ki = kmax – Dk
Emax is the highest energy for an image along the path, Eref is the energy of the higher energy endpoint, and kmax and Dk are parameters with units of force per length. Because the spring force applies only in directions along the path and because the potential of the energy surface is zeroed along the path, the calculation is relatively insensitive to the magnitude of the spring constants. Care must be taken, however, to select a spring constant that does not result in higher frequency motions than those found in the system of interest6. At each step, before calculating the spring forces that compose F║, the images, starting with the second image, are rotated and translated onto the previous image to find the RMSD minimum.Energy minimisation of the path is complicated by the fact that the forces are truncated according to the tangent direction, making it impossible to define a Lagrangian6. Conjugate gradient minimisation, therefore, cannot be used to find the minimum energy path. A velocity Verlet algorithm for quenched molecular dynamics has been used to find the minimum1.
With this method, the component of the velocity parallel to the force is kept, but perpendicular components are scaled:
If (v×f ≥ 0) then v = (v×f)f
If (v×f < 0) then v = x(v×f)f
where f is the 3N-dimensional unit force vector, v is the 3N-dimensional velocity vector, and x is a scaling factor less than one. Recently, a super-linear minimisation method was described using an adopted basis Newton-Raphson minimiser6.
The implementation of NEB in sander allows minimisation by simulated annealing. This requires no hypothesis for a starting path, but does require careful judgment of the temperature and length of time required to populate the minimum energy path. The initial coordinates can have multiple copies of the structure superimposed on the start and endpoints. When adjacent structures are superimposed, the tangent, t, is 0 in every direction. This case is explicitly handled so that the calculation is stable.
1) Jónsson H, Mills G, Jacobsen KW. 1998. Nudged elastic band method for finding minimum energy paths of transitions. In: Berne BJ, Ciccoti G, Coker DF, eds. Classical and Quantum Dynamics in Condensed Phase Simulations. Singapore: World Scientific. pp 385-404.
2) Mills G, Jonsson H. 1994. Quantum and thermal effects in H2 dissociative adsorption: Evaluation of free energy barriers in multidimensional quantum systems. Physical Rev Lett 72:1124-1127.
3) Elber R, Karplus M. 1987. A method for determining reaction paths in large molecules: Application to myoglobin. Chem Phys Lett 139:375-380.
4) Henkelman G, Jónsson H. 2000. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. J Chem Phys 113:9978-9985.
5) Henkelman G, Uberuaga BP, Jónsson H. 2000. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J Chem Phys 113:9901-9904.
6) Chu J, Trout BL, Brooks BR. 2003. A super-linear minimization scheme for the nudged elastic band method. J Chem Phys 119:12708-12717.