(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
6.6 Nudged Elastic Band (NEB) simulations - SECTION 1
By Christina Bergonzo, Carlos Simmerling & Ross Walker
1. Background
In the nudged elastic band method^{1,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 method^{3} 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^{║} = [(k_{i+1}|P_{i+1} - P_{i}| - k_{i}|P_{i} - P_{i-1}|)×t] t
where, if N is the number of atoms per image, F is the force on image i, P_{i} is the 3N dimensional position vector of image i, k_{i} 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(P_{i} - P_{i-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 path^{4} that depends upon the energies, E, of adjacent images:
If (E_{i+1} > E_{i} > E_{i-1}) then t = NORM(P_{i+1} - P_{i})
If (E_{i+1} < E_{i} < E_{i-1}) then t = NORM(P_{i} - P_{i-1})
Or, if E_{i} is a local minima or maxima, then:
If (E_{i+1} > E_{i-1}) then t = NORM[(P_{i+1 }- P_{i})E^{+} + (P_{i }- P_{i-1})E^{-}]
If (E_{i+1} < E_{i-1}) then t = NORM[(P_{i+1 }- P_{i})E^{-} + (P_{i} - P_{i-1})E^{+}]
where:
E^{+} = MAX[|E_{i+1} - E_{i}|,|E_{i-1} - E_{i}|]
E^{-} = MIN[|E_{i+1} - E_{i}|,|E_{i-1} - E_{i}|]
The spring constants can be constant or they can be scaled to move the images closer together in the regions of transition states^{5}:
If (E_{i} > E_{ref}) then k_{i} = k_{max} - Dk(E_{max} - E_{i})/(E_{max} - E_{ref})
If (E_{i} ≤ E_{ref}) then k_{i} = k_{max} - Dk
E_{max} is the highest energy for an image along the path, E_{ref} is the energy of the higher energy endpoint, and k_{max} 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 interest^{6}. 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 Lagrangian^{6}. 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 minimum^{1}.
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 minimiser^{6}.
Recently, a Partial Nudged Elastic Band implementation was discussed^{7} that allows the NEB method to be applied to a user defined subset of the system of interest. This new implementation, available only in Amber11, requires users to define the part of the system to apply NEB force decoupling to, as well as part of the system to RMS fit neighboring images to remove rotational and translational motion. This implementation allows NEB to be used in large systems where a local transition is desired, or in explicitly solvated systems.
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.
7) Bergonzo, C., Campbell, A.J., Walker, R.C., Simmerling, C. 2009. A partial nudged elastic band implementation for use with large or explicitly solvated systems. Intl J Quantum Chemistry 109:3781-3790.