Citation
Krishna Mohan, G. P. (2010, October 13). Dissociative chemisorption of methane on Ni(111). Retrieved from https://hdl.handle.net/1887/16033
Version: Corrected Publisher’s Version
License: Licence agreement concerning inclusion of doctoral thesis in the Institutional Repository of the University of Leiden
Downloaded from: https://hdl.handle.net/1887/16033
Note: To cite this publication please use the final published version (if applicable).
Chapter 2
Review of the theoretical methods
In this chapter a brief overview of all the major theoretical methods used in this thesis is presented. Section 2.1 illustrates the Born-Oppenheimer approximation.
In Section 2.2 the overview of plane-wave density functional theory is provided along with pseudopotential and projector augmented plane wave (PAW) approxi- mations. The nudged elastic band method to obtain the reaction path and tran- sition state along with different optimization methods are presented in Section 2.3. The theory of the vibrational normal modes is reviewed in Section 2.4. The interpolation scheme for the potential energy surface (PES) is described in Section 2.5. Finally, in Section 2.6 classical dynamics is reviewed and the last Section, 2.7, covers the aspects of quantum dynamics and MCTDH method.
2.1 The Born-Oppenheimer approximation
The Hamiltonian of an electronic system determines the energy of the system.
Although a complete Hamiltonian for a molecule is easily determined, the resulting Schrödinger equation is very difficult to solve, due to the multidimensional nature of the problem. A standard approximation to solve this problem is to separate the nuclear and electronic degrees of freedom. Since the mass of the nuclei is considerably heavier than that of the electrons, it can be assumed that the electrons will respond ‘instantaneously’ to the nuclear coordinates. This approximation is called the Born–Oppenheimer (BO) or adiabatic approximation [1]. This idea allows one to treat the nuclear coordinates as parameters. For many condensed matter systems, this presumption is a reasonable approximation.
The non-relativistic molecular Hamiltonian operator for a system of nuclei and electrons can be written with usual notations as,
H = −ˆ ~2 2
PN I=1
∇2I MI − ~2
2m
Pn
i=1∇2i −PNI=1Pni=1Ze2 ri,I +
PN I=1
PN0 J =1
"
ZIZJe2 rI,J
#
+Pni=1Pnj=10
"
e2 ri,j
#
,
(2.1)
where m is the mass of the electron, MI is the mass of the I0th nucleus, N is the number of nuclei and n is the number of electrons. We can write Eq. 2.1 into a simple form, viz.,
21
H = ˆˆ TN + ˆTe+ ˆVN e+ ˆVN N + ˆVee. (2.2) Here, ˆTN and ˆTe are the kinetic energies of the nuclei and of the electrons, respectively. ˆVN e and ˆVee denote the electron - nuclei attraction and electron - electron repulsion energies, and ˆVN N represents the nucleus - nucleus repulsion energy. The time-independent Schrödinger equation of the system is written as,
HΨ(rˆ i, RI) = EΨ(ri, RI). (2.3) Here, Ψ(ri, RI) represents the complete eigenfunction of ˆH and is a function of the electronic (ri) and nuclear (RI) positions.
As the mass of the proton is about 1840 times that of the electrons, a first approximation is to assume the nuclei to be stationary with respect to the electrons.
This assumption readily suggests that in solving for the electronic motion one can neglect the kinetic energy of the nuclei ( ˆTN). Additionally, it allows us to write the full eigenfunction in a seperable form as,
Ψ(ri, RI) = φ(ri, {RI})χ(RI). (2.4) In Eq.2.4, φ(ri, {RI}) is an electronic wavefunction that parametrically depends on the nuclear positions while χ(RI) is the nuclear wavefunction. Using the above expression, one can set up the Schrödinger equation for the nuclei and for the electrons, seperately. In other words, the Schrödinger equation for electrons is written as,
[ ˆTe+ ˆVN e+ ˆVee+ ˆVN N]φ(ri, {RI}) = Un({RI})φ(ri, {RI}) (2.5) and that for the nuclei is,
[ ˆTN + Un({RI})]χ(RI) = Eχχ(RI). (2.6) In the above equations, Un({RI}) is the electronic energy at fixed nuclear co- ordinates. Eχ is the eigenvalue of the nuclear wavefunction. It should be noted that the χ(RI) describes the nuclear motions like rotations, vibrations and trans- lations. In BO approximation, the changes in the energies of the electrons in molecules that are brought about by nuclear motions are essentially adiabatic.
This also shows that the interatomic potentials and forces are only functions of the nuclear positions.
2.2 The gas-surface interaction
In this study the gas-surface electronic interaction plays a major role. This interac- tion potential is needed to generate further derived quantities and/or observables to describe the system. This interaction potential is calculated by the density functional theory (DFT), within the Born-Oppenheimer approximation. DFT is currently one of the most successful electronic structure theories for the large systems. The main aspects of plane-wave DFT and its subsequent features are reviewed in the below subsection.
2.2. THE GAS-SURFACE INTERACTION 23
2.2.1 The outline of the density functional theory
DFT is formulated by P. Hohenberg and W. Kohn in 1964 [2], by proving that:
(1) The external potential and hence the ground state energy of the system can be uniquely determined by the ground state electron density, or this energy is a functional of the electron density. One immediate consequence of this theorem is that the familiar relation of the ground state density can be reversed i.e. ψ0 = ψ0[n0], where ψ0 and n0 are the ground state electronic wavefunction and density, respectively. (2) The electron density variationally minimizes the total energy of the system, in other words, the determination of the exact ground state is possible.
Thus an N-electron Schrödinger wave equation can be written within the Born- Oppenheimer approximation as,
Hψˆ 0 =hT + ˆˆ V + ˆUiψ0 =
N
X
i
−~2 2m∇2i +
N
X
i
V (ri) +X
i<j
U (ri, rj)
ψ0 = Eψ0. (2.7) Here, ˆH is the electronic molecular Hamiltonian, N is the number of electrons and ˆU is the electron-electron interaction. The operators ˆT and ˆU are ‘universal operators’ because they are independent of the system, while ˆV is generally an external (Coulombic) potential which is system dependent. Following Hohenberg and Kohn, we can represent any ground state property, O as
hOi[n0] = hψ0[n0] | ˆO | ψ0[n0]i, (2.8) where n0 is the ground state electron density. The ground state energy can be represented as,
E[n0] =Dψ0[n0] T + ˆˆ V + ˆU ψ0[n0]E. (2.9) Since V is the only functional which is system dependent, we can write a general expression of E[n],
E[n] = T [n] + U [n] + ˆ
V (r)n(r)d3r = FHK[n] + ˆ
V (r)n(r)d3r (2.10) where FHK is an universal functional.
This expression can be minimized to get n0 and hence all the ground state properties. More interestingly the above expression has an integral over three dimensions, so that the integration schemes became more simpler and cheaper compared to that of many-body wavefunction methods (where the equivalent in- tegrals are tedious due to the 3N dimensionality).
2.2.2 Kohn-Sham formulation
W. Kohn and L. J. Sham proposed a self-consistent minimization scheme to min- imize E[n], thus making DFT calculations possible for real systems [3]. It start with non-interacting system of electrons with an external potential. Furthermore in this system U = 0, and T = Ts[n], where Ts[n] is the non-interacting kinetic energy. Then Eq. 2.10 can be rewritten as,
Es[n] = ˆ
Vs(r)n(r)d3r + Ts[n] (2.11) and corresponding exchange and correlation (xc) energy is defined by
Exc = FHK[n] + Ts−
ˆ ρ(r1)ρ(r2)
|r1− r2| dr1dr2. (2.12) The xc potential can now be written as
Vxc(r) = δExc[n]
δn[r] . (2.13)
From the above three equations one can set up an expression for the non-interacting external potential or ‘effective’ potential as,
Vs(r1) = V (r1) +
ˆ ρ(r2)
|r1− r2|dr2 + Vxc(r1). (2.14) Eq. 2.14 forms N one particle Kohn-Sham (KS) set of coupled differential equations [4]:
−∇2
2 + Vs(r)
!
ψi(r) = εiψi(r). (2.15) These equations can be solved iteratively yielding eigenvalues or KS orbital energies {εi} and eigenfunctions or one electron KS orbital functions {ψi}. The density is calculated by summing up all the occupied KS orbitals,
n0 =
N
X
i=1
|ψi|2 (2.16)
An important aspect of the KS formulation is that it implicitly includes the treatment of the xc part. For many systems this makes a DFT calculation readily comparable with -say, an MP2 level of theory with a computational cost of HF method [4]. The choice of a suitable functional form is a very important step in any DFT calculation. In this work GGA (generalized gradient approximation) type functionals are used viz. PW91 [5, 6] and RPBE [7]. GGA approximates the xc energy via
Exc= ˆ
xc(n, ∇n)n(r)d3r, (2.17) where xc is the exchange-correlation energy density.
2.2.3 Plane-wave DFT.
A planewave basis set is a natural choice for the basis set to solve the Kohn- Sham equations for spatially extended systems like metals, polymers etc. Using the translational symmetry of the atomic arrangements, the form of the effective potential for an extended system can be written as,
2.2. THE GAS-SURFACE INTERACTION 25
Vef f(r + T) = Vef f(r), (2.18) where r and T are the arbitrary and (direct-) lattice vectors, respectively. By virtue of the Bloch’s theorem, eigenfunction of such system can be written as a product of a plane wave and a periodic function viz.,
ψ(r + T) = eikTψ(r). (2.19)
In the above equation, k is a wave vector and it is restricted to be within the first Brillouin zone (BZ). Eq. 2.19 can be stated in an alternative form such that all eigenfunctions of single-particle Schrödinger equation with a periodic potential can be written in terms of a periodic function ukj modulated by a plane wave with wave vector k. Equivalently, this can be expressed for each and every band index, j such that
ψkj(r) = eikrukj(r). (2.20) By assuming this eigenfunctions to be normalized with respect to the unit cell, the periodic function ukj can be expanded in a Fourier series:
ukj(r) =X
G
ckjGeiGr. (2.21)
Here, G is a lattice vector defined in the reciprocal space under the condition, G.T = 2πn, where n is any integer. Now the Kohn-Sham equation (Eq. 2.15) can be re-written as
−∇2
2 + Vef f(r)
!
ψkj(r) = εkjψkj(r), (2.22) where εkj is the Kohn-Sham eigenvalue of the band j.
Vef f(r) can be easily decomposed into the following functionals of the density as,
Vef f(r) = Vext(r) + VH[n(r)] + Vxc[n(r)],
in which Vext(r), VH[n(r)] and Vxc[n(r)] are the external potential of the nuclei, the Hartree and the exchange-correlation potentials, respectively. The n(r) is the ground state electron density defined as,
n(r) = 2 Ωc (2π)3
X
j
ˆ
BZ
|ψkj(r)|2Θ(EF − kj)d3k, (2.23) where, Ωcis the volume of a unit cell, EF is the Fermi energy and Θ is a Heaviside step function (which is one for positive and zero for negative arguments) and are excluded the eigenvectors corresponding to empty states from the sum.
For practical purpose, the above equation is approximated - i.e. the integrals are discretized and are written as,
2 Ωc (2π)3
X
j
ˆ
BZ
...Θ(EF − kj)d3k → 1 Nkpt
X
k
fkj... , (2.24)
where fkj are the occupation number taking values of either zero or one and Nkpt is defined as a regular mesh consisting of k-points. There are several schemes proposed to construct such meshes in the literature [8, 9]. The major advantage of this approximation is that only a finite number of k-points are needed to evaluate the charge density and hence the total energy of the system. Another consequence of this approach is that the larger the super cells (which defines the boundary condition), the fewer k-points are needed to perform the calculation. This can be interpreted, because the Brillouin cell becomes smaller in the reciprocal space when a larger supercell is used.
By virtue of the plane wave basis set, a Fourier representation of the Kohn- Sham equations (Eq. 2.15) is possible. By inserting Eq. 2.21 into Eq. 2.15 and multiplying the left-hand-side with e−i(k+G0)r and integrating over r, one gets a matrix eigenvalue equation:
X
G
(~2
2m) k k + G k2 δG0G+ Vef f(G0−G)ckjG = kjckjG0. (2.25) In practice the Fourier expansion (Eq. 2.21) of the wave function is truncated using a kinetic energy cut-off defined in terms of plane wave vectors (k + G). This cut-off energy, Epw is defined as below,
(~2
2m) k k + G k25 Epw. (2.26) Epw is one of the most important parameters and hence it should be converged reasonably well.
The next major approximation used in plane wave DFT is the use of pseu- dopotentials. Although core electrons do not participate in chemical bonding as they are strongly localized around the nucleus, the rigorous representation of these electrons in the PW scenario is a difficult task. For example, the oscillatory wave function of the core electron (see ΦAE in the Fig. 2.1) requires a large number of plane waves and hence its calculations becomes practically unfeasible. This situa- tion can be avoided by neglecting the orthogonal relationship of the core-valence electrons using a non-local potential, called a pseudopotential. The solutions of the resulting approximation (pseudo-orbitals) are very close to the valence orbitals in the outer regions of the systems but lack many oscillations (i.e. nodeless valence wavefunctions) in the core regions required for the core-valence orthogonality.
The main parameter to decide the accuracy of the pseudopotential approach is the cut-off radius (Rcut) where the Columb potential from the pseudopoten- tials and bare Coulomb potential of the system has the same value. Thus the pseudo orbital will match the corresponding all electron orbital outside this point as indicated in Fig. 2.1. By reducing Rcut one can make a more accurate pseu- dopotential (‘hard’ pseudopotential) but at the same time the DFT calculations becomes computationally more expensive. This technical difficulty can be handled by using Vanderbilt ultrasoft pseudo-potentials [10, 11]. In this approach, fewer plane waves are required in expansions of the electron valence states by relaxing the norm-conservation (of the valence pseudo orbitals) criteria.
In this work we have used two plane wave codes: DACAPO [12] and VASP [13, 14, 15, 16]. The DACAPO is an orthogonalized plane wave (OPW) code
2.3. REACTION PATH CALCULATION 27
Figure 2.1: A Schematic diagram of the PSP approach. The all electron wavefunc- tion, ΦAEis approximated with the pseudo wavefunction ΦPSand the correspoding all electron potential Z/r is replaced by the pseudo potential, VPS.
and it uses a plane wave basis for the valence electronic states and describes the core-electron interactions with Vanderbilt ultrasoft pseudo-potentials. The VASP code also uses plane wave basis set but the interaction between ions and electrons can be described by the projector-augmented wave (PAW) method [17, 18]. The PAW approach is superior to OPW methods [19] since it preserves the local node structure of orbitals in its scheme, in other words, it treats the inner- electronic wavefunction.
2.3 Reaction path calculation
The reaction path (RP) can be a miniature representation of the whole PES, by connecting initial, final and transition states with a suitable reaction coordinate.
Though there is no universal method to locate the TS or to obtain a reaction path, there are a variety of useful alternative methods [20]. Here an outline of the nudged elastic band (NEB) method and optimization schemes are presented.
2.3.1 The Nudged elastic band methods
The nudged elastic band (NEB) is one of the ‘iterative’ methods to find the reaction path and/or the transition state [21, 22, 23]. The NEB defines a minimum energy path (MEP) by connecting, generally, two minima on the PES.
In the NEB method, as a first approximation, an initial elastic band is made by a linear interpolation between the initial and final configurations. This band consists of a number of geometries, called ‘images’. Then the force acting upon an image i is,
FiN EB = −∇V (Ri) + Fis||, (2.27)
where
Fis|| = ki+1(Ri+1− Ri) − ki(Ri− Ri−1). (2.28) In the above equations, Ri is the Cartesian geometry of an image i and ki
is a ‘spring constant’ between Ri − Ri−1 and it defines the stiffness of the har- monic spring (or an ‘elastic band’) connecting these images. Since FiN EB has two independent terms it can be written as,
FiN EB = Fi⊥+ Fis||, (2.29) where Fi⊥ and Fis|| (See Fig. 2.2) are the perpendicular and parallel force compo- nents to the band, respectively. By defining a tangent vector of an image ˆτi along the MEP, the above components can be written as,
Fi⊥ = −∇V (Ri) + (∇V (Ri). ˆτi)ˆτi (2.30) and
Fis|| = (ki+1(Ri+1− Ri) − ki(Ri − Ri−1)) ˆτi. (2.31)
Figure 2.2: The initial (black dots) and final (circles) images of a NEB calculation.
The different force components of an image i are shown in the inset.
Since an image along the MEP satisfies a necessary condition, i.e., Fi⊥ =
−∇V (Ri) + (∇V (Ri).ˆτi) ˆτi = 0, a minimization algorithm can be applied to move the initial images in the direction defined by Fi⊥. This process can be iterated until all the images in the calculation lie on the MEP. The convergence rate of a NEB calculation is mainly dependent on the closeness between the true MEP and its initial estimate and the number of images.
Several modifications of the NEB have been developed, such as climbing-image NEB [24] and Adaptive-NEB [25], to give an accurate estimate of the saddle point.
2.3. REACTION PATH CALCULATION 29
2.3.2 Optimization methods
The methods described above use the optimization or minimization of the forces with or without constraints. Thus a short review of these methods is presented here. Fig. 2.3 shows an illustration of the common optimization methods applied to the two dimensional Rosenbrock function (f (x, y) = (1 − x)2+ 100(y − x2)2), which is widely used as a test function for optimizing problems [26, 27].
• Steepest descent:
The steepest descent (SD) method is the simplest and the slowest method among the optimization methods discussed. The gradient - or the force components - are used to find the directions of the displacement vector, which can be represented as,
xn+1= xn− hi∇V (xn), n ≥ 1.
where xnis the sequence of points which converged to a local minimum, hi is the step size, which is usually positive and an order of magnitude less than unity, and ∇V (xn) is the gradient at the point xn.
• Quasi-Newton:
The Quasi-Newton (QN) method is an improved version of Newton’s iter- ative method, which requires a second derivative matrix (Hessian). This method does not require the Hessian matrix explicitly, but an approximated inverse of the Hessian matrix, Q, and the derivative matrix. Several imple- mentations are available for QN to update the Q matrix, for example the BFGS (Broyden, Fletcher, Goldfarb, Shanno) scheme. This scheme may be written as,
1. dn= −Qn∇V (xn) 2. αn=−∇V (xdT n)T·dn
kQndn (Wolfe’s condition) 3. xn+1= xn+ αndn, n ≥ 1.
where dn is the direction vector and αn is the variable step length which is usually obtained by a line search procedure, by minimizing ∇V (xn)T · dn.
• Conjugate gradient:
In the conjugate gradient method, the direction of descent is obtained by using the previous ‘conjugate’ (or linearly independent) directions and the present SD direction. The attractive feature is that it does not require any knowledge of the Hessian matrix. The general description of the iterations is as follows:
1. β = [∇V (x[∇V (xn)]T∇V (xn)
n−1)]T∇V (xn−1) (Fletcher-Reeves formula) 2. dn= −∇V (xn) + βdn−1
(a) (b) (c)
Figure 2.3: Contour plots of Rosenbrock function of two variables with different minimizers. The paths are started at [-1,1] and terminated at its global minimum [1,1]. The steepest descent path (a) converges with 724 steps, the quasi-Newton method (b) with 38 and the conjugate gradient method (c) with 22 steps.
3. xn+1 = xn+ dn, n ≥ 1.
2.4 Vibrational normal modes
The vibrational energy of a molecule is one of the three main quantized internal energies, apart from the electronic and rotational energies. Many physically inter- esting and important properties like zero point energy are related to the vibrational motion of the molecule. In this section a small description of vibrational normal modes is provided.
A molecule can be considered as a N-coupled mass points and thus viewed as a coupled simple harmonic oscillator system - a dynamic system vibrating according to Hooke’s law. This N-atom oscillator can be then described by the Lagrangian, L as
L = T − V, (2.32)
where T and V are the kinetic and potential energy terms. To solve this, par- ticularly in a ‘normal mode’ approach, one can define a set of mass-weighted normal mode coordinates, qk, k = {1, 2, 3..., 3N } and qi = √
mj4xj, qi+1 =
√mj4yj, qi+2 = √
mj4zj..., where j labels the atoms and i the 3N Cartesian coordinates. Using these coordinates, Eq. 2.32 can easily be uncoupled, i.e., avoiding the cross terms and leaving T and V as diagonal matrices [28]. For this the above equation is first written as,
L = 1 2
3N
X
i=1
q.2i − 1 2
3N
X
i,j=1
∂V
∂qi∂qj
0
qiqj, (2.33)
2.4. VIBRATIONAL NORMAL MODES 31
Figure 2.4: The nine vibrational normal modes of methane are plotted with the degenracies and displayed in the descending order in the vibrational energy.
or in vector-matrix form as,
L = 1 2
q.t .q − 1
2qtF q. (2.34)
Here, the superscript, t indicates the transpose and F is the second derivative or the Hessian matrix. The above equation determines the dynamics of the vibrating system. The F can be computed by taking numerical derivatives instead of the analytical gradients. Since F is a symmetric matrix, an orthogonal or similarity transform is exists to diagonalize F :
WtF W = Ω, (2.35)
where Ω is the diagonal matrix, W is a matrix of order 3N, and WtW = I, where I is an identity matrix. Eq. 2.35 is nothing but an eigenvalue equation, where each element of Ω is an eigenvalue (harmonic frequency) associated with a specific column of W , the normal mode.
There are 9 vibrational normal modes present in theCH4 system. From the point group analysis of the Td character table of the methane one can see that the
symmetry species of the normal modes of methane:
ΓCH4 = A1+ E + 2T2. (2.36)
Here, A1, E and T2 are the irreducible representations in which the T2 modes are infrared and Raman active while the A1 and E modes are only Raman active.
The mode E is doubly degenerate and the T2 modes are triply degenerate. Using spectroscopic notations, we can represent these modes as ν1, ν2, ν3 (ν4) for A1, E and T2, respectively. The ν1 mode describes the symmetric stretch mode, ν2 represents symmetric bend modes, ν3 corresponds to the asymmetric stretches and ν4 describes the asymmetric bending modes (see Fig. 2.4). The below table lists the normal mode frequencies of these vibrations [29].
CH4 normal modes Frequencies (cm−1)
ν1 2916
ν2 1533
ν3 3019
ν4 1311
2.5 Growing the potential energy surface
The reaction dynamics of a system can be elucidated if a PES of the system is known, by performing classical and quantum scattering calculations. Obtaining and fitting the PES, especially for a high dimensional system CH4/Ni(111) is a computationally demanding task. In this work, we employed the GROW code [30, 31, 32] which is based on the modified Shepard (MS) interpolation method to fit the DFT potential energy surface. The computational cost is reduced by ap- plying the following approximations (1) the surface atoms are taken as frozen (the static-surface model), and (2) the use of the Born-Oppenheimer approximation (neglecting electron-hole pair excitations).
2.5.1 Electronic structure calculations
The DFT calculations are carried out using the DACAPO code [12]. The gen- eralized gradient approximation is used to describe the exchange-correlation en- ergy of the electrons. We used the RPBE functional [7] for the GGA and it was constructed from the PBE functional to satisfy a constraint upon approximate exchange-correlation functionals - the Leib-Oxford bound[33]. Many DFT studies in various gas-surface reactions show that the functional gives better chemisorp- tion energies than the PBE functional [7]. The ion cores are approximated using ultrasoft pseudo-potentials and a plane-wave basis set is used for the electronic orbitals.
2.5. GROWING THE POTENTIAL ENERGY SURFACE 33
2.5.2 Modified Shepard interpolation method
We used the modified Shepard (MS) interpolation method extended to gas-surface reactions by Crespos et al. to grow the PES [34]. In this method, interpolated data are calculated by a weighted series of Taylor expansions centered on DFT data points, sampled throughout the configuration space of the system. To monitor the growth of the PES, the convergence of the reaction probability of interest is calculated. The main feature of the MS method is that the new points are added only to the dynamically relevant region of configuration space. This thus avoids the use of regular grids and DFT optimizations.
To define the configuration space, a symmetry unique triangle is defined (see Fig 2.5). To construct this triangle, the approximation was made that the in- teraction of CH4 with the fcc site is equal to the interaction of CH4 with the hcp site. Inverse interatomic distances, Qi = 1/Ri, are used as coordinates in this method and the vector defining any configuration of the system is given by Q = 1/R1, . . . , 1/Rn(n−1)/2
, where n is the number of atoms used to describe the system. These inverse interatomic distances are used to get a physically more reasonable behavior when two atoms are close to each other. For each geometry Q, a set of 3n − 6 algebraically independent linear combinations of the n (n − 1) /2 interatomic distances, ζ (Q), can be defined [31] in terms of the inverse distance as,
ζm =
n(n−1)/2
X
k=1
UmkQk (m = 1, . . . , 3n − 6) , (2.37) where Umk is the transformation matrix from reciprocal bond lengths to the new coordinates. Now, the potential energy at a configuration Q can be expanded as a second-order Taylor expansion Ti(Q)
Ti(Q) = V [Q (i)] +P3n−6k=1 [ζk− ζk(i)] ∂ζ∂V
k
Q=Q(i)
+12P3n−6k=1 P3n−6j=1 [ζk− ζk(i)] [ζj − ζj(i)] ∂ζ∂2V
k∂ζj
Q=Q(i)
(2.38)
where V [Q (i)] is the potential value at Q (i), and the gradients at this point are computed analytically by DACAPO. The second derivatives are computed from the gradients using forward differencing using a lower level of theory, i.e. by applying a 2 layer Ni slab and a lower Epw and a k-point grid containing less points than the values used for the potential.
The total potential energy at any configuration Q is then taken as
V (Q) = X
g∈G Ndata
X
i=1
wg◦i(Q) Tg◦i(Q) (2.39) where Tg◦i(Q) represents a second-order Taylor expansion and wg◦i(Q) a normal- ized weight function (see Refs. [31] and [34] for more details). Ndata is the number of DFT data configurations in the interpolation, G is the symmetry group, and g ◦ i denotes the transformation of the i’th data point by the group element g. To take into account the full symmetry of the system, a sum is taken over both the DFT data points and their symmetry equivalent points.
Figure 2.5: The shaded triangle represents the symmetry unique triangle of the unit cell. The large circle indicates the Ni atoms and the dark and open circle represents fcc and hcp sites respectively.
As was said before, the MS interpolation scheme adds DFT points not uni- formly in the PES but only in the regions where the dynamics is apparently im- portant. This scheme is in contrast with other interpolation methods [35, 36]. The classical trajectory method is used to find these regions by spotting the presence of the molecules during dynamics calculations. The new points are added according to one of the following criteria [30, 31, 32]: (1) The h-weight criterion, which is based on the assumption that the best location for a new point would be in the region most frequently visited by the trajectories, (2) the variance criterium, which is based on the assumption that the accuracy of the PES will be best improved if the new point is added in the region where the determination of the energy is suspected to be the most inaccurate.
2.5.3 Implementation of the MS method in GROW
An iterative procedure is used to grow the PES [30, 31, 32]. A short outline is given below.
1. An initial PES is constructed mainly from the NEB calculations, where the data set consists of optimized initial and final geometries and the transi- tions state. In addition to this we used a sparse set of optimized DFT CH4 geometries distributed over the unit cell for the GROW calculation.
2. Using this initial PES a few classical trajectories (typically ten) are run. We have grown the PES using simultaneously several translational energies, to assure the accurate representation of the PES for the whole energy range of interest. Along each trajectory the molecular geometries explored by the simulation are periodically stored.
3. From the stored geometries a new point is chosen and added to the PES according to one of the two criteria described in the above section.
4. The continuation of the process: the program loops back to step 2
2.6. THE CLASSICAL TRAJECTORY METHOD 35 5. Convergence test: After about 100 points are added, the accuracy of the PES is checked by computing an accurate value of the desired observable (the reaction probability) from a more extensive classical trajectory simulation ( approximately, 10 000 trajectories). The growth process is stopped when the observable is considered to be converged, within a given tolerance.
2.6 The classical trajectory method
Classical mechanics based on Newton’s law of motion may provide an approximate picture of the collision dynamics of a chemical reaction. The numerical solution of Newton’s equations of motion yields the position of each and every atom as a function of time, in what are called trajectories. Statistical mechanics provides a link between the macroscopic observables and the microscopic world obtained from the classical trajectory (CT) simulations via the microcanonical ensemble and its ergodic property. These relations can be written for an observable, a, as
hai =
´ dx a(x)δ( ˆH(x) − E)
´ dx δ( ˆH(x) − E) = lim
τ →∞
ˆτ
0
dt a(xt) ≡ ˜a. (2.40)
The above equation describes that the expectation value, a, of a microcanonical ensemble is equal to the time average over the classical trajectories.
Technically, several algorithms exist to perform the numerical integration and in this work, a velocity-Verlet algorithm has been used [37]. This algorithm is a simple, robust and reversible method and allows large time steps. The algorithm can be represented in terms of the variables viz., position(q), force (f ), and velocity (v) as,
qn+1 = qn+ vnδt + f (q2mn)δt2 vn+1 = vn+f (qn+12m)+f (qn)δt
(2.41)
A main initial procedure in CT method is to choose a set of initial conditions where a Monte Carlo sampling method is used to simulate the molecular initial conditions for each set of initial parameters (such as translational, rotational energies etc).
In the quasi classical trajectory (QCT) method the initial zero point energy (ZPE) of the molecule is included in the dynamics [38]. This method allows us to include quantum effects in the classical simulation in an approximate way. The QCT is susceptible to the so-called ZPE violation problem but for an activated sys- tem this problem is expected to play a minor role at energies above the threshold to reaction [39]. The QCT method generally gives accurate results for activated molecule-surface reactions [40, 41, 42, 43, 44]. This is mainly achieved by treating the ‘vibrational softening’ (adiabatic transfer of energy from internal to transla- tional motion which takes place when a molecule approaches the surface) in the QCT method [45].
2.7 Quantum Dynamics and MCTDH method
2.7.1 Quantum dynamics: general concepts
To extract the details of the reaction dynamics at the quantum level, one needs to solve the Schrödinger equation of motion of the atomic nuclei:
i~∂Ψ(R, t)
∂t = ˆH(R, t) = ( ˆTn+ E(R))Ψ(R, t). (2.42) Here, ˆTnis the KE operator of the nuclei, E(R) is the electronic energy (i.e. Poten- tial energy from the PES), and Ψ(R, t) represents a time dependent wavefunction of all the nuclear coordinates, R.
Since the time dependent Schrödinger equation (TDSE) is a linear-hyperbolic differential equation it is an initial value problem. In other words, to calculate the time evolution of the system only the initial state specification is needed. If the Ψ(R, t0) is known, then the general solution to the above Eq. 2.42 is given by,
Ψ(R, t) = ˆU (t − t0)Ψ(R, t0), (2.43) where ˆU (t − t0) is the time evolution operator, which is defined as:
U (t, tˆ 0) = e−(i/~) ˆH(t−t0). (2.44) The basic operation (propagation) to solve the TDSE is to apply the Hamiltonian operator to the wavefunction (Ψ(R, t0)) recursively, to generate a new wavefunction [46]. From Eq. 2.43 and 2.44, it can be readily seen that,
Ψ(R, t) = ˆU (t − t0)Ψ(R, t0) = {1 − i ˆH(t − t0)
~
!
− 1 2!
H(t − tˆ 0)
~
!2
+i 3!
H(t − tˆ 0)
~
!3
+ ...}Ψ(R, t0). (2.45)
In general, the initial state, Ψ(R, t0) is constructed from the wave packet and its evolution (the translational motion) is then interrogated to compute the desired information, such as reaction probability.
2.7.2 Representation of the wave function
An initial step to solve the TDSE is to expand the time dependent wavefunction in a proper basis set. This is normally carried out by expanding Ψ(R, t0) in a suitable basis set and then by evaluating the operator action on basis functions. One can use the collocation method, the discrete variable representation (DVR) [47, 48, 49, 50], the finite basis representation (FBR) [47, 48, 49, 50], FFT techniques, and a variational basis (or close coupling) representation. The computational efficiency of any time dependent wave packet (TDWP) method is strongly dependent on the representation of the basis sets.
2.7. QUANTUM DYNAMICS AND MCTDH METHOD 37 A. Collocation method
Collocation refers to the numerical solution of a set of equations such that the result is exact at a discrete set of points [46]. Let Ψ(r) be a given wave function which can be represented by an expansion of N linearly independent analytic basis function (fn(r)) obeying the appropriate boundary conditions:
Ψ(r) =
N
X
n=1
Cnfn(r). (2.46)
Here, the coefficients, Cn are found by relating the wave functions at N points in r to the expansion and solving N linear equations by a matrix inversion, viz.,
Ψ(ri) =
N
X
n=1
Cnfn(ri) ≈ Ψ = fC. (2.47) In Eq. 2.47, f is defined as an N dimensional square matrix which transforms the coefficients of the vector C into the wave function vector,Ψ (an approximation to Ψ). The inverse matrix, f−1 can be used to obtain the vector C from the vector Ψ by C = f−1Ψ. This demonstrates that the action of a non-local operator, ˆO(r) on Ψ can be computed in a set of finite basis of fn , using:
O(r)Ψ(r) ≈ˆ
N
X
n=1
CnO(r)fˆ n(r) ≈
N
X
n=1
Cnfn0(r). (2.48)
B. FBR and DVR methods
The collocation method can be greatly simplified if the linearly independent basis functions fn(r) of Eq. 2.46, are chosen to be eigen functions of the operator, ˆO(r) and are taken to be mutually orthogonal over the N discrete points at ri. In other words, if we can write a relation of basis functions such as,
N
X
k=1
fn∗(rk)fm(rk) = δn,m, (2.49) the matrix f in Eq. 2.47 does not have to be inverted anymore and the following identity holds for Eq. 2.47:
N
X
k=1
fm∗(rk)Ψ(rk) =
N
X
n=1 N
X
k=1
Cnfn(rk)fm∗(rk) =
N
X
n=1
Cnδn,m = Cm. (2.50)
Since fn are the eigenfunctions of ˆO Eq. 2.48 reduces to
O(r)Ψ(r) ≈ˆ
N
X
n=1
CnO(r)fˆ n(r) ≈
N
X
n=1
Cnonfn(r) ≈
N
X
n=1
Cn0fn(r). (2.51)
Here the operator ˆO(r) is considered to be diagonal in the C representation.
The action of ˆO on Ψ can now be evaluated in three steps:
1. Transforming the DVR of Ψ (Ψ) into the finite basis representation (FBR) of Ψ i.e. C, by using f−1= f† because of Eq. 2.49.
2. Applying the operator ˆO(r) on the FBR of Ψ(r) to get C’ - the FBR of O(r)Ψ(r). At this step the operator ˆˆ O(r) can be represented by a diagonal matrix viz., Omn= δmnon by virtue of the FBR.
3. Transforming C’, the FBR of ˆO(r)Ψ(r), back to the DVR of ˆO(r)Ψ(r) i.e.
Ψ0, by using f.
The DVR and FBR are considered to be conjugate to each other since the matrices f anf f−1 are used to transform between these representation of the wave function.
Additionally, in the vector-matrix notation a complete procedure of applying ˆO(r) on Ψ(r) can now be expressed as:
Ψ0 = fOf−1Ψ (2.52)
This method of applying an operator to a grid based wave function is called the pseudo-spectral (PS) method. The DVRs are generally the time-independent bases and widely used in many wavepacket propagation calculations. By using a variety of basis functions, a number of different DVRs has been developed: the harmonic oscillator DVR used for vibrational motion, Legendre DVR for rotations, and exponential and sine DVRs used for free motion with or without periodic boundary conditions.
C. The Fourier representation
There are three properties of the plane wave functions that make the representation of a wave function much simpler. Firstly, they are eigenfunctions of the kinetic energy operator (i.e. ∇2reik.r = −k2eik.r), secondly, they are mutually orthogonal to a series of N equally spaced grids rn= nL/N over the range of [0,L]:
N −1
X
n=0
e−2πip.rnL e2πiq.rnL = N δp,q. (2.53)
and thirdly, they form a complete set for periodic or bandwidth limited functions of r [51]. Note that Eq. 2.53 is exact and is resembles Eq. 2.49. Using Eq. 2.53 the discrete space conjugate to rnis given by km= 2πmL where 0 ≤ m ≤ N − 1. In this representation the matrix elements of the matrix f in Eq. 2.47 , fmn has a simple form: fmn = e2πimn/N. It immediatly follows that the FBR or the coefficients of the discretized wave function is given by C = f−1Ψ. These transforms based on plane wave basis are called Fourier transforms and in its advanced versions (FFT) the problem scales approximately as N log2N [52].
2.7.3 The MCTDH method
In general, the numerical effort of conventional wave packet propagation schemes increases exponentially with the number of degrees of freedom. The multi configu- ration time-dependent Hartree (MCTDH) method is an in-principle exact method and its ability lies in the fact that it uses a variational principle to derive equations
2.7. QUANTUM DYNAMICS AND MCTDH METHOD 39 of motion [53, 54, 55]. This method uses optimized time-dependent one-particle basis functions and thus arrives at a compact representation of the state vector at the expense of more complicated equations of motion. The single-particle functions in MCTDH method are usually represented by standard DVR or FFT-schemes.
A. Wavefunction Ansatz and Equations of Motion
The solution of the TDSE (Eq. 2.42) is supposed to be distinct and and well- resolved in position. This immediately leads to a better approximation for the time dependent wavefunction: using linear combination of the single particle functions (SPFs). The MCTDH method uses the following wavefunction ansatz to solve the TDSE for a system with f degrees of freedom (DOFs) described by nuclear coordinates Q1, ..., Qf :
ψ (Q1, . . . Qf, t) =
n1
X
j1=1
. . .
nf
X
jf=1
Aj1...jf (t)
f
Y
κ=1
ϕ(κ)jκ (Qκ, t) . (2.54) Here, Aj1, . . . , Ajf are the MCTDH expansion coefficients and ϕ (Qκ, t) are the nκ single particle functions associated with the degree of freedom κ. The Ansatz is similar to the standard wavepacket expansion, except that the SPFs provide a time- dependent basis set. In general, the SPFs are represented on a time-independent primitive basis set with Nκ points for the κ’th degrees of freedom, viz.,
ϕ(κ)jκ (Qκ, t) =
Nκ
X
iκ=1
c(κ)iκ,jκ(t) χ(κ)iκ (Qκ) . (2.55) A DVR or Fourier representation is usually chosen to represent the basis func- tion, χ(κ)iκ (Qκ). Using a technique called mode-combination, SPFs can be used to describe the motion in more than one degree of freedom [56].
To formulate the working equations of MCTDH one needs to introduce the single-hole functions. The single-hole functions Ψ(κ)l are defined as the linear com- bination of Hartree products of (f − 1) SPFs that do not contain a SPF for the coordinate Qκ
Ψ(κ)l =
κ
X
J
AJκ
lϕ(1)j1 . . . ϕ(κ−1)jκ−1 ϕ(κ+1)jκ+1 . . . ϕ(f )jf . (2.56) Here Jlκ denotes a composite index J with the κ’th entry set at l, and PκJ is the sum over the indices for all degrees of freedom excluding the κ’th. By defining the composite index J and the configurations ΦJ : AJ = Aj1...jf and ΦJ = Qfκ=1ϕ(κ)jκ , the equations of motion can be derived using this ansatz in the Dirac–Frenkel variational method. The resulting coupled differential equations can be then written as:
i ˙AJ =X
J
h ΦJ| H| ΦLi AL (2.57)
i ˙ϕ(κ)j =X
lm
1 − P(κ) ρ(κ)−1
jl hHi(κ)lm ϕ(κ)m , (2.58)
