Nudged Elastic Band calculations

1   Self-diffusion on the Al(110) surface

The script neb1.py will find the minimum energy path (MEP) for diffusion along the rows. If you fill in 9 instead of ?? in the line:

for i in range(??):

you will get 9 images between initial and final states, which is more than enough. The barrier is 0.11 eV:

In order to find the barrier for diffusion across the rows, you must change this line defining the final state:

final[-1].SetCartesianPosition((a / 2, 3 * b / 2, h))

to:

final[-1].SetCartesianPosition((3 * a / 2, b / 2, h))

The barrier for this process is 0.56 eV:

The MEP for the exchange process is found by replacing the above line by:

final[24].SetCartesianPosition((3 * a / 2, 3 * b / 2, h))
final[-1].SetCartesianPosition((a, b, 0))

The exchange process is more complex than the two simple hopping events, and we need to relax the path using 15 x 10 steps - so change 5 to 15 in this line:

for i in range(5):

The barrier is 0.15 eV:

Notice that the path for the exchange process goes through a meta stable state. Plot the energy along the path like this:

plottrajectory jump3.traj

2   Analytical problem on transition state theory (TST)

The four basic assumptions of TST:

1. The dynamics of the atoms can be described by Newton's equation, i.e. classical dynamics.
2. The Born-Oppenheimer approximation is valid, i.e. the time scale of the motion of electrons and atoms is different enough that we can solve for the electronic degrees of freedom for fixed positions of the nuclei. This gives a potential energy surface (PES) for the motion of the nuclei.
3. The transitions are slow enough that a Boltzmann distribution of energy is established and maintained for each degree of freedom of the reactant(s).
4. A region, the so called transition state (TS), of the PES in between reactants and products can be defined in such a way that any reactive trajectory must go through this region and if the system makes it there and is on its way towards products it will continue to go to the product region of the PES and stay there for an extended time. The potential surface should have a simple and narrow barrier (without dips) in between the initial and final state. Then TST can be expected to be a good approximation.

Here is how to plot the potential using gnuplot:

gnuplot> Vs=0.2
gnuplot> b=3
gnuplot> a=2
gnuplot> V(x,y,z)=Vs*(exp(-cos(2*pi*x/b)-cos(2*pi*y/b)-2*a*z)-2*exp(-a*z))
gnuplot> splot [-b:b] [-1:0] V(x,0,y)

The minima are at:

(x, y, z) = (ib, jb, -2/α),

where i and j are integers. The saddle points are located at:

(x, y, z) = (ib/2, jb/2, 0),

where i + j is an odd integer.

The Taylor expansion close to the minima at (0, 0, -2/α) is:

V(x, y, z) = V_se²[-1 + 2(π/b)²(x² + y²) + α²(z +2/α)²],

and close to the saddle point at (b/2, 0, 0) we have:

V(x, y, z) = V_s[-1 - 2(π/b)²(x - b/2)² + 2(π/b)²y ² + α:sup:2z²]

The energy barrier is:

E_a = V_s(e² - 1) = 1.28 eV

The frequencies can be found from the force constants like this:

ν = (k/m)^1/2/(2π).

From the Taylor expansion close to the minimum, we get:

k_x = 4V_s(πe/b)²

k_y = 4V_s(πe/b)²

k_z = 2V_s(αe)²

Plugging in the numbers, we get the frequencies:

3.98 x 10¹³/sec

3.98 x 10¹³/sec

5.38 x 10¹³/sec

At the saddle point (transition state) there are two positive force constants:

k_y = 4V_s(π/b)²

k_z = 2V_sα²

corresponding to the frequencies:

1.46 x 10¹³/sec

1.98 x 10¹³/sec

The third frequency is imaginary, because the force constant k_x is negative. This mode is along the x-axis which is the reaction coordinate.

The prefactor in the Arrhenius expression is equal to the product of the three frequencies at the minimum divided by the product of the two frequencies at the transition state:

(V_s/m)^1/2e³/b = 2.95 x 10¹⁴/sec

The correction to the activation energy from zero point motion is -0.20 eV, so the final expression for the rate is

k = 2.95 x 10¹⁴/sec exp[-1.07 eV/(k_BT)]

The average length of time in between diffusion hops at room temperature (298 K) is 4805 sec and 0.11 msec at 400 K.