Baixe introduction to Ab Initio Molecular Dynamics e outras Notas de estudo em PDF para Engenharia de Produção, somente na Docsity! Lecture Notes Introduction to Ab Initio Molecular Dynamics Jürg Hutter Physical Chemistry Institute University of Zurich Winterthurerstrasse 190 8057 Zurich, Switzerland hutter@pci.unizh.ch September 14, 2002 Contents 1 Molecular Dynamics 4 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Microcanonical Ensemble . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Numerical Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.5 Extended System Approach . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.5.1 Barostats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.5.2 Thermostats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2 Ab initio Molecular Dynamics 12 2.1 Born–Oppenheimer Molecular Dynamics . . . . . . . . . . . . . . . . . . . 14 2.1.1 Forces in BOMD . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Car–Parrinello Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . 16 2.2.1 How to Control Adiabaticity ? . . . . . . . . . . . . . . . . . . . . . 18 2.2.2 Forces in CPMD . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.3 Velocity Verlet Equations for CPMD . . . . . . . . . . . . . . . . . 19 2.3 Comparing BOMD and CPMD . . . . . . . . . . . . . . . . . . . . . . . . 20 3 Plane Waves 24 3.1 Unit Cell and Plane Wave Basis . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Kinetic Energy and Local Potentials . . . . . . . . . . . . . . . . . . . . . 26 3.3 Electrostatic Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.4 Exchange and Correlation Energy . . . . . . . . . . . . . . . . . . . . . . . 29 3.5 Car–Parrinello Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4 Pseudopotentials 33 4.1 Why Pseudopotentials ? . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2 Norm–Conserving Pseudopotentials . . . . . . . . . . . . . . . . . . . . . . 35 4.2.1 Hamann–Schlüter–Chiang Conditions . . . . . . . . . . . . . . . . . 35 4.2.2 Bachelet-Hamann-Schlüter (BHS) form . . . . . . . . . . . . . . . . 37 4.2.3 Kerker Pseudopotentials . . . . . . . . . . . . . . . . . . . . . . . . 37 1 Lecture 1 Molecular Dynamics 1.1 Introduction The aim of molecular dynamics is to model the detailed microscopic dynamical behavior of many different types of systems as found in chemistry, physics or biology. The history of molecular dynamics goes back to the mid 1950’s when first computer simulations on simple systems were performed[1]. Excellent modern textbooks [2, 3] on this topic can be found and collections of articles from summer schools are a very good source for in depth information on more specialized aspects [4, 5, 6]. Molecular Dynamics is a technique to investigate equilibrium and transport properties of many–body systems. The nuclear motion of the particles is modeled using the laws of classical mechanics. This is a very good approximation for molecular systems as long as the properties studied are not related to the motion of light atoms (i.e. hydrogen) or vibrations with a frequency ν such that hν > kBT . Extensions to classical molecular dynamics to incorporate quantum effects or full quantum dynamics (see for example Refs. [7, 8] for a starting point) is beyond the scope of this lecture series. 1.2 Equations of Motion We consider a system of N particles moving under the influence of a potential function U [9, 10]. Particles are described by their positions R and momenta P = MV. The union of all positions (or momenta) {R1,R2, . . . ,RN} will be called RN (PN). The potential is assumed to be a function of the positions only; U(RN). The Hamiltonian H of this system is H(RN ,PN) = N∑ I=1 P2I 2MI + U(RN) . (1.1) 4 The forces on the particle are derived from the potential FI(R N) = −∂U(R N) ∂RI . (1.2) The equations of motion are according to Hamilton’s equation ṘI = ∂H ∂PI = PI MI (1.3) ṖI = − ∂H ∂RI = − ∂U ∂RI = FI(R N) , (1.4) from which we get Newton’s second law (using PI = MIṘI) MIR̈I = FI(R N) . (1.5) The equations of motion can also be derived using the Lagrange formalism. The Lagrange function is L(RN , ṘN) = N∑ I=1 1 2 MIṘ 2 I − U(RN) , (1.6) and the associated Euler–Lagrange equation d dt ∂L ∂Ṙi = ∂L ∂Ri (1.7) leads to the same final result. The two formulations are equivalent, but the ab initio molecular dynamics literature almost exclusively uses the Lagrangian techniques. 1.3 Microcanonical Ensemble The following section is slightly shortened from the very clear and concise feature article by Tuckerman and Martyna [10]. The equations of motion are time reversible (invariant to the transformation t→ −t) and the total energy is a constant of motion ∂E ∂t = ∂H(RN , ṘN) ∂t = 0 . (1.8) These properties are important to establish a link between molecular dynamics and sta- tistical mechanics. The latter connects the microscopic details of a system the physical observables such as equilibrium thermodynamic properties, transport coefficients, and spectra. Statistical mechanics is based on the Gibbs’ ensemble concept. That is, many individual microscopic configurations of a very large system lead to the same macroscopic 5 properties, implying that it is not necessary to know the precise detailed motion of every particle in a system in order to predict its properties. It is sufficient to simply average over a large number of identical systems, each in a different configuration; i.e. the macroscopic observables of a system are formulated in term of ensemble averages. Statistical ensembles are usually characterized by fixed values of thermodynamic variables such as energy, E; temperature, T ; pressure, P ; volume, V ; particle number, N ; or chemical potential µ. One fundamental ensemble is called the microcanonical ensemble and is characterized by constant particle number, N ; constant volume, V ; and constant total energy, E, and is denoted the NV E ensemble. Other examples include the canonical or NV T ensemble, the isothermal–isobaric or NPT ensemble, and the grand canonical or µV T ensemble. The thermodynamic variables that characterize an ensemble can be regarded as experimental control parameters that specify the conditions under which an experiment is performed. Now consider a system of N particles occupying a container of volume V and evolving under Hamilton’s equation of motion. The Hamiltonian will be constant and equal to the total energy E of the system. In addition, the number of particles and the volume are assumed to be fixed. Therefore, a dynamical trajectory (i.e. the positions and momenta of all particles over time) will generate a series of classical states having constant N , V , and E, corresponding to a microcanonical ensemble. If the dynamics generates all possible states then an average over this trajectory will yield the same result as an average in a microcanonical ensemble. The energy conservation condition, H(RN , ṘN) = E which imposes a restriction on the classical microscopic states accessible to the system, defines a hypersurface in the phase space called a constant energy surface. A system evolving according to Hamilton’s equation of motion will remain on this surface. The assumption that a system, given an infinite amount of time, will cover the entire constant energy hypersurface is known as the ergodic hypothesis. Thus, under the ergodic hypothesis, averages over a trajectory of a system obeying Hamilton’s equation are equivalent to averages over the microcanonical ensemble. In addition to equilibrium quantities also dynamical properties are defined through ensemble averages. Time correlation functions are important because of their relation to transport coefficients and spectra via linear response theory [11, 12]. The important points are: by integration Hamilton’s equation of motion for a number of particles in a fixed volume, we can create a trajectory; time averages and time correlation functions of the trajectory are directly related to ensemble averages of the microcanonical ensemble. 1.4 Numerical Integration In a computer experiment we will not be able to generate the true trajectory of a sys- tem with a given set of initial positions and velocities. For all potentials U used in real applications only numerical integration techniques can be applied. These techniques are 6 possibility of non–isotropic external stress; the additional fictitious degrees of freedom in the Parrinello–Rahman approach [17] are the lattice vectors of the supercell. These variable–cell approaches make it possible to study dynamically structural phase transitions in solids at finite temperatures. The basic idea to allow for changes in the cell shape consists in constructing an extended Lagrangian where the lattice vectors a1, a2 and a3 of the simulation cell are additional dynamical variables. Using the 3 × 3 matrix h = [a1, a2, a3] (which fully defines the cell with volume det h = Ω) the real–space position RI of a particle in this original cell can be expressed as RI = hSI (1.14) where SI is a scaled coordinate with components Si,u ∈ [0, 1] that defines the position of the ith particle in a unit cube (i.e. Ωunit = 1) which is the scaled cell [17]. The resulting metric tensor G = hth converts distances measured in scaled coordinates to distances as given by the original coordinates. The variable–cell extended Lagrangian can be postulated L = N∑ I 1 2 MI ( ṠtIGṠI ) − U(h,SN) + 1 2 W Tr ḣtḣ− p Ω , (1.15) with additional nine dynamical degrees of freedom that are associated to the lattice vectors of the supercell h. Here, p defines the externally applied hydrostatic pressure, W defines the fictitious mass or inertia parameter that controls the time–scale of the motion of the cell h. In particular, this Lagrangian allows for symmetry–breaking fluctuations – which might be necessary to drive a solid–state phase transformation – to take place spontaneously. The resulting equations of motion read MIS̈I,u = − 3∑ v=1 ∂U(h,SN) ∂RI,v ( ht )−1 vu −MI 3∑ v=1 3∑ s=1 G−1uv ĠvsṠI,s (1.16) W ḧuv = Ω 3∑ s=1 ( Πtotus − p δus ) ( ht )−1 sv , (1.17) where the total internal stress tensor Πtotus = 1 Ω ∑ I MI ( ṠtIGṠI ) us + Πus (1.18) is the sum of the thermal contribution due to nuclear motion at finite temperature and the internal stress tensor Π derived from the interaction potential. A modern formulation of barostats that combines the equation of motion also with thermostats (see next section) was given by Martyna et al. [19]. 9 1.5.2 Thermostats Standard molecular dynamics generates the microcanonical or NV E ensemble where in addition the total momentum is conserved [3]. The temperature is not a control variable and cannot be preselected and fixed. But it is evident that also within molecular dy- namics the possibility to control the average temperature (as obtained from the average kinetic energy of the nuclei and the energy equipartition theorem) is welcome for phys- ical reasons. A deterministic algorithm of achieving temperature control in the spirit of extended system dynamics [18] by a sort of dynamical friction mechanism was devised by Nosé and Hoover [16, 20], see e.g. Refs. [2, 16, 3] for reviews of this technique. Thereby, the canonical or NV T ensemble is generated in the case of ergodic dynamics. It is well–known that the standard Nosé–Hoover thermostat method suffers from non– ergodicity problems for certain classes of Hamiltonians, such as the harmonic oscilla- tor [20]. A closely related technique, the so–called Nosé–Hoover–chain thermostat [21], cures that problem and assures ergodic sampling of phase space even for the patholog- ical harmonic oscillator. This is achieved by thermostatting the original thermostat by another thermostat, which in turn is thermostatted and so on. In addition to restoring ergodicity even with only a few thermostats in the chain, this technique is found to be much more efficient in imposing the desired temperature. The underlying equations of motion read MIR̈I = −∇IEKS −MI ξ̇1ṘI (1.19) Qn1 ξ̈1 = [∑ I MIṘ 2 I − gkBT ] −Qn1 ξ̇1ξ̇2 Qnk ξ̈k = [ Qnk−1ξ̇ 2 k−1 − kBT ] −Qnk ξ̇kξ̇k+1 (1− δkK) where k = 2, . . . , K . By inspection of Eq. (1.19) it becomes intuitively clear how the thermostat works: ξ̇1 can be considered as a dynamical friction coefficient. The resulting ”dissipative dynamics” leads to non–Hamiltonian flow, but the friction term can aquire positive or negative sign according to its equation of motion. This leads to damping or acceleration of the nuclei and thus to cooling or heating if the instantaneous kinetic energy of the nuclei is higher or lower than kBT which is preset. As a result, this extended system dynamics can be shown to produce a canonical ensemble in the subspace of the nuclear coordinates and momenta. In spite of being non–Hamiltonian, Nosé–Hoover (–chain) dynamics is also distinguished by conserving an energy quantity of the extended system; see Eq. (1.21). The desired average physical temperature is given by T and g denotes the number of dynamical degrees of freedom to which the nuclear thermostat chain is coupled (i.e. con- straints imposed on the nuclei have to be subtracted). It is found that this choice requires a very accurate integration of the resulting equations of motion (for instance by using a high–order Suzuki–Yoshida integrator [22]). The integration of these equations of motion is discussed in detail in Ref. [22] using the velocity Verlet algorithm. One of the advan- tages of the velocity Verlet integrator is that it can be easily used together with higher 10 order schemes for the thermostats. Multiple time step techniques can be used and time reversibility of the overall algorithm is preserved. Integrate Thermostats for dt/2 V(:) := V(:) + dt/(2M(:))*F(:) R(:) := R(:) + dt*V(:) Calculate new forces F(:) V(:) := V(:) + dt/(2M(:))*F(:) Integrate Thermostats for dt/2 The choice of the ”mass parameters” assigned to the thermostat degrees of freedom should be made such that the overlap of their power spectra and the ones of the thermostatted subsystems is maximal [22]. The relations Qn1 = gkBT ω2n , Qnk = kBT ω2n , (1.20) assures this if ωn is a typical phonon or vibrational frequency of the nuclear subsystem (say of the order of 2000 to 4000 cm−1). There is a conserved energy quantity in the case of thermostatted molecular dynamics. This constant of motion reads ENVTcons = ∑ I 1 2 MiṘ 2 i + U(R N) + K∑ k=1 1 2 Qnk ξ̇ 2 k + K∑ k=2 kBTξk + gkBTξ1 (1.21) for Nosé–Hoover–chain thermostatted molecular dynamics. 11 A unitary transformation within the space of the occupied orbitals leads to the canonical form HKSφi(r) = iφi(r) (2.11) of the Kohn–Sham equations, with the eigenvalues {i}. This set of equations has to be solved self–consistently in order to yield the density, the orbitals and the Kohn–Sham potential for the electronic ground state. The functional derivative of the Kohn–Sham functional with respect to the orbitals, the Kohn–Sham force acting on the orbitals, can be expressed as δEKS δφ?i (r) = fiH KS e φi(r) . (2.12) Crucial to any application of density functional theory is the approximation of the un- known exchange–correlation functional. Investigations on the performance of different functionals for different type of properties and applications are abundant in the recent literature. A discussion focused on the framework of ab initio molecular dynamics is for instance given in Ref. [27]. Two important classes of functionals are the ”Generalized Gradient Approximation” (GGA) functionals EGGAxc [n] = ∫ dr n(r) εGGAxc (n(r);∇n(r)) , (2.13) where the functional depends only on the density and its gradient at a given point in space, and hybrid functionals, where the GGA type functionals are combined with a fraction of exact exchange energy from Hartree–Fock theory. 2.1 Born–Oppenheimer Molecular Dynamics The interaction energy U(RN) in the molecular dynamics method has the same physical meaning as the Kohn–Sham energy within the Born–Oppenheimer (BO) approximation. The Kohn–Sham energy depends only on the nuclear positions and defines the hypersur- face for the movement of the nuclei. The Lagrangian for BO dynamics is therefore LBO(RN , ṘN) = N∑ I=1 1 2 MIṘ 2 I −min{φi} EKS[{φi};RN ] , (2.14) and the minimization is constraint to orthogonal sets of {φi}. The equations of motions are MIR̈I = −∇I [ min {φi} EKS[{φi};RN ] ] . (2.15) The BOMD program in pseudocode is simply derived from the previous version. 14 V(:) := V(:) + dt/(2M(:))*F(:) R(:) := R(:) + dt*V(:) Optimize Kohn-Sham Orbitals (EKS) Calculate forces F(:) = dEKS/dR(:) V(:) := V(:) + dt/(2M(:))*F(:) Extensions to other ensembles along the ideas outlined in the last lecture are straight forward. In fact, a classical molecular dynamics program can easily be turned into a BOMD program by replacing the energy and force routines by the corresponding routines from a quantum chemistry program. 2.1.1 Forces in BOMD The forces needed in an implementation of BOMD are d dRI [ min {φi} EKS[{φi};RN ] ] . (2.16) They can be calculated from the extended energy functional EKS = EKS + ∑ ij Λij (〈φi | φj〉 − δij) (2.17) to be dEKS dRI = ∂EKS ∂RI + ∑ ij Λij ∂ ∂RI 〈φi | φj〉 + ∑ i ∂EKS ∂〈φi | + ∑ j Λij | φj〉  ∂〈φi | ∂RI . (2.18) The Kohn–Sham orbitals are assumed to be optimized, i.e. the term in brackets is (almost) zero and the forces simplify to FKS(RI) = − ∂EKS ∂RI + ∑ ij Λij ∂ ∂RI 〈φi | φj〉 . (2.19) The accuracy of the forces used in BOMD depends linearly on the accuracy of the min- imization (see Fig. 2.1) of the Kohn–Sham energy. This is an important point we will further investigate when we compare BOMD to the Car–Parrinello method. 15 Figure 2.1: Accuracy of nuclear forces for a system of 8 silicon atoms in a cubic unit cell at 10 Ry cutoff using norm–conserving pseudopotentials 2.2 Car–Parrinello Molecular Dynamics The basic idea of the Car–Parrinello approach can be viewed to exploit the time–scale separation of fast electronic and slow nuclear motion by transforming that into classical– mechanical adiabatic energy–scale separation in the framework of dynamical systems the- ory. In order to achieve this goal the two–component quantum / classical problem is mapped onto a two–component purely classical problem with two separate energy scales at the expense of loosing the explicit time–dependence of the quantum subsystem dynam- ics. This is achieved by considering the extended Kohn–Sham energy functional EKS to be dependent on {Φi} and RN . In classical mechanics the force on the nuclei is obtained from the derivative of a Lagrangian with respect to the nuclear positions. This suggests that a functional derivative with respect to the orbitals, which are interpreted as classical fields, might yield the force on the orbitals, given a suitable Lagrangian. Car and Parrinello postulated the following Lagrangian [28] using EKS. LCP[RN , ṘN , {Φi}, {Φ̇i}] = ∑ I 1 2 MIṘ 2 I + ∑ i 1 2 µ 〈 Φ̇i ∣∣∣Φ̇i〉− EKS [{Φi},RN] .(2.20) The corresponding Newtonian equations of motion are obtained from the associated Euler–Lagrange equations d dt ∂LCP ∂ṘI = ∂LCP ∂RI (2.21) d dt δLCP δ〈Φ̇i | = δLCP δ〈Φi | (2.22) like in classical mechanics, but here for both the nuclear positions and the orbitals. Note that the constraints contained in EKS are holonomic [9]. Following this route of ideas, Car–Parrinello equations of motion are found to be of the form MIR̈I(t) = − ∂EKS ∂RI + ∑ ij Λij ∂ ∂RI 〈Φi | Φj〉 (2.23) µΦ̈i(t) = − δEKS δ〈Φi | + ∑ j Λij | Φj〉 (2.24) where µ is the ”fictitious mass” or inertia parameter assigned to the orbital degrees of freedom; the units of the mass parameter µ are energy times a squared time for reasons of dimensionality. Note that the constraints within EKS lead to ”constraint forces” in 16 and therefore they are only correct up to the accuracy achieved in the wavefunction optimization. In CPMD these are the correct forces and calculated from analytic energy expressions are correct to machine precision. Constraint forces are Fc(Φi) = ∑ j Λij | Φj〉 (2.32) Fc(RI) = ∑ ij Λij ∂ ∂RI 〈Φi | Φj〉 , (2.33) where the second force only appears for basis sets (or metrics) with a nuclear position dependent overlap of wavefunctions. 2.2.3 Velocity Verlet Equations for CPMD The appearance of the constraint terms complicates the velocity Verlet method slightly for CPMD. The Lagrange parameters Λij have to be calculated to be consistent with the discretization method employed. How to include constraints in the velocity Verlet algorithm has been explained by Andersen [15]. In the following we will assume that the overlap is not position dependent, as it is the case for plane wave basis sets. The more general case will be explained when needed in a later lecture. These basic equations and generalizations thereof can be found in a series of papers by Tuckerman et al. [22]. For the case of overlap matrices that are not position dependent the constraint term only appears in the equations for the orbitals. µ | Φ̈i(t)〉 =| ϕi(t)〉+ ∑ j Λij | Φj(t)〉 (2.34) where the definition | ϕi(t)〉 = −fiHKS | Φi(t)〉 (2.35) has been used. The velocity Verlet scheme for the wavefunctions has to incorporate the constraints by using the rattle algorithm. The explicit formulas will be derived in the framework of plane waves in a later lecture. The structure of the algorithm for the wavefunctions is given below CV(:) := CV(:) + dt/(2m)*CF(:) CRP(:) := CR(:) + dt*CV(:) Calculte Lagrange multiplier L CR(:) := CRP(:) + L*CR(:) Calculate forces CF(:) = HKS*CR(:) CV(:) := CV(:) + dt/(2m)*CF(:) Calculte Lagrange multiplier L CV(:) := CV(:) + L*CR(:) 19 Figure 2.2: Conserved energy Econs from Car–Parrinello (CP) and Born–Oppenheimer (BO) molecular dynamics simulations of a model system for various time steps and con- vergence criteria; see text for further details and Table 2.1 for the corresponding timings. Top: solid line: CP, 5 a.u.; open circles: CP, 10 a.u.; filled squares: BO, 10 a.u., 10−6. Middle: open circles: CP, 10 a.u.; filled squares: BO, 10 a.u., 10−6; filled triangles: BO, 100 a.u., 10−6; open diamonds: BO, 100 a.u., 10−5. Bottom: open circles: CP, 10 a.u.; open diamonds: BO, 100 a.u., 10−5; dashed line: BO, 100 a.u., 10−4. 2.3 Comparing BOMD and CPMD The comparison of the overall performance of Car–Parrinello and Born–Oppenheimer molecular dynamics in terms of computer time is a delicate issue. For instance it depends crucially on the choice made concerning the accuracy of the conservation of the energy Econs. Thus, this is to some extend subject of “personal taste” as to what is considered to be a ”sufficiently accurate” energy conservation. In addition, this comparison might lead to different conclusions as a function of type and size of the system investigated. Never- theless, in order to shed light on this point, microcanonical simulations of 8 silicon atoms were performed with various parameters using Car–Parrinello and Born–Oppenheimer molecular dynamics. This large–gap system was initially extremely well equilibrated and the runs were extended to 8 ps (and a few to 12 ps with no noticeable difference) at a temperature of about 360–370 K (with ±80 K root–mean–square fluctuations). The wavefunction was expanded up to Ecut = 10 Ry at the Γ–point of a simple cubic supercell and LDA was used to describe the interactions. In both cases the velocity Verlet scheme was used to integrate the equations of motion. In Car–Parrinello molecular dynamics two different time steps were used, 5 a.u. and 10 a.u. (corresponding to about 0.24 fs), in conjunction with a fictitious electron mass of µ = 400 a.u.. Within Born–Oppenheimer molecular dynamics the minimization of the energy functional was done using the highly efficient DIIS (direct inversion in the iterative subspace) scheme [31]. In this case, the time step was either 10 a.u. or 100 a.u. and three convergence criteria were used; note that the large time step corresponding to 2.4 fs is already at the limit to be used to in- vestigate typical molecular systems (with frequencies up to 4000 cm−1). The convergence criterion is based on the largest element of the wavefunction gradient which was required to be smaller than 10−6, 10−5 or 10−4 a.u.. The outcome of this comparison is shown in Fig. 2.2 in terms of the time evolution of the conserved energy Econs on scales that cover more than three orders of magnitude in absolute accuracy. Within the present comparison ultimate energy stability was obtained using Car–Parrinello molecular dynamics with the shortest time step of 5 a.u., which con- serves the energy of the total system to about 6×10−8 a.u. per picosecond, see solid line in Fig. 2.2(top). Increasing the time step to 10 a.u. leads to an energy conservation of about 3×10−7 a.u./ps and much larger energy fluctuations, see open circles in Fig. 2.2(top). The 20 Table 2.1: Timings in cpu seconds and energy conservation in a.u. / ps for Car–Parrinello (CP) and Born–Oppenheimer (BO) molecular dynamics simulations of a model system for 1 ps of trajectory on an IBM RS6000 / model 390 (Power2) workstation using the CPMD package [32]; see Fig. 2.2 for corresponding energy plots. Method Time step (a.u.) Convergence (a.u.) Conservation (a.u./ps) Time (s) CP 5 — 6×10−8 3230 CP 7 — 1×10−7 2310 CP 10 — 3×10−7 1610 BO 10 10−6 1×10−6 16590 BO 50 10−6 1×10−6 4130 BO 100 10−6 6×10−6 2250 BO 100 10−5 1×10−5 1660 BO 100 10−4 1×10−3 1060 computer time needed in order to generate one picosecond of Car–Parrinello trajectory increases linearly with the time step, see Table 2.1. The most stable Born–Oppenheimer run was performed with a time step of 10 a.u. and a convergence of 10−6. This leads to an energy conservation of about 1×10−6 a.u./ps, see filled squares in Fig. 2.2(top). As the maximum time step in Born–Oppenheimer dynamics is only related to the time scale associated to nuclear motion it could be increased from 10 to 100 a.u. while keeping the convergence at the same tight limit of 10−6. This worsens the energy conservation slightly (to about 6×10−6 a.u./ps), whereas the energy fluctuations increase dramati- cally, see filled triangles in Fig. 2.2(middle) and note the change of scale compared to Fig. 2.2(top). The overall gain is an acceleration of the Born–Oppenheimer simulation by a factor of about seven to eight, see Table 2.1. In the Born–Oppenheimer scheme, the computer time needed for a fixed amount of simulated physical time decreases only sub linearly with increasing time step since the initial guess for the iterative minimization degrades in quality as the time step is made larger. Further savings of computer time can be easily achieved by decreasing the quality of the wavefunction convergence from 10−6 to 10−5 and finally to 10−4, see Table 2.1. This is unfortunately tied to a signif- icant decrease of the energy conservation from 6×10−6 a.u./ps at 10−6 (filled triangles) to about 1×10−3 a.u./ps at 10−4 (dashed line) using the same 100 a.u. time step, see Fig. 2.2(bottom) but note the change of scale compared to Fig. 2.2(middle). In conclusion, Born–Oppenheimer molecular dynamics can be made as fast as (or even faster than) Car–Parrinello molecular dynamics (as measured by the amount of cpu time spent per picosecond) at the expense of sacrificing accuracy in terms of energy conservation. 21 where the sum over G vectors in Eq. (3.15) expands over double the range given by the wavefunction expansion. This is one of the main advantages of the plane wave basis. Whereas for atomic orbital basis sets the number of functions needed to describe the density grows quadratically with the size of the system, there is only a linear dependence for plane waves. In actual calculations the infinite sums over G vectors and cells has to be truncated. Furthermore, we have to approximate the integral over the Brillouin zone by a finite sum over special k–points ∫ dk→ ∑ k wk , (3.16) where wk are the weights of the integration points. From now on we will assume that the Brillouin zone integration can be done efficiently by a single point at k = 0, the so called Γ–point. The truncation of the plane wave basis rests on the fact that the Kohn–Sham potential V KS(G) converges rapidly with increasing modulus of G. For this reason only G vectors with a kinetic energy lower than a given maximum cutoff 1 2 |G|2 ≤ Ecut (3.17) are included in the basis. With this choice of the basis the precision of the calculation within the approximations of density functional theory is controlled by one parameter Ecut only. The number of plane waves for a given cutoff depends on the unit cell. A good estimate for the size of the basis is NPW = 1 2π2 Ω E 3/2 cut , (3.18) where Ecut is in Hartree units. The basis set needed to describe the density calculated from the Kohn-Sham orbitals has a corresponding cutoff that is four times the cutoff of the orbitals. The number of plane waves needed at a given density cutoff is therefore eight times the number of plane waves needed for the orbitals. 3.2 Kinetic Energy and Local Potentials Plane waves are eigenfunctions of the kinetic energy operator 1 2 ∇2 eiG·r = −1 2 | G |2 eiG·r . (3.19) The kinetic energy is therefore easily calculated in Fourier space Ekin = ∑ i ∑ G 1 2 fi |G|2 |ci(G)|2, , (3.20) 24 and the same is true for the wavefunction forces Fkin = 1 2 |G|2 ci(G) . (3.21) The plane waves do not depend on the atomic positions, therefore there are no Pulay forces and no contribution of the kinetic energy to the forces on the nuclei. Local operators act multiplicatively on wavefunctions in real space∫ dr′ V (r, r′)Φ(r′) = Vloc(r)Φ(r) . (3.22) The matrix elements of local operators can be calculated from the plane wave expansion of the operator in real space 〈G1 | Vloc(r) | G2〉 = 1 Ω ∑ G Vloc(G) ∫ dr e−iG1·reiG·reiG2·r (3.23) = 1 Ω ∑ G Vloc(G) ∫ dr ei(G−G1+G2)·r (3.24) = 1 Ω Vloc(G1 −G2) . (3.25) The expectation value only depends on the density Eloc = ∑ i fi〈Φi | Vloc | Φi〉 (3.26) = ∫ dr Vloc(r) n(r) (3.27) = 1 Ω ∑ G V ?loc(G) n(G) . (3.28) Expectation values are calculated in Fourier space as a sum over G–vectors. The local potential is multiplied by the density and therefore only those components of the local potential that are non–zero in the density have to be calculated. Forces are calculated in real space by multiplying the wavefunctions with the potential on the real space grid. 3.3 Electrostatic Energy The electrostatic energy of a system of nuclear charges ZI at positions RI and an electronic charge distribution n(r) consists of three parts: the Hartree energy of the electrons, the interaction energy of the electrons with the nuclei and the internuclear interactions EES = 1 2 ∫ ∫ dr dr′ n(r)n(r′) |r− r′| + ∑ I ∫ drV Icore(r)n(r) + 1 2 ∑ I 6=J ZIZJ |RI −RJ | . (3.29) 25 The Ewald method (see e.g. Ref. [2]) can be used to avoid singularities in the individual terms when the system size is infinite. In order to achieve this a Gaussian core charge distribution associated with each nuclei is defined nIc(r) = − ZI (RcI) 3π −3/2 exp −(r−RI RcI )2 . (3.30) It is convenient at this point to use a special definition for the core potential and define it to be the potential of the Gaussian charge distribution of Eq. (3.30) V Icore(r) = ∫ dr′ nIc(r ′) |r− r′| = − ZI |r−RI | erf [ |r−RI | RcI ] , (3.31) where erf is the error function. This potential has the correct long range behavior but we will have to add a correction potential for the short range part. The interaction energy of this Gaussian charge distributions is now added and subtracted from the total electrostatic energy EES = 1 2 ∫ ∫ dr dr′ n(r)n(r′) |r− r′| + 1 2 ∫ ∫ dr dr′ nc(r)nc(r ′) |r− r′| + ∫ ∫ dr dr′ nc(r)n(r ′) |r− r′| + 1 2 ∑ I 6=J ZIZJ |RI −RJ | − 1 2 ∫ ∫ dr dr′ nc(r)nc(r ′) |r− r′| ,(3.32) where nc(r) = ∑ I n I c(r). The first three terms can be combined to the electrostatic energy of a total charge distribution ntot(r) = n(r) + nc(r). The remaining terms are rewritten as a double sum over nuclei and a sum over self–energy terms of the Gaussian charge distributions EES = 1 2 ∫ ∫ dr dr′ ntot(r)ntot(r ′) |r− r′| + 1 2 ∑ I 6=J ZIZJ |RI −RJ | erfc  |RI −RJ |√ RcI 2 + RcJ 2 −∑ I 1√ 2π Z2I RcI , (3.33) where erfc denotes the complementary error function. For a periodically repeated system the total energy per unit cell is derived from the above expression by using the solution to Poisson’s equation in Fourier space for the first term and make use of the rapid convergence of the second term in real space. The total charge is expanded in plane waves with expansion coefficients ntot(G) = n(G) + ∑ I nIc(G)SI(G) (3.34) = n(G)− 1 Ω ∑ I ZI√ 4π exp [ −1 2 G2RcI 2 ] SI(G) . (3.35) 26 RI(t+ δt) = RI(t) + δt ˙̃ RI(t+ δt) ˙̃cI(t+ δt) = ċI(t) + δt 2µ fi(t) c̃i(t+ δt) = ci(t) + δt ˙̃ci(t+ δt) ci(t+ δt) = c̃i(t+ δt) + ∑ j Xij cj(t) calculate FI(t+ δt) calculate fi(t+ δt) ṘI(t+ δt) = ˙̃ RI(t+ δt) + δt 2MI FI(t+ δt) ċ′i(t+ δt) = ˙̃ci(t+ δt) + δt 2µ fi(t+ δt) ċi(t+ δt) = ċ′i(t+ δt) + ∑ j Yij cj(t+ δt) , where RI(t) and ci(t) are the atomic positions of particle I and the Kohn–Sham orbital i at time t respectively. Here, FI are the forces on atom I, and fi are the forces on Kohn– Sham orbital i. The matrices X and Y are directly related to the Lagrange multipliers by Xij = δt2 2µ Λpij (3.46) Yij = δt 2µ Λvij . (3.47) Notice that in the rattle algorithm the Lagrange multipliers to enforce the orthonormal- ity for the positions Λp and velocities Λv are treated as independent variables. Denoting with C the matrix of wavefunction coefficients ci(G), the orthonormality constraint can be written as C†(t+ δt)C(t+ δt)− I = 0 (3.48)[ C̃ + XC ]† [ C̃ + XC ] − I = 0 (3.49) C̃†C̃ + XC̃†C + C†C̃X† + XX† − I = 0 (3.50) XX† + XB + B†X† = I−A , (3.51) where the new matrices Aij = c̃ † i (t + δt)c̃j(t + δt) and Bij = c † i (t)c̃j(t + δt) have been introduced in Eq. (3.51). The unit matrix is denoted by the symbol I. By noting that A = I +O(δt2) and B = I +O(δt), Eq. (3.51) can be solved iteratively using X(n+1) = 1 2 [ I−A + X(n) (I−B) + (I−B)X(n) − ( X(n) )2] (3.52) 29 and starting from the initial guess X(0) = 1 2 (I−A) . (3.53) In Eq. (3.52) it has been made use of the fact that the matrices X and B are real and symmetric, which follows directly from their definitions. Eq. (3.52) can usually be iterated to a tolerance of 10−6 within a few iterations. The rotation matrix Y is calculated from the orthogonality condition on the orbital ve- locities ċ†i (t+ δt)cj(t+ δt) + c † i (t+ δt)ċj(t+ δt) = 0. (3.54) Applying Eq. (3.54) to the trial states Ċ′ + YC yields a simple equation for Y Y = −1 2 (Q + Q†), (3.55) where Qij = c † i (t + δt)ċ ′† i (t + δt). The fact that Y can be obtained without iteration means that the velocity constraint condition Eq. (3.54) is satisfied exactly at each time step. 30 Lecture 4 Pseudopotentials The norm–conserving pseudopotential approach provides an effective and reliable means for performing calculations on complex molecular, liquid and solid state systems using plane wave basis sets. In this approach only the chemically active valence electrons are dealt with explicitely. The inert core electrons are eliminated within the frozen– core approximation, being considered together with the nuclei as rigid non–polarizable ion cores. In turn, all electrostatic and quantum–mechanical interactions of the valence electrons with the cores, as the nuclear Coulomb attraction screened by the core electrons, Pauli repulsion and exchange and correlation between core and valence electrons, are accounted for by angular momentum dependent pseudopotentials. These reproduce the true potential and valence orbitals outside a chosen core region but remain much weaker and smoother inside. The valence electrons are described by smooth pseudo orbitals which play the same role as the true orbitals, but avoid the nodal structure near the nuclei that keeps the core and valence states orthogonal in an all–electron framework. The respective Pauli repulsion largely cancels the attractive parts of the true potential in the core region, and is built into the therefore rather weak pseudopotential. This pseudoization of the valence wavefunctions along with the removal of the core states eminently facilitates a numerically accurate solution of the Kohn–Sham equations and the Poisson equation, and enables the use of plane waves as an expedient basis set in electronic structure calculations. By virtue of the norm–conservation property and when constructed carefully pseudopotentials present a rather marginal approximation, and indeed allow for an adequate description of the valence electrons over the entire chemically relevant range of systems. 4.1 Why Pseudopotentials ? • Pseudopotentials should be additive and transferable. Additivity can most easily be achieved by building pseudopotentials for atoms in reference states. Transferability means that one and the same pseudopotential should be adequate for an atom in 31 1. V (1) l (r) = VAE(r)[1− f1 ( r rcl ) ] rcl : core radius ≈ 0.4 - 0.6 Rmax, where Rmax is the outermost maximum of the real wave function. 2. V (2) l (r) = V (1) l (r) + clf2 ( r rcl ) determine cl so that ̂l = l in (T + V (2) l (r))w (2) l (r) = ̂lw (2) l (r) 3. Φl(r) = γl [ w (2) l (r) + δlr l+1f3 ( r rcl )] where γl and δl are chosen such that Φl(r)→ Ψl(r) for r ≥ rcl and γ2l ∫ |w(2)l (r) + δlrl+1f3 ( r rcl ) |2dr = 1 4. Invert the Schrödinger equation for ̂l and Φl(r) to get V l val(r). 5. Unscreen Vlval(r) to get V l ps(r). V lps(r) = V l val(r)− VH(nv)− Vxc(nv) where VH(ρv) and Vxc(ρv) are the Hartree and exchange and correlation potentials of the pseudo valence density. Hamann, Schlüter and Chiang chose the following cutoff functions f1(x) = f2(x) = f3(x) = exp(−x4). These pseudopotentials are angular momentum dependent. Each angular momentum state has its own potential that can be determined independently from the other poten- tials. It is therefore possible to have a different reference configuration for each angular momentum. This allows it for example to use excited or ionic states to construct the pseudopotential for l states that are not occupied in the atomic ground state. The total pseudopotential in a solid state calculation then takes the form Vps(r) = ∑ L V Lps(r)PL where L is a combined index {l,m} and PL is the projector on the angular momentum state {l,m}. 34 4.2.2 Bachelet-Hamann-Schlüter (BHS) form Bachelet et al. [43] proposed an analytic fit to the pseudopotentials generated by the HSC recipe of the following form Vps(r) = Vcore(r) + ∑ L ∆V ionL (r) Vcore(r) = − Zv r [ 2∑ i=1 ccorei erf (√ αcorei r )] ∆V ionL (r) = 3∑ i=1 ( Ai + r 2Ai+3 ) exp(−αir2) The cutoff functions were slightly modified to be f1(x) = f2(x) = f3(x) = exp(−x3.5). They generated pseudopotentials for almost the entire periodic table (for the local density approximation), where generalizations of the original scheme to include spin–orbit effects for heavy atoms were made. Useful is also their list of atomic reference states. BHS did not tabulate the Ai coefficients as they are often very big numbers but another set of numbers Ci, where Ci = − 6∑ l=1 AlQil and Ai = − 6∑ l=1 ClQ −1 il with Qil =  0 for i > l[ Sil − ∑i−1 k=1Q 2 ki ]1/2 for i = l 1 Qii [ Sil − ∑i−1 k=1QkiQkl ]1/2 for i < l where Sil = ∫∞ 0 r 2ϕi(r)ϕl(r)dr, and ϕi(r) = { e−αir 2 for i = 1, 2, 3 r2e−αir 2 for i = 4, 5, 6 . 4.2.3 Kerker Pseudopotentials Also in this approach [44] pseudopotentials with the HSC properties are constructed. But instead of using cutoff functions (f1, f2, f3) the pseudo wavefunctions are directly constructed from the all–electron wavefunctions by replacing the all–electron wavefunc- tion inside some cutoff radius by a smooth analytic function that is matched to the all– electron wavefunction at the cutoff radius. The HSC properties then translate into a set of equations for the parameters of the analytic form. After having determined the pseudo 35 wavefunction the Schrödinger equation is inverted and the resulting potential unscreened. Note that the cutoff radius of this type of pseudopotential construction scheme is consid- erably larger than the one used in the HSC scheme. Typically the cutoff radius is chosen slightly smaller than Rmax, the outermost maximum of the all–electron wavefunction. The analytic form proposed by Kerker is Φl(r) = r l+1ep(r) with p(r) = αr4 + βr3 + γr2 + δ . The term linear in r is missing to avoid a singularity of the potential at r = 0. The HSC conditions can be translated into a set of equations for the parameters α, β, γ, δ. 4.2.4 Trouiller–Martins Pseudopotentials The Kerker method was generalized by Trouiller and Martins [45] to polynomials of higher order. The rational behind this was to use the additional parameters (the coefficients of the higher terms in the polynomial) to construct smoother pseudopotentials. The Trouiller–Martins wavefunctions has the following form Φl(r) = r l+1ep(r) with p(r) = c0 + c2r 2 + c4r 4 + c6r 6 + c8r 8 + c10r 10 + c12r 12 and the coefficients cn are determined from • norm–conservation • For n = 0 . . . 4 dnΦ drn ∣∣∣∣∣ r=rc = dnΨ drn ∣∣∣∣∣ r=rc • dΦ dr ∣∣∣∣∣ r=0 = 0 4.2.5 Kinetic Energy Optimized Pseudopotentials This scheme is based on the observation that the total energy and the kinetic energy have similar convergence properties when expanded in plane waves. Therefore, the kinetic 36 = ∑ L 1 Ω ∑ i wi∆V L(ri)P L? i (G)P L i (G ′) , (4.9) where the definition for the projectors P PLi (G) = 〈YL | G〉riω (4.10) has been introduced. The number of projectors per atom is the number of integration points (5 - 20 for low to high accuracy) multiplied by the number of angular momenta. For the case of s and p non–local components and 15 integration points this accounts to 60 projectors per atom. The integration of the projectors can be done analytically PLi (G) = ∫ ω Y ?L (ω)e iGridω (4.11) = ∫ ω Y ?L (ω)4π ∞∑ l=0 iljl(Gri) l∑ m′=−l Y ?lm′(ω)Ylm′(G)dω (4.12) = 4πiljl(Gri)YL(Ĝ) , (4.13) where the expansion of a plane wave in spherical harmonics has been used. jl are the spherical Bessel functions and Ĝ the angular components of the Fourier vector G. 4.3.2 Kleinman–Bylander Scheme The other method is based on the resolution of the identity in a local basis set∑ α | χα〉〈χα |= 1 , (4.14) where {χα} are orthonormal functions. This identity can now be introduced in the inte- grals for the non–local part V nl(G,G′) = ∑ L ∫ ∞ 0 dr〈G | YL〉ωr2∆V L(r)〈YL | G′〉ω = ∑ α,β ∑ L ∫ ∞ 0 dr〈G | χα〉〈χα | YL〉ωr2∆V L(r)〈YL | χβ〉ω〈χβ | G′〉 ,(4.15) and the angular integrations are easily performed using the decomposition of the basis in spherical harmonics χα(r) = χ lm α (r)Ylm(ω) . (4.16) This leads to V nl(G,G′) = ∑ α,β ∑ L 〈G | χα〉 ∫ ∞ 0 drχlmα (r)r 2∆V L(r)χlmβ (r)〈χβ | G′〉 (4.17) = ∑ α,β ∑ L 〈G | χα〉∆V lαβ〈χβ | G′〉 (4.18) 39 which is the non–local pseudopotential in fully separable form. The coupling elements of the pseudopotential ∆V lαβ = ∫ ∞ 0 drχlmα (r)r 2∆V L(r)χlmβ (r) (4.19) are independent of the plane wave basis and can be calculated for each type of pseudopo- tential once the expansion functions χ are known. The final question is now what is an optimal set of basis function χ. Kleinman and Bylander[47] proposed to use the eigenfunctions of the pseudo atom, i.e. the solutions to the calculations of the atomic reference state using the pseudopotential Hamiltonian. This choice of a single reference function per angular momenta guarantees nevertheless the correct result for the reference state. Now assuming that in the molecular environment only small perturbations of the wavefunctions close to the atoms occur, this minimal basis should still be adequate. The Kleinman–Bylander form of the projectors is∑ L | χL〉〈∆V LχL | 〈χL∆V LχL〉 = 1 , (4.20) where χL are the atomic pseudo wavefunctions. The plane wave matrix elements of the non–local pseudopotential in Kleinman–Bylander form is V KB(G,G′) = 〈G | ∆V LχL〉〈∆V LχL | G′〉〈χL∆V LχL〉 . (4.21) Generalizations of the Kleinman–Bylander scheme to more than one reference function were introduced by Blöchl [48] and Vanderbilt [49]. They make use of several reference functions, calculated at a set of reference energies. In transforming a semilocal to the corresponding Kleinman–Bylander(KB) pseudopoten- tial one needs to make sure that the KB–form does not lead to unphysical ”ghost” states at energies below or near those of the physical valence states as these would undermine its transferability. Such spurious states can occur for specific (unfavorable) choices of the underlying semilocal and local pseudopotentials. They are an artefact of the KB–form nonlocality by which the nodeless reference pseudo wavefunctions need to be the lowest eigenstate, unlike for the semilocal form [50]. Ghost states can be avoided by using more than one reference state or by a proper choice of the local component and the cutoff radii in the basic semilocal pseudopotentials. The appearance of ghost states can be analyzed by investigating the following properties: • Deviations of the logarithmic derivatives of the energy of the KB–pseudopotential from those of the respective semilocal pseudopotential or all–electron potential. • Comparison of the atomic bound state spectra for the semilocal and KB–pseudopotentials. • Ghost states below the valence states are identified by a rigorous criteria by Gonze et al. [50]. 40 4.4 Dual–Space Gaussian Pseudopotentials Pseudopotentials in the Kleinman–Bylander form have the advantage of requiring minimal amount of work in a plane wave calculation by still keeping most of the transferability and general accuracy of the underlying semilocal pseudopotential. However, one wonders if it would not be possible to generate directly pseudopotentials in the separable form fulfilling the Hamann–Schlüter–Chiang conditions. It was found [51] that indeed it is possible to optimize a small set of parameters defining an analytical form for the local and non–local form of a pseudopotential that fulfills those conditions and reproduces even additional properties leading to highly transferable pseudopotentials. The local part of the pseudopotential is given by Vloc(r) = −Zion r erf ( r̄/ √ 2 + exp [ −1 2 r̄2 ]) × [ C1 + C2r̄ 2 + C3r̄ 4 + C6r̄ 6 ] , (4.22) where erf denotes the error function and r̄ = r/rloc. Zion is the ionic charge of the atomic core, i.e. the total charge minus the charge of the valence electrons. The non–local contribution to the pseudopotential is a sum of separable terms Vl(r, r ′) = 3∑ i=1 3∑ j=1 l∑ m=−l Ylm(r̂)p l i(r) h l ij p l j(r)Y ? lm(r̂ ′) , (4.23) where the projectors pli(r) are Gaussians of the form pli(r) = √ 2rl+2(i−1) exp [ − r2 2r2 l ] r l+(4i−1)/2 l √ Γ [ l + 4i−1 2 ] , (4.24) where Γ is the gamma function. The projectors are normalized∫ ∞ 0 r2pli(r)p l i(r)dr = 1 . (4.25) This pseudopotential also has an analytical form in Fourier space. In both real and Fourier space, the projectors have the form of a Gaussian multiplied by a polynomial. Due to this property the dual–space Gaussian pseudopotential is the optimal compromise between good convergence properties in real and Fourier space. The multiplication of the wavefunction with the non–local pseudopotential arising from an atom can be limited to a small region around the atom as the radial projectors asymptotically tend to zero outside the covalent radius of the atom. In addition, a very dense integration grid is not required, as the projector is reasonably smooth because of its good decay properties in Fourier space. The parameters of the pseudopotential are found by minimizing a target function. This function is build up as the sum of the differences of properties calculated from the all– electron atom and the pseudo–atom. Properties included are the integrated charge and the eigenvalues of occupied and the lowest unoccupied states. 41 Lecture 5 Implementation 5.1 Total Energy and Gradients 5.1.1 Plane Wave Expansion The plane wave expansions introduced in the last lecture were for local potentials V local(r) = ∑ G V local(G) exp[iG · r] , (5.1) Kohn–Sham orbitals Φ(r) = 1√ Ω ∑ G ci(G) exp[iG · r] , (5.2) and the electron density n(r) = ∑ G n(G) exp[iG · r] . (5.3) 5.1.2 Total Energy Molecular dynamics calculations with interaction potentials derived from density func- tional theory require the evaluation of the total energy and derivatives with respect to the parameters of the Lagrangian. The total energy can be calculated as a sum of kinetic, external (local and non-local pseudopotential), exchange and correlation, and electrostatic energy Etotal = Ekin + E PP local + E PP nonlocal + Exc + EES . (5.4) The individual terms are defined by Ekin = ∑ i ∑ G 1 2 fi |G|2 |ci(G)|2 (5.5) 44 EPPlocal = ∑ I ∑ G ∆V Ilocal(G)SI(G)n ?(G) (5.6) EPPnonlocal = ∑ i fi ∑ I ∑ α,β∈I ( FαI,i )? hIαβF β I,i (5.7) Exc = Ω ∑ G xc(G)n ?(G) (5.8) EES = 2πΩ ∑ G 6=0 |ntot(G)|2 G2 + Eovrl − Eself . (5.9) The overlap between the projectors of the non-local pseudopotential and the Kohn–Sham orbitals has been introduced in the equation above FαI,i = 1√ Ω ∑ G P Iα(G)SI(G) c ? i (G) . (5.10) 5.1.3 Wavefunction Gradient Analytic derivatives of the total energy with respect to the parameters of the calculation are needed for stable molecular dynamics calculations. All derivatives needed are easily accessible in the plane wave pseudopotential approach. In the following Fourier space formulas are presented 1 fi ∂Etotal ∂c?i (G) = 1 2 G2 ci(G) + ∑ G′ V ?loc(G−G′)ci(G′) + ∑ I ∑ α,β ( FαI,i )? hIαβP I β (G)SI(G) , (5.11) where Vloc is the total local potential Vloc(G) = ∑ I ∆V Ilocal(G)SI(G) + Vxc(G) + 4π ntot(G) G2 . (5.12) Wavefunction gradients are needed in optimization calculations and in the Car-Parrinello molecular dynamics approach. 5.1.4 Nuclear Gradient The derivative of the total energy with respect to nuclear positions is needed for structure optimization and in molecular dynamics, that is ∂Etotal ∂RI,s = ∂EPPlocal ∂RI,s + ∂EPPnonlocal ∂RI,s + ∂EES ∂RI,s , (5.13) 45 as the kinetic energy Ekin and the exchange and correlation energy Exc do not depend directly on the atomic positions, the relevant parts are ∂EPPlocal ∂RI,s = −Ω ∑ G iGs ∆V I local(G)SI(G)n ?(G) (5.14) ∂EPPnonlocal ∂RI,s = ∑ i fi ∑ α,β∈I (FαI,i)? hIαβ  ∂F βI,i ∂RI,s + ( ∂FαI,i ∂RI,s )? hIα,βF β I,i  (5.15) ∂EES ∂RI,s = −Ω ∑ G 6=0 iGs n?tot(G) G2 nIc(G)SI(G) + ∂Eovrl ∂RI,s . (5.16) The contribution of the projectors of the non-local pseudopotentials is calculated from ∂FαI,i ∂RI,s = − 1√ Ω ∑ G iGs P I α(G)SI(G) c ? i (G,k) . (5.17) Finally, the real space part contribution of the Ewald sum is ∂Eovrl ∂RI,s = ∑′ J ∑ L  ZIZJ|RI −RJ − L|3 erfc  |RI −RJ − L|√ RcI 2 + RcJ 2  + 2√ π 1√ RcI 2 + RcJ 2 ZIZJ |RI −RJ − L|2 exp −|RI −RJ − L|2√ RcI 2 + RcJ 2  ×(RI,s −RJ,s − Ls) . (5.18) The self energy Eself is independent of the atomic positions and does not contribute to the forces. 5.2 Fast Fourier Transforms A function given as a finite linear combination of plane waves can also be defined as a set of functional values on a equally spaced grid in real space. The sampling theorem (see e.g. Ref. [13]) gives the maximal grid spacing that still allows to hold the same information as the expansion coefficients of the plane waves. The real space sampling points R are defined R = h Nq , (5.19) where N is a diagonal matrix with the entries 1/Ns and q is a vector of integers ranging from 0 to Ns − 1 (s = x, y, z). To fulfill the sampling theorem Ns has to be bigger than 2 max(gs) + 1. To be able to use fast Fourier techniques, Ns must be decomposable into small prime numbers (typically 2, 3, and 5). In applications the smallest number Ns that fulfills the above requirements is chosen. 46 orbitals are only calculated at the Γ–point then the wavefunctions can be taken as real quantities. The plane wave expansion coefficients of real functions have to following symmetry property c(−G) = c?(G) , (5.31) and c(0) is real. Therefore it is possible to store only half of the coefficients and recalculate the others whenever needed. In addition the symmetry can be used in calculating overlap integrals 〈Φi | Φj〉 = ∑ G c?i (G)cj(G) . (5.32) The sum can be restricted to half of the G–vectors 〈Φi | Φj〉 = ci(0)cj(0) + ∑ G 2 Re(c?i (G)cj(G)) . (5.33) This sum can be implemented efficiently using real arithmetic avoiding multiplication of complex numbers. Another direct use of the symmetry of the wavefunctions can be made when using Fourier transforms. The Fourier transform pair is defined by F (ω) = ∫ ∞ −∞ f(t)eiωtdt f(t) = 1 2π ∫ ∞ −∞ F (ω)e−iωtdω , where in our case t is the direct (real) space and ω is reciprocal space (G - space). We want to make use of the special structure of our wavefunctions f(t)isreal⇒ F (ω) = F (−ω)∗ , (5.34) that allows to perform two transforms together. First we investigate a real to complex transform. We define a new function g(t) = f1(t) + if2(t) then we get for the transformed function G(ω) = F1(ω) + iF2(ω) G(−ω) = F1(−ω) + iF2(−ω) = F1(ω) ∗ + iF2(ω) ∗ We can calculate the two new functions G(ω) +G(−ω) and G(ω)−G(−ω). G(ω) +G(−ω) = F1(ω) + F1(ω)∗ + i (F2(ω) + F2(ω)∗) = 2Re [F1(ω)] + 2iRe [F2(ω)] G(ω)−G(−ω) = 2iIm [F1(ω)]− 2Im [F2(ω)] 49 Table 5.1: Comparison of number of one dimensional FFT’s needed in a full transform and a transform making use of the sparsity of the data set in reciprocal space Reciprocal space to direct space transform Direct space to reciprocal space transform full transform sparse transform full transform sparse transform N2 π 16 N2 N2 N2 N2 1 2 N2 N2 1 2 N2 N2 N2 N2 π 16 N2 3 N2 ( 3 2 + π 16 ) N2 3 N2 ( 3 2 + π 16 ) N2 and we find F1(ω) = 1 2 (Re [G(ω) +G(−ω)] + iIm [G(ω)−G(−ω)]) F2(ω) = 1 2 (Im [G(ω) +G(−ω)] + iRe [G(ω)−G(−ω)]) For the complex to real transform we define G(ω) = F1(ω) + iF2(ω) then we get for the functions in direct (time) space g(t) = f1(t) + f2(t) f1(t) = Re [g(t)] f2(t) = Im [g(t)] Finally, we can take advantage of the fact that the wavefunction cutoff is only 1/4 of the density cutoff. Therefore in reciprocal space only the values of grid points inside a sphere of radius N/4 are non-zero, where for simplicity we assume a simple cubic box with N3 grid points. In a three dimensional FFT there will be many one dimensional transforms that can be avoided. In table 5.4 the amount of work for a full transform and a transform that makes use of the sparsity of the data set are compared. In both transforms the savings amount to about a factor of two. 5.5 Exchange and Correlation Functionals Gradient corrected exchange and correlation functionals and their potentials are defined as Exc = ∫ F (n,∇n)dr (5.35) 50 Vxc = δExc δn = ∂F ∂n − 3∑ α=1 d drα [ ∂F ∂(∇αn) ] = ∂F ∂n − 3∑ α=1 ∂ ∂rα [ 1 |∇n| ∂F ∂|∇n| ∂n ∂rα ] (5.36) The function F and its derivatives with respect to n and |∇n| have to be calculated in real space. The derivatives with respect to r can be calculated most easily in reciprocal space. The following scheme outlines the steps necessary to perform the calculation 1. n(R) FFT−→ n(G) 2. ∂n ∂rα = iGαn(G) 3. iGαn(G) INV FFT−→ ∂n ∂rα (R) 4. |∇n(R)| 5. F (n(R),∇n(R)) 6. ∂F (R) ∂n 7. 1 |∇n(R)| ∂F (R) ∂|∇n(R)| 8. H1(R) = ∂F (R) ∂n FFT−→ H1(G) H2α(R) = 1 |∇n(R)| ∂F (R) ∂|∇n(R)| ∂n ∂rα (R) FFT−→ H2α(G) 9. Vxc(G) = H 1(G)− 3∑ α=1 iGαH 2 α(G) 10. Vxc(G) INV FFT−→ Vxc(R) 51 the nearest neighbor distance; independent of the position of the maximum of the AE wavefunction. Only for the augmentation charges a small cutoff radius must be used to restore the moments and the charge distribution of the AE wavefunction accurately. The pseudized augmentation charges are usually treated on a regular grid in real space, which is not necessarily the same as the one used for the representation of the wavefunctions. The relation between the ultrasoft pseudopotential method and other plane wave based methods was discussed by Singh [35]. A closely related method to Vanderbilts ultrasoft pseudopotentials was introduced by Blöchl [112]. In the projector augmented–wave method (PAW) a linear transformation is defined that connects the PS and AE wavefunctions. Already Blöchl did mention the similarities of his approach to the ultrasoft pseudopotentials. A formal derivation of the ultrasoft pseudopotentials from the PAW equations was done by Kresse and Joubert [54]. 6.2.1 The PAW Transformation We start the derivation of the ultrasoft pseudopotentials and the PAW method by intro- ducing a linear transformation [112, 55]. The PAW method is based on a formal division of the whole space Ω in distinct regions: a collection of non overlapping spherical regions around each atom: atomic spheres region ⋃ a Ωa, and the remainder, the interstitial region ΩI: Ω = ΩI + ⋃ a Ωa . (6.7) It is clear that the plane wave basis, being the ideal choice in the interstitial region ΩI will have great difficulties describing the wavefunctions in the atomic spheres region. In the PAW method this problem is circumvented by introducing auxiliary wavefunctions which satisfies the following requirements. First, the auxiliary wavefunction Φ̃i(r) can be obtained from the AE wavefunction Φi(r) via a invertible linear transformation T | Φ̃i〉 = T | Φi〉 (6.8) | Φi〉 = T −1 | Φ̃i〉 (6.9) Second, Φ̃i(r) is smooth, i.e. can be represented by plane wave basis set of a practicle size, everywhere, including the atomic spheres region Φ̃i(r) = 1√ Ω ∑ G ci(G)e iG·r . (6.10) The first requirement ensures that the task of solving the Kohn–Sham equations can be equivalently reformulated in terms of Φ̃i(r), whereas the second requirement allows the entire process to be performed using the plane wave basis set. The actual construction of Φ̃i(r) from a given Φi(r) proceeds as follows. For each atom, we define a finite set of local basis functions {χaα} that is expected to accurately describe the 54 oscillating behavior of the relevant wavefunction Φi(r) within the corresponding atomic sphere. Associated with {χaα} we introduce a set of localized projector functions {paα} such that 〈paβ | χaα〉 = δαβ (6.11) paα(r) = 0 , ∀r ⊂ ΩI . (6.12) Using {χaα} and {paα}, the wavefunction Φi(r) in the atomic sphere region can be repre- sented as Φi(r) = ∑ α cai,αχ a α(r) + ∆ a i (r) , ∀r ⊂ Ωa . (6.13) The coefficients cai,α in the expansion (6.13) are given by cai,α = 〈paα | Φi〉 . (6.14) The correction | ∆ai 〉 = (1− ∑ α | χaα〉〈paα |) | Φi〉 , (6.15) reflects the incompleteness of the set {χaα}. As the size of the basis {χaα} gets larger, the local basis representation of Φi(r) becomes more accurate, and ∆ a i (r) goes to zero. To define a mapping into Φ̃i(r) an auxiliary smooth basis set {χ̃aα} is formed, subject to the following conditions. First, the basis functions χ̃aα(r) are smooth, i.e. expandable in terms of the plane wave basis of a practical size, everywhere including the atomic sphere region. Second, χ̃aα(r) merges differentiable into χ a α(r) outside the atomic sphere: χ̃aα(r) = χ a α(r) , ∀r ⊂ ΩI . (6.16) Third, both χ̃aα(r) and differences χ̃ a α(r) − χaα(r) form linearly independent sets. The smooth wavefunction Φ̃i(r) can be obtained based on the following prescription. Inside the atomic sphere region it is generated by replacing each occurrence of χaα(r) with χ̃ a α(r) in the expansion (6.13) Φ̃i(r) = ∑ α cai,αχ a α(r) + ∆ a i (r) , ∀r ⊂ Ωa , (6.17) whereas in the interstitial region it simply coincides with Φi(r): Φ̃i(r) = Φi(r) ∀r ⊂ ΩI . (6.18) The transformation can therefore be expressed as T = 1 + ∑ a ∑ α (| χ̃aα〉− | χaα〉)〈paα | . (6.19) 55 Its inverse can be obtained as T −1 = 1 + ∑ a ∑ α (| χaα〉− | χ̃aα〉)〈p̃aα | , (6.20) where a set of smooth projector functions {p̃aα} is defined as 〈p̃aα |= ∑ β (pa | χ̃a〉)−1αβ 〈p a β | . (6.21) It can be shown that similar to {paα}, the smooth projector functions {p̃aα} have the following properties 〈p̃aβ | χ̃aα〉 = δαβ (6.22) p̃aα(r) = 0 , ∀r ⊂ ΩI . (6.23) Furthermore, it is straightforward to prove that 〈p̃aα |= 〈paα | T −1 (6.24) and therefore the local basis expansion coefficients and the remainder can be alternatively represented as cai,α = 〈p̃aα | Φ̃i〉 . | ∆ai 〉 = (1− ∑ α | χ̃aα〉〈p̃aα |) | Φ̃i〉 . (6.25) The above two expressions show that if the basis {χaα} provides an accurate local repre- sentation for Φi(r), then the smooth basis {χ̃aα} provides an accurate local representation for Φ̃i(r) and vice versa. This is an important observation, since it is our objective to completely eliminate Φi(r)and seek for Φ̃i(r) directly. From a practical point of view, it is the inverse transformation T −1 that plays a major role in all the applications. The expression for T −1 involves basis sets Φi(r) and Φ̃i(r) and smooth projector functions {p̃aα}. If desired the projector functions {paα} can be found from 〈paα |= ∑ β (〈p̃a | χa〉)−1αβ 〈p̃ a β | . (6.26) In the following we will make the approximation that the local expansions are accurate and the remainder terms can be neglected. 6.2.2 Expectation Values Consider the expectation value of the general local or quasilocal operator O with respect to the Kohn–Sham orbital Φ 〈O〉 = 〈Φ | O | Φ〉 . (6.27) 56 〈Φ̃i | S | Φ̃j〉 = δij (6.51) S = 1 + ∑ a ∑ αβ qaαβ | p̃aα〉〈p̃aβ | (6.52) EKS = ∑ i fi〈Φ̃i | − 1 2 ∇2 | Φ̃i〉+ ∫ dr Vloc(r)ñ(r) + ∑ i fi ∑ a ∑ αβ 〈Φ̃i | p̃aα〉 ( Daαβ + ∫ dr Vloc(r)Qaαβ(r) ) 〈p̃aβ | Φ̃i〉 +EHxc[n(r)] + Eion . (6.53) These are the working equations for the Kohn–Sham method using ultrasoft pseudopo- tentials. The pseudopotentials are specified through the functions Vloc(r), the local pseu- dopotential; the augmentation charges Qaαβ(r) and their integrated values q a αβ; the non- local matrix elements Daαβ and the projector functions {p̃aα(r)}. Note that we made the transition from a formally all–electron method to a pseudopotential treatment by as- suming that the external potential is given as a norm–conserving pseudopotential. By doing so we of course also introduced the approximations coming with pseudopotentials, e.g. the frozen core approximation and the linearization of the exchange and correlation functional. Methods to calculate the parameters and functions needed for an ultrasoft pseudopotential are described in the literature [49, 56, 57]. The only difference in the energy expression to the form with fully nonlocal pseudopoten- tials is that the total charge density includes the augmentation charges Qaαβ(r). However, for the calculation of derivatives the special form of the metric S introduces many new terms. Starting from the extended energy functional Euspp({Φi},Λ,RN) = EKS({Φi},RN) + ∑ ij Λij ( 〈Φ̃i | S | Φ̃j〉 − δij ) , (6.54) one arrives at the force expressions δEuspp δ〈Φi | = fi H | Φi〉+ ∑ j Λij S | Φ̃j〉 (6.55) H = −1 2 ∇2 + Veff(r) + ∑ a ∑ αβ D̄aαβ | p̃aα〉〈p̃aβ | (6.56) Veff(r) = δEKS δn(r) = Vloc(r) + ∫ dr′ n(r′) | r− r′ | + Vxc(r ′) , (6.57) where Vxc(r ′) = δExc[n]/δn(r) and all the terms arising from the augmentation part of the electron density have been grouped together with the nonlocal part of the pseudopotential, by defining new coefficients D̄aαβ = D a αβ + ∫ dr Veff(r) Qaαβ(r) . (6.58) 59 The forces on the ions are calculated from the derivative of Euspp wrt. the nuclei positions ∂Euspp ∂RI = ∂EKS ∂RI + ∑ ij Λij〈Φ̃i | ∂S ∂RI | Φ̃j〉 (6.59) ∂S ∂RI = ∑ αβ qIαβ  ∣∣∣∣∣ ∂p̃Iα∂RI 〉 〈p̃Iβ | + | p̃Iα〉 〈 ˜∂pIβ ∂RI ∣∣∣∣∣∣  (6.60) ∂EKS ∂RI = ∂Eion ∂RI + ∫ dr ∂Vloc(r) ∂RI ñ(r) + ∫ dr Veff(r) ∑ αβ ∂QIαβ(r) ∂RI [∑ i 〈Φ̃i | p̃Iα〉〈p̃Iβ | Φ̃i〉 ] + ∑ αβ D̄Iαβ ∑ i [ 〈Φ̃i ∣∣∣∣∣ ∂p̃Iα∂RI 〉 〈p̃Iβ | Φ̃i〉+ 〈Φ̃i | p̃Iα〉 〈 ∂p̃Iβ ∂RI ∣∣∣∣∣ Φ̃i〉 ] (6.61) 6.2.4 PAW Energy Expression In contrast to the ultrasoft pseudopotential derivation the PAW approach avoids the intro- duction of a pseudopotential but rather works in the frozen core approximation directly. The sum over states is therfore restricted to valence electrons but the electronic densities always include contributions from the core electrons. In addition, the expansion functions χaα(r) have to be orthogonal to the core state on the atom (see Refs. ( [112, 54, 55]) for details). In the following we will assume a all–electron treatment although this is in practice never done. The PAW method works directly with the three sets of functions {χai }, {χ̃ai }, and {p̃ai }. These functions together with a local potential Vloc fully define the PAW energy expres- sion. Similar to the expectation values the total energy is divided into individual terms EKS = Ẽ + ∑ a Ea − ∑ a Ẽa . (6.62) The smooth part Ẽ, which is evaluated on regular grids in Fourier or real space, and the one–center contributions Ea and Ẽa, which are evaluated on radial grids in an angular momentum representation. The three contributions to EKS are Ẽ = ∑ i fi〈Φ̃i | − 1 2 ∇2 | Φ̃i〉+ ∫ dr Vloc(r) ñ(r) + EH[ñ+ n̂] + Exc[ñ] (6.63) Ea = ∑ i fi ∑ αβ 〈Φ̃i | p̃aα〉〈χα | − 1 2 ∇2 | χβ〉〈p̃aβ | Φ̃i〉+ EH[n a + nZa] + Exc[n a] (6.64) 60 Ea = ∑ i fi ∑ αβ 〈Φ̃i | p̃aα〉〈χ̃α | − 1 2 ∇2 | χ̃β〉〈p̃aβ | Φ̃i〉+ ∫ dr Vloc(r) ña(r) + EH[ñ a + n̂] + Exc[ñ a] . (6.65) The potential Vloc is an arbitrary potential localized in the augmentation regions. Its contribution to the total energy vanishes exactly because ñ(r) = ña(r) within the atomic spheres. Since the potential contributes only if the partial wave expansion is not complete, it is used to minimize truncation errors. In addition the point charge density nZa of the nuclei and a compensation charge density n̂ was introduced. The compensation charge density has the same multipole moments as the density na + nZa − ña and is localized within the atomic regions. Whereas the division of the kinetic energy and the exchange– correlation energy are straightforward from the formulas derived for the expectation value of local operators, the derivation of the electrostatic terms is rather involved [112, 54, 55] and will be omitted here. In fact, another compensation charge [112] is needed in order to be able to calculate the electrostatic term in the first energy contribution solely within an energy cutoff dictated by Φ̃i(r). 6.2.5 Integrating the CP Equations We will now derive the velocity Verlet equations for the Car–Parrinello molecular dy- namics method using ultrasoft pseudopotentials are the PAW method. The equations of motion are µ|φ̈i〉 = |ϕi〉+ ∑ j ΛijS({RI})|φj〉 MIR̈I = FI + ∑ i,j Λij〈φi|∇IS({RI})|φj〉 (6.66) where the electronic equation is cast in the abstract Dirac notation, and the forces are defined to be ϕi(r) = − δE δφ∗i (r) FI = − ∂E ∂RI (6.67) The electronic force is often written in the form |ϕi〉 = −fiH|φi〉 (6.68) where H was defined before. In the velocity Verlet scheme, the positions and veloc- ities are treated explicitly. That is, one carries the information {RI(t), ṘI(t)} and 61 Note that Eq.(6.80), although linear in Y is most easily solved iteratively. However, because the ionic velocities in Eq.(6.78) have been used explicitly in Eq.(6.79), there is no self-consistency problem. Once Eq.(6.80) has been solved, then the new velocities can be determined via Eq.(6.78) straightforwardly. The solution for Y(= Y(t + ∆t)) can be used to obtain an initial guess for the matrix X in the next step via X(0)(t+ ∆t) = X(t) + ∆tY(t+ ∆t) . (6.83) In the velocity Verlet algorithms, the iteration of the position step requires several ap- plications of the matrix ∆I . This matrix can either be stored or recalculated as needed. However, due to the size of this matrix (3 × NA × N2s , where NA is the number of atoms and Ns is the number of electronic states in the system), storing it for large systems is often not possible and recalculating it is expensive. To circumvent these problems the constraint nonorthogonal orbital method has been devised [22]. 6.3 Metals; Free Energy Functional In the free energy approach [58], the excited states are populated according to the Fermi– Dirac (finite–temperature equilibrium) distribution which is based on the assumption that the electrons ”equilibrate” more rapidly than the timescale of the nuclear motion. This means that the set of electronic states evolves at a given temperature ”isothermally” (rather than adiabatically) under the inclusion of incoherent electronic transitions at the nuclei move. Thus, instead of computing the force acting on the nuclei from the elec- tronic ground–state energy it is obtained from the electronic free energy as defined in the canonical ensemble. By allowing such electronic transitions to occur the free energy approach transcends the usual Born–Oppenheimer approximation. However, the approx- imation of an instantaneous equilibration of the electronic subsystem implies that the electronic structure at a given nuclear configuration {RI} is completely independent from previous configurations along a molecular dynamics trajectory. Due to this assumption the notion ”free energy Born–Oppenheimer approximation” was coined in Ref. [59] in a similar context. Certain non–equilibrium situations can also be modeled within the free energy approach by starting off with an initial orbital occupation pattern that does not correspond to any temperature in its thermodynamic meaning, see e.g. Refs. [60] for such applications. The free energy functional as defined in Refs. [58] is introduced most elegantly by starting the discussion for the special case of non–interacting Fermions Hs = − 1 2 ∇2 − ∑ I ZI |RI − r| (6.84) in a fixed external potential due to a collection of nuclei at positions {RI}. The associated grand partition function and its thermodynamic potential (”grand free energy”) are given 64 by Ξs(µV T ) = det 2 (1 + exp [−β (Hs − µ)]) (6.85) Ωs(µV T ) = −kBT ln Ξs(µV T ) , (6.86) where µ is the chemical potential acting on the electrons and the square of the determinant stems from considering the spin–unpolarized special case only. This reduces to the well– known grand potential expression Ωs(µV T ) = −2kBT ln det (1 + exp [−β (Hs − µ)]) = −2kBT ∑ i ln ( 1 + exp [ −β ( (i)s − µ )]) (6.87) for non–interacting spin–1/2 Fermions where {(i)s } are the eigenvalues of a one–particle Hamiltonian such as Eq. (6.84); here the standard identity ln detM = Tr lnM was invoked for positive definite M. According to thermodynamics the Helmholtz free energy F(NV T ) associated to Eq. (6.86) can be obtained from a Legendre transformation of the grand free energy Ω(µV T ) Fs(NV T ) = Ωs(µV T ) + µN + ∑ I<J ZIZJ |RI −RJ | (6.88) by fixing the average number of electrons N and determining µ from the conventional thermodynamic condition N = − ( ∂Ω ∂µ ) V T . (6.89) In addition, the internuclear Coulomb interactions between the classical nuclei were in- cluded at this stage in Eq. (6.88). Thus, derivatives of the free energy Eq. (6.88) with respect to ionic positions −∇IFs define forces on the nuclei that could be used in a (hypothetical) molecular dynamics scheme using non–interacting electrons. The interactions between the electrons can be ”switched on” by resorting to Kohn–Sham density functional theory and the concept of a non–interacting reference system. Thus, instead of using the simple one–particle Hamiltonian Eq. (6.84) the effective Kohn–Sham Hamiltonian Eq. (2.9) has to be utilized. As a result, the grand free energy Eq. (6.85) can be written as ΩKS(µV T ) = −2kBT ln [ det ( 1 + exp [ −β ( HKS − µ )])] (6.90) HKS = −1 2 ∇2 − ∑ I ZI |RI − r| + VH(r) + δΩxc[n] δn(r) (6.91) HKSφi = iφi (6.92) 65 where Ωxc is the exchange–correlation functional at finite temperature. By virtue of Eq. (6.87) one can immediately see that ΩKS is nothing else than the “Fermi–Dirac weighted sum” of the bare Kohn–Sham eigenvalues {i}. Whence, this term is the ex- tension to finite temperatures of the ”band–structure energy” contribution to the total electronic energy. In order to obtain the correct total electronic free energy of the interacting electrons the corresponding extra terms (properly generalized to finite temperatures) have to be included in ΩKS. This finally allows one to write down the generalization of the Helmholtz free energy of the interacting many–electron case FKS(NV T ) = ΩKS(µV T ) + µ ∫ dr n(r) + ∑ I<J ZIZJ |RI −RJ | −1 2 ∫ dr VH(r) n(r) + Ωxc − ∫ dr δΩxc[n] δn(r) n(r) (6.93) in the framework of a Kohn–Sham–like formulation. The corresponding one–particle density at the Γ–point is given by n(r) = ∑ i fi(β) |φi(r)|2 (6.94) fi(β) = (1 + exp [β (i − µ)])−1 , (6.95) where the fractional occupation numbers {fi} are obtained from the Fermi–Dirac distri- bution at temperature T in terms of the Kohn–Sham eigenvalues {i}. Finally, ab initio forces can be obtained as usual from the nuclear gradient of FKS, which makes molecular dynamics possible. By construction, the total free energy Eq. (6.93) reduces to that of the non–interacting toy model Eq. (6.88) once the electron–electron interaction is switched off. Another use- ful limit is the ground–state limit β → ∞ where the free energy FKS(NV T ) yields the standard Kohn–Sham total energy expression EKS after invoking the appropriate limit Ωxc → Exc as T → 0. Most importantly, stability analysis [58] of Eq. (6.93) shows that this functional shares the same stationary point as the exact finite–temperature functional due to Mermin [62], see e.g. the textbooks [23, 24] for introductions to density functional formalisms at finite temperatures. This implies that the self–consistent density, which defines the stationary point of FKS, is identical to the exact one. This analysis reveals furthermore that, unfortunately, this stationary point is not an extremum but a saddle point so that no variational principle and, numerically speaking, no direct minimization algorithms can be applied. For the same reason a Car–Parrinello fictitious dynamics ap- proach to molecular dynamics is not a straightforward option, whereas Born–Oppenheimer dynamics based on diagonalization can be used directly. The band–structure energy term can be evaluated by diagonalizing the Kohn–Sham Hamiltonian after a suitable ”preconditioning” [58]. Specifically, a second–order Trotter 66 Figure 6.1: Schematic view of the Hockney method to decouple images in the electrostatic energy. The first line shows the artificially replicated system consisting of the original computational cell and an empty duplicated cell. The Green’s function of this new periodic system is shown in the lower part of the figure. dimensional case is presented. The charge density is assumed to be non-zero only within an interval L and sampled on N equidistant points. These points are denoted by xp. The potential can then be written VH(xp) = L N ∞∑ p′=−∞ G(xp − xp′)n(xp′) (6.100) = L N N∑ p′=0 G(xp − xp′)n(xp′) (6.101) for p = 0, 1, 2, . . . N , whereG(xp−xp′) is the corresponding Green’s function. In Hockney’s algorithm this equation is replaced by the cyclic convolution ṼH(xp) = L N 2N+1∑ p′=0 G̃(xp − xp′)ñ(xp′) (6.102) where p = 0, 1, 2, . . . 2N + 1, and ñ(xp) = { n(xp) 0 ≤ p ≤ N 0 N ≤ p ≤ 2N + 1 (6.103) G̃(xp) = G(xp) − (N + 1) ≤ p ≤ N (6.104) ñ(xp) = ñ(xp + L) (6.105) G̃(xp) = G̃(xp + L) (6.106) The solution ṼH(xp) can be obtained by a series of fast Fourier transforms and has the desired property ṼH(xp) = VH(xp) for 0 ≤ p ≤ N . (6.107) To remove the singularity of the Green’s function at x = 0, G(x) is modified for small x and the error is corrected by using the identity G(x) = 1 x erf [ x rc ] + 1 x erfc [ x rc ] , (6.108) where rc is chosen such, that the short-ranged part can be accurately described by a plane wave expansion with the density cutoff. In an optimized implementation Hockney’s method requires the double amount of memory and two additional fast Fourier transforms 69 Table 6.1: Fourier space formulas for the Hartree energy, see text for definitions. Dim. periodic (G2/4π)VH(G) VH(0) 0 – (1− cos [RG])n(G) 2πR2n(0) 1 z (1 +R (GxyJ1(RGxy)K0(Rgz) −gzJ0(RGxy)K1(Rgz)))n(G) 0 2 x, y (1− (−1)gz exp [−GZ/2])n(G) 0 3 x, y, z n(G) 0 on the box of double size. Hockney’s method can be generalized to systems with period- icity in one (wires) and two (slabs) dimensions. It was pointed out [77] that Hockney’s method gives the exact solution to Poisson’s equation for isolated systems if the boundary condition (zero density at the edges of the box) are fulfilled. A different, fully reciprocal space based method, that can be seen as an approximation to Hockney’s method, was recently proposed [78]. The final expression for the Hartree energy is also based on the splitting of the Green’s function in Eq. (6.108) EES = 2πΩ ∑ G V MTH (G)n ? tot(G) + Eovrl − Eself . (6.109) The potential function is calculated from two parts, V MTH (G) = V̄H(G) + ṼH(G) , (6.110) where ṼH(G) is the analytic part, calculated from a Fourier transform of erfc ṼH(G) = 4π G2 ( 1− exp [ −G 2r2c 4 ]) n(G) (6.111) and V̄H(G) is calculated from a discrete Fourier transform of the Green’s function on an appropriate grid. The calculation of the Green’s function can be done at the beginning of the calculation and has not to be repeated again. It is reported [78] that a cutoff of ten to twenty percent higher than the one employed for the charge density gives converged results. The same technique can also be applied for systems that are periodic in one and two dimensions. If the boundary conditions are appropriately chosen, the discrete Fourier transforms for the calculation of V̄H(G) can be performed analytically [79]. This is possible for the limiting case where rc = 0 and the boundary conditions are on a sphere of radius R for the cluster. For a one-dimensional system we choose a torus of radius R and for the two- dimensional system a slab of thickness Z. The electrostatic potential for these systems are listed in Table 6.1, where Gxy = [ g2x + g 2 y ]1/2 and Jn and Kn are the Bessel functions of the first and second kind of integer order n. 70 Hockney’s method requires a computational box such that the charge density is negligible at the edges. This is equivalent to the supercell approach [80]. Practical experience tells that a minimum distance of about 3 Å of all atoms to the edges of the box is sufficient for most systems. The Green’s function is then applied to the charge density in a box double this size. The Green’s function has to be calculated only once at the beginning of the calculation. The other methods presented in this chapter require a computational box of double the size of the Hockney method as they are applying the artificially periodic Green’s function within the computational box. This can only be equivalent to the exact Hockney method if the box is enlarged to double the size. In plane wave calculations computational costs grow linearly with the volume of the box. Therefore Hockney’s method will prevail over the others in accuracy, speed, and memory requirements in the limit of large systems. The direct Fourier space methods have advantages through their easy implementation and for small systems, if not full accuracy is required, i.e. if they are used with smaller computational boxes. In addition, they can be of great use in calculations with classical potentials. 71 where for the small overlap matrix the notation Sm,m′(qs, qs+1) = ∫ L 0 dx Φ?qs,m(x)e −i 2π L xΦqs+1,m′(x) , (7.10) and ΦqM ,m(x) = Φq0,m(x) is implicitly understood. Finally we get for the electric polar- ization Pel = − e 2π lim L→∞ Im ln M−1∏ s=0 detS(qs, qs+1) . (7.11) The total dipole of a system is calculated as the sum of nuclear and electronic contributions Ptot = Pnuc + Pel . (7.12) Only the total dipole will be independent of the gauge and a proper choice of reference point will be discussed later. First we will give the formulas for the electronic contribution in three dimension and for the case where only the Γ point of the supercell Brillouin zone is used Pαel = − 2e 2π | Gα | Im ln detSα . (7.13) The additional factor of 2 comes from the assumed spin degeneracy. The matrix S is defined as Sαmn = 〈Φm | e−iGαx | Φn〉 . (7.14) The index α = 1, 2, 3 labels the reciprocal–lattice basis vectors {Gα}, Pαel is the projection of the electronic dipole moment along the direction defined by Gα, and Φm are the Γ point Kohn–Sham orbitals. The IR absorption coefficient α(ω) can be calculated from the formula [84] α(ω) = 4π ω tanh(βh̄ω/2) 3h̄n(ω)cV ∫ ∞ −∞ dt e−iωt〈P (t) · P (0)〉 , (7.15) where V is the volume of the supercell, T the temperature, β = 1/kBT , n(ω) is the refractive index, c is the speed of light in vacuum. The angular brackets indicate a statistical average. The correlation function 〈P (t) · P (0)〉 is calculated classically and quantum effect corrections are taken into account through the factor tanh(βh̄ω/2). More sophisticated treatments of quantum effects are also available. They consist of replacing the classical correlation function with a more involved procedure [85]. As mentioned above only the total dipole moment is independent of the reference point. This can cause some problems during a molecular dynamics simulation. The electronic contribution can only be calculated modulo the supercell size and a unfortunate choice of reference might lead to frequent changes of the dipole moment by amounts of the cell size. Therefore it is most convenient to use a dynamic reference calculated from the nuclear 74 positions and charges. If the reference point is chosen to be the center of charge of the nuclei, the nuclear contribution to the dipole will always be zero Q = 1∑ I ZI ∑ I ZIRI . (7.16) ZI are the nuclear charges and RI the nuclear positions within the supercell. The elec- tronic contribution is then Pel = − 2e 2π hd (7.17) d = tan−1 [ImDα/ReDα] (7.18) Dα = exp [ i(hT )−1Q ] det Sα , (7.19) where h is the matrix defining the supercell. 7.3 Localized Orbitals, Wannier Functions The representation of the electronic ground state in terms of localized Wannier orbitals [86] provides a powerful tool in the study of periodic solids. Recent advances in the formulation of a theory of electronic polarization [81] and the development of linear-scaling methods [63] have rejuvenated the use of Wannier functions as an analysis tool. Namely, Wannier functions afford an insightful picture to the nature of chemical bonding and aid in the understanding of classical chemical concepts (e.g. nonbonding electron pairs or valency) in terms of quantum mechanics. Wannier functions (WF) are defined in terms of a unitary transformation performed on the occupied Bloch orbitals (BO) [86]. One major problem in a practical calculation is their non-uniqueness. This is a result of the indeterminacy of the BO’s, which are, in the case of a single band, only determined up to a phase factor, in the multi-band case, up to an arbitrary unitary transformation among all occupied orbitals at every point in the Brillouin zone. As proposed recently by Marzari and Vanderbilt [87], one can resolve this non-uniqueness by requiring that the total spread of the localized function be minimal. This criterion is in close analogy with the Boys-Foster method [88] for finite systems, here one uses the spread defined through the conventional position operator. The new technique has been successfully applied to crystal systems and to small molecules within a general k-point scheme[87, 89]. An extension to disordered systems within the Γ-point approximation was recently performed[90]. This is of particular interest when one would like a localized orbital picture within the framework of Car-Parrinello molecular dynamics (CPMD). Here we examine the problem focusing on the Γ-point approximation only. Upon minimization of the spread functional the appropriate unitary transformation to the localized orbitals can be calculated. With explicit knowledge of the spread functional 75 we can derive the complete expressions required to implement the iterative minimization procedure. We begin by reviewing the work of Resta [93]. In his treatment, the fundamental object for studying localization of an electronic state within Born-Von Karman boundary conditions is the dimensionless complex number, z = ∫ L dx exp(i2πx/L) |ψ(x)|2 . (7.20) Here, L is the linear dimension, and ψ(x) denotes the wavefunction. By considering the definition of the spread of the wavefunction to be Ω = 〈x2〉 − 〈x〉2, where 〈· · ·〉 denotes an expectation value, Resta has shown that to O(1/L2) the functional for the spread in one-dimension to be, Ω = 1 (2π)2 ln |z|2. (7.21) One goal of this study is to generalize Eq. (7.20) to three-dimensions and obtain the appro- priate generalization of Eq. (7.21). Thus, we choose to study the following dimensionless complex number within Born-Von Karman boundary conditions, zI = ∫ V dr exp(iGI · r) |ψ(r)|2 . (7.22) Here, I labels a general reciprocal lattice vector, GI = lIb1 +mIb2 + nIb3, where bα are the primitive reciprocal lattice vectors, the integers l, m, and n are the Miller indices, V is the volume of the supercell, and ψ(r) denotes the wavefunction. We must find an appropriate function of the zI ’s that gives the three dimensional spread in the case of an arbitrary simulation cell. We proceed by noting that in a molecular dynamics simulation the cell parameters (primitive lattice vectors) to describe systems of general symmetry are given by a1, a2 and a3. It is convenient to form a matrix of these cell parameters, ↔ h= (a1, a2, a3) where the volume V of the simulation cell is given by the determinant of ↔ h. It is also very useful to define scaled coordinates, s = ↔ h −1 · r that lie in the unit cube. In molecular dynamics simulations, this allows one to perform periodic boundary conditions for systems with general symmetry by first transforming to the unit cube, performing cubic periodic boundary conditions, and transforming back to the general cell with the action of ↔ h One can also compute the reciprocal space vectors for systems of general symmetry with knowledge of the matrix of cell parameters. Thus, the I-th reciprocal lattice vector, GI = 2π (↔ h −1)T · ĝI . (7.23) Here, the superscript T denotes transposition, and ĝI = (lI ,mI , nI) is the I-th Miller index. We then substitute this expression into Eq. (7.22) and use the definition of r to obtain, zI = det ↔ h ∫ 1 0 ds exp ( i2πĝTI · s ) |ψ( ↔ h · s)|2 . (7.24) 76 where ĝ1 = (1, 0, 0), ĝ2 = (0, 1, 0), and ĝ3 = (0, 0, 1), and Im denotes the imaginary part. Again, the salient feature of Eq.(7.34) is that the expectation value of the exponential is invariant with respect to the choice of cell. Thus, a general equation for the expectation value of the position operator in supercells of arbitrary symmetry is, rα,n = − ∑ β ↔ hαβ 2π Im log zα,n . (7.35) We now proceed to determine the weights ωI as defined in the sum rule Eq. (7.30) for supercells of general symmetry. Recall that the metric, ↔ g will contain at most six in- dependent entries as defined by the case of least symmetry, triclinic. Thus, Eq. (7.30) is a linear set of six equations with six unknowns. We have freedom to choose the six Miller indices, ĝI of which we are to take the linear combinations of. For computational convenience of computing zi we choose the first six indices that take you from one to the next point in the Brillouin zone. Namely, ĝ1 = (1, 0, 0), ĝ2 = (0, 1, 0), ĝ3 = (0, 0, 1), ĝ4 = (1, 1, 0), ĝ5 = (1, 0, 1), ĝ6 = (0, 1, 1). With this choice of ĝi the explicit system of equations based on Eq. (7.30) takes the following simple form, 1 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 1 0 1 0 1 0 0 0 0 0 1 0 0 1 0 1 1   ω1 ω2 ω3 ω4 ω5 ω6  =  g11 g12 g13 g22 g23 g33  (7.36) Thus, the solution to Eq. (7.36) yields the following set of general weights, ω1 = g11 − g12 − g13 ω2 = g22 − g12 − g23 ω3 = g33 − g13 − g23 ω4 = g12 ω5 = g13 ω6 = g23 (7.37) Eq. (7.37) indeed reduces to the specific cases computed in Ref.[91]. However, here, the case for triclinic symmetry is also included. Thus, with knowledge of the cell parameters, in conjunction with Eq. (7.31) allows one to compute the maximally localized WF. 79 Lecture 8 Chemical Reactions with CPMD 8.1 Simulated Annealing By adding a friction term, Car–Parrinello molecular dynamics can be turned into a damped second order dynamics scheme. The friction can be applied both to the nuclear degrees of freedom and the electronic coordinates. The resulting dynamics equation are a powerful method to simultaneously optimize the atomic structure and the Kohn–Sham orbitals [37, 109]. Harmonic reference system integration and plane wave dependent electron masses, are especially helpful in this context, as the derived dynamics does not have a direct physical relevance. Introducing a friction force proportional to the constants γn and γe the equations of motion can readily be integrated using the velocity Verlet algorithm. The friction terms translate into a simple rescaling of the velocities at the beginning and end of the time step according to ṘI(t) = γnṘI(t) ċi(t) = γeċi(t) VELOCITY VERLET UPDATE ṘI(t+ δt) = γnṘI(t+ δt) ċi(t+ δt) = γeċi(t+ δt) . It was shown [37, 109] that this scheme leads to efficient optimizations of all degrees of freedom and a combination with geometrical constraints is straightforward. One advantage of the velocity Verlet integrator is that it can be easily combined with multiple time scale algorithms [22] and still results in reversible dynamics. The most successful implementation of a multiple time scale scheme in connection with the plane wave–pseudopotential method is the harmonic reference system idea [110, 22]. The high frequency motion of the plane waves with large kinetic energy is used as a reference system for the integration. The dynamics of this reference system is harmonic and can 80 be integrated analytically. In addition, this can be combined with the basic notion of a preconditioner already introduced in the section on optimizations. The electronic mass used in the Car–Parrinello scheme is a fictitious construct and it is allowed to generalize the idea by introducing different masses for different ”classical” degrees of freedom [111, 109, 22], µ(G) = { µ H(G,G) ≤ α (µ/α) (1 2 G2 + V(G,G)) H(G,G) ≥ α , (8.1) where H and V are the matrix elements of the Kohn–Sham matrix and the potential respectively. The reference electron mass is µ. With the preconditioned masses and the harmonic reference system, the equations of motion of the system are µ(G)c̈i(G) = −λ(G)ci(G) + δΦi(G) + ∑ j Λijcj(G) . (8.2) where δΦi(G) is the force on orbital i minus −λ(G). From Eq. (8.2) it is easy to see that the frequencies ω(G) = √ λ(G)/µ(G) are independent of G and that there is only one harmonic frequency equal to √ α/µ. The revised formulas for the integration of the equations of motion for the velocity Verlet algorithm can be found in the literature [22]. The implications of the G vector dependent masses can be seen by revisiting the formulas for the characteristic frequencies of the electronic system. The masses µ are chosen such that all frequencies ωij are approximately the same, thus optimizing both, adiabaticity and maximal time step. The disadvantage of this method is that the average electron mass seen by the nuclei is drastically enhanced, leading to renormalization corrections [112] on the masses MI that are significantly higher than in the standard approach and not as simple to estimate by an analytical expression. 8.2 Rare Events A problem that is frequently encountered in simulations of chemical reactions is that of computing the free energy profile along a reaction path on a potential energy surface characterized by a reaction coordinate q. Not only does the free energy profile provide a thermodynamic picture along the reaction path, but it also permits determination of activation energies and estimation of associated rate constants via classical or quantum transition state theory. In many instances, a multidimensional free energy profile or sur- face characterized by a set of reaction coordinates may be of particular interest. However, the ability to characterize a process by a reaction coordinate or coordinates generally assumes some prior knowledge of the process under consideration. This is not always the case and one has to choose alternative methods like transition path sampling [95] or bias potential methods [96, 97] discussed in the subsequent sections. 81 calculate FI(t+ δt) Ṙ′I = ˙̃ RI + δt 2MI FI(t+ δt) ṘI(t+ δt) = Ṙ′I + δt 2MI gv(t+ δt) , where the constraint forces are defined by gp(t) = − ∑ i λip ∂σ(i)({RI(t)}) ∂RI (8.11) gv(t) = − ∑ i λiv ∂σ(i)({RI(t)}) ∂RI . (8.12) The Lagrange multiplier have to be determined to ensure that the constraint on the positions and velocities are exactly fulfilled at the end of the time step. For the position, the constraint condition is σ(i)({RI(t+ δt)}) = 0 . (8.13) Eq. (8.13) is in general a system of nonlinear equations in the Lagrange multipliers λip. These equations can be solved using a generalized Newton algorithm [13] that can be com- bined with a convergence acceleration scheme based on the direct inversion in the iterative subspace method [31, 108]. The error vectors for a given set of Lagrange multipliers λ are calculated from ei(λ) = − ∑ j J−1ij (λ)σ (j)(λ) . (8.14) The Jacobian J is defined by Jij(λ) = ∂σ(i)(λ) ∂λj (8.15) = ∑ I ∂σ(i)(λ) ∂RI(λ) ∂RI(λ) ∂λj (8.16) = − ∑ I δt2 2MI f cI (λ)f c I (0) , (8.17) where f cI (λ) = ∑ i λ i∂σ(i)/∂RI . Typically only a few iterations are needed to converge the Lagrange multipliers to an accuracy of 1× 10−8. The constraint condition for the velocities can be cast into a system of linear equations. Again, as in the case of the orthonormality constraints in the Car–Parrinello method, the Lagrange multiplier for the velocity update can be calculated exactly without making use of an iterative scheme. Defining the derivative matrix AiI = ∂σ(i) ∂RI , (8.18) 84 the velocity constraints are σ̇(i)(t+ δt) = 0 (8.19)∑ I ∂σ(i) ∂RI ṘI = 0 (8.20) − ∑ j (∑ I δt2 2MI AiIAjI ) λvj = ∑ I AiIṘ′I . (8.21) The only information needed to implement a new type of constraint are the formulas for the functional value and its derivative with respect to the nuclear coordinates involved in the constraint. 8.2.2 Adiabatic Free Energy Sampling This method [104, 105] is based on the creation of a dynamical adiabatic separation between the reaction coordinate and remaining degrees of freedom. In particular, a dy- namical scheme is constructed in which the reaction coordinate evolves slowly relative to the other degrees of freedom and is simultaneously maintained at a high temperature. The latter condition ensures that all activation barriers along the reaction path can be easily crossed and can be enforced by coupling the reaction coordinate to its own heat bath or thermostat. The former condition permits the remaining degrees of freedom to fully relax in response to the motion of the reaction coordinate and ,thereby, sample a large portion of their configuration space as the reaction coordinate slowly evolves. An analysis of the resulting dynamics reveals that the free energy profile will be given by Eq. 8.3 with β replaced by βs = 1/kTs, where Ts is the temperature of the reaction coor- dinate. This approach eliminates biasing from the simulation procedure and the need for post processing of the output data. Let us consider the case where a straightforward partitioning of the system into a reactive subsystem (S) and the environment (E) is possible. A general separation of degrees of freedom will require a slightly more elaborate formulation. For each subsystem, positions qSi ,q E j , momenta p S i ,p E j , and masses M S i ,M E j are defined and the number of atoms in the subsystems is NS and NE. The Hamiltonian of the system is then H(q,p) = NS∑ i ( pSi )2 2MSi + NE∑ j ( pEj )2 2MEj + +V ( qS,qE ) . (8.22) Each subsystem is then coupled to its own thermostat with a target temperature T S and TE. The crucial step, which is needed to obtain the dynamics of the reactive system on the physical free energy surface at the specific temperature, is the rescaling of the atomic masses so that MS/ME 1. This ensures that the frequency spectra of the subsystems do not overlap. In particular the reactive system is varying slowly, as compared to the 85 environment. This has two consequences: 1) The systems are decoupled and are in equilibrium although they are at different temperatures; 2) the slowly moving system feels the average force of the environment, which is at the physical temperature, and its dynamics is performed on this free energy surface. Note that this separation of systems is analogous to the two systems used in the Car–Parrinello method. The fast degrees of freedom qEj ,p E j sample the canonical distribution at temperature T E for every configuration qSj ,p S j of the reactive system. The dynamics of the slow degrees of freedom approximates the one resulting from a Hamiltonian with the effectively averaged potential. This potential V̄ ( qS ) is, for a full decoupling given by exp [ −βEV̄ ( qS )] = ∫ exp [ −βEV ( qS,qE )] dqE∫ exp [−βEV (qS,qE)] dqEdqS . (8.23) This potential of mean force equals the physical free energy surface of the reactive system. An enhanced canonical sampling of this free energy surface can be obtained by choosing T S > TE. Because we have assumed adiabatic evolution, ergodic behavior, and hence equilibrium sampling, the probability distribution of the full system in configuration space is given by ρ ( qS,qE ) = exp [ −(βS − βE)V̄ ( qS )] ∫ exp [ −βSV̄ (qS) ] dqS × exp [ −βEV ( qS,qE )] ∫ exp [−βEV (qS,qE)] dqEdqS . (8.24) This probability distribution can be used to calculate all thermodynamic properties of the physical system from a trajectory generated within the adiabatically decoupled system. Efficiency and accuracy of this method depends on the choice of subsystems and the free parameters MS/ME and the temperature T S. 8.2.3 Bias Potential Methods A long–standing problem of molecular dynamics is that its utility is limited to processes that occur on a time scale of nanoseconds or less. This is especially true for ab initio molecular dynamics where the time scale is typically only a few tens of picoseconds. For many systems, the dynamics can be characterized as a sequence of infrequent tran- sitions from one potential basin or state to another. In these cases, longer times scales can be accessed using transition state theory (TST), an elegant approach with a long history [3]. In TST, one takes the transition rate between states as the flux through a dividing surface separating the states. This flux is an equilibrium property of the system, and so does not require that actual dynamics be performed. TST assumes that each crossing of the dividing surface corresponds to a true reactive event, in which the system passes from one state to another and then loses all memory of this transition before the next event. In actuality, some surface crossings can be dynamically connected, i.e. if 86 Lecture 9 Linear Response 9.1 Adiabatic Density–Functional Perturbation The- ory The total energy and charge density are basic quantities of DFT, and give access to a wide number of experimental observables. Most of these quantities can be derived from deriva- tives of the energy or density with respect to external perturbations. As an example, the force exerted on a nucleus is equal to minus the derivative of the total energy of the system. The calculation of such energy derivatives can be done by finite–difference methods: the total energy is obtained at slightly different values of the external perturbation, then the derivative of the total energy curve with respect to the small disturbance is calculated nu- merically. Although this is a very convenient method, the recent practice has shown that perturbative techniques within DFT are more powerful. Such techniques are very similar to the treatment of perturbations within Hartree–Fock theory (the coupled–perturbed Hartree–Fock formalism) [113]. Such techniques were discovered and rediscovered within DFT many times. They are based either on the Sternheimer equation, Green’s functions, sum–over–states techniques, or on the Hylleras variational technique. Generalizations to arbitrary order of perturbations, based on the 2n+1 theorem of perturbation theory and on generalized Sternheimer equations were also proposed. 9.2 Coupled Perturbed Kohn–Sham Equations Let us consider the extended Kohn–Sham energy functional EKS[{Φi}] = EKS[{Φi}] + λEpert.[{Φi}] + ∑ ij Λij (〈Φi | Φj〉 − δij) . (9.1) All quantities are expanded in powers of λ in the form X(λ) = X(0) + λX(1) + λ2X(2) + · · · , (9.2) 89 Table 9.1: List of some properties related to derivatives of the total energy. α is an external electrical field, xi a displacement of a nuclei, Bα a magnetic field, mi the nuclear magnetic moment. Derivative Observable d2E dαdβ polarizability d2E dxidxj harmonic force constants d2E dxidα dipole derivatives, infrared intensities d2E dxidαdβ polarizability derivative, Raman intensities d2E dBαdBβ magnetazibility d2E dmidBβ nuclear magnetic shielding tensor where X can be the Kohn–Sham energy EKS, the Kohn–Sham orbitals Φi(r), the electron density n(r), the Lagrange multipliers Λ, or the Kohn–Sham Hamiltonian HKS. Because the Kohn–Sham energy satisfies a variational principle under constraints, it is possible to derive a constraint variational principle for the 2nth order derivative of the energy with respect to the nth order derivative of the wavefunctions Φi(r) [115, 114]: when the expansion of the wavefunction up to an order of n− 1 is known, the variational principle for the 2nth order derivative of the energy is as follows: E(2n) = min Φ (n) i ( E [∑ k λkφ (k) i ])(2n) (9.3) under the constraint, in the parallel transport gauge n∑ k=0 〈Φ(n−k)i | Φ (k) j 〉 = 0 (9.4) for all occupied states i and j. The requirement that the constraint condition (orthogonal- ity in our case) is fulfilled at every order in the perturbation was used to derive the above equations. It turns out that there is a certain degree of arbitrariness in the definition of the Lagrange multipliers. By adding additional conditions that leave energy invariant a special gauge for the Lagrange multipliers is chosen. The best known of such a gauge is the canonical gauge for the zeroth order wavefunction in Hartree–Fock and Kohn–Sham theory. In this gauge the matrix of Lagrange multipliers is diagonal. For the perturbation treatment at hand the parallel transport gauge defined above is most convenient. The explicit expressions for E(2n) can be derived by introducing Eq. 9.1 into Eq. 9.3. For zeroth order we get back the Kohn–Sham equations H(0) | Φ(0)i 〉 = ∑ j Λij | Φ(0)j 〉 , (9.5) 90 and for second order the following expression is obtained E(2) [ {Φ(0)i }, {Φ (1) i } ] = ∑ ij 〈Φ(1)i | H(0)δij + Λij | Φ (1) j 〉 + ∑ i { 〈Φ(1)i | δEpert δ〈Φi | + δEpert δ | Φi〉 | Φ(1)i 〉 } + 1 2 ∫ ∫ K(r, r′) n(1)(r) n(1)(r′) dr dr′ + ∫ d dλ δEHxc δn(r) ∣∣∣∣∣ n(0) n(1)(r) dr + 1 2 d2EHxc dλ2 ∣∣∣∣∣ n(0) , (9.6) where the first–order wavefunctions are varied under the constraints 〈Φ(0)i | Φ (1) j 〉 = 0 , (9.7) for all occupied states i and j. The density and the first order density are given by n(0)(r) = ∑ i fi | Φ(0)i (r) |2 (9.8) n(1)(r) = ∑ i fi[Φ (0)? i (r)Φ (1) i (r) + Φ (1)? i (r)Φ (0) i (r)] . (9.9) Further, we have introduced the zeroth order Kohn–Sham Hamiltonian H(0)[{Φ(0)i }] = 1 2 ∇2 + V(0)ext(r) + ∫ n(0)(r′) | r− r′ | dr′ + δExc δn(r) ∣∣∣∣∣ n(0) , (9.10) the combined Hartree and exchange–correlation energy EHxc[n] = 1 2 ∫ ∫ n(0)(r)n(0)(r′) | r− r′ | dr′dr + Exc[n] , (9.11) and the second order energy kernel K(r, r′) = δ 2EHxc[n] δn(r)δn(r′) ∣∣∣∣∣ n(0) . (9.12) Since E(2) is variational with respect to Φ (1) i the Euler–Lagrange, or in this case self– consistent Sternheimer equations [116, 117] can be deduced Pc ∑ j ( H(0)δij + Λ(0)ij ) Pc | Φ(1)j 〉 = −PcH(1) | Φ (0) i 〉 , (9.13) where Pc is the projector upon the unoccupied states Pc = 1− ∑ i | Φ(0)i 〉〈Φ (0) i | (9.14) 91 and taking the second derivative: ∂2EKS ∂RIα∂RJβ = ∑ k (〈 Φ (0) k ∣∣∣ ∂2Vext ∂RIα∂RJβ ∣∣∣Φ(0)k 〉+ 〈Φ(1)k,Iα∣∣∣ ∂Vext∂RJβ ∣∣∣Φ(0)k 〉+ 〈Φ(0)k ∣∣∣ ∂Vext∂RJβ ∣∣∣Φ(1)k,Iα〉 ) = ∑ k [〈 Φ (0) k ∣∣∣ ∂2Vloc ∂RIα∂RJβ ∣∣∣Φ(0)k 〉+∑ L ( 〈Φ(0)k ∣∣∣∣∣ ∂2pL∂RIα∂RJβ 〉 ωL〈pL | Φ(0)k 〉 +〈Φ(0)k | pL〉ωL 〈 ∂2pL ∂RIα∂RJβ ∣∣∣∣∣Φ(0)k 〉+ 2〈Φ(0)k ∣∣∣∣∣ ∂pL∂RIα 〉 ωL 〈 ∂pL ∂Rβj ∣∣∣∣∣Φ(0)k 〉 ) + 〈 Φ (1) k,Iα ∣∣∣ ∂Vloc ∂RJβ ∣∣∣Φ(0)k 〉+ 〈Φ(0)k ∣∣∣ ∂Vloc∂RJβ ∣∣∣Φ(1)k,Iα〉 + ∑ L ( 〈Φ(1)k,Iα ∣∣∣∣∣ ∂pL∂RJβ 〉 ωL〈pL | Φ(0)k 〉+ 〈Φ (1) k,Iα | pL〉ωL 〈 ∂pL ∂RJβ ∣∣∣∣∣Φ(0)k 〉 +〈Φ(0)k ∣∣∣∣∣ ∂pL∂RJβ 〉 ωL〈pL | Φ(1)k,Iα〉+ 〈Φ (0) k | pL〉ωL 〈 ∂pL ∂Rβj ∣∣∣∣∣Φ(1)k,Iα〉 )] (9.25) 9.3.1 Selected Eigenmodes of the Hessian Harmonic vibrational frequencies are calculated from the eigenvalues of the dynamical matrix. If only part of the eigenvalues are needed an iterative diagonalization scheme, e.g. Lanczos or Davidson diagonalization could be used. This might be useful in geom- etry optimization or in the location of saddle points, as well as when only certain type of vibrations are of interest, for example surface modes or hydrogen stretching modes. Iterative diagonalization scheme are based on matrix vector multiplication a = Hq . (9.26) This is of special interest when the calculation of the matrix H can be avoided and only an algorithm for the application of a vector to the matrix has to be available. Let us consider a collective displacement [119] q = ∑ I,α qIαeIα , (9.27) where the indices I, α run over all atoms and coordinates x, y, z respectively, and eIα is a 3Natom Cartesian unit vector. The matrix vector multiplication then becomes aJβ = ∑ Iα ∂2E ∂RJβ∂RIα qIα . (9.28) This can be calculated using density functional perturbation theory where we again as- sumed only one projector per angular momentum channel. Therefore the perturbation 94 functional is Epert[n] = ∫ dr n(r) ∑ Iα qIα ∂Vext(r−RI) ∂RIα = ∫ dr n(r) ∑ Iα qIα [ ∂Vloc(r−RI) ∂RIα + ∑ k ∑ L ( 〈Φk ∣∣∣∣∣ ∂pL∂RIα 〉 ωL〈pL|Φk〉+ 〈Φk|pL〉ωL 〈 ∂pL ∂RIα ∣∣∣∣∣Φk〉 )] . (9.29) Using this perturbation expression the response Φqi of the Kohn–Sham wavefunction to the collective displacement qI,α can be calculated. And with the response function the matrix vector multiplication becomes aJβ = ∑ Iα ∂2E ∂RJβ∂RIα qIα = ∂2E ∂RJβ∂q = ∑ k (∑ Iα qIα 〈 Φ (0) k ∣∣∣ ∂2Vext ∂RIα∂RJβ ∣∣∣Φ(0)k 〉 + 〈 Φ (1) k,q ∣∣∣ ∂Vext ∂RJβ ∣∣∣Φ(0)k 〉+ 〈Φ(0)k ∣∣∣ ∂Vext∂RJβ ∣∣∣Φ(1)k,q〉 ) = ∑ k [∑ Iα qIα 〈 Φ (0) k ∣∣∣ ∂2Vloc ∂RIα∂RJβ ∣∣∣Φ(0)k 〉 + ∑ L ∑ Iα qIα ( 〈Φ(0)k ∣∣∣∣∣ ∂2pL∂RIα∂RJβ 〉 ωL〈pL | Φ(0)k 〉 +〈Φ(0)k | pL〉ωL 〈 ∂2pL ∂RIα∂RJβ ∣∣∣∣∣Φ(0)k 〉+ 2〈Φ(0)k ∣∣∣∣∣ ∂pL∂RIα 〉 ωL 〈 ∂pL ∂Rβj ∣∣∣∣∣Φ(0)k 〉 ) + 〈 Φ (1) k,q ∣∣∣ ∂Vloc ∂RJβ ∣∣∣Φ(0)k 〉+ 〈Φ(0)k ∣∣∣ ∂Vloc∂RJβ ∣∣∣Φ(1)k,q〉 + ∑ L ( 〈Φ(1)k,q ∣∣∣∣∣ ∂pL∂RJβ 〉 ωL〈pL | Φ(0)k 〉+ 〈Φ (1) k,q | pL〉ωL 〈 ∂pL ∂RJβ ∣∣∣∣∣Φ(0)k 〉 +〈Φ(0)k ∣∣∣∣∣ ∂pL∂RJβ 〉 ωL〈pL | Φ(1)k,q〉+ 〈Φ (0) k | pL〉ωL 〈 ∂pL ∂Rβj ∣∣∣∣∣Φ(1)k,q〉 )] (9.30) 9.4 Polarizability As discussed in the introduction, a case in which the perturbation cannot be expressed in a Hamiltonian form is that of an external electric field, which couples with the electric polarization Pele = e〈r〉 in a periodic system. The latter using the modern theory of 95 polarization in the Γ only sampling of the Brillouin zone can be written in terms of Berry phase ϕµ[83]: ϕµ = Im log detQ (µ), (9.31) where the matrix Q(µ) is defined as Q (µ) i,j = 〈 Φi ∣∣∣eiGµ·r∣∣∣Φj〉 , (9.32) Gµ is the smallest vector in a periodically repeated cubic cell in the direction µ. Eq.(9.31) can be generalized to cells of arbitrary shape. This formula is in principle valid in the limit of an infinite dimension of the cell, but in a non conducting material this is a good approximation even with relatively small supercells. P eleµ is then given by: P eleµ = 2 |e| |Gµ| ϕµ. (9.33) This induces a perturbation in the Kohn–Sham functional of the type: λEpert [{|Φi〉}] = −E ·Pele = − ∑ ν Eν 2 |e| |Gν | Im log det 〈 Φi ∣∣∣eiGν ·r∣∣∣Φj〉 . (9.34) In this case the perturbative parameter is Eν and |Φi〉 ' |ϕ(0)i 〉−Eν |Φν (1) i 〉 The derivative δEpert /δ〈ϕ(0)i | and its ket conjugate can be evaluated using the formula for the derivative of a matrix A with respect to a generic variable x: d dx ln detA = Tr dAij dx A−1ji . (9.35) The perturbative term becomes: 2 |e| |Gν | Im ∑ i,j (〈 Φν (1) i ∣∣∣eiGν ·r∣∣∣Φ(0)j 〉+ 〈Φ(0)i ∣∣∣eiGν ·r∣∣∣φν(1)j 〉)Q(ν)−1j,i  . (9.36) Using this perturbative term we can calculate the first order correction to the wavefunc- tions {Φ(1)i }. This allows to evaluate the induced dipole moment in the µ direction: δP eleµ = 2 |e| |Gµ| δΦµ = − ∑ ν 2 |e| |Gν | Im ∑ i,j (〈 Φν (1) i ∣∣∣eiGµ·r∣∣∣Φ(0)j 〉+ 〈Φ(0)i ∣∣∣eiGµ·r∣∣∣Φν(1)j 〉)Q(ν)−1j,i Eν(9.37) and the polarizability αµ,ν = −∂Pµ/∂Eν : αµ,ν = 2 |e| |Gν | Im ∑ i,j (〈 Φν (1) i ∣∣∣eiGµ·r∣∣∣Φ(0)j 〉+ 〈Φ(0)i ∣∣∣eiGµ·r∣∣∣Φν(1)j 〉)Q(ν)−1j,i  . (9.38) 96