# Computational Methods for Ab Initio Molecular Dynamics

Authors: Eric Paquet; Herna L. Viktor
Journal: Advances in Chemistry
 (2018)

---

## Abstract

Ab initio molecular dynamics is an irreplaceable technique for the realistic simulation of complex molecular systems and processes from first principles. This paper proposes a comprehensive and self-contained review of ab initio molecular dynamics from a computational perspective and from first principles. Quantum mechanics is presented from a molecular dynamics perspective. Various approximations and formulations are proposed, including the Ehrenfest, Born–Oppenheimer, and Hartree–Fock molecular dynamics. Subsequently, the Kohn–Sham formulation of molecular dynamics is introduced as well as the afferent concept of density functional. As a result, Car–Parrinello molecular dynamics is discussed, together with its extension to isothermal and isobaric processes. Car–Parrinello molecular dynamics is then reformulated in terms of path integrals. Finally, some implementation issues are analysed, namely, the pseudopotential, the orbital functional basis, and hybrid molecular dynamics.


---

## Body

## 1. Introduction

Ab initio molecular dynamics (AIMD) is an irreplaceable technique for the realistic simulation of complex molecular systems and processes associated with biological organisms [1, 2] such as monoclonal antibodies as illustrated in Figure 1. With ab initio molecular dynamics, it is possible to predict in silico phenomena for which an in vivo experiment is either too difficult or expensive, or even currently deemed infeasible [3, 4]. Ab initio molecular dynamics essentially differs from molecular dynamics (MD) in two ways. Firstly, AIMD is based on the quantum Schrödinger equation while its classical counterpart relies on Newton’s equation. Secondly, MD relies on semiempirical effective potentials which approximate quantum effects, while AIMD is based on the real physical potentials.Figure 1
Monoclonal antibody 1F1 isolated from a 1918 influenza pandemic (Spanish Flu) survivor.In this paper, we present a review of ab initio molecular dynamics from a computational perspective and from first principles. Our main aim is to create a document which is self-contained, coherent, and concise but still as complete as possible. As the theoretical details are essential for real-life implementation, we have provided the equations for the relevant physical aspect and approximations presented. Most important, we expose how these equations and concepts are related to one another.The paper is organised as follows. In Section2, quantum mechanics is reviewed from a molecular dynamics perspective. Two important approximations are introduced, namely, the adiabatic and the Born–Oppenheimer approximations. Subsequently, in Section 3, the Ehrenfest molecular dynamics is presented, which allows for a semiclassical treatment of the nuclei. This opens the door, in Section 4, to the Born–Oppenheimer molecular dynamics formulation. This is followed, in Section 5, by the Hartree–Fock formulation which approximates the atomic antisymmetric wave function by a single determinant of the orbitals. Sections 6 and 7 introduce the Kohn–Sham formulation whose primary objective is to replace the interacting electrons by a fictitious, but equivalent, system of noninteracting particles. Electrons are reintroduced as dynamical quantities, in Section 8, through the Car–Parrinello formulation. Isothermal and isobaric processes are addressed in this section. In Section 9, a path integral formulation of molecular dynamics is presented, which makes allowance for a quantum treatment of the nuclear degrees of freedom. Finally, in Section 10, some important computational issues are addressed such as a simplification of the electronic interaction with the pseudopotential, the representation of orbitals in terms of a functional basis, the use of the Fourier and wavelet transform in order to reduce the computational complexity, and the simulation of larger systems with hybrid molecular dynamics. This is followed by a conclusion in Section 11.


## 2. Quantum Molecular Dynamics and the Schrödinger Equation: A Molecular Perspective

We review some important notions of quantum mechanics from a molecular perspective. In quantum mechanics [5–7], the Hamiltonian [5–7] is the operator corresponding to the total energy of the molecular system associated with the electrons and the nuclei. The total Hamiltonian is the sum of the kinetic energies of all the particles plus the potential energy of the constituent particles; in occurrence the electrons and the nuclei(1)Hr,R=-∑Iħ22MIΔI+Her,R,where ħ is the Planck constant and MI is the mass of a given nucleus. The electronic and nuclear Cartesian coordinates for a given particle are defined, respectively, as(2)ri≙ri,x,ri,y,ri,zT,RI≙RI,x,RI,y,RI,zT.The ensemble of all electronic coordinates is represented by r while the ensemble of all nuclear coordinates is denoted by R. In addition, the gradient and the Laplacian operators for a given electron i and nucleus I are defined, respectively, as(3)∇i≙∂∂ri,Δi≙∂∂ri·∂∂ri,∇I≙∂∂RI,ΔI≙∂∂RI·∂∂RI.The gradient operator, a vector, is related to the momentum while the Laplacian operator, a scalar quantity, is associated with the kinetic energy [5–7]. The electronic Hamiltonian is defined as the sum of the kinetic energy of the electrons, the potential energy associated with each pair of electrons, the potential energy associated with each pair of electron-nucleus, and the potential energy concomitant to each pair of nuclei:(4)Her,R≙-∑iħ22meΔi+14πε0∑i<je2ri-rj-14πε0∑I,ie2ZIRI-rj+14πε0∑I<Je2ZIZJRI-RJ,where me is the electronic mass, e is the electronic charge, Z is the atomic number (the number of protons in a given nucleus), and ε0 is the vacuum permittivity. The quantum molecular system is characterised by a wave function Φr,R;t whose evolution is determined by the time-dependent Schrödinger equation [5–7]:(5)iħ∂∂tΦr,R;t=Hr,RΦr,R;t,where i=-1. The Schrödinger equation admits a stationary or time-independent solution:(6)Her,RΨkr,R=EkRΨkr,R,where EkR is the energy associated with the electronic wave functionΨkr,R and k is a set of quantum numbers that labelled the Eigenstates or wave functions as well as the Eigenvalues or energies associated with the stationary Schrödinger equation. The total wave function Φr,R;t may be expended in terms of time-dependent nuclear wave functions χkR;t and stationary electronic wave functions Ψkr,R:(7)Φr,R;t=∑k=0∞χkR;tΨkr,R.It should be noticed that this expansion is exact and does not involve any approximation. The electronic wave functions are solution of the stationary Schrödinger equation, an Eigen equation which involves the electronic Hamiltonian. If one substitutes this expansion in the time-dependent Schrödinger equation, one obtains a system of equations for the evolution of the nuclear wave functions:(8)iħ∂χk∂t=-∑Iħ22MIΔI+EkRχk+∑lCklχl.As may be seen, the nuclear wave functions are coupled to the electronic wave functions. The strength of the coupling is determined by the coupling coefficients which form a matrix:(9)Ckl=∫drΨk∗-∑Iħ22MIΔIΨl+∑I1MI∫drΨk∗-iħ∇IΨl-iħ∇I.This system of equations is too complex to be solved directly. Consequently, in the next subsection, we consider two important approximations of the Schrödinger equation, namely, the adiabatic and the Born–Oppenheimer approximations, which aim to reduce such a complexity. These approximations, as well as those that later follow, reduce substantially the duration of the calculations allowing for larger molecular systems to be simulated and longer time-scales to be explored [8, 9].


### 2.1. Adiabatic and Born–Oppenheimer Approximations

The above-mentioned system of differential equations is too complex to be solved directly. Some approximations must be performed in order to reduce the computational complexity, while maintaining the predictive accuracy of the calculations. The first approximation, which is called the adiabatic approximation [5–7], assumes that the electronic wave functions adapt quasi-instantaneously to a variation of the nuclear configuration. This approximation is justified by the fact that the electrons are much lighter than the nucleus. Consequently, such an approximation is valid, unless the motion of the electron becomes relativistic, which may happen if the electrons are too close to the nucleus. Mathematically, this approximation implies that(10)∇IΨl=0.As a result, the coupling matrix becomes diagonal:(11)Ckl=-∑Iħ22MI∫drΨk∗ΔIΨkδkl,where δkl is the well-known Kronecker delta.The second approximation, which is called the Born–Oppenheimer approximation [6, 10, 11], assumes that the electronic and the nuclear motions are separable as a result of the difference between nuclear and electronic masses. As a result, the expansion reduces to a single Eigen state:(12)Φr,R;t≈Ψkr,RχkR;t.In addition, as the energy associated with the wave function is usually much larger than the corresponding coupling constant:(13)Ckk≪Ek,the time-dependent nuclear Schrödinger further reduces to(14)iħ∂χk∂t=-∑Iħ22MIΔI+EkRχk.Furthermore, it is assumed that the electronic wave function is in its ground or nonexcited state:(15)k=0⟹HeΨ0=E0Ψ0.This occurs if the thermal energy is lower than the energy difference between the ground state and the first excited state, that is, if the temperature is low.In the next section, we introduce another approximation, in which the motion of heavy nuclei is described by a semiclassical equation.


## 2.1. Adiabatic and Born–Oppenheimer Approximations

The above-mentioned system of differential equations is too complex to be solved directly. Some approximations must be performed in order to reduce the computational complexity, while maintaining the predictive accuracy of the calculations. The first approximation, which is called the adiabatic approximation [5–7], assumes that the electronic wave functions adapt quasi-instantaneously to a variation of the nuclear configuration. This approximation is justified by the fact that the electrons are much lighter than the nucleus. Consequently, such an approximation is valid, unless the motion of the electron becomes relativistic, which may happen if the electrons are too close to the nucleus. Mathematically, this approximation implies that(10)∇IΨl=0.As a result, the coupling matrix becomes diagonal:(11)Ckl=-∑Iħ22MI∫drΨk∗ΔIΨkδkl,where δkl is the well-known Kronecker delta.The second approximation, which is called the Born–Oppenheimer approximation [6, 10, 11], assumes that the electronic and the nuclear motions are separable as a result of the difference between nuclear and electronic masses. As a result, the expansion reduces to a single Eigen state:(12)Φr,R;t≈Ψkr,RχkR;t.In addition, as the energy associated with the wave function is usually much larger than the corresponding coupling constant:(13)Ckk≪Ek,the time-dependent nuclear Schrödinger further reduces to(14)iħ∂χk∂t=-∑Iħ22MIΔI+EkRχk.Furthermore, it is assumed that the electronic wave function is in its ground or nonexcited state:(15)k=0⟹HeΨ0=E0Ψ0.This occurs if the thermal energy is lower than the energy difference between the ground state and the first excited state, that is, if the temperature is low.In the next section, we introduce another approximation, in which the motion of heavy nuclei is described by a semiclassical equation.


## 3. Ehrenfest Molecular Dynamics

Another possible approximation is to assume that the nuclear motion is semiclassical as it is the case when the nuclei are relatively heavy. This implies that, instead of being determined by the Schrödinger equation, theaverage nuclear motion is determined by Newton’s equation. The right form for this equation is given by the Ehrenfest theorem [6, 12–14] which states that the potential, which governs the classical motion of the nucleons, is equal to the expectation or average, in the quantum sense, of the electronic Hamiltonian with respect to the electronic wave function:(16)MIR¨It=-∇IΨHeΨ=-∇IHeΨ,where the bra–ket notation is to be understood as(17)ΨOΨ≡∫drΨ∗rOΨrfor any operator O and where R¨I≡d2RI/dt2. This equation may be further simplified if the adiabatic and the Born–Oppenheimer approximations are enforced:(18)∇IΨHeΨ⟹∇IΨ0HeΨ0⟹Ψ0∇IHeΨ0.In the specific context of Ehrenfest molecular dynamics, the electrons follow a time-dependent Schrödinger equation:(19)iħ∂∂tΨ=HeΨ.The fact that this equation is time-dependent ensures that the orthogonality of the orbitals is maintained at all times. The motivation is to be found in the fact that the electronic Hamiltonian is a Hermitian operator [6]. As stated earlier, this is not the case in the Born–Oppenheimer approximation in which the electronic wave function is governed by the stationary Schrödinger equation which does not maintain the orthonormality over time.In the next section, we apply the adiabatic and Born–Oppenheimer approximation in the context of ab initio molecular dynamics.


## 4. Born–Oppenheimer Molecular Dynamics

In Born–Oppenheimer molecular dynamics [6, 10, 11], it is assumed that the adiabatic and the Born–Oppenheimer approximations are valid and that the nuclei follow a semiclassical Newton equation whose potential is determined by the Ehrenfest theorem. It is further assumed that the electronic wave function is in its ground state (lowest energy). Since the stationary Schrödinger equation is used for the electronic degrees of freedom, the orthogonality condition must be enforced with a Lagrange multiplier as this condition is not preserved by the stationary Schrödinger equation over time. The orthogonality of the orbitals ψi is a physical requirement which must be enforced at all time in order to be able to compute real observable quantities. The Born–Oppenheimer molecular dynamics is characterised by a Lagrangian [6] which is defined as the difference between the kinetic and the potential energy:(20)LBO=12∑IMIR˙I2-Ψ0HeΨ0+∑i,jΛijψi∣ψj-δij.The first term of the Lagrangian corresponds to the kinetic energy of the nuclei, the second term corresponds to the nuclear potential as obtained from the Ehrenfest theorem, and the last term is a Lagrange function which ensures that the orbitals remain orthonormal at all time. The Lagrange multipliers are denoted by Λij. Since the electrons follow a Fermi–Dirac statistics [6], they obey the Pauli Exclusion Principle (two electrons cannot be in the same quantum state) which means that the electronic wave function must be an antisymmetric function of its orbitals, namely, the wave functions associated with the individual electrons.The equations of motion associated with this Lagrangian are obtained from the Euler–Lagrange equations [6]. There is one system of Euler–Lagrange equations for the nuclear degrees of freedom:(21)ddt∂L∂R˙I-∂L∂RI=0,and one system of Euler–Lagrange equations for their electronic counterpart:(22)ddtδLδψ˙i-δLδψi=0⟹ddtδLδψ˙i∗r-δLδψi∗r=0.From the Euler–Lagrange equations, one obtains(23)MIR¨It=-∇Iminψ⁡Ψ0HeΨ0for the nuclear equations of motion, and(24)Heψir=∑jΛijψirfor the electronic equations of motion. One should notice the presence of the Lagrange multiplier in the electronic equations of motion which ensures that the orbital remains orthonormal at all time.In the next section, we introduce another approximation in order to further reduce the complexity of the electronic Hamiltonian.


## 5. Hartree–Fock Molecular Dynamics

From now on, we shall adopt the atomic unit system:(25)me=e=ħ=14πε0=1in order to alleviate the notation. In Hartree–Fock molecular dynamics [15], the antisymmetric electronic wave function is approximated by a single determinant of the electronic orbitals ψiri:(26)Ψ=1N!det⁡ψ1r1ψ1r2⋯ψ1rNψ2r1ψ2r2⋯ψ2rN⋮⋮⋱⋮ψNr1ψNrN⋯ψNrN.Two operators are associated with the electronic interaction. The first one, the Coulomb operator, corresponds to the electrostatic interaction between the orbitals:(27)Jjrψir≡∫dr′ψj∗r′1r-r′ψjr′ψir.The second operator, the exchange operator [6, 15], is associated with the exchange energy, a quantum mechanical effect that occurs as a result of the Pauli Exclusion Principle:(28)Kjrψir≡∫dr′ψj∗r′1r-r′ψir′ψjr.From the Coulomb and the exchange operators, an electronic Hamiltonian may be inferred:(29)HeHF=-12Δ+Vextr+2∑jJjr-∑jKjr≡12Δ+VHFr,where the external potential is defined as(30)Vextr≙-∑IZIRI-r+∑I<JZIZJRI-RJ.This Hamiltonian determines, in turn, the electronic Lagrangian:(31)LeHF=-Ψ0HeHFΨ0+∑i,jΛijψi∣ψj-δij.A Lagrange multiplier term is added to the Lagrangian in order to ensure that the orbitals remain orthonormal at all time: a quantum mechanical requirement [6]. From the Euler–Lagrange equations, one obtains the equations of motion for the electrons:(32)HeHFψi=∑jΛijψiwhich differ from (24) by the Hamiltonian.In the next subsection, we seek to replace the interacting electrons by a fictitious but equivalent system of noninteracting particles in order to further reduce the computational complexity.


## 6. Kohn–Sham Molecular Dynamics

The Kohn–Sham formulation [16, 17] aims to replace the interacting electrons by a fictitious system of noninteracting particles that generate the same electronic density as the physical system of interacting particles. The electronic density is defined as(33)n≙∑ifiψi∣ψi,where fi is the occupation number of a given orbital.In this formulation of AIMD, the Ehrenfest term is replaced by the Kohn–Sham energy:(34)minΨ0⁡Ψ0HeΨ0⁡=minψEKSψ.The latter is define as(35)EKSn≙-12∑i=1NfiψiΔψi+∫drVextrnr+12∫drVHrnr+EXCn,where the first term is the total electronic kinetic energy associated with the electrons, the second term is the electrostatics interaction energy between the electronic density and the external potential, and the third term is the self-electrostatic interaction energy associated with the electronic density. The latter involves the interaction of the electronic density with it self-created electrostatic potential. This potential, called the Hartree potential, is obtained by solving the Poisson equation [18]:(36)ΔVH=-4πn⟹VHr=∫dr′nr′r-r′.The last term is the celebrated exchange-correlation energy or density functional [16, 19–21] which takes into account the residual electronic interaction, that is, the self-interaction of the electrons. Unfortunately, the density functional has no closed-form solution but many approximations are known. Some of these approximations are considered in the next section.From the exchange-correlation energy, one may define the exchange-correlation potential:(37)VXCr≙δEXCnδnrwhich is the functional derivative [22] of the exchange-correlation energy with respect to the electronic density. In addition, one may define the Kohn–Sham potential:(38)VKS≙VH+Vext+VXC.As for the other approaches, a Lagrangian is defined in which the orthonormality conditions are enforced with Lagrange multipliers. The electronic equations of motion, which are determined by the Euler–Lagrange equations are(39)HeKSψi=∑jΛijψj,where the fictitious one-particle Hamiltonian is given by(40)HeKS≡-12Δ+VKS.A canonical form may be obtained if a unitary transformation is applied to the previous equation:(41)HeKSψi=εiψi.In the next section, we introduce some density functionals in order to replace the interacting electrons by an equivalent but yet simpler system of noninteracting particles.


## 7. Exchange-Correlation Energy

The detailed analysis of density functionals, also known as exchange-correlation energies [23], is beyond the scope of this review. We refer the interested reader to [21] for an exhaustive analysis. In this section, we briefly introduce a few common density functionals. As mentioned earlier, the aim of the density functional is to express the electronic interaction in terms of the sole electronic density.In the simplest case, the exchange-correlation energy is a functional of the electronic density alone for which the most important representative is the local density approximation:(42)EXCLDAn=-343π1/3∫drn4/3.One may take into consideration, in addition to the electronic density, the gradient of the electronic density, in which case the approach is called the generalised gradient approximation:(43)EXCGGAn=∫drn4/3FXCζ,where the dimensionless reduced gradient is defined as(44)ζ≙∇nn4/3≡23π21/3s.Among these approximations is the B88 approximation:(45)FXCB88ζ=346π1/3+βζσ21+6βζσsinh-1⁡ζσ,where σ refers to the spin and the Perdew, Burke, and Ernzerhof (PBE) approximation:(46)FXCPBEs=343π1/3+μs21+μs2κ-1.The exchange energy, which is associated with the Pauli Exclusion Principle, may be characterised by the Hartree–Fock energy:(47)EXCHF=-12∑i,j∫drdr′ψi∗rψjrψj∗r′ψir′r-r′which was introduced earlier. These density functionals may be linearly combined in order to increase their precision such as in the case of the B3 approximation:(48)EXCB3=aEXCHF+1-aEXCLDA+bΔEXCB88+cEXCGGA+1-cEXCLDA,where a, b, and c are empirical parameters. Obviously, these parameters are application dependent.Recently, it has been proposed to evaluate the functional density with machine learning techniques. The functional density is learned by examples instead of directly solving the Kohn–Sham equations [17, 24, 25]. As a result, substantially less time is required to complete the calculations allowing for larger system to be simulated and longer time-scales to be explored.In the next section, we further improve the precision of the calculations by reintroducing the dynamic electronic degrees of freedom which, until now, have been absent from the equations of motion.


## 8. Car–Parrinello Molecular Dynamics

### 8.1. Equations of Motion

The Kohn–Sham dynamics, as formulated in the previous sections, does not take the dynamics of the electrons into account despite the fact that it is present: an unattractive unphysical feature. Indeed, one obtains from the Lagrangian and the Euler–Lagrange equations(49)ddt∂L∂ψ˙i∗r≡0which means the orbitals are not dynamical fields of the molecular system. Nevertheless, the Lagrangian may be modified in order to also include the dynamic nature of electronic degrees of freedom through the introduction of a fictitious electronic kinetic energy term. The extended Lagrangian, which is called the Car–Parrinello Lagrangian [11, 26, 27](50)LCP=12∑IMIR˙I2+∑iμψ˙i∣ψ˙i+Ψ0HeKSΨ0+∑i,jΛijψi∣ψj-δijdiffers from the original Kohn–Sham Lagrangian by the second term which associates a fictitious kinetic energy to the electronic orbitals. The parameter μ acts as a fictitious electronic mass or inertia. As in the previous sections, the nuclear equations of motion(51)MIR¨I=-∇IΨ0HeKSΨ0+∑i,jΛij∇Iψi∣ψjand the electronic equations of motion(52)μψ¨it=-HeKSψit+∑jΛijψjtmay be inferred from the Euler–Lagrange equations. These equations involve a second order time derivative of the orbitals which means that the latter are now proper dynamical quantities.In the next two subsections, we consider ab initio molecular dynamics at constant temperature and at constant pressure.


### 8.2. Massive Thermostating and Isothermal Processes

In the previous sections, we have implicitly assumed that the molecular system under consideration was isolated. Of course, this is not compatible with most experimental conditions [27–29] in which the system is either kept at constant temperature (isothermal) in a heat bath or at constant pressure (isobaric) as a consequence, for instance, of the atmospheric pressure. In this subsection and in the next, we explain how to address these important issues.As we have recently explained in a computational review on molecular dynamics [30], an isothermal process cannot be obtained by adding an extra term to the Lagrangian. Rather, the equations of motions must be modified directly as a result of the non-Hamiltonian nature of the isothermal process [31]. In order to obtain an isothermal molecular system, a friction term must be recursively added to each degree of freedom [28, 29]. Therefore, the electronic equations of motion must be modified as follows:(53)μψ¨i=-HeKSψi+∑jΛijψj-μη˙1ψ˙i,Qe1η¨1=2μ∑iψ˙i∣ψ˙i-Ne2βe-Qe1η˙1η˙2,Qelη¨l=Qel-1η˙l-12-1βe-Qelη˙lη˙l+11-δlL,l=2,…,Lwhile the nuclear equations of motion must take a similar form:(54)MIR¨I=-∇IEKS-MIξ˙1R˙I,Qn1ξ¨1=∑IMIR˙I2-gβ-Qn1ξ˙1ξ˙2,Qnkξ¨k=Qnk-1ξ˙k-12-1β-Qnkξ˙kξ˙k+11-δkKk=2,…,K,where β≡1/kBT, T is the temperature of the heat bath, kB is the Boltzmann constant, the ηl are frictional electronic degrees of freedom, the Qel are the friction coefficients associated with the electronic orbitals, the ξk are frictional nuclear degrees of freedom, the Qnk are the friction coefficients associated with the nuclei, and g denotes the number of dynamical degrees of freedom to which the nuclear thermostat chain is coupled. The friction terms are proportional to the velocity of the corresponding degrees of freedom as it is customary. Not only must friction terms be added to the physical degrees of freedom, but each nonphysical friction term, in turn, must be thermalised by another nonphysical friction term until an isothermal state is properly achieved. The terms η˙1 and ξ˙1 may be considered as dynamical friction coefficients for the physical degrees of freedom.In the next section, we present an alternative approach based on the generalised Langevin equation.


### 8.3. Generalised Langevin Equation and Isothermal Processes

Massive thermostating is not the sole approach to simulate an isothermal process. Indeed, the latter may be achieved by means of the generalised Langevin equation [32–34]. In this subsection, we restrict ourselves to one generalised coordinate in order not to clutter the notation. The generalisation to more than one coordinate is immediate. The generalised Langevin equation, which is a differential stochastic equation, may be written as(55)q˙=pp˙p˙=-V′q0-Appp+Bpξ,where Ap is the drift matrix:(56)Ap≙appapTa-pA,Bp is the diffusion matrix, q is a generalised coordinate associated with an atom or a nucleus (position), p is the corresponding generalised momentum, and p is a set of np hidden nonphysical momenta. The structure of the diffusion matrix Bp is similar to the one of the corresponding drift matrix. The random process is characterised by an uncorrelated Gaussian noise:(57)ξitξi0=δijδt.The hidden momenta, associated with the generalised Langevin equation, may be marginalised because the evolution of the variables p,pT, in the free particle limit, is described by a linear Markovian stochastic differential equation. As a result, the generalised Langevin equation becomes equivalent to a Langevin equation with memory kernel and noise correlation [32–34]:(58)q˙=p,p˙=-V′q-∫-∞tKt-spsds+ζ.The memory kernel, which is a dissipative term associated with friction, is given by the expression(59)Kt=2appδt-apTexp-tAa-pwhereas the noise correlation, which characterised the fluctuations associated with the random noise, is provided by(60)Ht≙ζtζ0=dppδt+apTexp⁡-tAZap-dp,where the matrixes Z and Dp are defined as(61)Z≙∫0∞e-AtDe-ATtdt,(62)Dp≙BpBpT,respectively. The canonical, isothermal ensemble is obtained by applying the fluctuation-dissipation theorem [35]. The fluctuation-dissipation theorem states that the Fourier transforms F of the noise correlation and the memory kernel must be related by(63)F∘Hω=kBTF∘Kω.The fluctuation-dissipation theorem implies in turn that(64)Dp=BpBpT=kBTAp+ApT.Therefore, the fluctuation-dissipation theorem fixes the diffusion matrix once the drift matrix is determined. As a result, an isothermal process follows immediately.In the next subsection, we address isobaric molecular processes, which are very common as most experiments are performed at atmospheric pressure.


### 8.4. Isobaric Processes

As opposed to the isothermal process, the isobaric process is a Hamiltonian process which means that it may be obtained by adding extra terms to the Lagrangian [28, 36].In order to model an isobaric molecular process, the volume occupied by the molecular system is divided into congruent parallelepipeds called Bravais cells. These cells are characterised by three oriented vectors which correspond to the three primitive edges spanning their volume. These vectors are concatenated in order to form a Bravais matrix:(65)h=a1∣a2∣a3.The Bravais cell is characterised by its volume and by a local metric:(66)Ω=det⁡h,G=hTh.In order to introduce the volume as a dynamic quantity into the Lagrangian, the nuclear and the electronic coordinates are reformulated in terms of the Bravais matrix and of the scaled coordinates:(67)RI=hSI,ri=hsi.The Bravais matrix and the scale coordinates constitute distinct degrees of freedom. In order to obtain an isobaric system, the Car–Parrinello Lagrangian must be supplemented with three additional terms which appear on the third line of the Lagrangian:(68)L=∑iμψ˙ihs∣ψ˙ihs+EKSψ,hS+∑i,jΛijψihs∣ψjhs-δij+12∑IMIS˙ITGS˙I+12Wtrh˙Th˙-pΩ.The first term represents the kinetic energy associated with the scale coordinates. One should notice the presence of the metric or quadratic form in the inner product. The second term represents the kinetic energy associated with the Bravais matrix. The matrix norm is the square of the Frobenius norm while W is a fictitious mass or inertia attributed to the Bravais cells. The last term is associated with the pressure pof the barostat (ambient pressure). The equations of motion are obtained, as usual, from the Euler–Lagrange equations. In particular, the Euler–Lagrange equations for the Bravais cells take the form:(69)ddt∂L∂h˙uv-∂L∂huv=0.In the next section, we present a path integral formulation of ab initio molecular dynamics which allows a quantum formulation of both the electronic and nuclear degrees of freedom.


## 8.1. Equations of Motion

The Kohn–Sham dynamics, as formulated in the previous sections, does not take the dynamics of the electrons into account despite the fact that it is present: an unattractive unphysical feature. Indeed, one obtains from the Lagrangian and the Euler–Lagrange equations(49)ddt∂L∂ψ˙i∗r≡0which means the orbitals are not dynamical fields of the molecular system. Nevertheless, the Lagrangian may be modified in order to also include the dynamic nature of electronic degrees of freedom through the introduction of a fictitious electronic kinetic energy term. The extended Lagrangian, which is called the Car–Parrinello Lagrangian [11, 26, 27](50)LCP=12∑IMIR˙I2+∑iμψ˙i∣ψ˙i+Ψ0HeKSΨ0+∑i,jΛijψi∣ψj-δijdiffers from the original Kohn–Sham Lagrangian by the second term which associates a fictitious kinetic energy to the electronic orbitals. The parameter μ acts as a fictitious electronic mass or inertia. As in the previous sections, the nuclear equations of motion(51)MIR¨I=-∇IΨ0HeKSΨ0+∑i,jΛij∇Iψi∣ψjand the electronic equations of motion(52)μψ¨it=-HeKSψit+∑jΛijψjtmay be inferred from the Euler–Lagrange equations. These equations involve a second order time derivative of the orbitals which means that the latter are now proper dynamical quantities.In the next two subsections, we consider ab initio molecular dynamics at constant temperature and at constant pressure.


## 8.2. Massive Thermostating and Isothermal Processes

In the previous sections, we have implicitly assumed that the molecular system under consideration was isolated. Of course, this is not compatible with most experimental conditions [27–29] in which the system is either kept at constant temperature (isothermal) in a heat bath or at constant pressure (isobaric) as a consequence, for instance, of the atmospheric pressure. In this subsection and in the next, we explain how to address these important issues.As we have recently explained in a computational review on molecular dynamics [30], an isothermal process cannot be obtained by adding an extra term to the Lagrangian. Rather, the equations of motions must be modified directly as a result of the non-Hamiltonian nature of the isothermal process [31]. In order to obtain an isothermal molecular system, a friction term must be recursively added to each degree of freedom [28, 29]. Therefore, the electronic equations of motion must be modified as follows:(53)μψ¨i=-HeKSψi+∑jΛijψj-μη˙1ψ˙i,Qe1η¨1=2μ∑iψ˙i∣ψ˙i-Ne2βe-Qe1η˙1η˙2,Qelη¨l=Qel-1η˙l-12-1βe-Qelη˙lη˙l+11-δlL,l=2,…,Lwhile the nuclear equations of motion must take a similar form:(54)MIR¨I=-∇IEKS-MIξ˙1R˙I,Qn1ξ¨1=∑IMIR˙I2-gβ-Qn1ξ˙1ξ˙2,Qnkξ¨k=Qnk-1ξ˙k-12-1β-Qnkξ˙kξ˙k+11-δkKk=2,…,K,where β≡1/kBT, T is the temperature of the heat bath, kB is the Boltzmann constant, the ηl are frictional electronic degrees of freedom, the Qel are the friction coefficients associated with the electronic orbitals, the ξk are frictional nuclear degrees of freedom, the Qnk are the friction coefficients associated with the nuclei, and g denotes the number of dynamical degrees of freedom to which the nuclear thermostat chain is coupled. The friction terms are proportional to the velocity of the corresponding degrees of freedom as it is customary. Not only must friction terms be added to the physical degrees of freedom, but each nonphysical friction term, in turn, must be thermalised by another nonphysical friction term until an isothermal state is properly achieved. The terms η˙1 and ξ˙1 may be considered as dynamical friction coefficients for the physical degrees of freedom.In the next section, we present an alternative approach based on the generalised Langevin equation.


## 8.3. Generalised Langevin Equation and Isothermal Processes

Massive thermostating is not the sole approach to simulate an isothermal process. Indeed, the latter may be achieved by means of the generalised Langevin equation [32–34]. In this subsection, we restrict ourselves to one generalised coordinate in order not to clutter the notation. The generalisation to more than one coordinate is immediate. The generalised Langevin equation, which is a differential stochastic equation, may be written as(55)q˙=pp˙p˙=-V′q0-Appp+Bpξ,where Ap is the drift matrix:(56)Ap≙appapTa-pA,Bp is the diffusion matrix, q is a generalised coordinate associated with an atom or a nucleus (position), p is the corresponding generalised momentum, and p is a set of np hidden nonphysical momenta. The structure of the diffusion matrix Bp is similar to the one of the corresponding drift matrix. The random process is characterised by an uncorrelated Gaussian noise:(57)ξitξi0=δijδt.The hidden momenta, associated with the generalised Langevin equation, may be marginalised because the evolution of the variables p,pT, in the free particle limit, is described by a linear Markovian stochastic differential equation. As a result, the generalised Langevin equation becomes equivalent to a Langevin equation with memory kernel and noise correlation [32–34]:(58)q˙=p,p˙=-V′q-∫-∞tKt-spsds+ζ.The memory kernel, which is a dissipative term associated with friction, is given by the expression(59)Kt=2appδt-apTexp-tAa-pwhereas the noise correlation, which characterised the fluctuations associated with the random noise, is provided by(60)Ht≙ζtζ0=dppδt+apTexp⁡-tAZap-dp,where the matrixes Z and Dp are defined as(61)Z≙∫0∞e-AtDe-ATtdt,(62)Dp≙BpBpT,respectively. The canonical, isothermal ensemble is obtained by applying the fluctuation-dissipation theorem [35]. The fluctuation-dissipation theorem states that the Fourier transforms F of the noise correlation and the memory kernel must be related by(63)F∘Hω=kBTF∘Kω.The fluctuation-dissipation theorem implies in turn that(64)Dp=BpBpT=kBTAp+ApT.Therefore, the fluctuation-dissipation theorem fixes the diffusion matrix once the drift matrix is determined. As a result, an isothermal process follows immediately.In the next subsection, we address isobaric molecular processes, which are very common as most experiments are performed at atmospheric pressure.


## 8.4. Isobaric Processes

As opposed to the isothermal process, the isobaric process is a Hamiltonian process which means that it may be obtained by adding extra terms to the Lagrangian [28, 36].In order to model an isobaric molecular process, the volume occupied by the molecular system is divided into congruent parallelepipeds called Bravais cells. These cells are characterised by three oriented vectors which correspond to the three primitive edges spanning their volume. These vectors are concatenated in order to form a Bravais matrix:(65)h=a1∣a2∣a3.The Bravais cell is characterised by its volume and by a local metric:(66)Ω=det⁡h,G=hTh.In order to introduce the volume as a dynamic quantity into the Lagrangian, the nuclear and the electronic coordinates are reformulated in terms of the Bravais matrix and of the scaled coordinates:(67)RI=hSI,ri=hsi.The Bravais matrix and the scale coordinates constitute distinct degrees of freedom. In order to obtain an isobaric system, the Car–Parrinello Lagrangian must be supplemented with three additional terms which appear on the third line of the Lagrangian:(68)L=∑iμψ˙ihs∣ψ˙ihs+EKSψ,hS+∑i,jΛijψihs∣ψjhs-δij+12∑IMIS˙ITGS˙I+12Wtrh˙Th˙-pΩ.The first term represents the kinetic energy associated with the scale coordinates. One should notice the presence of the metric or quadratic form in the inner product. The second term represents the kinetic energy associated with the Bravais matrix. The matrix norm is the square of the Frobenius norm while W is a fictitious mass or inertia attributed to the Bravais cells. The last term is associated with the pressure pof the barostat (ambient pressure). The equations of motion are obtained, as usual, from the Euler–Lagrange equations. In particular, the Euler–Lagrange equations for the Bravais cells take the form:(69)ddt∂L∂h˙uv-∂L∂huv=0.In the next section, we present a path integral formulation of ab initio molecular dynamics which allows a quantum formulation of both the electronic and nuclear degrees of freedom.


## 9. Path Integral Molecular Dynamics: Toward a Quantum Formulation of Nuclei

The ab initio path integral technique is based on the formulation of quantum statistical mechanics in terms of Feynman path integrals. Contrarily to the previous approaches, the path integral method allows for a quantum formulation which includes, in addition to the electronic degrees of freedom, their nuclear counterpart. Such a formulation is essential for systems containing light nuclei [37, 38].


### 9.1. Path Integral Formulation

The quantum path integral formulation [36, 39–42] is based on the partition function of the quantum statistical canonical ensemble which is defined as the trace of the exponential of the Hamiltonian operator:(70)Z=trexp⁡-βH.The partition function describes the statistical properties of the molecular system. Since the operators associated with the electrons and the nuclei do not commute, the exponential of the Hamiltonian operator must be evaluated with the Trotter factorization:(71)exp⁡-β-∑Iħ22MIΔI+He=limP→∞exp⁡βP∑Iħ22MIΔIexp⁡-βPHeP.If the completeness relation(72)∫dR∑kΨkR,RR,ΨkR=I,which is an identity operator, is recursively inserted into the Trotter formula and if the adiabatic approximation is performed:(73)∇IΨk=0⟹Ψks+1Rs+1∣ΨksRs≈δks+1ks,the partition function becomes(74)Z=∑klimP→∞⁡∏s=1P∏I=1N∫×exp⁡-β∑s=1P∑I=1N12MIωP2RIs-RIs+12+1PEkRs×∏s=1P∏I=1NAPdRIs.The quantities(75)AP≡MIP2πβħ23/2,ωP2=Pħ2β2are the amplitude and the angular frequency associated with the quantum harmonic oscillators. The latter appear as a consequence of the path integral formulation. The computational time s is a discrete evolution parameter associated with the evolution of the molecular system. As such, it represents a particular time-slice. The path integral associated with the partition function is a weighted sum over all possible nuclear trajectories. The nuclear configuration at a particular time-slice s is provided by the ensemble of all the individual nuclear configurations RIs at this specific instant. The weighting function corresponds to the exponential factor which consists of two terms: the first one is related to the harmonic potential energy associated with the nuclei while the second is the energy associated with the electrons.As in the previous sections, the Born–Oppenheimer approximation is legitimate if the thermal energy is much smaller than the difference between the electronic ground state and the excited states:(76)Ek-E0≫kBT.It follows that the partition function becomes(77)Z=limP→∞⁡∏s=1P∏I=1N∫×exp⁡-β∑s=1P∑I=1N12MIωP2RIs-RIs+12+1PE0Rs×∏s=1P∏I=1NAPdRIs.From this partition function, an extended Lagrangian may be defined:(78)LPI=∑s=1P∑I=1N12MI′PIs2-12MIωP2RIs-RIs+12-1PE0Rs,where MI′ are fictitious masses or unphysical auxiliary parameters associated with the nuclei whereas their physical masses are identified by MI. One should notice that N×P fictitious momenta PIs=MI′R˙Is have been formally introduced. These momenta are required in order to evaluate numerically the path integral with Monte Carlo techniques [30, 43]. The reader should notice that the momenta affect neither the partition function (up to an irrelevant normalisation factor) nor the physical observables. The ground state energy E0Rs must be evaluated concurrently with the nuclear energy for each time-slice.The nuclear coordinates are not linearly independent. Nevertheless, they may become linearly independent if they are expressed in terms of their normal modes. The normal decomposition [39, 41, 44] is obtained by representing each coordinate in terms of complex Fourier series:(79)RIs⟶∑k=1PaIke2πis-1k-1/P.The complex Fourier coefficients are further expanded in terms of their real and imaginary parts:(80)uI1=aI1,uIP=aIP+2/2,uI2s-2=ReaIs,uI2s-1=ImaIs.As for the electronic structure, it may be obtained from one of the previously introduced approaches. For instance, for a Car–Parrinello technique formulated in terms of a Kohn–Sham Hamiltonian, the path integral Lagrangian in normal coordinates becomes(81)LPICP=1P∑s=1P∑iμψ˙is∣ψ˙is-EKSψs,Rus+∑i,jΛijsψis∣ψjs-δij+∑s=1P12∑IMI′su˙Is2-12∑IMIsωP2uIs2,where the normal mode masses are defined as(82)MI1=0,MIs=2P1-cos⁡2πs-1PMI,s=2,…,Pto which the fictitious normal mode masses are closely related:(83)MI′1=MI,MI′s=γMIs,s=2,…,P,where the centroid adiabaticity parameter γ belongs to the interval 0<γ≤1. These masses are functions of the computational time. All the other quantities were defined earlier.In the next subsection, we address the calculation of physical observables in the path integral framework.


### 9.2. Physical Observables

An observable [6, 7, 45] is a dynamic quantity that may be measured experimentally such as the average energy or the heat capacity. Numerous observables may be obtained directly from the partition function, for instance,(i)
the expectation of the energy (average energy):(84)E=-∂ln⁡Z∂β(ii)
the variance of the energy or energy fluctuation:(85)ΔE2=∂2Z∂β2(iii)
the heat or thermal capacity:(86)Cv=1kBTΔE2(iv)
the entropy:(87)S=∂∂TkBTln⁡Z(v)
the Helmholtz free energy:(88)A=-kBTln⁡Zamong others. In addition, if a small perturbation is applied to a molecular system, the expectation of the energy associated with this perturbation is(89)E=E0+λW⟹W=-1β∂∂λlnZ,where λ is a small running parameter.In the next subsection, we present an isothermal formulation of path integrals.


### 9.3. Car–Parrinello Path Integral Molecular Dynamics and Massive Thermostating

Isothermal processes are essential for two reasons. The first one, which is rather evident, is that the simulation of realistic physical conditions often involves their simulation [34]. The second reason must be found in a rather subtle shortcoming of the formalism.The massive thermostating of the path integral [34] is a “sine qua non” condition in order to obtain physically realistic results. Indeed, the harmonic potential leads to inefficient and nonergodic dynamics when microcanonical trajectories are used to generate ensemble averages.In contrast, the thermostats generate ergodic, canonical averages at the expense of introducing sets of auxiliary chain variables which add to the complexity of the calculation, but which are nevertheless required for its precision and physical trustworthiness.As we saw in Section8.2, it is not possible to generate an isothermal process simply by adding extra terms to the Lagrangian. Rather, the equations of motion, which are obtained from the Euler–Lagrange equations, must be modified accordingly. Therefore, friction terms must be added recursively to the various degrees of freedom. As a result, the electronic equations of motion become(90)μψ¨is=-δEKSψs,Rsδψis∗+∑jΛijsψis-μη˙1sψisQe1η¨1s=2μ∑iψis∣ψis-Ne2βe-Qe1η˙1sη˙2s,Qelη¨ls=Qel-1η˙l-1s2-1βe-Qelη˙lsη˙l+1s1-δ2L,l=2,…,L;s=1,…,Pwhereas the nuclear equations of motion transform into(91)MIR¨Is=-1P∂EKSψs,Rs∂RIs-MIωP22RIs-RIs+1-RIs-1-MI′ξ˙1sR˙Is,Qn1ξ¨1s=∑IMI′R˙Is2-gβ-Qn1ξ˙1sξ˙2s,Qnkξ¨ks=Qnk-1ξ˙k-1s2-1β-Qnkξ˙ksξ˙k+1s1-δkK,k=2,…,K;s=1,…,P.The quantities appearing in these equations are all defined in Section 8.2. One should notice that number of equations is quite large due to the evolution parameter swhich is absent from the standard Car–Parrinello formulation as reported in Section 9.3.In the next section, we address some important implementational issues.


## 9.1. Path Integral Formulation

The quantum path integral formulation [36, 39–42] is based on the partition function of the quantum statistical canonical ensemble which is defined as the trace of the exponential of the Hamiltonian operator:(70)Z=trexp⁡-βH.The partition function describes the statistical properties of the molecular system. Since the operators associated with the electrons and the nuclei do not commute, the exponential of the Hamiltonian operator must be evaluated with the Trotter factorization:(71)exp⁡-β-∑Iħ22MIΔI+He=limP→∞exp⁡βP∑Iħ22MIΔIexp⁡-βPHeP.If the completeness relation(72)∫dR∑kΨkR,RR,ΨkR=I,which is an identity operator, is recursively inserted into the Trotter formula and if the adiabatic approximation is performed:(73)∇IΨk=0⟹Ψks+1Rs+1∣ΨksRs≈δks+1ks,the partition function becomes(74)Z=∑klimP→∞⁡∏s=1P∏I=1N∫×exp⁡-β∑s=1P∑I=1N12MIωP2RIs-RIs+12+1PEkRs×∏s=1P∏I=1NAPdRIs.The quantities(75)AP≡MIP2πβħ23/2,ωP2=Pħ2β2are the amplitude and the angular frequency associated with the quantum harmonic oscillators. The latter appear as a consequence of the path integral formulation. The computational time s is a discrete evolution parameter associated with the evolution of the molecular system. As such, it represents a particular time-slice. The path integral associated with the partition function is a weighted sum over all possible nuclear trajectories. The nuclear configuration at a particular time-slice s is provided by the ensemble of all the individual nuclear configurations RIs at this specific instant. The weighting function corresponds to the exponential factor which consists of two terms: the first one is related to the harmonic potential energy associated with the nuclei while the second is the energy associated with the electrons.As in the previous sections, the Born–Oppenheimer approximation is legitimate if the thermal energy is much smaller than the difference between the electronic ground state and the excited states:(76)Ek-E0≫kBT.It follows that the partition function becomes(77)Z=limP→∞⁡∏s=1P∏I=1N∫×exp⁡-β∑s=1P∑I=1N12MIωP2RIs-RIs+12+1PE0Rs×∏s=1P∏I=1NAPdRIs.From this partition function, an extended Lagrangian may be defined:(78)LPI=∑s=1P∑I=1N12MI′PIs2-12MIωP2RIs-RIs+12-1PE0Rs,where MI′ are fictitious masses or unphysical auxiliary parameters associated with the nuclei whereas their physical masses are identified by MI. One should notice that N×P fictitious momenta PIs=MI′R˙Is have been formally introduced. These momenta are required in order to evaluate numerically the path integral with Monte Carlo techniques [30, 43]. The reader should notice that the momenta affect neither the partition function (up to an irrelevant normalisation factor) nor the physical observables. The ground state energy E0Rs must be evaluated concurrently with the nuclear energy for each time-slice.The nuclear coordinates are not linearly independent. Nevertheless, they may become linearly independent if they are expressed in terms of their normal modes. The normal decomposition [39, 41, 44] is obtained by representing each coordinate in terms of complex Fourier series:(79)RIs⟶∑k=1PaIke2πis-1k-1/P.The complex Fourier coefficients are further expanded in terms of their real and imaginary parts:(80)uI1=aI1,uIP=aIP+2/2,uI2s-2=ReaIs,uI2s-1=ImaIs.As for the electronic structure, it may be obtained from one of the previously introduced approaches. For instance, for a Car–Parrinello technique formulated in terms of a Kohn–Sham Hamiltonian, the path integral Lagrangian in normal coordinates becomes(81)LPICP=1P∑s=1P∑iμψ˙is∣ψ˙is-EKSψs,Rus+∑i,jΛijsψis∣ψjs-δij+∑s=1P12∑IMI′su˙Is2-12∑IMIsωP2uIs2,where the normal mode masses are defined as(82)MI1=0,MIs=2P1-cos⁡2πs-1PMI,s=2,…,Pto which the fictitious normal mode masses are closely related:(83)MI′1=MI,MI′s=γMIs,s=2,…,P,where the centroid adiabaticity parameter γ belongs to the interval 0<γ≤1. These masses are functions of the computational time. All the other quantities were defined earlier.In the next subsection, we address the calculation of physical observables in the path integral framework.


## 9.2. Physical Observables

An observable [6, 7, 45] is a dynamic quantity that may be measured experimentally such as the average energy or the heat capacity. Numerous observables may be obtained directly from the partition function, for instance,(i)
the expectation of the energy (average energy):(84)E=-∂ln⁡Z∂β(ii)
the variance of the energy or energy fluctuation:(85)ΔE2=∂2Z∂β2(iii)
the heat or thermal capacity:(86)Cv=1kBTΔE2(iv)
the entropy:(87)S=∂∂TkBTln⁡Z(v)
the Helmholtz free energy:(88)A=-kBTln⁡Zamong others. In addition, if a small perturbation is applied to a molecular system, the expectation of the energy associated with this perturbation is(89)E=E0+λW⟹W=-1β∂∂λlnZ,where λ is a small running parameter.In the next subsection, we present an isothermal formulation of path integrals.


## 9.3. Car–Parrinello Path Integral Molecular Dynamics and Massive Thermostating

Isothermal processes are essential for two reasons. The first one, which is rather evident, is that the simulation of realistic physical conditions often involves their simulation [34]. The second reason must be found in a rather subtle shortcoming of the formalism.The massive thermostating of the path integral [34] is a “sine qua non” condition in order to obtain physically realistic results. Indeed, the harmonic potential leads to inefficient and nonergodic dynamics when microcanonical trajectories are used to generate ensemble averages.In contrast, the thermostats generate ergodic, canonical averages at the expense of introducing sets of auxiliary chain variables which add to the complexity of the calculation, but which are nevertheless required for its precision and physical trustworthiness.As we saw in Section8.2, it is not possible to generate an isothermal process simply by adding extra terms to the Lagrangian. Rather, the equations of motion, which are obtained from the Euler–Lagrange equations, must be modified accordingly. Therefore, friction terms must be added recursively to the various degrees of freedom. As a result, the electronic equations of motion become(90)μψ¨is=-δEKSψs,Rsδψis∗+∑jΛijsψis-μη˙1sψisQe1η¨1s=2μ∑iψis∣ψis-Ne2βe-Qe1η˙1sη˙2s,Qelη¨ls=Qel-1η˙l-1s2-1βe-Qelη˙lsη˙l+1s1-δ2L,l=2,…,L;s=1,…,Pwhereas the nuclear equations of motion transform into(91)MIR¨Is=-1P∂EKSψs,Rs∂RIs-MIωP22RIs-RIs+1-RIs-1-MI′ξ˙1sR˙Is,Qn1ξ¨1s=∑IMI′R˙Is2-gβ-Qn1ξ˙1sξ˙2s,Qnkξ¨ks=Qnk-1ξ˙k-1s2-1β-Qnkξ˙ksξ˙k+1s1-δkK,k=2,…,K;s=1,…,P.The quantities appearing in these equations are all defined in Section 8.2. One should notice that number of equations is quite large due to the evolution parameter swhich is absent from the standard Car–Parrinello formulation as reported in Section 9.3.In the next section, we address some important implementational issues.


## 10. Computational Implementation

In this section, we present some important implementational considerations. In particular, we approximate the electronic interaction with an effective pseudopotential. In addition, we demonstrate that the computational complexity may be reduced if the orbitals are expressed in terms of a suitable functional basis. While ab initio molecular dynamics is restricted to small molecules, because of its computational complexity, the method may be extended to larger systems if a hybrid approach is followed. The latter is introduced in the last subsection.


### 10.1. Pseudopotential

For many calculations, the complete knowledge of the electronic interaction is not essential for the required precision. Consequently, for the sake of computational efficiency, the exact electronic potential may be approximated by means of an effective potential, called the pseudopotential [46–48]. Moreover, the relativistic effects associated with the core electrons may be implicitly incorporated into the pseudopotential, without recourse to explicit and intricate approaches.A commonly employed pseudopotential, the dual-space Gaussian pseudopotential [46, 47], is composed of two parts, a local potential and a nonlocal potential:(92)VPPr,r′=VLr+VNLr,r′.The local potential, which described the local interaction, is defined as(93)VLr≙-Zionrerfr2rL+exp⁡-12rrL2×C1+C2rrL2+C3rrL4+C4rrL6,where erf is the error function while rL, C1, C2, C3, and C4 are adjustable empirical parameters. The nonlocal potential(94)VNLr,r′=∑i∑l∑m=-llYlmθ,φpilrVilpilrYlm∗θ′,φ′that describes the nonlocal interaction is defined in terms of the complex spherical harmonics functions Ylmθ,φ and of the Gaussian projectors:(95)pilr=2rl+2i-1×e-1/2r/rl2rll+4i-1/2Γl+4i-1/2.In these equations, l is the azimuthal quantum number and m is the magnetic quantum number while Γ is the well-known gamma function.In the next subsection, we explore the advantages of projecting the orbitals on a functional basis.


### 10.2. Orbitals and Basis Functions

Any continuous function, such as an orbital, may be represented as a linear combination of basis functions:(96)ψi=∑jcijfj,where cij are the projection coefficients and fj are the basis functions. Such decomposition may be either physically motivated, computationally motivated, or both. For instance, if physical solutions of the Schrödinger equation are known, it is possible to project the orbitals on these solutions [49–51]. The most common Schrödinger basis functions are the Slater-type basis functions [49]:(97)fmSr=NmSrxmxrymyrzmzexp⁡-ζmrand the Gaussian-type basis functions [10, 50]:(98)fmGr=NmGrxmxrymyrzmzexp⁡-αmr2.where NmS and NmG are normalisation constants while m represents the magnetic quantum numbers. Once an orbital is projected on a basis, it is entirely determined by its projection coefficients. Thus, the resulting representation is parametric. As a result, the determination of the projection coefficients is equivalent to the determination of the orbitals: the closer the basis functions are to the real solution, the more efficient the calculation is. This is due to the fact that fewer coefficients are required to adequately represent the underlying orbital. An even more realistic representation may be obtained if the basis functions are centred upon their respective nuclei [13]:(99)fir=∑μ,IciμIϕμr-RI,where ϕμ is an atomic basis function and RI is the location of a given nucleus.The orbitals may also be represented by plane waves [14]:(100)fkPWr=1Ωexp⁡ik·r,where Ω is the volume of the cell associated with the underlying discrete grid and k is the wave vector associated with the plane wave. From a physical point of view, plane waves form an appropriate basis when the orbitals are smooth functions. Otherwise, a large number of plane waves are required for an accurate reconstruction of the orbitals. Consequently, the computational burden associated with the parametric representation becomes rapidly prohibitive. Nevertheless, if the orbitals are smooth functions, one may take advantage of the Fourier transform in order to reduce the complexity involved in the calculations of the derivatives. For instance, the Laplacian of a Fourier transformed function is obtained by multiplying this function by the square of the module of the wave vector:(101)Δ⟹Fk2,k2⟹F-1Δ, where F is the Fourier transform. If the Fourier transform is calculated with the Fast Fourier Transform algorithm, the complexity of such a calculation, for a N×N×N discrete, grid, reduces to(102)ON6FT⟹ONlog⁡N3FFT,which explains why plane wave basis and the Fourier transforms are prevalent in ab initio molecular dynamics simulations. The computational efficiency may be further improved if the Fourier transformations are evaluated with graphics processing units (GPU), as their architecture makes them particularly suited for these calculations, especially regarding speed [12].The orbital functions are not always smooth as a result of the strength of the electronic interaction. In this particular case, the plane wave basis constitutes an inappropriate choice as a very large number of basis elements are required in order to describe the quasi-discontinuities. Nevertheless, one would be inclined to retain the low computational complexity associated with plane waves and the Fourier transform while being able to efficiently describe the rougher parts of the orbitals. This may be achieved with a wavelet basis and the wavelet transforms [20, 48, 52–54]. As opposed to the plane wave functions, which span the whole space, the wavelets are spatially localised. Therefore, they are particularly suited for the description of discontinuities and fast-varying functions. In addition, the wavelet basis provides a multiresolution representation of the orbitals which means that the calculations may be performed efficiently at the required resolution level. The wavelet functions are filter banks [55]. In one dimension, they involve two functions:(i) the scaling function:(103)ϕx=2∑i=1-mmhjϕ2x-iwhich is a high-pass filter responsible for the multiresolution aspect of the transform;(ii) the mother wavelet:(104)κx=2∑i=1-mmgjκ2x-iwhich is a band-pass filter responsible for the basis localisation. The coefficients of the scaling and the mother wavelet function are related by(105)gi=-1ih-i+1.In three dimensions, at the lowest resolution, the scaling function is given by(106)ϕijkr=ϕxa-iϕya-jϕza-kwhile the mother wavelets become(107)κijk1r=κxa-iϕya-jϕza-k,κijk2r=ϕxa-iκya-jϕza-k,κijk3r=κxa-iκya-jϕza-k,κijk4r=ϕxa-iϕya-jκza-k,κijk5r=κxa-iϕya-jκza-k,κijk6r=ϕxa-iκya-jκza-k,κijk7r=κxa-iκya-jκza-k,where a is the resolution of the underlying computational grid. As a result, any orbital function may be approximated by a truncated expansion:(108)ψr=∑i,j,ksijkϕijkr+∑i,j,k∑ν=17dijkνκijkνr,where sijk and dijkν are the wavelet coefficients. Naturally, higher resolutions may be achieved; the maximum resolution is determined by the resolution of the computational grid. As for the Fourier transform, the wavelet transform admits Fast Wavelet Transform implementation which has the same complexity as its Fourier counterpart.In the next subsection, we briefly present a hybrid approach which involves both ab initio molecular dynamics and molecular dynamics.


### 10.3. Hybrid Quantum Mechanics–Molecular Dynamics Methods

Because of its inherent complexity, the applicability of ab initio molecular dynamics is essentially limited to small molecular domains. Molecular dynamics, which we recently reviewed in [30], is suitable for larger molecular system. A drawback is that it is not adapted for the simulation of chemical complexes, since it is not based on first principles as its quantum mechanics counterpart. Rather, molecular dynamics relies on empirical potentials [56] and classical mechanics which allow for the simulation of much larger molecular complexes.Therefore, in order to simulate larger systems, hybrid approaches [12, 57–59] must be followed. In such a method, a small region of interest is simulated with ab initio molecular dynamics, while the rest of the molecular system is approximated with molecular mechanics:(109)LQM/MM=LCPQM+LMM+LQM-MM.An extra Lagrangian term is required in order to ensure a proper coupling between the quantum degrees of freedom and the classical degrees of freedom. This Lagrangian consists of a bounded and an unbounded part. The bounded part consists of the stretching, bending, and torsional terms which are characterised by their distances, angles, and dihedral angles, respectively. The unbounded part contains the electrostatic interaction between the molecular mechanics atoms and the quantum density as well as the steric interaction which may be approximated, for instance, by a van der Waals potential.


## 10.1. Pseudopotential

For many calculations, the complete knowledge of the electronic interaction is not essential for the required precision. Consequently, for the sake of computational efficiency, the exact electronic potential may be approximated by means of an effective potential, called the pseudopotential [46–48]. Moreover, the relativistic effects associated with the core electrons may be implicitly incorporated into the pseudopotential, without recourse to explicit and intricate approaches.A commonly employed pseudopotential, the dual-space Gaussian pseudopotential [46, 47], is composed of two parts, a local potential and a nonlocal potential:(92)VPPr,r′=VLr+VNLr,r′.The local potential, which described the local interaction, is defined as(93)VLr≙-Zionrerfr2rL+exp⁡-12rrL2×C1+C2rrL2+C3rrL4+C4rrL6,where erf is the error function while rL, C1, C2, C3, and C4 are adjustable empirical parameters. The nonlocal potential(94)VNLr,r′=∑i∑l∑m=-llYlmθ,φpilrVilpilrYlm∗θ′,φ′that describes the nonlocal interaction is defined in terms of the complex spherical harmonics functions Ylmθ,φ and of the Gaussian projectors:(95)pilr=2rl+2i-1×e-1/2r/rl2rll+4i-1/2Γl+4i-1/2.In these equations, l is the azimuthal quantum number and m is the magnetic quantum number while Γ is the well-known gamma function.In the next subsection, we explore the advantages of projecting the orbitals on a functional basis.


## 10.2. Orbitals and Basis Functions

Any continuous function, such as an orbital, may be represented as a linear combination of basis functions:(96)ψi=∑jcijfj,where cij are the projection coefficients and fj are the basis functions. Such decomposition may be either physically motivated, computationally motivated, or both. For instance, if physical solutions of the Schrödinger equation are known, it is possible to project the orbitals on these solutions [49–51]. The most common Schrödinger basis functions are the Slater-type basis functions [49]:(97)fmSr=NmSrxmxrymyrzmzexp⁡-ζmrand the Gaussian-type basis functions [10, 50]:(98)fmGr=NmGrxmxrymyrzmzexp⁡-αmr2.where NmS and NmG are normalisation constants while m represents the magnetic quantum numbers. Once an orbital is projected on a basis, it is entirely determined by its projection coefficients. Thus, the resulting representation is parametric. As a result, the determination of the projection coefficients is equivalent to the determination of the orbitals: the closer the basis functions are to the real solution, the more efficient the calculation is. This is due to the fact that fewer coefficients are required to adequately represent the underlying orbital. An even more realistic representation may be obtained if the basis functions are centred upon their respective nuclei [13]:(99)fir=∑μ,IciμIϕμr-RI,where ϕμ is an atomic basis function and RI is the location of a given nucleus.The orbitals may also be represented by plane waves [14]:(100)fkPWr=1Ωexp⁡ik·r,where Ω is the volume of the cell associated with the underlying discrete grid and k is the wave vector associated with the plane wave. From a physical point of view, plane waves form an appropriate basis when the orbitals are smooth functions. Otherwise, a large number of plane waves are required for an accurate reconstruction of the orbitals. Consequently, the computational burden associated with the parametric representation becomes rapidly prohibitive. Nevertheless, if the orbitals are smooth functions, one may take advantage of the Fourier transform in order to reduce the complexity involved in the calculations of the derivatives. For instance, the Laplacian of a Fourier transformed function is obtained by multiplying this function by the square of the module of the wave vector:(101)Δ⟹Fk2,k2⟹F-1Δ, where F is the Fourier transform. If the Fourier transform is calculated with the Fast Fourier Transform algorithm, the complexity of such a calculation, for a N×N×N discrete, grid, reduces to(102)ON6FT⟹ONlog⁡N3FFT,which explains why plane wave basis and the Fourier transforms are prevalent in ab initio molecular dynamics simulations. The computational efficiency may be further improved if the Fourier transformations are evaluated with graphics processing units (GPU), as their architecture makes them particularly suited for these calculations, especially regarding speed [12].The orbital functions are not always smooth as a result of the strength of the electronic interaction. In this particular case, the plane wave basis constitutes an inappropriate choice as a very large number of basis elements are required in order to describe the quasi-discontinuities. Nevertheless, one would be inclined to retain the low computational complexity associated with plane waves and the Fourier transform while being able to efficiently describe the rougher parts of the orbitals. This may be achieved with a wavelet basis and the wavelet transforms [20, 48, 52–54]. As opposed to the plane wave functions, which span the whole space, the wavelets are spatially localised. Therefore, they are particularly suited for the description of discontinuities and fast-varying functions. In addition, the wavelet basis provides a multiresolution representation of the orbitals which means that the calculations may be performed efficiently at the required resolution level. The wavelet functions are filter banks [55]. In one dimension, they involve two functions:(i) the scaling function:(103)ϕx=2∑i=1-mmhjϕ2x-iwhich is a high-pass filter responsible for the multiresolution aspect of the transform;(ii) the mother wavelet:(104)κx=2∑i=1-mmgjκ2x-iwhich is a band-pass filter responsible for the basis localisation. The coefficients of the scaling and the mother wavelet function are related by(105)gi=-1ih-i+1.In three dimensions, at the lowest resolution, the scaling function is given by(106)ϕijkr=ϕxa-iϕya-jϕza-kwhile the mother wavelets become(107)κijk1r=κxa-iϕya-jϕza-k,κijk2r=ϕxa-iκya-jϕza-k,κijk3r=κxa-iκya-jϕza-k,κijk4r=ϕxa-iϕya-jκza-k,κijk5r=κxa-iϕya-jκza-k,κijk6r=ϕxa-iκya-jκza-k,κijk7r=κxa-iκya-jκza-k,where a is the resolution of the underlying computational grid. As a result, any orbital function may be approximated by a truncated expansion:(108)ψr=∑i,j,ksijkϕijkr+∑i,j,k∑ν=17dijkνκijkνr,where sijk and dijkν are the wavelet coefficients. Naturally, higher resolutions may be achieved; the maximum resolution is determined by the resolution of the computational grid. As for the Fourier transform, the wavelet transform admits Fast Wavelet Transform implementation which has the same complexity as its Fourier counterpart.In the next subsection, we briefly present a hybrid approach which involves both ab initio molecular dynamics and molecular dynamics.


## 10.3. Hybrid Quantum Mechanics–Molecular Dynamics Methods

Because of its inherent complexity, the applicability of ab initio molecular dynamics is essentially limited to small molecular domains. Molecular dynamics, which we recently reviewed in [30], is suitable for larger molecular system. A drawback is that it is not adapted for the simulation of chemical complexes, since it is not based on first principles as its quantum mechanics counterpart. Rather, molecular dynamics relies on empirical potentials [56] and classical mechanics which allow for the simulation of much larger molecular complexes.Therefore, in order to simulate larger systems, hybrid approaches [12, 57–59] must be followed. In such a method, a small region of interest is simulated with ab initio molecular dynamics, while the rest of the molecular system is approximated with molecular mechanics:(109)LQM/MM=LCPQM+LMM+LQM-MM.An extra Lagrangian term is required in order to ensure a proper coupling between the quantum degrees of freedom and the classical degrees of freedom. This Lagrangian consists of a bounded and an unbounded part. The bounded part consists of the stretching, bending, and torsional terms which are characterised by their distances, angles, and dihedral angles, respectively. The unbounded part contains the electrostatic interaction between the molecular mechanics atoms and the quantum density as well as the steric interaction which may be approximated, for instance, by a van der Waals potential.


## 11. Conclusions

In this paper, we have presented a comprehensive, but yet concise, review of ab initio molecular dynamics from a computational perspective and from first principles. Although it is always hazardous to speculate about the future, we foresee two important breakthroughs. Fourier transforms, which constitute an essential part of many molecular simulations, may be evaluated with high performance on graphical processing units or GPU [54]. The same remark applies to the wavelet transform whose role is essential in describing discontinuous orbitals [52]. This opens the door to high performance simulations, while tempering the limitations associated with computational complexity.For decades, one of the main bottlenecks of molecular simulations has been the calculation of the density functionals. As in numerous fields, machine learning constitutes a promising avenue for the fast and efficient evaluation of these functionals without recourse to explicit calculations [17, 24, 25, 60]. Here, the density functional is simply learned from existing examples with machine learning techniques. This paves the way for the simulation for larger and more complex molecular systems. We plan to review machine learning-based approaches in a future publication.



---

*Source: 9839641-2018-04-29.xml*