Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Computational Methods for Ab Initio Molecular Dynamics

Computational Methods for Ab Initio Molecular Dynamics Hindawi Advances in Chemistry Volume 2018, Article ID 9839641, 14 pages https://doi.org/10.1155/2018/9839641 Review Article 1 2 Eric Paquet and Herna L. Viktor National Research Council, 1200 Montreal Road, Ottawa, ON, Canada University of Ottawa, 800 King Edward Avenue, Ottawa, ON, Canada Correspondence should be addressed to Eric Paquet; eric.paquet@nrc-cnrc.gc.ca Received 23 February 2018; Accepted 20 March 2018; Published 29 April 2018 Academic Editor: Yusuf Atalay Copyright © 2018 Eric Paquet and Herna L. Viktor. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 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 aer ff ent 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. 1. Introduction important, we expose how these equations and concepts are related to one another. Ab initio molecular dynamics (AIMD) is an irreplaceable eTh paper is organised as follows. In Section 2, quantum technique for the realistic simulation of complex molecular mechanics is reviewed from a molecular dynamics perspec- systems and processes associated with biological organisms tive. Two important approximations are introduced, namely, [1, 2] such as monoclonal antibodies as illustrated in Figure 1. the adiabatic and the Born–Oppenheimer approximations. With ab initio molecular dynamics, it is possible to predict in Subsequently, in Section 3, the Ehrenfest molecular dynamics silico phenomena for which an in vivo experiment is either is presented, which allows for a semiclassical treatment of too difficult or expensive, or even currently deemed infeasible the nuclei. This opens the door, in Section 4, to the Born– [3, 4]. Ab initio molecular dynamics essentially differs from Oppenheimer molecular dynamics formulation. This is fol- molecular dynamics (MD) in two ways. Firstly, AIMD is lowed, in Section 5, by the Hartree–Fock formulation which based on the quantum Schrodin ¨ ger equation while its clas- approximates the atomic antisymmetric wave function by a single determinant of the orbitals. Sections 6 and 7 introduce sical counterpart relies on Newton’s equation. Secondly, MD relies on semiempirical eeff ctive potentials which approxi- the Kohn–Sham formulation whose primary objective is to mate quantum effects, while AIMD is based on the real physi- replace the interacting electrons by a cfi titious, but equiva- lent, system of noninteracting particles. Electrons are reintro- cal potentials. In this paper, we present a review of ab initio molecular duced as dynamical quantities, in Section 8, through the Car– dynamics from a computational perspective and from rfi st Parrinello formulation. Isothermal and isobaric processes are principles. Our main aim is to create a document which is addressed in this section. In Section 9, a path integral for- self-contained, coherent, and concise but still as complete mulation of molecular dynamics is presented, which makes as possible. As the theoretical details are essential for real- allowance for a quantum treatment of the nuclear degrees life implementation, we have provided the equations for the of freedom. Finally, in Section 10, some important compu- relevant physical aspect and approximations presented. Most tational issues are addressed such as a simplification of the 2 Advances in Chemistry ∇ Š , 𝜕R 𝜕 𝜕 Δ Š ⋅ . 𝜕R 𝜕R 𝐼 𝐼 (3) eTh gradient operator, a vector, is related to the momentum while the Laplacian operator, a scalar quantity, is associated with the kinetic energy [5–7]. eTh electronic Hamiltonian is defined as the sum of the kinetic energy of the electrons, the Figure 1: Monoclonal antibody 1F1 isolated from a 1918 influenza potential energy associated with each pair of electrons, the pandemic (Spanish Flu) survivor. potential energy associated with each pair of electron-nu- cleus, and the potential energy concomitant to each pair of nuclei: 2 2 electronic interaction with the pseudopotential, the repre- ℏ 1 𝑒 H (r, R) Š − ∑ Δ + ∑ sentation of orbitals in terms of a functional basis, the use 𝑒 𝑖 󵄩 󵄩 󵄩 󵄩 2𝑚 4𝜋𝜀 󵄩 󵄩 r − r 𝑒 0 𝑖 𝑖<𝑗 󵄩 󵄩 𝑖 𝑗 of the Fourier and wavelet transform in order to reduce the 󵄩 󵄩 (4) computational complexity, and the simulation of larger sys- 𝑒 𝑍 𝑍 1 𝑒 𝑍 1 𝐼 𝐽 tems withhybrid moleculardynamics.Thisisfollowedbya − ∑ + ∑ , 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 4𝜋𝜀 4𝜋𝜀 󵄩 󵄩 󵄩 R − R 󵄩 conclusion in Section 11. R − r 0 0 𝐼 𝐽 𝐼,𝑖 󵄩 𝐼 𝑗 󵄩 𝐼<𝐽 󵄩 󵄩 󵄩 󵄩 where 𝑚 is the electronic mass, 𝑒 is the electronic charge, 2. Quantum Molecular Dynamics 𝑍 istheatomicnumber(thenumberofprotons inagiven and the Schrödinger Equation: nucleus), and 𝜀 is the vacuum permittivity. eTh quantum A Molecular Perspective molecular system is characterised by a wave function Φ(r, R;𝑡) whose evolution is determined by the time-dependent We review some important notions of quantum mechanics Schrodin ¨ ger equation [5–7]: from a molecular perspective. In quantum mechanics [5– 7], the Hamiltonian [5–7] is the operator corresponding to 𝑖ℏ Φ (r, R;𝑡 )= H (r, R)Φ (r, R;𝑡 ), (5) 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 where𝑖= −1.TheSchr odin ¨ ger equation admits a stationary energy of the constituent particles; in occurrence the elec- or time-independent solution: trons and the nuclei H r, R Ψ r, R =𝐸 R Ψ r, R , ( ) ( ) ( ) ( ) (6) 𝑒 𝑘 𝑘 𝑘 H (r, R)=−∑ Δ + H (r, R), (1) where 𝐸 (R)is the energy associated with the electronic 𝐼 𝑒 𝑘 2𝑀 wave functionΨ (r,R)and 𝑘 is a set of quantum numbers that labelled the Eigenstates or wave functions as well as where ℏ is the Planck constant and 𝑀 is the mass of a given the Eigenvalues or energies associated with the stationary nucleus. The electronic and nuclear Cartesian coordinates for Schrodin ¨ ger equation. The total wave function Φ(r,R;𝑡) may a given particle are defined, respectively, as be expended in terms of time-dependent nuclear wave functions 𝜒 (R;𝑡) and stationary electronic wave functions Ψ (r,R): r Š [𝑟 ,𝑟 ,𝑟 ] , 𝑖 𝑖,𝑥 𝑖,𝑦 𝑖,𝑧 (2) ∞ Φ (r, R;𝑡 )= ∑𝜒 (R;𝑡 )Ψ (r, R). (7) R Š [𝑅 ,𝑅 ,𝑅 ] . 𝑘 𝑘 𝐼 𝐼,𝑥 𝐼,𝑦 𝐼,𝑧 𝑘=0 The ensemble of all electronic coordinates is represented by It should be noticed that this expansion is exact and does r while the ensemble of all nuclear coordinates is denoted by not involve any approximation. eTh electronic wave functions R. In addition, the gradient and the Laplacian operators for a are solution of the stationary Schrodin ¨ ger equation, an Eigen given electron 𝑖 and nucleus 𝐼 are defined, respectively, as equation which involves the electronic Hamiltonian. If one substitutes this expansion in the time-dependent Schrodin ¨ ger equation, one obtains a system of equations for the evolution ∇ Š , of the nuclear wave functions: 𝜕r 𝜕 𝜕 𝑖ℏ =[−∑ Δ +𝐸 (R)]𝜒 + ∑C 𝜒 . (8) Δ Š ⋅ , 𝐼 𝑘 𝑘 𝑘𝑙 𝑙 𝜕𝑡 2𝑀 𝜕r 𝜕r 𝐼 𝑖 𝑖 𝜕𝜒 Advances in Chemistry 3 As may be seen, the nuclear wave functions are coupled to Furthermore, it is assumed that the electronic wave function the electronic wave functions. eTh strength of the coupling is is in its ground or nonexcited state: determined by the coupling coefficients which form a matrix: 𝑘=0󳨐⇒ (15) C = ∫ 𝑑rΨ [−∑ Δ ]Ψ 𝑘𝑙 𝐼 𝑙 𝑘 H Ψ =𝐸 Ψ . 𝑒 0 0 0 2𝑀 (9) Thisoccursifthethermalenergyislower thantheenergy + ∑ [∫ 𝑑r Ψ [−𝑖ℏ∇ ]Ψ] [−𝑖ℏ∇ ]. 𝑘 𝐼 𝑙 𝐼 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, This system of equations is too complex to be solved directly. in which the motion of heavy nuclei is described by a semi- Consequently, in the next subsection, we consider two impor- classical equation. tant approximations of the Schrodin ¨ ger equation, namely, the adiabatic and the Born–Oppenheimer approximations, which aim to reduce such a complexity. es Th e approxima- 3. Ehrenfest Molecular Dynamics tions, as well as those that later follow, reduce substantially the duration of the calculations allowing for larger molecular Another possible approximation is to assume that the nuclear motion is semiclassicalasitisthe casewhen thenuclei systems to be simulated and longer time-scales to be explored [8, 9]. are relatively heavy. This implies that, instead of being determined by the Schrodin ¨ ger equation, the average nuclear motion is determined by Newton’s equation. The right form 2.1. Adiabatic and Born–Oppenheimer Approximations. The for this equation is given by the Ehrenfest theorem [6, 12– above-mentioned system of differential equations is too 14] which states that the potential, which governs the classical complex to be solved directly. Some approximations must be motion of the nucleons, is equal to the expectation or average, performed in order to reduce the computational complexity, in the quantum sense, of the electronic Hamiltonian with while maintaining the predictive accuracy of the calculations. respect to the electronic wave function: The first approximation, which is called the adiabatic approx- imation [5–7], assumes that the electronic wave functions 󵄨 󵄨 ̈ 󵄨 󵄨 𝑀 R (𝑡 )=−∇ ⟨Ψ󵄨 H 󵄨 Ψ⟩ = −∇ ⟨H ⟩ , (16) 𝐼 𝐼 𝐼 𝑒 𝐼 𝑒 󵄨 󵄨 Ψ adapt quasi-instantaneously to a variation of the nuclear configuration. This approximation is justified by the fact where the bra–ket notation is to be understood as that the electrons are much lighter than the nucleus. Conse- quently, such an approximation is valid, unless the motion ⟨Ψ |O| Ψ⟩ ≡ ∫ 𝑑r Ψ (r)OΨ (r) (17) of the electron becomes relativistic, which may happen if the electrons are too close to the nucleus. Mathematically, this 2 2 approximation implies that for any operator O and where R ≡𝑑 R /𝑑𝑡 .This 𝐼 𝐼 equation may be further simplified if the adiabatic and the ∇ Ψ =0. (10) 𝐼 𝑙 Born–Oppenheimer approximations are enforced: As a result, the coupling matrix becomes diagonal: 󵄨 󵄨 󵄨 󵄨 ∇ ⟨Ψ󵄨 H 󵄨 Ψ⟩⇒ 󳨐 𝐼 𝑒 󵄨 󵄨 󵄨 󵄨 󵄨 󵄨 ∇ ⟨Ψ H Ψ ⟩󳨐⇒ C =[−∑ ∫ 𝑑rΨ Δ Ψ ]𝛿 , (11) 󵄨 󵄨 (18) 𝑘𝑙 𝐼 𝑘 𝑘𝑙 𝐼 0 𝑒 0 𝑘 󵄨 󵄨 2𝑀 󵄨 󵄨 󵄨 󵄨 ⟨Ψ 󵄨 ∇ H 󵄨 Ψ ⟩. 0 𝐼 𝑒 0 󵄨 󵄨 where 𝛿 is the well-known Kronecker delta. 𝑘𝑙 The second approximation, which is called the Born– In the specific context of Ehrenfest molecular dynamics, the Oppenheimer approximation [6, 10, 11], assumes that the electrons follow a time-dependent Schrodin ¨ ger equation: electronic and the nuclear motions are separable as a result of the difference between nuclear and electronic masses. As a 𝜕 (19) 𝑖ℏ Ψ= H Ψ. result, the expansion reduces to a single Eigen state: 𝜕𝑡 Φ (r, R;𝑡 )≈Ψ (r, R)𝜒 (R;𝑡 ). (12) The fact that this equation is time-dependent ensures that the 𝑘 𝑘 orthogonality of the orbitals is maintained at all times. eTh In addition, as the energy associated with the wave function motivation is to be found in the fact that the electronic Hamil- is usually much larger than the corresponding coupling con- tonian is a Hermitian operator [6]. As stated earlier, this is not stant: thecaseintheBorn–Oppenheimer approximationinwhich the electronic wave function is governed by the stationary C ≪𝐸 , (13) 𝑘𝑘 𝑘 Schrodinger equation which does not maintain the orthonor- the time-dependent nuclear Schrodinger further reduces to mality over time. In thenext section,we apply theadiabaticandBorn– Oppenheimer approximation in the context of ab initio mo- 𝑖ℏ =[−∑ Δ +𝐸 R ]𝜒 . ( ) (14) 𝐼 𝑘 𝑘 𝜕𝑡 2𝑀 lecular dynamics. 𝜕𝜒 4 Advances in Chemistry 4. Born–Oppenheimer Molecular Dynamics for the electronic equations of motion. One should notice the presence of the Lagrange multiplier in the electronic equa- In Born–Oppenheimer molecular dynamics [6, 10, 11], it tions of motion which ensures that the orbital remains ortho- is assumed that the adiabatic and the Born–Oppenheimer normal at all time. approximations are valid and that the nuclei follow a semi- In the next section, we introduce another approximation classical Newton equation whose potential is determined in order to further reduce the complexity of the electronic by the Ehrenfest theorem. It is further assumed that the Hamiltonian. electronic wave function is in its ground state (lowest energy). Since the stationary Schrodinger equation is used for the elec- 5. Hartree–Fock Molecular Dynamics tronic degrees of freedom, the orthogonality condition must be enforced with a Lagrange multiplier as this condition is not Fromnowon,we shalladopttheatomicunitsystem: preserved by the stationary Schrodin ¨ ger equation over time. eTh orthogonality of the orbitals 𝜓 is a physical requirement 𝑖 1 𝑚 =𝑒=ℏ= =1 (25) whichmustbeenforcedatalltime in ordertobeable to 𝑒 4𝜋𝜀 compute real observable quantities. The Born–Oppenheimer molecular dynamics is characterised by a Lagrangian [6] in order to alleviate the notation. In Hartree–Fock molecular which is defined as the difference between the kinetic and the dynamics [15], the antisymmetric electronic wave function potential energy: is approximated by a single determinant of the electronic orbitals 𝜓 (r): 𝑖 𝑖 BO 2 󵄨 󵄨 ̇ 󵄨 󵄨 L = ∑𝑀 R −⟨Ψ 󵄨 H 󵄨 Ψ ⟩ 𝐼 0 𝑒 0 𝐼 󵄨 󵄨 𝐼 𝜓 (r)𝜓 (r) ⋅⋅⋅ 𝜓 (r ) 1 1 1 2 1 𝑁 (20) [ ] 𝜓 (r)𝜓 (r) ⋅⋅⋅ 𝜓 (r ) [ ] 2 1 2 2 2 𝑁 + ∑Λ (⟨𝜓 |𝜓 ⟩−𝛿 ). 𝑖 𝑗 1 [ ] Ψ= det . (26) [ ] 𝑖,𝑗 . . . [ ] 𝑁! . . . [ . . d . ] eTh rfi st term of the Lagrangian corresponds to the kinetic 𝜓 (r)𝜓 (r ) ⋅⋅⋅ 𝜓 (r ) [ ] 𝑁 1 𝑁 𝑁 𝑁 𝑁 energy of the nuclei, the second term corresponds to the nuclear potential as obtained from the Ehrenfest theorem, Two operators are associated with the electronic interaction. and the last term is a Lagrange function which ensures that eTh rfi st one, the Coulomb operator, corresponds to the the orbitals remain orthonormal at all time. eTh Lagrange electrostatic interaction between the orbitals: multipliers are denoted by Λ .Since theelectrons follow a Fermi–Dirac statistics [6], they obey the Pauli Exclusion 󸀠 ∗ 󸀠 󸀠 J r 𝜓 r ≡ [∫ 𝑑r 𝜓 (r ) 𝜓 (r )]𝜓 r . (27) ( ) ( ) ( ) Principle (two electrons cannot be in the same quantum 𝑗 𝑖 𝑗 󵄩 󵄩 𝑗 𝑖 󵄩 󸀠 󵄩 󵄩 r − r 󵄩 󵄩 󵄩 state) which means that the electronic wave function must be an antisymmetric function of its orbitals, namely, the wave The second operator, the exchange operator [6, 15], is associ- functions associated with the individual electrons. ated with the exchange energy, a quantum mechanical eect ff The equations of motion associated with this Lagrangian thatoccursasaresultofthePauliExclusionPrinciple: are obtained from the Euler–Lagrange equations [6]. eTh re isonesystemofEuler–Lagrangeequations forthe nuclear 󸀠 ∗ 󸀠 󸀠 K (r)𝜓 (r)≡[∫ 𝑑r 𝜓 (r ) 𝜓 (r )] 𝜓 (r). (28) 𝑗 𝑖 󵄩 󵄩 𝑖 𝑗 degrees of freedom: 𝑗 󵄩 󵄩 󵄩 r − r 󵄩 󵄩 󵄩 𝑑 𝜕L 𝜕L − =0, From the Coulomb and the exchange operators, an electronic (21) 𝜕 R 𝜕R 𝐼 Hamiltonian may be inferred: and one system of Euler–Lagrange equations for their elec- HF H =− Δ+𝑉 (r)+2∑J (r)− ∑K (r) tronic counterpart: ext 𝑗 𝑗 𝑗 𝑗 (29) 𝑑 𝛿L 𝛿L − =0󳨐⇒ 󵄨 󵄨 1 HF 󵄨 󵄨 𝛿⟨ 𝜓 󵄨 𝛿⟨𝜓 󵄨 ≡ Δ+𝑉 (r), 𝑖 𝑖 󵄨 󵄨 (22) 𝑑 𝛿L 𝛿L − =0. where the external potential is defined as ∗ ∗ 𝛿 𝜓 (r) (r) 𝑖 𝑖 𝑍 𝑍 𝐼 𝐽 From the Euler–Lagrange equations, one obtains 𝑉 (r) Š − ∑ + ∑ . (30) ext 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 R − r󵄩 󵄩 R − R 󵄩 𝐼 𝐼 𝐽 𝐼 󵄩 󵄩 𝐼<𝐽 󵄩 󵄩 󵄨 󵄨 ̈ 󵄨 󵄨 𝑀 R (𝑡 )=−∇ min ⟨Ψ 󵄨 H 󵄨 Ψ ⟩ 𝐼 𝐼 𝐼 0 𝑒 0 (23) 󵄨 󵄨 This Hamiltonian determines, in turn, the electronic Lagrangian: for the nuclear equations of motion, and 󵄨 󵄨 HF HF 󵄨 󵄨 󵄨 󵄨 H 𝜓 (r)= ∑Λ 𝜓 (r) L =−⟨Ψ H Ψ ⟩+ ∑Λ (⟨𝜓 |𝜓 ⟩−𝛿 ). 󵄨 󵄨 𝑒 𝑖 𝑖 0 0 𝑖 𝑗 𝑒 𝑒 (24) 󵄨 󵄨 (31) 𝑗 𝑖,𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝛿𝜓 𝑑𝑡 𝑑𝑡 𝑑𝑡 𝑖𝑗 𝑖𝑗 𝑖𝑗 Advances in Chemistry 5 A Lagrange multiplier term is added to the Lagrangian in Some of these approximations are considered in the next sec- order to ensure that the orbitals remain orthonormal at all tion. time: a quantum mechanical requirement [6]. From the Euler– From the exchange-correlation energy, one may define Lagrange equations, one obtains the equations of motion for the exchange-correlation potential: the electrons: [𝑛 ] XC 𝑉 (r) Š (37) HF XC H 𝜓 = ∑Λ 𝜓 (r) 𝑒 𝑖 𝑖 (32) which is the functional derivative [22] of the exchange-corre- lation energy with respect to the electronic density. In addi- which differ from (24) by the Hamiltonian. tion, one may define the Kohn–Sham potential: In the next subsection, we seek to replace the interacting electrons by a cfi titious but equivalent system of noninter- 𝑉 Š 𝑉 +𝑉 +𝑉 . (38) KS 𝐻 ext XC acting particles in order to further reduce the computational complexity. 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 6. Kohn–Sham Molecular Dynamics determined by the Euler–Lagrange equations are eTh Kohn–Sham formulation [16, 17] aims to replace the KS interacting electrons by a fictitious system of noninteracting H 𝜓 = ∑Λ 𝜓 , 𝑒 𝑖 𝑗 (39) particles that generate the same electronic density as the physical system of interacting particles. eTh electronic density where the cfi titious one-particle Hamiltonian is given by is defined as KS H ≡− Δ+𝑉 . (40) 𝑛 Š ∑𝑓 ⟨𝜓 |𝜓 ⟩, 𝑖 𝑖 𝑖 𝑒 KS (33) A canonical form may be obtained if a unitary transformation where 𝑓 is the occupation number of a given orbital. is appliedtothe previous equation: In this formulation of AIMD, the Ehrenfest term is KS replaced by the Kohn–Sham energy: (41) H 𝜓 =𝜀 𝜓 . 𝑖 𝑖 𝑖 󵄨 󵄨 KS 󵄨 󵄨 In the next section, we introduce some density functionals in min ⟨Ψ 󵄨 H 󵄨 Ψ ⟩= min𝐸 [𝜓] . 0 𝑒 0 (34) 󵄨 󵄨 order to replace the interacting electrons by an equivalent but yet simpler system of noninteracting particles. The latter is define as 1 7. Exchange-Correlation Energy KS 𝐸 [𝑛 ] Š − ∑𝑓 ⟨𝜓 |Δ| 𝜓 ⟩+ ∫ 𝑑r𝑉 (r)𝑛 (r) 𝑖 𝑖 𝑖 ext 𝑖=1 The detailed analysis of density functionals, also known as (35) exchange-correlation energies [23], is beyond the scope of + ∫ 𝑑r𝑉 (r)𝑛 (r)+𝐸 [𝑛 ], 𝐻 XC 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 wherethefirst term isthetotal electronic kineticenergyasso- the density functional is to express the electronic interaction ciated with the electrons, the second term is the electrostatics interaction energy between the electronic density and the in terms of the sole electronic density. In thesimplestcase, theexchange-correlationenergyisa external potential, and the third term is the self-electrostatic interaction energy associated with the electronic density. eTh functional of the electronic density alone for which the most latter involves the interaction of the electronic density with it important representative is the local density approximation: self-created electrostatic potential. This potential, called the 1/3 3 3 LDA 4/3 Hartree potential, is obtained by solving the Poisson equation (42) 𝐸 [𝑛 ]=− ( ) ∫ 𝑑r 𝑛 . XC 4 𝜋 [18]: One may take into consideration, in addition to the electronic Δ𝑉 =−4𝑛󳨐𝜋⇒ density, the gradient of the electronic density, in which case the approach is called the generalised gradient approxima- (36) 𝑛( r ) 𝑉 r = ∫ 𝑑r . tion: ( ) 𝐻 󵄩 󵄩 󵄩 󸀠 󵄩 󵄩 r − r 󵄩 󵄩 󵄩 GGA 4/3 𝐸 [𝑛 ]= ∫ 𝑑r 𝑛 𝐹 (𝜁 ), (43) XC XC The last term is the celebrated exchange-correlation energy or density functional [16, 19–21] which takes into account the where the dimensionless reduced gradient is defined as residual electronic interaction, that is, the self-interaction of theelectrons.Unfortunately,thedensity functional hasno 1/3 ‖∇𝑛 ‖ (44) 𝜁 Š ≡ 2 (3𝜋 ) 𝑠. closed-form solution but many approximations are known. 4/3 𝑖𝑗 𝑖𝑗 𝛿𝑛 𝛿𝐸 6 Advances in Chemistry Among these approximations is the B88 approximation: Lagrangian, which is called the Car–Parrinello Lagrangian [11, 26, 27] 1/3 3 6 B88 𝜎 1 󵄨 󵄨 𝐹 (𝜁 )= ( ) + , (45) CP 2 KS 󵄨 󵄨 XC −1 󵄨 󵄨 L = ∑𝑀 R + ∑𝜇⟨ 𝜓 ̇ | 𝜓 ̇ ⟩+⟨Ψ H Ψ ⟩ 󵄨 󵄨 4 𝜋 𝐼 𝐼 𝑖 𝑖 0 𝑒 0 1+6𝜁𝛽 sinh 𝜁 󵄨 󵄨 𝜎 𝜎 𝐼 𝑖 (50) where 𝜎 refers to thespinand thePerdew, Burke, andErnzer- + ∑Λ (⟨𝜓 |𝜓 ⟩−𝛿 ) 𝑖 𝑗 hof (PBE) approximation: 𝑖,𝑗 1/3 differs from the original Kohn–Sham Lagrangian by the 3 3 PBE 𝐹 (𝑠 )= ( ) + . (46) XC second term which associates a fictitious kinetic energy to the 2 −1 4 𝜋 1+𝜇𝑠 𝜅 electronic orbitals. eTh parameter 𝜇 acts as a fictitious elec- tronic mass or inertia. As in the previous sections, the nuclear The exchange energy, which is associated with the Pauli equations of motion Exclusion Principle, may be characterised by the Hartree– Fock energy: 󵄨 󵄨 KS 󵄨 󵄨 󵄨 󵄨 𝑀 R =−∇ ⟨Ψ H Ψ ⟩+ ∑Λ ∇ ⟨𝜓 |𝜓 ⟩ 󵄨 󵄨 𝐼 𝐼 𝐼 0 𝑒 0 𝐼 𝑖 𝑗 󵄨 󵄨 (51) 𝑖,𝑗 ∗ ∗ 󸀠 󸀠 𝜓 (r)𝜓 (r)𝜓 (r)𝜓 (r ) 1 𝑗 𝑖 𝑖 𝑗 HF 󸀠 𝐸 =− ∑ ∫ 𝑑r 𝑑r (47) 󵄩 󵄩 XC 󸀠 and the electronic equations of motion 󵄩 󵄩 2 r − r 󵄩 󵄩 𝑖,𝑗 󵄩 󵄩 KS 𝜇 𝜓 ̈ 𝑡 =−H 𝜓 𝑡 + ∑Λ 𝜓 𝑡 () () () 𝑖 𝑒 𝑖 𝑗 (52) 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: may be inferred from the Euler–Lagrange equations. eTh se equations involve a second order time derivative of the B3 HF LDA B88 GGA orbitals which means that the latter are now proper dynamical 𝐸 =𝑎𝐸 + (1−𝑎 )𝐸 +𝑏Δ𝐸 +𝑐𝐸 XC XC XC XC XC quantities. (48) LDA In the next two subsections, we consider ab initio + (1−𝑐 )𝐸 , XC molecular dynamics at constant temperature and at constant pressure. where 𝑎, 𝑏,and 𝑐 are empirical parameters. Obviously, these parameters are application dependent. Recently, it has been proposed to evaluate the functional 8.2. Massive Thermostating and Isothermal Processes. In the density with machine learning techniques. The functional previous sections, we have implicitly assumed that the molec- density is learned by examples instead of directly solving the ular system under consideration was isolated. Of course, this Kohn–Sham equations [17, 24, 25]. As a result, substantially is not compatible with most experimental conditions [27–29] less time is required to complete the calculations allowing for in which the system is either kept at constant temperature larger system to be simulated and longer time-scales to be (isothermal) in a heat bath or at constant pressure (isobaric) explored. as a consequence, for instance, of the atmospheric pressure. In the next section, we further improve the precision In this subsection and in the next, we explain how to address of the calculations by reintroducing the dynamic electronic these important issues. degrees of freedom which, until now, have been absent from As we have recently explained in a computational review the equations of motion. 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 modiefi d directly as a result 8. Car–Parrinello Molecular Dynamics of the non-Hamiltonian nature of the isothermal process [31]. In order to obtain an isothermal molecular system, a friction 8.1. Equations of Motion. The Kohn–Sham dynamics, as for- term must be recursively added to each degree of freedom mulated in the previous sections, does not take the dynamics [28, 29]. Therefore, the electronic equations of motion must of the electrons into account despite the fact that it is present: be modified as follows: an unattractive unphysical feature. Indeed, one obtains from the Lagrangian and the Euler–Lagrange equations 󵄨 󵄨 󵄨 󵄨 KS 󵄨 󵄨 󵄨 ̈ 󵄨 ̇ ̇ 𝜇 󵄨 𝜓 ⟩= −H 󵄨 𝜓 ⟩+ ∑Λ 𝜓 ⟩− 𝜇 𝜂 󵄨 𝜓 ⟩, 𝑖 𝑒 𝑖 𝑗 1 𝑖 󵄨 󵄨 󵄨 𝑑 𝜕L ≡0 (49) 𝜕 𝜓 (r) 1 𝑒 1 ̈ ̇ ̇ ̇ ̇ 𝑄 𝜂 =2[𝜇 ∑ ⟨𝜓 | 𝜓 ⟩− ]−𝑄 𝜂 𝜂 , 𝑒 1 𝑖 𝑖 𝑒 1 2 2𝛽 𝑖 (53) which means the orbitals are not dynamical fields of the molecular system. Nevertheless, the Lagrangian may be 𝑙 𝑙−1 2 𝑙 𝑄 𝜂 ̈ =[𝑄 𝜂 ̇ − ]−𝑄 𝜂 ̇ 𝜂 ̇ (1 − 𝛿 ), 𝑙 𝑙 𝑙+1 𝑙𝐿 modified in order to also include the dynamic nature of 𝑒 𝑒 𝑙−1 𝑒 electronic degrees of freedom through the introduction of a cfi titious electronic kinetic energy term. The extended 𝑙 = 2, ...,𝐿 𝑑𝑡 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝜇𝑠 𝑖𝑗 𝑖𝑗 𝛽𝜁 Advances in Chemistry 7 while the nuclear equations of motion must take a similar The hidden momenta, associated with the generalised Lan- form: gevin equation, may be marginalised because the evolution of the variables [𝑝, p] , in the free particle limit, is described KS ̈ ̇ ̇ 𝑀 R =−∇ 𝐸 −𝑀 𝜉 R , 𝐼 𝐼 𝐼 𝐼 1 𝐼 by a linear Markovian stochastic differential equation. As a result,thegeneralisedLangevinequationbecomesequivalent 1 2 1 ̈ ̇ ̇ to aLangevinequationwithmemorykerneland noisecorre- 𝑄 𝜉 =[∑𝑀 R − ]−𝑄 𝜉 𝜉 , 𝑛 1 𝐼 𝐼 𝑛 1 2 lation [32–34]: (54) 𝑘 𝑘−1 2 𝑘 𝑞=𝑝 ̇ , ̈ ̇ ̇ ̇ 𝑄 𝜉 =[𝑄 𝜉 − ]−𝑄 𝜉 𝜉 (1 − 𝛿 ) 𝑘 𝑘 𝑘+1 𝑛 𝑛 𝑘−1 𝑛 (58) 𝑝=−𝑉 (𝑞) − ∫ 𝐾 (𝑡−𝑠 )𝑝 (𝑠 )+ 𝜁. 𝑘 = 2,...,𝐾, −∞ where 𝛽≡1/𝑘 𝑇, 𝑇 is the temperature of the heat bath, 𝐵 The memory kernel, which is a dissipative term associated 𝑘 is the Boltzmann constant, the 𝜂 are frictional electronic with friction, is given by the expression 𝐵 𝑙 degrees of freedom, the 𝑄 are the friction coefficients associated with the electronic orbitals, the 𝜉 are frictional 𝐾 (𝑡 )=2𝑎 𝛿 (𝑡 )− a exp [− |𝑡 | A]a (59) nuclear degrees of freedom, the 𝑄 are the friction coefficients associated with the nuclei, and 𝑔 denotes the number of whereas the noise correlation, which characterised the fluc- dynamical degrees of freedom to which the nuclear thermo- tuations associated with the random noise, is provided by statchainiscoupled.Thefrictiontermsareproportional to the velocity of the corresponding degrees of freedom as it 𝐻 (𝑡 ) Š ⟨𝜁 (𝑡 )𝜁 (0)⟩ is customary. Not only must friction terms be added to the (60) =𝑑 𝛿 (𝑡 )+ a exp [− |𝑡 | A](Za − d ), 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. where the matrixes Z and D are defined as The terms 𝜂 ̇ and 𝜉 may be considered as dynamical friction 1 1 coefficients for the physical degrees of freedom. −A𝑡 −A 𝑡 Z Š ∫ 𝑒 D𝑒 𝑑𝑡, (61) In the next section, we present an alternative approach basedonthegeneralisedLangevinequation. D Š B B , (62) 𝑝 𝑝 8.3. Generalised Langevin Equation and Isothermal Processes. respectively. eTh canonical, isothermal ensemble is obtained Massive thermostating is not the sole approach to simulate by applying the u fl ctuation-dissipation theorem [35]. eTh an isothermal process. Indeed, the latter may be achieved by u fl ctuation-dissipation theorem states that the Fourier trans- means of the generalised Langevin equation [32–34]. In this formsF of the noise correlation and the memory kernel must subsection, we restrict ourselves to one generalised coordi- be related by nate in order not to clutter the notation. The generalisation to more than one coordinate is immediate. eTh generalised (F ∘𝐻 )( 𝜔 )=𝑘 𝑇 (F ∘𝐾 )( 𝜔 ). (63) Langevin equation, which is a differential stochastic equation, may be written as eTh u fl ctuation-dissipation theorem implies in turn that 𝑝 𝑝 −𝑉 (𝑞) 𝑇 𝑇 D = B B =𝑘 𝑇 (A + A ). (64) 𝑞=𝑝 ̇ [ ]=[ ]− A [ ]+ B [𝜉 ], (55) 𝑝 𝑝 𝑝 𝐵 𝑝 𝑝 𝑝 𝑝 p 0 p eTh refore, the fluctuation-dissipation theorem xe fi s the dif- where A is the drift matrix: 𝑝 fusion matrix once the drift matrix is determined. As a result, an isothermal process follows immediately. 𝑎 a In the next subsection, we address isobaric molecular A Š [ ], (56) processes, which are very common as most experiments are a A performed at atmospheric pressure. B is the diffusion matrix, 𝑞 is a generalised coordinate 8.4. Isobaric Processes. As opposed to the isothermal process, associated with an atom or a nucleus (position), 𝑝 is the corre- the isobaric process is a Hamiltonian process which means sponding generalised momentum, and p is a set of 𝑛 hidden that it may be obtained by adding extra terms to the Lagran- nonphysical momenta. eTh structure of the diffusion matrix gian [28, 36]. B is similar to the one of the corresponding drift matrix. eTh In order to model an isobaric molecular process, the random process is characterised by an uncorrelated Gaussian volume occupied by the molecular system is divided into noise: congruent parallelepipeds called Bravais cells. eTh se cells are ⟨𝜉 𝑡 𝜉 0 ⟩=𝛿 𝛿 𝑡 . () ( ) () (57) characterised by three oriented vectors which correspond to 𝑖 𝑖 𝑖𝑗 𝑝𝑝 𝑝𝑝 𝑝𝑝 𝑑𝑠 𝑘𝐾 8 Advances in Chemistry the three primitive edges spanning their volume. These vec- path integrals. Contrarily to the previous approaches, the tors are concatenated in order to form a Bravais matrix: path integral method allows for a quantum formulation which includes, in addition to the electronic degrees of h =[a | a | a ]. (65) 1 2 3 freedom, their nuclear counterpart. Such a formulation is essential for systems containing light nuclei [37, 38]. The Bravais cell is characterised by its volume and by a local metric: 9.1. Path Integral Formulation. eTh quantum path integral formulation [36, 39–42] is based on the partition function of Ω= det h, the quantum statistical canonical ensemble which is defined (66) G = h h. as the trace of the exponential of the Hamiltonian operator: Z = tr exp [−𝛽H]. (70) In order to introduce the volume as a dynamic quantity into theLagrangian,thenuclearandtheelectroniccoordinatesare The partition function describes the statistical properties of reformulated in terms of the Bravais matrix and of the scaled the molecular system. Since the operators associated with the coordinates: electrons and the nuclei do not commute, the exponential of R = hS , the Hamiltonian operator must be evaluated with the Trotter 𝐼 𝐼 (67) factorization: r = hs . 𝑖 𝑖 exp [−𝛽 (− Δ + H )] The Bravais matrix and the scale coordinates constitute 𝐼 𝑒 2𝑀 distinct degrees of freedom. In order to obtain an isobaric sys- (71) tem, the Car–Parrinello Lagrangian must be supplemented 𝑃 𝛽 ℏ 𝛽 with three additional terms which appear on the third line of = lim (exp [ ∑ Δ ] exp [− H ]) . 𝐼 𝑒 𝑃→∞ 𝑃 2𝑀 𝑃 the Lagrangian: 𝐼 KS If the completeness relation L = ∑𝜇⟨ 𝜓 ̇ hs | 𝜓 ̇ hs ⟩+𝐸 [𝜓, hS] ( ) ( ) 𝑖 𝑖 󵄨 󵄨 󵄨 󵄨 ∫ 𝑑R ∑ 󵄨 Ψ (R), R⟩⟨R,Ψ (R)󵄨 = I, 𝑘 𝑘 (72) 󵄨 󵄨 + ∑Λ (⟨𝜓 hs |𝜓 hs ⟩−𝛿 ) ( ) ( ) 𝑖 𝑗 𝑘 (68) 𝑖,𝑗 which is an identity operator, is recursively inserted into the 1 1 𝑇 𝑇 ̇ ̇ ̇ ̇ Trotter formula and if the adiabatic approximation is per- + ∑𝑀 (S GS )+ 𝑊 tr (h h)−𝑝Ω. 𝐼 𝐼 𝐼 2 2 formed: The first term represents the kinetic energy associated with ∇ Ψ =0󳨐⇒ (73) the scale coordinates. One should notice the presence of the (𝑠+1) (𝑠) ⟨Ψ (𝑠+1) (R )|Ψ(𝑠) (R )⟩ ≈ 𝛿 (𝑠+1) (𝑠) , metric or quadratic form in the inner product. eTh second 𝑘 𝑘 𝑘 𝑘 term represents the kinetic energy associated with the Bravais the partition function becomes matrix. eTh matrix norm is the square of the Frobenius norm while 𝑊 is a cfi titious mass or inertia attributed to the Bravais 𝑃 𝑁 cells. The last term is associated with the pressure pof the Z = ∑ lim ∏ [∏ ∫ ] 𝑃→∞ barostat (ambient pressure). eTh equations of motion are 𝑠=1 𝐼=1 obtained,asusual,fromtheEuler–Lagrange equations.In 𝑃 𝑁 1 1 2 (𝑠) (𝑠+1) (𝑠) particular, the Euler–Lagrange equations for the Bravais cells × exp [−𝛽 ∑ {∑ 𝑀 𝜔 (R − R ) + 𝐸 (R )}] (74) 𝐼 𝑃 𝐼 𝐼 𝑘 2 𝑃 take the form: 𝑠=1 𝐼=1 𝑃 𝑁 𝑑 𝜕L 𝜕L (𝑠) − =0. × ∏∏𝐴 𝑑R . 𝑃 𝐼 (69) 𝜕 (h) 𝜕( h) 𝑠=1 𝐼=1 𝑢V 𝑢V eTh quantities In the next section, we present a path integral formulation of ab initio molecular dynamics which allows a quantum 3/2 𝑀 𝑃 formulation of both the electronic and nuclear degrees of 𝐴 ≡( ) , 2𝜋𝛽ℏ freedom. (75) 𝜔 = 2 2 9. Path Integral Molecular Dynamics: ℏ 𝛽 Toward a Quantum Formulation of Nuclei are the amplitude and the angular frequency associated with The ab initio path integral technique is based on the formula- the quantum harmonic oscillators. eTh latter appear as a con- tion of quantum statistical mechanics in terms of Feynman sequence of the path integral formulation. eTh computational 𝑑𝑡 𝑖𝑗 𝑖𝑗 Advances in Chemistry 9 time 𝑠 is a discrete evolution parameter associated with the The complex Fourier coefficients are further expanded in evolutionofthe molecularsystem. As such,itrepresentsa terms of their real and imaginary parts: particular time-slice. eTh path integral associated with the (1) (1) u = a , partition function is a weighted sum over all possible nuclear 𝐼 𝐼 trajectories. eTh nuclear congfi uration at a particular time- (𝑃) ((𝑃+2)/2) u = a , slice 𝑠 isprovidedbytheensembleofalltheindividualnuclear 𝐼 𝐼 (80) (𝑠) configurations R at this specicfi instant. eTh weighting (2𝑠−2) (𝑠) u = Re (a ), 𝐼 𝐼 function corresponds to the exponential factor which consists of two terms: the rfi st one is related to the harmonic potential (2𝑠−1) (𝑠) u = Im (a ). 𝐼 𝐼 energy associated with the nuclei while the second is the energy associated with the electrons. As for the electronic structure, it may be obtained from As in the previous sections, the Born–Oppenheimer one of the previously introduced approaches. For instance, approximation is legitimate if the thermal energy is much for a Car–Parrinello technique formulated in terms of a smaller than the difference between the electronic ground Kohn–Sham Hamiltonian, the path integral Lagrangian in state and the excited states: normal coordinates becomes 𝐸 −𝐸 ≫𝑘 𝑇. (76) 1 PICP (𝑠) (𝑠) 𝑘 0 𝐵 ̇ ̇ L = ∑ {∑𝜇⟨ 𝜓 | 𝜓 ⟩ 𝑖 𝑖 𝑠=1 𝑖 It follows that the partition function becomes KS (𝑠) (𝑠) −𝐸 [𝜓 , (R (u)) ] 𝑃 𝑁 Z = lim ∏ [∏ ∫ ] (81) 𝑃→∞ (𝑠) (𝑠) (𝑠) 𝑠=1 𝐼=1 + ∑Λ (⟨𝜓 |𝜓 ⟩−𝛿 ) 𝑖 𝑗 𝑖,𝑗 𝑃 𝑁 1 1 2 (𝑠) (𝑠+1) (𝑠) × exp [−𝛽 ∑ {∑ 𝑀 𝜔 (R − R ) + 𝐸 (R )}] (77) 𝐼 𝑃 𝐼 𝐼 0 2 𝑃 𝑃 𝑠=1 𝐼=1 1 (𝑠) 2 1 2 󸀠 (𝑠) (𝑠) 2 (𝑠) + ∑ { ∑𝑀 (u̇ ) − ∑𝑀 𝜔 (u ) }, 𝐼 𝐼 𝐼 𝑃 𝐼 𝑃 𝑁 2 2 𝑠=1 𝐼 𝐼 (𝑠) × ∏∏𝐴 𝑑R . 𝑠=1 𝐼=1 where the normal mode masses are defined as (1) From this partition function, an extended Lagrangian may be 𝑀 =0, defined: (82) 2𝜋 (𝑠−1 ) (𝑠) 𝑀 =2𝑃(1 − cos )𝑀 , 𝑠 = 2,...,𝑃 PI 𝐼 𝐼 𝑃 𝑁 to which the cti fi tious normal mode masses are closely re- 2 2 1 1 (𝑠) 2 (𝑠) (𝑠+1) = ∑ {∑ ( (P ) − 𝑀 𝜔 (R − R ) ) lated: 𝐼 𝐼 𝑃 𝐼 𝐼 2𝑀 2 (78) 𝑠=1 𝐼=1 𝐼 (1) 𝑀 =𝑀 , 1 (83) − 𝐸 (R)} , (𝑠) 0 󸀠 (𝑠) 𝑀 =𝛾𝑀 , 𝑠 = 2,...,𝑃, 𝐼 𝐼 󸀠 where the centroid adiabaticity parameter 𝛾 belongs to the where 𝑀 are cfi titious masses or unphysical auxiliary param- interval 0<𝛾≤1. eTh se masses arefunctions ofthecompu- eters associated with the nuclei whereas their physical masses tational time. All the other quantities were defined earlier. areidentiefi dby 𝑀 . One should notice that 𝑁×𝑃 ficti- In the next subsection, we address the calculation of (𝑠) 󸀠 (𝑠) tious momenta P =𝑀 R have been formally introduced. 𝐼 𝐼 𝐼 physical observables in the path integral framework. These momenta are required in order to evaluate numerically thepathintegralwithMonteCarlotechniques[30,43].The 9.2. Physical Observables. Anobservable [6,7,45]is ady- reader should notice that the momenta aeff ct neither the namic quantity that may be measured experimentally such as partition function (up to an irrelevant normalisation factor) 𝑠 the average energy or the heat capacity. Numerous observ- nor the physical observables. The ground state energy 𝐸 (R) ables may be obtained directly from the partition func- must be evaluated concurrently with the nuclear energy for tion, for instance, each time-slice. The nuclear coordinates are not linearly independent. (i) the expectation of the energy (average energy): Nevertheless, they may become linearly independent if they 𝜕 lnZ are expressed in terms of their normal modes. eTh normal ⟨𝐸 ⟩=− (84) decomposition [39, 41, 44] is obtained by representing each coordinate in terms of complex Fourier series: (ii) the variance of the energy or energy uc fl tuation: (𝑠) (𝑘) 2𝜋𝑖(𝑠−1)(𝑘−1)/𝑃 𝜕 Z R 󳨀→ ∑a 𝑒 . (79) ⟨(Δ𝐸 )⟩= (85) 𝐼 𝐼 𝑘=1 𝜕𝛽 𝜕𝛽 𝑖𝑗 𝑖𝑗 10 Advances in Chemistry (iii) the heat or thermal capacity: 𝑁 1 (𝑠) (𝑠) (𝑠) 1 (𝑠) (𝑠) 𝑄 𝜂 ̈ =2[𝜇 ∑ ⟨𝜓 |𝜓 ⟩− ]−𝑄 𝜂 ̇ 𝜂 ̇ , 𝑒 1 𝑖 𝑖 𝑒 1 2 2𝛽 𝐶 = ⟨(Δ𝐸 )⟩ (86) 𝑘 𝑇 𝑙 (𝑠) 𝑙−1 (𝑠) ̈ ̇ 𝑄 𝜂 =[𝑄 (𝜂 ) − ] 𝑒 𝑙 𝑒 𝑙−1 𝑙 (𝑠) (𝑠) (iv) the entropy: −𝑄 𝜂 ̇ 𝜂 ̇ (1 − 𝛿 ), 2𝐿 𝑒 𝑙 𝑙+1 𝑙 = 2,...,𝐿; 𝑠 = 1,...,𝑃 (87) 𝑆= (𝑘 𝑇 lnZ) 𝜕𝑇 (90) whereas the nuclear equations of motion transform into (v) the Helmholtz free energy: KS (𝑠) (𝑠) [𝜓 , R ] (𝑠) 𝐴=−𝑘 𝑇 lnZ (88) 𝑀 R =− 𝐼 𝐼 (𝑠) 𝜕R 2 (𝑠) (𝑠+1) (𝑠−1) among others. In addition, if a small perturbation is applied to −𝑀 𝜔 (2R − R − R ) 𝑃 𝐼 𝐼 𝐼 a molecular system, the expectation of the energy associated 󸀠 (𝑠) (𝑠) with this perturbation is ̇ −𝑀 𝜉 R , 𝐼 1 𝐼 (91) 𝐸=𝐸 +𝜆𝑊⇒ 󳨐 1 (𝑠) 󸀠 (𝑠) 1 (𝑠) (𝑠) 0 ̈ ̇ ̇ 𝑄 𝜉 =[∑𝑀 (R ) − ]−𝑄 𝜉 𝜉 , 𝑛 1 𝐼 𝐼 𝑛 1 2 (89) 𝐼 1 𝜕 ⟨𝑊 ⟩=− lnZ, 𝛽 2 1 𝑘 (𝑠) 𝑘−1 (𝑠) 𝑘 (𝑠) (𝑠) ̈ ̇ ̇ ̇ 𝑄 𝜉 =[𝑄 (𝜉 ) − ]−𝑄 𝜉 𝜉 (1 − 𝛿 ), 𝑛 𝑘 𝑛 𝑘−1 𝑛 𝑘 𝑘+1 where 𝜆 is a small running parameter. 𝑘 = 2,...,𝐾 ; 𝑠 = 1,...,𝑃. In the next subsection, we present an isothermal formu- lation of path integrals. eTh quantities appearing in these equations are all den fi ed in Section 8.2. One should notice that number of equations is 9.3. Car–Parrinello Path Integral Molecular Dynamics and quite large due to the evolution parameter swhich is absent Massive Thermostating. Isothermal processes are essential for from the standard Car–Parrinello formulation as reported in two reasons. eTh rfi st one, which is rather evident, is that the Section 9.3. simulation of realistic physical conditions often involves their In the next section, we address some important imple- simulation [34]. eTh second reason must be found in a rather mentational issues. subtle shortcoming of the formalism. eTh massive thermostating of the path integral [34] is 10. Computational Implementation a “sine qua non” condition in order to obtain physically realistic results. Indeed, the harmonic potential leads to In this section, we present some important implementational inefficient and nonergodic dynamics when microcanonical considerations. In particular, we approximate the electronic trajectories are used to generate ensemble averages. interaction with an effective pseudopotential. In addition, In contrast, the thermostats generate ergodic, canonical we demonstrate that the computational complexity may be averages at the expense of introducing sets of auxiliary chain reduced if the orbitals are expressed in terms of a suitable variables which add to the complexity of the calculation, but functional basis. While ab initio molecular dynamics is which are nevertheless required for its precision and physical restricted to small molecules, because of its computational trustworthiness. complexity, the method may be extended to larger systems if As we saw in Section 8.2, it is not possible to generate an a hybrid approach is followed. The latter is introduced in the isothermal process simply by adding extra terms to the La- last subsection. grangian. Rather, the equations of motion, which are obtained from the Euler–Lagrange equations, must be modified 10.1. Pseudopotential. For many calculations, the complete accordingly. eTh refore, friction terms must be added recur- knowledge of the electronic interaction is not essential for sively to the various degrees of freedom. As a result, the electronic equations of motion become the required precision. Consequently, for the sake of com- putational efficiency, the exact electronic potential may be KS (𝑠) (𝑠) approximated by means of an eeff ctive potential, called the [𝜓 , R ] 󵄨 󵄨 (𝑠) (𝑠) (𝑠) 󵄨 󵄨 󵄨 ̈ 󵄨 pseudopotential [46–48]. Moreover, the relativistic effects 𝜇 𝜓 ⟩= − + ∑Λ 𝜓 ⟩ 󵄨 ∗ 󵄨 𝑖 𝑖 󵄨 󵄨 (𝑠) 𝛿(𝜓 ) 𝑗 associated with the core electrons may be implicitly incorpo- rated into the pseudopotential, without recourse to explicit (𝑠) (𝑠) ̇ 󵄨 −𝜇 𝜂 𝜓 ⟩ and intricate approaches. 1 𝑖 𝑖𝑗 𝛿𝐸 𝑘𝐾 𝜕𝜆 𝜕𝐸 Advances in Chemistry 11 A commonly employed pseudopotential, the dual-space projected on a basis, it is entirely determined by its projection Gaussian pseudopotential [46, 47], is composed of two parts, coefficients. us, Th the resulting representation is parametric. a local potential and a nonlocal potential: As a result, the determination of the projection coefficients is equivalent to the determination of the orbitals: the closer the PP 󸀠 𝐿 NL 󸀠 𝑉 (r,r )=𝑉 (𝑟 )+𝑉 (r,r ). (92) basis functions are to the real solution, the more efficient the calculation is. This is due to the fact that fewer coefficients The local potential, which described the local interaction, is are required to adequately represent the underlying orbital. defined as An even more realistic representation may be obtained if the basis functions are centred upon their respective nuclei [13]: 𝑍 𝑟 1 𝑟 𝐿 ion 𝑉 (𝑟 ) Š − erf ( )+ exp (− ( ) ) 𝑟 √ 2 𝑟 2𝑟 𝑓 (r)= ∑𝑐 𝜙 (r− R ), 𝑖 𝑖𝜇 𝜇 𝐼 (99) (93) 𝜇,𝐼 2 4 6 𝑟 𝑟 𝑟 ×(𝐶 +𝐶 ( ) +𝐶 ( ) +𝐶 ( ) ), 1 2 3 4 where 𝜙 is an atomic basis function and R is the location of 𝑟 𝑟 𝑟 𝜇 𝐼 𝐿 𝐿 𝐿 a given nucleus. The orbitals may also be represented by plane waves [14]: where erf is the error function while 𝑟 , 𝐶 , 𝐶 , 𝐶 ,and 𝐶 𝐿 1 2 3 4 are adjustable empirical parameters. eTh nonlocal potential PW 𝑓 (r)= exp [𝑖k ⋅ r], (100) NL 󸀠 k V (r,r ) Ω ℓ (94) where Ω is the volume of the cell associated with the underly- 𝑚 ℓ ℓ ℓ 𝑚 ∗ 󸀠 󸀠 = ∑∑ ∑ 𝑌 (𝜃, 𝜑) 𝑝 (𝑟 )𝑉 𝑝 (𝑟 )(𝑌 ) (𝜃 ,𝜑 ) ing 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. that describes the nonlocal interaction is defined in terms of Otherwise, a large number of plane waves are required for the complex spherical harmonics functions 𝑌 (𝜃, 𝜑) and of an accurate reconstruction of the orbitals. Consequently, the the Gaussian projectors: computational burden associated with the parametric rep- 2 resentation becomes rapidly prohibitive. Nevertheless, if the −(1/2)(𝑟/𝑟 ) ℓ ℓ+2(𝑖−1) orbitals are smooth functions, one may take advantage of the 𝑝 (𝑟 )= 2𝑟 × . (95) ℓ+4𝑖−1/2 √ Fourier transform in order to reduce the complexity involved 𝑟 Γ (ℓ+ (4𝑖 − 1 )/2) in the calculations of the derivatives. For instance, the Lapla- In these equations, ℓ is the azimuthal quantum number and 𝑚 cian of a Fourier transformed function is obtained by multi- is the magnetic quantum number while Γ is the well-known plying this function by the square of the module of the wave gamma function. vector: In the next subsection, we explore the advantages of projecting theorbitalsonafunctional basis. Δ 󳨐⇒ 󳨐 𝑘 , (101) −1 10.2. Orbitals and Basis Functions. Any continuous function, 2 𝑘 󳨐 󳨐 󳨐⇒ Δ, such as an orbital, may be represented as a linear combination of basis functions: whereF is the Fourier transform. If the Fourier transform is 󵄨 calculated with the Fast Fourier Transform algorithm, the 𝜓 ⟩= ∑𝑐 𝑓 ⟩, 𝑖 𝑗 󵄨 (96) complexity of such a calculation, for a 𝑁×𝑁×𝑁 discrete, grid, reduces to where 𝑐 are the projection coefficients and |𝑓 ⟩ are the 󵄨 󵄨 󵄨 󵄨 󵄨 󵄨 𝑂 (𝑁 ) 󳨐⇒ 𝑂 ((𝑁 log 𝑁 ) ) , (102) basis functions. Such decomposition may be either physically 󵄨 󵄨 󵄨 FT 󵄨 FFT motivated, computationally motivated, or both. For instance, ¨ which explains why plane wave basis and the Fourier if physical solutions of the Schrodinger equation are known, transforms are prevalent in ab initio molecular dynamics it is possible to project the orbitals on these solutions [49– simulations. The computational efficiency may be further 51]. The most common Schr odin ¨ ger basis functions are the improved if the Fourier transformations are evaluated with Slater-type basis functions [49]: graphics processing units (GPU), as their architecture makes 𝑆 𝑆 𝑚 𝑦 𝑚 𝑥 𝑧 them particularly suited for these calculations, especially 𝑓 (r)=𝑁 𝑟 𝑟 𝑟 exp [−𝜁 r ] (97) ‖ ‖ m m 𝑥 𝑧 m regarding speed [12]. and the Gaussian-type basis functions [10, 50]: The orbital functions are not always smooth as a result of the strength of the electronic interaction. In this particular 𝐺 𝐺 𝑚 𝑦 𝑚 2 𝑥 𝑧 𝑓 (r)=𝑁 𝑟 𝑟 𝑟 exp [−𝛼 ‖r‖ ]. (98) 𝑦 case,the planewavebasis constitutesaninappropriate choice m m 𝑥 𝑧 m as a very large number of basis elements are required in 𝑆 𝐺 where 𝑁 and 𝑁 are normalisation constants while m order to describe the quasi-discontinuities. Nevertheless, one m m represents the magnetic quantum numbers. Once an orbital is would be inclined to retain the low computational complexity 𝑖𝑗 𝑖𝑗 12 Advances in Chemistry associated with plane waves and the Fourier transform while where 𝑠 and 𝑑 are the wavelet coefficients. Naturally, 𝑘 𝑘 being able to efficiently describe the rougher parts of the higher resolutions may be achieved; the maximum resolution orbitals. This may be achieved with a wavelet basis and the is determined by the resolution of the computational grid. wavelet transforms [20, 48, 52–54]. As opposed to the plane As for the Fourier transform, the wavelet transform admits wave functions, which span the whole space, the wavelets are Fast Wavelet Transform implementation which has the same spatially localised. er Th efore, they are particularly suited for complexity as its Fourier counterpart. the description of discontinuities and fast-varying functions. In the next subsection, we briefly present a hybrid ap- In addition, the wavelet basis provides a multiresolution rep- proachwhichinvolvesbothabinitiomoleculardynamicsand resentation of the orbitals which means that the calculations molecular dynamics. may be performed efficiently at the required resolution level. The wavelet functions are filter banks [55]. In one dimension, 10.3. Hybrid Quantum Mechanics–Molecular Dynamics Meth- they involve two functions: ods. Because of its inherent complexity, the applicability of ab initio molecular dynamics is essentially limited to small (i) the scaling function: molecular domains. Molecular dynamics, which we recently reviewed in [30], is suitable for larger molecular system. 𝜙 (𝑥 )= 2 ∑ ℎ 𝜙 (2𝑥−𝑖 ) (103) A drawback is that it is not adapted for the simulation of 𝑖=1−𝑚 chemical complexes, since it is not based on first principles which is a high-pass filter responsible for the multiresolution as its quantum mechanics counterpart. Rather, molecular dy- aspect of the transform; namics relies on empirical potentials [56] and classical mechanics which allow for the simulation of much larger mo- (ii) the mother wavelet: lecular complexes. Therefore, in order to simulate larger systems, hybrid 𝜅 (𝑥 )= 2 ∑ 𝑔 𝜅 (2𝑥 − 𝑖 ) (104) approaches [12, 57–59] must be followed. In such a method, a 𝑖=1−𝑚 small region of interest is simulated with ab initio molecular whichisaband-passfilterresponsibleforthebasis localisa- dynamics, while the rest of the molecular system is approxi- tion. eTh coefficients of the scaling and the mother wavelet mated with molecular mechanics: function are related by QM/MM QM MM QM-MM L = L + L + L . (109) CP 𝑔 =−1 ℎ . (105) 𝑖 −𝑖+1 An extra Lagrangian term is required in order to ensure a In three dimensions, at the lowest resolution, the scaling proper coupling between the quantum degrees of freedom function is given by and the classical degrees of freedom. This Lagrangian consists 𝑥 𝑦 𝑧 of a bounded and an unbounded part. eTh bounded part con- 𝜙 (r)=𝜙( −𝑖)𝜙( −𝑗)𝜙( −𝑘) (106) sists of the stretching, bending, and torsional terms which are 𝑎 𝑎 𝑎 characterised by their distances, angles, and dihedral angles, while the mother wavelets become respectively. eTh unbounded part contains the electrostatic 𝑥 𝑧 𝜅 (r)=𝜅( −𝑖)𝜙( −𝑗)𝜙( −𝑘), 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. 𝜅 (r)=𝜙( −𝑖)𝜅( −𝑗)𝜙( −𝑘), 𝑎 𝑎 𝑎 11. Conclusions 𝑥 𝑦 𝑧 𝜅 (r)=𝜅( −𝑖)𝜅( −𝑗)𝜙( −𝑘), 𝑎 𝑎 𝑎 In this paper, we have presented a comprehensive, but yet 𝑥 𝑧 4 concise, review of ab initio molecular dynamics from a com- 𝜅 (r)=𝜙( −𝑖)𝜙( −𝑗)𝜅( −𝑘), (107) 𝑎 𝑎 𝑎 putational perspective and from rfi st principles. Although it is always hazardous to speculate about the future, we 𝑥 𝑦 𝑧 𝜅 (r)=𝜅( −𝑖)𝜙( −𝑗)𝜅( −𝑘), foresee two important breakthroughs. Fourier transforms, 𝑎 𝑎 𝑎 which constitute an essential part of many molecular simula- 𝑥 𝑦 𝑧 tions, may be evaluated with high performance on graphical 𝜅 r =𝜙( −𝑖)𝜅( −𝑗)𝜅( −𝑘), ( ) processing units or GPU [54]. eTh same remark applies to 𝑎 𝑎 𝑎 the wavelet transform whose role is essential in describing 𝑥 𝑧 𝜅 (r)=𝜅( −𝑖)𝜅( −𝑗)𝜅( −𝑘), discontinuous orbitals [52]. This opens the door to high per- 𝑎 𝑎 𝑎 formance simulations, while tempering the limitations asso- where 𝑎 is the resolution of the underlying computational ciated with computational complexity. grid. As a result, any orbital function may be approximated For decades, one of the main bottlenecks of molecular by a truncated expansion: simulations has been the calculation of the density func- tionals. As in numerous fields, machine learning constitutes ] ] a promising avenue for the fast and efficient evaluation of 𝜓 (r)= ∑𝑠 𝜙 (r)+ ∑ ∑𝑑 𝜅 (r), (108) 𝑘 𝑘 𝑘 𝑘 these functionals without recourse to explicit calculations 𝑖,𝑗,𝑘 𝑖,𝑗,𝑘 ]=1 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 Advances in Chemistry 13 [17, 24, 25, 60]. Here, the density functional is simply learned [15] K. Kitaura and K. Morokuma, “A new energy decomposition scheme for molecular interactions within the Hartree-Fock from existing examples with machine learning techniques. approximation,” International Journal of Quantum Chemistry, This paves the way for the simulation for larger and more vol. 10, no. 2, pp. 325–340, 1976. complex molecular systems. We plan to review machine [16] H. S. Yu, X. He, and D. G. Truhlar, “MN15-L: A New Local Ex- learning-based approaches in a future publication. change-Correlation Functional for Kohn-Sham Density Func- tional eTh ory with Broad Accuracy for Atoms, Molecules, and Conflicts of Interest Solids,”JournalofChemicalTheoryandComputation ,vol.12,no. 3,pp.1280–1293,2016. eTh authors declare that there are no conflicts of interest re- [17] F.Brockherde,L.Vogt,L.Li,M.E.Tuckerman,K.Burke,andK.- garding the publication of this paper. R. Muller ¨ , “By-passing the Kohn-Sham equations with machine learning,” Nature Communications,vol.8,no. 1, articleno. 872, References [18] A. P. Selvadurai, Partial Dier ff ential Equations in Mechanics 2 , [1] P.Carloni,U.Rothlisberger, andM.Parrinello,“er Th oleand Springer Berlin Heidelberg, Berlin, Heidelberg, 2000. perspective of ab initio molecular dynamics in the study of [19] E. Engel and R. M. Dreizle, DensityFunctionalTheory:An biological systems,” Accounts of Chemical Research,vol.35, no. Advanced Course, Springer, London, UK, 2011. 6, pp. 455–464, 2002. [20] T. D. Engeness and T. A. Arias, “Multiresolution analysis for [2] L. Monticelli and E. Salonen, Eds., Biomolecular Simulations– efficient, high precision all-electron density-functional calcu- Methods and Protocols, Humana Press, Springer, London, UK, lations,” Physical Review B: Condensed Matter and Materials Physics,vol.65, no.16, pp.1–10,2002. [3] T.D.Kuhne ¨ , Ab-Initio Molecular Dynamics, 2013, https://arxiv [21] A. J. Cohen, P. Mori-Sanc ´ hez, and W. Yang, “Challenges for den- .org/abs/1201.5945. sity functional theory,” Chemical Reviews,vol.112,no.1,pp.289– [4] M. E. Tuckerman, “Ab initio molecular dynamics: Basic con- 320, 2012. cepts, current trends and novel applications,” Journal of Physics: [22] H. Ou-Yang and M. Levy, “Theorem for functional derivatives Condensed Matter, vol. 14, no. 50, pp. R1297–R1355, 2002. in density-functional theory,” Physical Review A: Atomic, Molec- [5] J.Broeckhoveand L.Lathouwers,Eds., Time-Dependent Quan- ular and Optical Physics,vol.44, no.1,pp. 54–58, 1991. tum Molecular Dynamics, NATO Advanced Research Workshop [23] P. Mori-Sanc ´ hez, A. J. Cohen, and W. Yang, “Self-interaction- on Time Dependent Quantum Molecular Dynamics: Theory and free exchange-correlation functional for thermochemistry and Experiment,Springer, Snowbird,USA,1992. kinetics,” The Journal of Chemical Physics ,vol.124,no.9,Article [6] F. Jensen, Introduction to Computational Chemistry,Wiley, ID 091102, 2006. Chichester, Uk, 2017. [24] L. Li,T.E.Baker,S.R.White,andK. Burke,“Puredensityfunc- [7] M.E.Tuckerman, Tuckerman, Statistical Mechanics: eor Th y and tional for strong correlation and the thermodynamic limit from Molecular Simulation, Oxford University Press, Oxford, UK, machine learning,” Physical Review B: Condensed Matter and Materials Physics,vol.94, no.24, ArticleID245129, 2016. [8] R. P. Steele, “Multiple-Timestep ab Initio Molecular Dynamics [25] F. Pereira,K. Xiao,D.A.R.S.Latino, C. Wu,Q.Zhang,andJ. Using an Atomic Basis Set Partitioning,” The Journal of Physical Aires-De-Sousa, “Machine Learning Methods to Predict Den- Chemistry A, vol. 119, no. 50, pp. 12119–12130, 2015. sity Functional eTh ory B3LYP Energies of HOMO and LUMO [9] N.Luehr,T.E.Markland,andT.J.Mart´ınez, “Multiple time Orbitals,” Journal of Chemical Information and Modeling,vol.57, step integrators in ab initio molecular dynamics,” The Journal no. 1, pp. 11–21, 2017. of Chemical Physics,vol.140,no.8, ArticleID084116,2014. [26] M. E. Tuckerman and M. Parrinello, “Integrating the Car-Par- [10] H.B. Schlegel,Li. X.,J.M.Millam,G.A.Voth,G.E. Scuseria, rinello equations. I. Basic integration techniques,” The Journal and M. J. Frisch, “Ab initio molecular dynamics: Propagating of Chemical Physics,vol.101,no.2, pp.1302–1315,1994. the density matrix with Gaussian orbitals. III. Comparison with [27] J. Lach, J. Goclon, and P. Rodziewicz, “Structural flexibility of BornOppenheimer dynamics,” Journal of Chemical Physics,vol. the sulfur mustard molecule at finite temperature from Car- 117,no.19,pp.9758–9763,2001. Parrinello molecular dynamics simulations,” Journal of Haz- [11] T. D. Kuhne ¨ , M. Krack, F. R. Mohamed, and M. Parrinello, “Effi- ardous Materials,vol.306,pp.269–277, 2016. cient and accurate car-parrinello-like approach to born-oppen- [28] G. Bussi, T. Zykova-Timan, and M. Parrinello, “Isothermal- heimer molecular dynamics,” Physical Review Letters, vol. 98, isobaric molecular dynamics using stochastic velocity rescal- no.6,ArticleID 066401,2007. ing,” The Journal of Chemical Physics ,vol.130,no.7, ArticleID [12] X. Andrade, A. Castro, D. Zueco et al., “Modified ehrenfest for- 074101, 2009. malism for efficient large-scale ab initio molecular dynamics,” [29] J. Lahnsteiner, G. Kresse, A. Kumar, D. D. Sarma, C. Franchini, Journal of Chemical eor Th y and Computation ,vol.5,no.4,pp. and M. Bokdam, “Room-temperature dynamic correlation be- 728–742, 2009. tween methylammonium molecules in lead-iodine based per- [13] M. Vacher, D. Mendive-Tapia, M. J. Bearpark, and M. A. Robb, ovskites: An ab initio molecular dynamics perspective,” Physical “eTh second-order Ehrenfest method: A practical CASSCF Review B: Condensed Matter and Materials Physics,vol.94, no. approach to coupled electron-nuclear dynamics,” Theoretical 21, Article ID 214114, 2016. Chemistry Accounts,vol.133,no. 7,pp.1–12, 2014. [30] E. Paquet and H. L. Viktor, “Molecular dynamics, monte carlo [14] F.Ding,J.J.Goings,H. Liu,D.B. Lingerfelt,andX. Li,“Ab initio simulations, and langevin dynamics: A computational review,” two-component Ehrenfest dynamics,” The Journal of Chemical BioMed Research International,vol.2015, ArticleID183918, Physics, vol. 143, no. 11, Article ID 114105, 2015. 2015. 14 Advances in Chemistry [31] M. E. Tuckerman, Y. Liu, G. Ciccotti, and G. J. Martyna, “Non- [48] L. Genovese, A. Neelov, S. Goedecker et al., “Daubechies wave- Hamiltonian molecular dynamics: Generalizing Hamiltonian lets as a basis set for density functional pseudopotential calcu- phase space principles to non-Hamiltonian systems,” The Jour- lations,” The Journal of Chemical Physics ,vol.129,no. 1, Article nal of Chemical Physics, vol. 115, no. 4, pp. 1678–1702, 2001. ID 014109, 2008. [49] M. A. Watson, N. C. Handy, and A. J. Cohen, “Density func- [32] M. Ceriotti, G. Bussi, and M. Parrinello, “Colored-noise ther- mostats ala ` Carte,” Journal of Chemical eTh ory and Computa- tional calculations, using Slater basis sets, with exact exchange,” The Journal of Chemical Physics ,vol.119,no. 13,pp. 6475–6481, tion,vol.6,no.4, pp.1170–1180,2010. [33] M. Ceriotti, G. Bussi, and M. Parrinello, “Langevin equation [50] X. Andrade and A. Aspuru-Guzik, “Real-space density func- with colored noise for constant-temperature molecular dynam- tional theory on graphical processing units: Computational ap- ics simulations,” Physical Review Letters,vol.102,no. 2,Article proach and comparison to Gaussian basis set methods,” Journal ID 020601, 2009. of Chemical Theory and Computation ,vol.9,no. 10,pp. 4360– [34] M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Mano- 4373, 2013. lopoulos, “Efficient stochastic thermostatting of path integral [51] S. V. Levchenko, X. Ren, J. Wieferink et al., “Hybrid functionals molecular dynamics,” eTh JournalofChemicalPhysics ,vol.133, for large periodic systems in an all-electron, numeric atom- no. 12, Article ID 124104, 2010. centered basis framework,” Computer Physics Communications, [35] R. Kubo, “eTh fluctuation-dissipation theorem,” Reports on Pro- vol. 192, pp. 60–69, 2015. gress in Physics,vol.29, no.1,articleno.306,pp. 255–284, 1966. [52] T. Deutsch and L. Genovese, “Wavelets for electronic structure [36] G. J. Martyna, A. Hughes, and M. E. Tuckerman, “Molecular dy- calculations,” Collection SFN,vol.12,pp.33–76,2011. namics algorithms for path integrals at constant pressure,” The [53] A. H. Prociuk and S. S. Iyengar, “A multiwavelet treatment of the Journal of Chemical Physics, vol. 110, no. 7, pp. 3275–3290, 1999. quantum subsystem in quantum wavepacket ab initio molecular [37] T. Spura, C. John, S. Habershon, and T. D. Kuhne ¨ , “Nuclear dynamics through an hierarchical partitioning of momentum quantum eeff cts in liquid water from path-integral simulations space,” JournalofChemicalTheoryandComputation ,vol.10,no. using an ab initio force-matching approach,” Molecular Physics, 8, pp. 2950–2963, 2014. vol.113,no.8, pp.808–822,2015. [54] L. Genovese, B. Videau, D. Caliste, J.-F. Meh ´ aut, S. Goedecker, [38] O. Marsalek and T. E. Markland, “Ab initio molecular dynamics andT.Deutsch,“Wavelet-BasedDensity FunctionalTheoryon with nuclear quantum eeff cts at classical cost: Ring polymer MassivelyParallelHybridArchitectures,” Electronic Structure contraction for density functional theory,” The Journal of Chem- Calculations on Graphics Processing Units: From Quantum ical Physics,vol.144,no. 5,ArticleID054112, 2016. Chemistry to Condensed Matter Physics, pp. 115–134, 2016. [39] M. E. Tuckerman, D. Marx, M. L. Klein et al., “Efficient and gen- [55] R. C. Gonzalez and R. E. Woods, Digital Image Processing,Pear- eral algorithms for path integral Car-Parrinello molecular dy- son Prentice Hall, Upper Saddle River, NJ, USA, 2008. namics,” The Journal of Chemical Physics ,vol.104,no.14,pp. [56] K. Vanommeslaeghe, E. Hatcher, and C. Acharya, “CHARMM 5579–5588, 1996. general force field: a force field for drug-like molecules compati- [40] P.Minary,G.J.Martyna,and M. E.Tuckerman, “Algorithms ble with the CHARMM all-atom additive biological force fields,” and novel applications based on the isokinetic ensemble. I. Bio- Journal of Computational Chemistry,vol.31,no.4,pp.671–690, physical andpathintegralmolecular dynamics,” The Journal of Chemical Physics,vol.118,no. 6,pp.2510–2526,2003. [57] J. R. Perilla, B. C. Goh, C. K. Cassidy et al., “Molecular dy- [41] M. Ceriotti, J. More, and D. E. Manolopoulos, “I-PI: A Python namics simulations of large macromolecular complexes,” Cur- interface for ab initio path integral molecular dynamics simu- rent Opinion in Structural Biology,vol.31, pp.64–74,2015. lations,” Computer Physics Communications,vol.185,no.3,pp. [58] N. Sahu and S. R. Gadre, “Molecular tailoring approach: A route 1019–1026, 2014. for ab initio treatment of large clusters,” Chemical Reviews,vol. [42] H. Y. Geng, “Accelerating ab initio path integral molecular dy- 47, no. 9, pp. 2739–2747, 2014. namics with multilevel sampling of potential surface,” Journal of [59] J. Dreyer,B.Giuseppe, E. Ippoliti,V.Genna,M.DeVivo, and Computational Physics, vol. 283, pp. 299–311, 2015. P. Carloni, “First principles methods in biology: continuum [43] K. A. Fichthorn and W. H. Weinberg, “eo Th retical foundations models to hybrid ab initio quantum mechanics/molecular of dynamical Monte Carlo simulations,” The Journal of Chemical mechanics,” in Simulating Enzyme Reactivity: Computational Physics,vol.95,no.2,pp. 1090–1096,1991. Methods in Enzyme Catalysis,I.Tun˜o´nandV.Moliner,Eds., chapter 9, Royal Society of Chemistry, London, UK, 2016. [44] D. Marx and M. Parrinello, “Ab initio path integral molecular dynamics: Basic ideas,” The Journal of Chemical Physics ,vol.104, [60] V. Botu and R. Ramprasad, “Adaptive machine learning frame- no. 11, pp. 4077–4082, 1995. work to accelerate ab initio molecular dynamics,” International JournalofQuantumChemistry,vol.115,no.16,pp.1074–1083, [45] O. Marsalek,P.-Y. Chen,R.Dupuisetal.,“Efficientcalculation of free energy differences associated with isotopic substitution using path-integral molecular dynamics,” Journal of Chemical Theory and Computation ,vol.10,no.4,pp. 1440–1453,2014. [46] S. Goedecker and M. Teter, “Separable dual-space Gaussian pseudopotentials,” Physical Review B: Condensed Matter and Materials Physics,vol.54,no.3,pp. 1703–1710,1996. [47] C. Hartwigsen, S. Goedecker, and J. Hutter, “Relativistic sep- arable dual-space Gaussian pseudopotentials from H to Rn,” Physical Review B: Condensed Matter and Materials Physics,vol. 58,no. 7,pp.3641–3662,1998. Journal of Nanomaterials Journal of International Journal of International Journal of The Scientific Analytical Methods Journal of Photoenergy in Chemistry World Journal Applied Chemistry Hindawi Hindawi Hindawi Publishing Corporation Hindawi Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 http://www www.hindawi.com .hindawi.com V Volume 2018 olume 2013 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Advances in International Journal of Physical Chemistry Medicinal Chemistry Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Submit your manuscripts at www.hindawi.com Journal of Bioinorganic Chemistry and Applications Materials Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Advances in Journal of International Journal of International Journal of BioMed Tribology Chemistry Research International Spectroscopy Electrochemistry Hindawi Hindawi Hindawi Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Enzyme Journal of Journal of International Journal of Biochemistry Analytical Chemistry Spectroscopy Nanotechnology Research Research International Hindawi Hindawi Hindawi Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Nanomaterials http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Advances in Chemistry Hindawi Publishing Corporation

Computational Methods for Ab Initio Molecular Dynamics

Advances in Chemistry , Volume 2018: 14 – Apr 29, 2018

Loading next page...
 
/lp/hindawi-publishing-corporation/computational-methods-for-ab-initio-molecular-dynamics-z8t5iLt7Xa

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Hindawi Publishing Corporation
Copyright
Copyright © 2018 Eric Paquet and Herna L. Viktor. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
ISSN
2356-6612
eISSN
2314-7571
DOI
10.1155/2018/9839641
Publisher site
See Article on Publisher Site

Abstract

Hindawi Advances in Chemistry Volume 2018, Article ID 9839641, 14 pages https://doi.org/10.1155/2018/9839641 Review Article 1 2 Eric Paquet and Herna L. Viktor National Research Council, 1200 Montreal Road, Ottawa, ON, Canada University of Ottawa, 800 King Edward Avenue, Ottawa, ON, Canada Correspondence should be addressed to Eric Paquet; eric.paquet@nrc-cnrc.gc.ca Received 23 February 2018; Accepted 20 March 2018; Published 29 April 2018 Academic Editor: Yusuf Atalay Copyright © 2018 Eric Paquet and Herna L. Viktor. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 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 aer ff ent 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. 1. Introduction important, we expose how these equations and concepts are related to one another. Ab initio molecular dynamics (AIMD) is an irreplaceable eTh paper is organised as follows. In Section 2, quantum technique for the realistic simulation of complex molecular mechanics is reviewed from a molecular dynamics perspec- systems and processes associated with biological organisms tive. Two important approximations are introduced, namely, [1, 2] such as monoclonal antibodies as illustrated in Figure 1. the adiabatic and the Born–Oppenheimer approximations. With ab initio molecular dynamics, it is possible to predict in Subsequently, in Section 3, the Ehrenfest molecular dynamics silico phenomena for which an in vivo experiment is either is presented, which allows for a semiclassical treatment of too difficult or expensive, or even currently deemed infeasible the nuclei. This opens the door, in Section 4, to the Born– [3, 4]. Ab initio molecular dynamics essentially differs from Oppenheimer molecular dynamics formulation. This is fol- molecular dynamics (MD) in two ways. Firstly, AIMD is lowed, in Section 5, by the Hartree–Fock formulation which based on the quantum Schrodin ¨ ger equation while its clas- approximates the atomic antisymmetric wave function by a single determinant of the orbitals. Sections 6 and 7 introduce sical counterpart relies on Newton’s equation. Secondly, MD relies on semiempirical eeff ctive potentials which approxi- the Kohn–Sham formulation whose primary objective is to mate quantum effects, while AIMD is based on the real physi- replace the interacting electrons by a cfi titious, but equiva- lent, system of noninteracting particles. Electrons are reintro- cal potentials. In this paper, we present a review of ab initio molecular duced as dynamical quantities, in Section 8, through the Car– dynamics from a computational perspective and from rfi st Parrinello formulation. Isothermal and isobaric processes are principles. Our main aim is to create a document which is addressed in this section. In Section 9, a path integral for- self-contained, coherent, and concise but still as complete mulation of molecular dynamics is presented, which makes as possible. As the theoretical details are essential for real- allowance for a quantum treatment of the nuclear degrees life implementation, we have provided the equations for the of freedom. Finally, in Section 10, some important compu- relevant physical aspect and approximations presented. Most tational issues are addressed such as a simplification of the 2 Advances in Chemistry ∇ Š , 𝜕R 𝜕 𝜕 Δ Š ⋅ . 𝜕R 𝜕R 𝐼 𝐼 (3) eTh gradient operator, a vector, is related to the momentum while the Laplacian operator, a scalar quantity, is associated with the kinetic energy [5–7]. eTh electronic Hamiltonian is defined as the sum of the kinetic energy of the electrons, the Figure 1: Monoclonal antibody 1F1 isolated from a 1918 influenza potential energy associated with each pair of electrons, the pandemic (Spanish Flu) survivor. potential energy associated with each pair of electron-nu- cleus, and the potential energy concomitant to each pair of nuclei: 2 2 electronic interaction with the pseudopotential, the repre- ℏ 1 𝑒 H (r, R) Š − ∑ Δ + ∑ sentation of orbitals in terms of a functional basis, the use 𝑒 𝑖 󵄩 󵄩 󵄩 󵄩 2𝑚 4𝜋𝜀 󵄩 󵄩 r − r 𝑒 0 𝑖 𝑖<𝑗 󵄩 󵄩 𝑖 𝑗 of the Fourier and wavelet transform in order to reduce the 󵄩 󵄩 (4) computational complexity, and the simulation of larger sys- 𝑒 𝑍 𝑍 1 𝑒 𝑍 1 𝐼 𝐽 tems withhybrid moleculardynamics.Thisisfollowedbya − ∑ + ∑ , 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 4𝜋𝜀 4𝜋𝜀 󵄩 󵄩 󵄩 R − R 󵄩 conclusion in Section 11. R − r 0 0 𝐼 𝐽 𝐼,𝑖 󵄩 𝐼 𝑗 󵄩 𝐼<𝐽 󵄩 󵄩 󵄩 󵄩 where 𝑚 is the electronic mass, 𝑒 is the electronic charge, 2. Quantum Molecular Dynamics 𝑍 istheatomicnumber(thenumberofprotons inagiven and the Schrödinger Equation: nucleus), and 𝜀 is the vacuum permittivity. eTh quantum A Molecular Perspective molecular system is characterised by a wave function Φ(r, R;𝑡) whose evolution is determined by the time-dependent We review some important notions of quantum mechanics Schrodin ¨ ger equation [5–7]: from a molecular perspective. In quantum mechanics [5– 7], the Hamiltonian [5–7] is the operator corresponding to 𝑖ℏ Φ (r, R;𝑡 )= H (r, R)Φ (r, R;𝑡 ), (5) 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 where𝑖= −1.TheSchr odin ¨ ger equation admits a stationary energy of the constituent particles; in occurrence the elec- or time-independent solution: trons and the nuclei H r, R Ψ r, R =𝐸 R Ψ r, R , ( ) ( ) ( ) ( ) (6) 𝑒 𝑘 𝑘 𝑘 H (r, R)=−∑ Δ + H (r, R), (1) where 𝐸 (R)is the energy associated with the electronic 𝐼 𝑒 𝑘 2𝑀 wave functionΨ (r,R)and 𝑘 is a set of quantum numbers that labelled the Eigenstates or wave functions as well as where ℏ is the Planck constant and 𝑀 is the mass of a given the Eigenvalues or energies associated with the stationary nucleus. The electronic and nuclear Cartesian coordinates for Schrodin ¨ ger equation. The total wave function Φ(r,R;𝑡) may a given particle are defined, respectively, as be expended in terms of time-dependent nuclear wave functions 𝜒 (R;𝑡) and stationary electronic wave functions Ψ (r,R): r Š [𝑟 ,𝑟 ,𝑟 ] , 𝑖 𝑖,𝑥 𝑖,𝑦 𝑖,𝑧 (2) ∞ Φ (r, R;𝑡 )= ∑𝜒 (R;𝑡 )Ψ (r, R). (7) R Š [𝑅 ,𝑅 ,𝑅 ] . 𝑘 𝑘 𝐼 𝐼,𝑥 𝐼,𝑦 𝐼,𝑧 𝑘=0 The ensemble of all electronic coordinates is represented by It should be noticed that this expansion is exact and does r while the ensemble of all nuclear coordinates is denoted by not involve any approximation. eTh electronic wave functions R. In addition, the gradient and the Laplacian operators for a are solution of the stationary Schrodin ¨ ger equation, an Eigen given electron 𝑖 and nucleus 𝐼 are defined, respectively, as equation which involves the electronic Hamiltonian. If one substitutes this expansion in the time-dependent Schrodin ¨ ger equation, one obtains a system of equations for the evolution ∇ Š , of the nuclear wave functions: 𝜕r 𝜕 𝜕 𝑖ℏ =[−∑ Δ +𝐸 (R)]𝜒 + ∑C 𝜒 . (8) Δ Š ⋅ , 𝐼 𝑘 𝑘 𝑘𝑙 𝑙 𝜕𝑡 2𝑀 𝜕r 𝜕r 𝐼 𝑖 𝑖 𝜕𝜒 Advances in Chemistry 3 As may be seen, the nuclear wave functions are coupled to Furthermore, it is assumed that the electronic wave function the electronic wave functions. eTh strength of the coupling is is in its ground or nonexcited state: determined by the coupling coefficients which form a matrix: 𝑘=0󳨐⇒ (15) C = ∫ 𝑑rΨ [−∑ Δ ]Ψ 𝑘𝑙 𝐼 𝑙 𝑘 H Ψ =𝐸 Ψ . 𝑒 0 0 0 2𝑀 (9) Thisoccursifthethermalenergyislower thantheenergy + ∑ [∫ 𝑑r Ψ [−𝑖ℏ∇ ]Ψ] [−𝑖ℏ∇ ]. 𝑘 𝐼 𝑙 𝐼 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, This system of equations is too complex to be solved directly. in which the motion of heavy nuclei is described by a semi- Consequently, in the next subsection, we consider two impor- classical equation. tant approximations of the Schrodin ¨ ger equation, namely, the adiabatic and the Born–Oppenheimer approximations, which aim to reduce such a complexity. es Th e approxima- 3. Ehrenfest Molecular Dynamics tions, as well as those that later follow, reduce substantially the duration of the calculations allowing for larger molecular Another possible approximation is to assume that the nuclear motion is semiclassicalasitisthe casewhen thenuclei systems to be simulated and longer time-scales to be explored [8, 9]. are relatively heavy. This implies that, instead of being determined by the Schrodin ¨ ger equation, the average nuclear motion is determined by Newton’s equation. The right form 2.1. Adiabatic and Born–Oppenheimer Approximations. The for this equation is given by the Ehrenfest theorem [6, 12– above-mentioned system of differential equations is too 14] which states that the potential, which governs the classical complex to be solved directly. Some approximations must be motion of the nucleons, is equal to the expectation or average, performed in order to reduce the computational complexity, in the quantum sense, of the electronic Hamiltonian with while maintaining the predictive accuracy of the calculations. respect to the electronic wave function: The first approximation, which is called the adiabatic approx- imation [5–7], assumes that the electronic wave functions 󵄨 󵄨 ̈ 󵄨 󵄨 𝑀 R (𝑡 )=−∇ ⟨Ψ󵄨 H 󵄨 Ψ⟩ = −∇ ⟨H ⟩ , (16) 𝐼 𝐼 𝐼 𝑒 𝐼 𝑒 󵄨 󵄨 Ψ adapt quasi-instantaneously to a variation of the nuclear configuration. This approximation is justified by the fact where the bra–ket notation is to be understood as that the electrons are much lighter than the nucleus. Conse- quently, such an approximation is valid, unless the motion ⟨Ψ |O| Ψ⟩ ≡ ∫ 𝑑r Ψ (r)OΨ (r) (17) of the electron becomes relativistic, which may happen if the electrons are too close to the nucleus. Mathematically, this 2 2 approximation implies that for any operator O and where R ≡𝑑 R /𝑑𝑡 .This 𝐼 𝐼 equation may be further simplified if the adiabatic and the ∇ Ψ =0. (10) 𝐼 𝑙 Born–Oppenheimer approximations are enforced: As a result, the coupling matrix becomes diagonal: 󵄨 󵄨 󵄨 󵄨 ∇ ⟨Ψ󵄨 H 󵄨 Ψ⟩⇒ 󳨐 𝐼 𝑒 󵄨 󵄨 󵄨 󵄨 󵄨 󵄨 ∇ ⟨Ψ H Ψ ⟩󳨐⇒ C =[−∑ ∫ 𝑑rΨ Δ Ψ ]𝛿 , (11) 󵄨 󵄨 (18) 𝑘𝑙 𝐼 𝑘 𝑘𝑙 𝐼 0 𝑒 0 𝑘 󵄨 󵄨 2𝑀 󵄨 󵄨 󵄨 󵄨 ⟨Ψ 󵄨 ∇ H 󵄨 Ψ ⟩. 0 𝐼 𝑒 0 󵄨 󵄨 where 𝛿 is the well-known Kronecker delta. 𝑘𝑙 The second approximation, which is called the Born– In the specific context of Ehrenfest molecular dynamics, the Oppenheimer approximation [6, 10, 11], assumes that the electrons follow a time-dependent Schrodin ¨ ger equation: electronic and the nuclear motions are separable as a result of the difference between nuclear and electronic masses. As a 𝜕 (19) 𝑖ℏ Ψ= H Ψ. result, the expansion reduces to a single Eigen state: 𝜕𝑡 Φ (r, R;𝑡 )≈Ψ (r, R)𝜒 (R;𝑡 ). (12) The fact that this equation is time-dependent ensures that the 𝑘 𝑘 orthogonality of the orbitals is maintained at all times. eTh In addition, as the energy associated with the wave function motivation is to be found in the fact that the electronic Hamil- is usually much larger than the corresponding coupling con- tonian is a Hermitian operator [6]. As stated earlier, this is not stant: thecaseintheBorn–Oppenheimer approximationinwhich the electronic wave function is governed by the stationary C ≪𝐸 , (13) 𝑘𝑘 𝑘 Schrodinger equation which does not maintain the orthonor- the time-dependent nuclear Schrodinger further reduces to mality over time. In thenext section,we apply theadiabaticandBorn– Oppenheimer approximation in the context of ab initio mo- 𝑖ℏ =[−∑ Δ +𝐸 R ]𝜒 . ( ) (14) 𝐼 𝑘 𝑘 𝜕𝑡 2𝑀 lecular dynamics. 𝜕𝜒 4 Advances in Chemistry 4. Born–Oppenheimer Molecular Dynamics for the electronic equations of motion. One should notice the presence of the Lagrange multiplier in the electronic equa- In Born–Oppenheimer molecular dynamics [6, 10, 11], it tions of motion which ensures that the orbital remains ortho- is assumed that the adiabatic and the Born–Oppenheimer normal at all time. approximations are valid and that the nuclei follow a semi- In the next section, we introduce another approximation classical Newton equation whose potential is determined in order to further reduce the complexity of the electronic by the Ehrenfest theorem. It is further assumed that the Hamiltonian. electronic wave function is in its ground state (lowest energy). Since the stationary Schrodinger equation is used for the elec- 5. Hartree–Fock Molecular Dynamics tronic degrees of freedom, the orthogonality condition must be enforced with a Lagrange multiplier as this condition is not Fromnowon,we shalladopttheatomicunitsystem: preserved by the stationary Schrodin ¨ ger equation over time. eTh orthogonality of the orbitals 𝜓 is a physical requirement 𝑖 1 𝑚 =𝑒=ℏ= =1 (25) whichmustbeenforcedatalltime in ordertobeable to 𝑒 4𝜋𝜀 compute real observable quantities. The Born–Oppenheimer molecular dynamics is characterised by a Lagrangian [6] in order to alleviate the notation. In Hartree–Fock molecular which is defined as the difference between the kinetic and the dynamics [15], the antisymmetric electronic wave function potential energy: is approximated by a single determinant of the electronic orbitals 𝜓 (r): 𝑖 𝑖 BO 2 󵄨 󵄨 ̇ 󵄨 󵄨 L = ∑𝑀 R −⟨Ψ 󵄨 H 󵄨 Ψ ⟩ 𝐼 0 𝑒 0 𝐼 󵄨 󵄨 𝐼 𝜓 (r)𝜓 (r) ⋅⋅⋅ 𝜓 (r ) 1 1 1 2 1 𝑁 (20) [ ] 𝜓 (r)𝜓 (r) ⋅⋅⋅ 𝜓 (r ) [ ] 2 1 2 2 2 𝑁 + ∑Λ (⟨𝜓 |𝜓 ⟩−𝛿 ). 𝑖 𝑗 1 [ ] Ψ= det . (26) [ ] 𝑖,𝑗 . . . [ ] 𝑁! . . . [ . . d . ] eTh rfi st term of the Lagrangian corresponds to the kinetic 𝜓 (r)𝜓 (r ) ⋅⋅⋅ 𝜓 (r ) [ ] 𝑁 1 𝑁 𝑁 𝑁 𝑁 energy of the nuclei, the second term corresponds to the nuclear potential as obtained from the Ehrenfest theorem, Two operators are associated with the electronic interaction. and the last term is a Lagrange function which ensures that eTh rfi st one, the Coulomb operator, corresponds to the the orbitals remain orthonormal at all time. eTh Lagrange electrostatic interaction between the orbitals: multipliers are denoted by Λ .Since theelectrons follow a Fermi–Dirac statistics [6], they obey the Pauli Exclusion 󸀠 ∗ 󸀠 󸀠 J r 𝜓 r ≡ [∫ 𝑑r 𝜓 (r ) 𝜓 (r )]𝜓 r . (27) ( ) ( ) ( ) Principle (two electrons cannot be in the same quantum 𝑗 𝑖 𝑗 󵄩 󵄩 𝑗 𝑖 󵄩 󸀠 󵄩 󵄩 r − r 󵄩 󵄩 󵄩 state) which means that the electronic wave function must be an antisymmetric function of its orbitals, namely, the wave The second operator, the exchange operator [6, 15], is associ- functions associated with the individual electrons. ated with the exchange energy, a quantum mechanical eect ff The equations of motion associated with this Lagrangian thatoccursasaresultofthePauliExclusionPrinciple: are obtained from the Euler–Lagrange equations [6]. eTh re isonesystemofEuler–Lagrangeequations forthe nuclear 󸀠 ∗ 󸀠 󸀠 K (r)𝜓 (r)≡[∫ 𝑑r 𝜓 (r ) 𝜓 (r )] 𝜓 (r). (28) 𝑗 𝑖 󵄩 󵄩 𝑖 𝑗 degrees of freedom: 𝑗 󵄩 󵄩 󵄩 r − r 󵄩 󵄩 󵄩 𝑑 𝜕L 𝜕L − =0, From the Coulomb and the exchange operators, an electronic (21) 𝜕 R 𝜕R 𝐼 Hamiltonian may be inferred: and one system of Euler–Lagrange equations for their elec- HF H =− Δ+𝑉 (r)+2∑J (r)− ∑K (r) tronic counterpart: ext 𝑗 𝑗 𝑗 𝑗 (29) 𝑑 𝛿L 𝛿L − =0󳨐⇒ 󵄨 󵄨 1 HF 󵄨 󵄨 𝛿⟨ 𝜓 󵄨 𝛿⟨𝜓 󵄨 ≡ Δ+𝑉 (r), 𝑖 𝑖 󵄨 󵄨 (22) 𝑑 𝛿L 𝛿L − =0. where the external potential is defined as ∗ ∗ 𝛿 𝜓 (r) (r) 𝑖 𝑖 𝑍 𝑍 𝐼 𝐽 From the Euler–Lagrange equations, one obtains 𝑉 (r) Š − ∑ + ∑ . (30) ext 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 󵄩 R − r󵄩 󵄩 R − R 󵄩 𝐼 𝐼 𝐽 𝐼 󵄩 󵄩 𝐼<𝐽 󵄩 󵄩 󵄨 󵄨 ̈ 󵄨 󵄨 𝑀 R (𝑡 )=−∇ min ⟨Ψ 󵄨 H 󵄨 Ψ ⟩ 𝐼 𝐼 𝐼 0 𝑒 0 (23) 󵄨 󵄨 This Hamiltonian determines, in turn, the electronic Lagrangian: for the nuclear equations of motion, and 󵄨 󵄨 HF HF 󵄨 󵄨 󵄨 󵄨 H 𝜓 (r)= ∑Λ 𝜓 (r) L =−⟨Ψ H Ψ ⟩+ ∑Λ (⟨𝜓 |𝜓 ⟩−𝛿 ). 󵄨 󵄨 𝑒 𝑖 𝑖 0 0 𝑖 𝑗 𝑒 𝑒 (24) 󵄨 󵄨 (31) 𝑗 𝑖,𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝛿𝜓 𝑑𝑡 𝑑𝑡 𝑑𝑡 𝑖𝑗 𝑖𝑗 𝑖𝑗 Advances in Chemistry 5 A Lagrange multiplier term is added to the Lagrangian in Some of these approximations are considered in the next sec- order to ensure that the orbitals remain orthonormal at all tion. time: a quantum mechanical requirement [6]. From the Euler– From the exchange-correlation energy, one may define Lagrange equations, one obtains the equations of motion for the exchange-correlation potential: the electrons: [𝑛 ] XC 𝑉 (r) Š (37) HF XC H 𝜓 = ∑Λ 𝜓 (r) 𝑒 𝑖 𝑖 (32) which is the functional derivative [22] of the exchange-corre- lation energy with respect to the electronic density. In addi- which differ from (24) by the Hamiltonian. tion, one may define the Kohn–Sham potential: In the next subsection, we seek to replace the interacting electrons by a cfi titious but equivalent system of noninter- 𝑉 Š 𝑉 +𝑉 +𝑉 . (38) KS 𝐻 ext XC acting particles in order to further reduce the computational complexity. 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 6. Kohn–Sham Molecular Dynamics determined by the Euler–Lagrange equations are eTh Kohn–Sham formulation [16, 17] aims to replace the KS interacting electrons by a fictitious system of noninteracting H 𝜓 = ∑Λ 𝜓 , 𝑒 𝑖 𝑗 (39) particles that generate the same electronic density as the physical system of interacting particles. eTh electronic density where the cfi titious one-particle Hamiltonian is given by is defined as KS H ≡− Δ+𝑉 . (40) 𝑛 Š ∑𝑓 ⟨𝜓 |𝜓 ⟩, 𝑖 𝑖 𝑖 𝑒 KS (33) A canonical form may be obtained if a unitary transformation where 𝑓 is the occupation number of a given orbital. is appliedtothe previous equation: In this formulation of AIMD, the Ehrenfest term is KS replaced by the Kohn–Sham energy: (41) H 𝜓 =𝜀 𝜓 . 𝑖 𝑖 𝑖 󵄨 󵄨 KS 󵄨 󵄨 In the next section, we introduce some density functionals in min ⟨Ψ 󵄨 H 󵄨 Ψ ⟩= min𝐸 [𝜓] . 0 𝑒 0 (34) 󵄨 󵄨 order to replace the interacting electrons by an equivalent but yet simpler system of noninteracting particles. The latter is define as 1 7. Exchange-Correlation Energy KS 𝐸 [𝑛 ] Š − ∑𝑓 ⟨𝜓 |Δ| 𝜓 ⟩+ ∫ 𝑑r𝑉 (r)𝑛 (r) 𝑖 𝑖 𝑖 ext 𝑖=1 The detailed analysis of density functionals, also known as (35) exchange-correlation energies [23], is beyond the scope of + ∫ 𝑑r𝑉 (r)𝑛 (r)+𝐸 [𝑛 ], 𝐻 XC 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 wherethefirst term isthetotal electronic kineticenergyasso- the density functional is to express the electronic interaction ciated with the electrons, the second term is the electrostatics interaction energy between the electronic density and the in terms of the sole electronic density. In thesimplestcase, theexchange-correlationenergyisa external potential, and the third term is the self-electrostatic interaction energy associated with the electronic density. eTh functional of the electronic density alone for which the most latter involves the interaction of the electronic density with it important representative is the local density approximation: self-created electrostatic potential. This potential, called the 1/3 3 3 LDA 4/3 Hartree potential, is obtained by solving the Poisson equation (42) 𝐸 [𝑛 ]=− ( ) ∫ 𝑑r 𝑛 . XC 4 𝜋 [18]: One may take into consideration, in addition to the electronic Δ𝑉 =−4𝑛󳨐𝜋⇒ density, the gradient of the electronic density, in which case the approach is called the generalised gradient approxima- (36) 𝑛( r ) 𝑉 r = ∫ 𝑑r . tion: ( ) 𝐻 󵄩 󵄩 󵄩 󸀠 󵄩 󵄩 r − r 󵄩 󵄩 󵄩 GGA 4/3 𝐸 [𝑛 ]= ∫ 𝑑r 𝑛 𝐹 (𝜁 ), (43) XC XC The last term is the celebrated exchange-correlation energy or density functional [16, 19–21] which takes into account the where the dimensionless reduced gradient is defined as residual electronic interaction, that is, the self-interaction of theelectrons.Unfortunately,thedensity functional hasno 1/3 ‖∇𝑛 ‖ (44) 𝜁 Š ≡ 2 (3𝜋 ) 𝑠. closed-form solution but many approximations are known. 4/3 𝑖𝑗 𝑖𝑗 𝛿𝑛 𝛿𝐸 6 Advances in Chemistry Among these approximations is the B88 approximation: Lagrangian, which is called the Car–Parrinello Lagrangian [11, 26, 27] 1/3 3 6 B88 𝜎 1 󵄨 󵄨 𝐹 (𝜁 )= ( ) + , (45) CP 2 KS 󵄨 󵄨 XC −1 󵄨 󵄨 L = ∑𝑀 R + ∑𝜇⟨ 𝜓 ̇ | 𝜓 ̇ ⟩+⟨Ψ H Ψ ⟩ 󵄨 󵄨 4 𝜋 𝐼 𝐼 𝑖 𝑖 0 𝑒 0 1+6𝜁𝛽 sinh 𝜁 󵄨 󵄨 𝜎 𝜎 𝐼 𝑖 (50) where 𝜎 refers to thespinand thePerdew, Burke, andErnzer- + ∑Λ (⟨𝜓 |𝜓 ⟩−𝛿 ) 𝑖 𝑗 hof (PBE) approximation: 𝑖,𝑗 1/3 differs from the original Kohn–Sham Lagrangian by the 3 3 PBE 𝐹 (𝑠 )= ( ) + . (46) XC second term which associates a fictitious kinetic energy to the 2 −1 4 𝜋 1+𝜇𝑠 𝜅 electronic orbitals. eTh parameter 𝜇 acts as a fictitious elec- tronic mass or inertia. As in the previous sections, the nuclear The exchange energy, which is associated with the Pauli equations of motion Exclusion Principle, may be characterised by the Hartree– Fock energy: 󵄨 󵄨 KS 󵄨 󵄨 󵄨 󵄨 𝑀 R =−∇ ⟨Ψ H Ψ ⟩+ ∑Λ ∇ ⟨𝜓 |𝜓 ⟩ 󵄨 󵄨 𝐼 𝐼 𝐼 0 𝑒 0 𝐼 𝑖 𝑗 󵄨 󵄨 (51) 𝑖,𝑗 ∗ ∗ 󸀠 󸀠 𝜓 (r)𝜓 (r)𝜓 (r)𝜓 (r ) 1 𝑗 𝑖 𝑖 𝑗 HF 󸀠 𝐸 =− ∑ ∫ 𝑑r 𝑑r (47) 󵄩 󵄩 XC 󸀠 and the electronic equations of motion 󵄩 󵄩 2 r − r 󵄩 󵄩 𝑖,𝑗 󵄩 󵄩 KS 𝜇 𝜓 ̈ 𝑡 =−H 𝜓 𝑡 + ∑Λ 𝜓 𝑡 () () () 𝑖 𝑒 𝑖 𝑗 (52) 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: may be inferred from the Euler–Lagrange equations. eTh se equations involve a second order time derivative of the B3 HF LDA B88 GGA orbitals which means that the latter are now proper dynamical 𝐸 =𝑎𝐸 + (1−𝑎 )𝐸 +𝑏Δ𝐸 +𝑐𝐸 XC XC XC XC XC quantities. (48) LDA In the next two subsections, we consider ab initio + (1−𝑐 )𝐸 , XC molecular dynamics at constant temperature and at constant pressure. where 𝑎, 𝑏,and 𝑐 are empirical parameters. Obviously, these parameters are application dependent. Recently, it has been proposed to evaluate the functional 8.2. Massive Thermostating and Isothermal Processes. In the density with machine learning techniques. The functional previous sections, we have implicitly assumed that the molec- density is learned by examples instead of directly solving the ular system under consideration was isolated. Of course, this Kohn–Sham equations [17, 24, 25]. As a result, substantially is not compatible with most experimental conditions [27–29] less time is required to complete the calculations allowing for in which the system is either kept at constant temperature larger system to be simulated and longer time-scales to be (isothermal) in a heat bath or at constant pressure (isobaric) explored. as a consequence, for instance, of the atmospheric pressure. In the next section, we further improve the precision In this subsection and in the next, we explain how to address of the calculations by reintroducing the dynamic electronic these important issues. degrees of freedom which, until now, have been absent from As we have recently explained in a computational review the equations of motion. 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 modiefi d directly as a result 8. Car–Parrinello Molecular Dynamics of the non-Hamiltonian nature of the isothermal process [31]. In order to obtain an isothermal molecular system, a friction 8.1. Equations of Motion. The Kohn–Sham dynamics, as for- term must be recursively added to each degree of freedom mulated in the previous sections, does not take the dynamics [28, 29]. Therefore, the electronic equations of motion must of the electrons into account despite the fact that it is present: be modified as follows: an unattractive unphysical feature. Indeed, one obtains from the Lagrangian and the Euler–Lagrange equations 󵄨 󵄨 󵄨 󵄨 KS 󵄨 󵄨 󵄨 ̈ 󵄨 ̇ ̇ 𝜇 󵄨 𝜓 ⟩= −H 󵄨 𝜓 ⟩+ ∑Λ 𝜓 ⟩− 𝜇 𝜂 󵄨 𝜓 ⟩, 𝑖 𝑒 𝑖 𝑗 1 𝑖 󵄨 󵄨 󵄨 𝑑 𝜕L ≡0 (49) 𝜕 𝜓 (r) 1 𝑒 1 ̈ ̇ ̇ ̇ ̇ 𝑄 𝜂 =2[𝜇 ∑ ⟨𝜓 | 𝜓 ⟩− ]−𝑄 𝜂 𝜂 , 𝑒 1 𝑖 𝑖 𝑒 1 2 2𝛽 𝑖 (53) which means the orbitals are not dynamical fields of the molecular system. Nevertheless, the Lagrangian may be 𝑙 𝑙−1 2 𝑙 𝑄 𝜂 ̈ =[𝑄 𝜂 ̇ − ]−𝑄 𝜂 ̇ 𝜂 ̇ (1 − 𝛿 ), 𝑙 𝑙 𝑙+1 𝑙𝐿 modified in order to also include the dynamic nature of 𝑒 𝑒 𝑙−1 𝑒 electronic degrees of freedom through the introduction of a cfi titious electronic kinetic energy term. The extended 𝑙 = 2, ...,𝐿 𝑑𝑡 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝜇𝑠 𝑖𝑗 𝑖𝑗 𝛽𝜁 Advances in Chemistry 7 while the nuclear equations of motion must take a similar The hidden momenta, associated with the generalised Lan- form: gevin equation, may be marginalised because the evolution of the variables [𝑝, p] , in the free particle limit, is described KS ̈ ̇ ̇ 𝑀 R =−∇ 𝐸 −𝑀 𝜉 R , 𝐼 𝐼 𝐼 𝐼 1 𝐼 by a linear Markovian stochastic differential equation. As a result,thegeneralisedLangevinequationbecomesequivalent 1 2 1 ̈ ̇ ̇ to aLangevinequationwithmemorykerneland noisecorre- 𝑄 𝜉 =[∑𝑀 R − ]−𝑄 𝜉 𝜉 , 𝑛 1 𝐼 𝐼 𝑛 1 2 lation [32–34]: (54) 𝑘 𝑘−1 2 𝑘 𝑞=𝑝 ̇ , ̈ ̇ ̇ ̇ 𝑄 𝜉 =[𝑄 𝜉 − ]−𝑄 𝜉 𝜉 (1 − 𝛿 ) 𝑘 𝑘 𝑘+1 𝑛 𝑛 𝑘−1 𝑛 (58) 𝑝=−𝑉 (𝑞) − ∫ 𝐾 (𝑡−𝑠 )𝑝 (𝑠 )+ 𝜁. 𝑘 = 2,...,𝐾, −∞ where 𝛽≡1/𝑘 𝑇, 𝑇 is the temperature of the heat bath, 𝐵 The memory kernel, which is a dissipative term associated 𝑘 is the Boltzmann constant, the 𝜂 are frictional electronic with friction, is given by the expression 𝐵 𝑙 degrees of freedom, the 𝑄 are the friction coefficients associated with the electronic orbitals, the 𝜉 are frictional 𝐾 (𝑡 )=2𝑎 𝛿 (𝑡 )− a exp [− |𝑡 | A]a (59) nuclear degrees of freedom, the 𝑄 are the friction coefficients associated with the nuclei, and 𝑔 denotes the number of whereas the noise correlation, which characterised the fluc- dynamical degrees of freedom to which the nuclear thermo- tuations associated with the random noise, is provided by statchainiscoupled.Thefrictiontermsareproportional to the velocity of the corresponding degrees of freedom as it 𝐻 (𝑡 ) Š ⟨𝜁 (𝑡 )𝜁 (0)⟩ is customary. Not only must friction terms be added to the (60) =𝑑 𝛿 (𝑡 )+ a exp [− |𝑡 | A](Za − d ), 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. where the matrixes Z and D are defined as The terms 𝜂 ̇ and 𝜉 may be considered as dynamical friction 1 1 coefficients for the physical degrees of freedom. −A𝑡 −A 𝑡 Z Š ∫ 𝑒 D𝑒 𝑑𝑡, (61) In the next section, we present an alternative approach basedonthegeneralisedLangevinequation. D Š B B , (62) 𝑝 𝑝 8.3. Generalised Langevin Equation and Isothermal Processes. respectively. eTh canonical, isothermal ensemble is obtained Massive thermostating is not the sole approach to simulate by applying the u fl ctuation-dissipation theorem [35]. eTh an isothermal process. Indeed, the latter may be achieved by u fl ctuation-dissipation theorem states that the Fourier trans- means of the generalised Langevin equation [32–34]. In this formsF of the noise correlation and the memory kernel must subsection, we restrict ourselves to one generalised coordi- be related by nate in order not to clutter the notation. The generalisation to more than one coordinate is immediate. eTh generalised (F ∘𝐻 )( 𝜔 )=𝑘 𝑇 (F ∘𝐾 )( 𝜔 ). (63) Langevin equation, which is a differential stochastic equation, may be written as eTh u fl ctuation-dissipation theorem implies in turn that 𝑝 𝑝 −𝑉 (𝑞) 𝑇 𝑇 D = B B =𝑘 𝑇 (A + A ). (64) 𝑞=𝑝 ̇ [ ]=[ ]− A [ ]+ B [𝜉 ], (55) 𝑝 𝑝 𝑝 𝐵 𝑝 𝑝 𝑝 𝑝 p 0 p eTh refore, the fluctuation-dissipation theorem xe fi s the dif- where A is the drift matrix: 𝑝 fusion matrix once the drift matrix is determined. As a result, an isothermal process follows immediately. 𝑎 a In the next subsection, we address isobaric molecular A Š [ ], (56) processes, which are very common as most experiments are a A performed at atmospheric pressure. B is the diffusion matrix, 𝑞 is a generalised coordinate 8.4. Isobaric Processes. As opposed to the isothermal process, associated with an atom or a nucleus (position), 𝑝 is the corre- the isobaric process is a Hamiltonian process which means sponding generalised momentum, and p is a set of 𝑛 hidden that it may be obtained by adding extra terms to the Lagran- nonphysical momenta. eTh structure of the diffusion matrix gian [28, 36]. B is similar to the one of the corresponding drift matrix. eTh In order to model an isobaric molecular process, the random process is characterised by an uncorrelated Gaussian volume occupied by the molecular system is divided into noise: congruent parallelepipeds called Bravais cells. eTh se cells are ⟨𝜉 𝑡 𝜉 0 ⟩=𝛿 𝛿 𝑡 . () ( ) () (57) characterised by three oriented vectors which correspond to 𝑖 𝑖 𝑖𝑗 𝑝𝑝 𝑝𝑝 𝑝𝑝 𝑑𝑠 𝑘𝐾 8 Advances in Chemistry the three primitive edges spanning their volume. These vec- path integrals. Contrarily to the previous approaches, the tors are concatenated in order to form a Bravais matrix: path integral method allows for a quantum formulation which includes, in addition to the electronic degrees of h =[a | a | a ]. (65) 1 2 3 freedom, their nuclear counterpart. Such a formulation is essential for systems containing light nuclei [37, 38]. The Bravais cell is characterised by its volume and by a local metric: 9.1. Path Integral Formulation. eTh quantum path integral formulation [36, 39–42] is based on the partition function of Ω= det h, the quantum statistical canonical ensemble which is defined (66) G = h h. as the trace of the exponential of the Hamiltonian operator: Z = tr exp [−𝛽H]. (70) In order to introduce the volume as a dynamic quantity into theLagrangian,thenuclearandtheelectroniccoordinatesare The partition function describes the statistical properties of reformulated in terms of the Bravais matrix and of the scaled the molecular system. Since the operators associated with the coordinates: electrons and the nuclei do not commute, the exponential of R = hS , the Hamiltonian operator must be evaluated with the Trotter 𝐼 𝐼 (67) factorization: r = hs . 𝑖 𝑖 exp [−𝛽 (− Δ + H )] The Bravais matrix and the scale coordinates constitute 𝐼 𝑒 2𝑀 distinct degrees of freedom. In order to obtain an isobaric sys- (71) tem, the Car–Parrinello Lagrangian must be supplemented 𝑃 𝛽 ℏ 𝛽 with three additional terms which appear on the third line of = lim (exp [ ∑ Δ ] exp [− H ]) . 𝐼 𝑒 𝑃→∞ 𝑃 2𝑀 𝑃 the Lagrangian: 𝐼 KS If the completeness relation L = ∑𝜇⟨ 𝜓 ̇ hs | 𝜓 ̇ hs ⟩+𝐸 [𝜓, hS] ( ) ( ) 𝑖 𝑖 󵄨 󵄨 󵄨 󵄨 ∫ 𝑑R ∑ 󵄨 Ψ (R), R⟩⟨R,Ψ (R)󵄨 = I, 𝑘 𝑘 (72) 󵄨 󵄨 + ∑Λ (⟨𝜓 hs |𝜓 hs ⟩−𝛿 ) ( ) ( ) 𝑖 𝑗 𝑘 (68) 𝑖,𝑗 which is an identity operator, is recursively inserted into the 1 1 𝑇 𝑇 ̇ ̇ ̇ ̇ Trotter formula and if the adiabatic approximation is per- + ∑𝑀 (S GS )+ 𝑊 tr (h h)−𝑝Ω. 𝐼 𝐼 𝐼 2 2 formed: The first term represents the kinetic energy associated with ∇ Ψ =0󳨐⇒ (73) the scale coordinates. One should notice the presence of the (𝑠+1) (𝑠) ⟨Ψ (𝑠+1) (R )|Ψ(𝑠) (R )⟩ ≈ 𝛿 (𝑠+1) (𝑠) , metric or quadratic form in the inner product. eTh second 𝑘 𝑘 𝑘 𝑘 term represents the kinetic energy associated with the Bravais the partition function becomes matrix. eTh matrix norm is the square of the Frobenius norm while 𝑊 is a cfi titious mass or inertia attributed to the Bravais 𝑃 𝑁 cells. The last term is associated with the pressure pof the Z = ∑ lim ∏ [∏ ∫ ] 𝑃→∞ barostat (ambient pressure). eTh equations of motion are 𝑠=1 𝐼=1 obtained,asusual,fromtheEuler–Lagrange equations.In 𝑃 𝑁 1 1 2 (𝑠) (𝑠+1) (𝑠) particular, the Euler–Lagrange equations for the Bravais cells × exp [−𝛽 ∑ {∑ 𝑀 𝜔 (R − R ) + 𝐸 (R )}] (74) 𝐼 𝑃 𝐼 𝐼 𝑘 2 𝑃 take the form: 𝑠=1 𝐼=1 𝑃 𝑁 𝑑 𝜕L 𝜕L (𝑠) − =0. × ∏∏𝐴 𝑑R . 𝑃 𝐼 (69) 𝜕 (h) 𝜕( h) 𝑠=1 𝐼=1 𝑢V 𝑢V eTh quantities In the next section, we present a path integral formulation of ab initio molecular dynamics which allows a quantum 3/2 𝑀 𝑃 formulation of both the electronic and nuclear degrees of 𝐴 ≡( ) , 2𝜋𝛽ℏ freedom. (75) 𝜔 = 2 2 9. Path Integral Molecular Dynamics: ℏ 𝛽 Toward a Quantum Formulation of Nuclei are the amplitude and the angular frequency associated with The ab initio path integral technique is based on the formula- the quantum harmonic oscillators. eTh latter appear as a con- tion of quantum statistical mechanics in terms of Feynman sequence of the path integral formulation. eTh computational 𝑑𝑡 𝑖𝑗 𝑖𝑗 Advances in Chemistry 9 time 𝑠 is a discrete evolution parameter associated with the The complex Fourier coefficients are further expanded in evolutionofthe molecularsystem. As such,itrepresentsa terms of their real and imaginary parts: particular time-slice. eTh path integral associated with the (1) (1) u = a , partition function is a weighted sum over all possible nuclear 𝐼 𝐼 trajectories. eTh nuclear congfi uration at a particular time- (𝑃) ((𝑃+2)/2) u = a , slice 𝑠 isprovidedbytheensembleofalltheindividualnuclear 𝐼 𝐼 (80) (𝑠) configurations R at this specicfi instant. eTh weighting (2𝑠−2) (𝑠) u = Re (a ), 𝐼 𝐼 function corresponds to the exponential factor which consists of two terms: the rfi st one is related to the harmonic potential (2𝑠−1) (𝑠) u = Im (a ). 𝐼 𝐼 energy associated with the nuclei while the second is the energy associated with the electrons. As for the electronic structure, it may be obtained from As in the previous sections, the Born–Oppenheimer one of the previously introduced approaches. For instance, approximation is legitimate if the thermal energy is much for a Car–Parrinello technique formulated in terms of a smaller than the difference between the electronic ground Kohn–Sham Hamiltonian, the path integral Lagrangian in state and the excited states: normal coordinates becomes 𝐸 −𝐸 ≫𝑘 𝑇. (76) 1 PICP (𝑠) (𝑠) 𝑘 0 𝐵 ̇ ̇ L = ∑ {∑𝜇⟨ 𝜓 | 𝜓 ⟩ 𝑖 𝑖 𝑠=1 𝑖 It follows that the partition function becomes KS (𝑠) (𝑠) −𝐸 [𝜓 , (R (u)) ] 𝑃 𝑁 Z = lim ∏ [∏ ∫ ] (81) 𝑃→∞ (𝑠) (𝑠) (𝑠) 𝑠=1 𝐼=1 + ∑Λ (⟨𝜓 |𝜓 ⟩−𝛿 ) 𝑖 𝑗 𝑖,𝑗 𝑃 𝑁 1 1 2 (𝑠) (𝑠+1) (𝑠) × exp [−𝛽 ∑ {∑ 𝑀 𝜔 (R − R ) + 𝐸 (R )}] (77) 𝐼 𝑃 𝐼 𝐼 0 2 𝑃 𝑃 𝑠=1 𝐼=1 1 (𝑠) 2 1 2 󸀠 (𝑠) (𝑠) 2 (𝑠) + ∑ { ∑𝑀 (u̇ ) − ∑𝑀 𝜔 (u ) }, 𝐼 𝐼 𝐼 𝑃 𝐼 𝑃 𝑁 2 2 𝑠=1 𝐼 𝐼 (𝑠) × ∏∏𝐴 𝑑R . 𝑠=1 𝐼=1 where the normal mode masses are defined as (1) From this partition function, an extended Lagrangian may be 𝑀 =0, defined: (82) 2𝜋 (𝑠−1 ) (𝑠) 𝑀 =2𝑃(1 − cos )𝑀 , 𝑠 = 2,...,𝑃 PI 𝐼 𝐼 𝑃 𝑁 to which the cti fi tious normal mode masses are closely re- 2 2 1 1 (𝑠) 2 (𝑠) (𝑠+1) = ∑ {∑ ( (P ) − 𝑀 𝜔 (R − R ) ) lated: 𝐼 𝐼 𝑃 𝐼 𝐼 2𝑀 2 (78) 𝑠=1 𝐼=1 𝐼 (1) 𝑀 =𝑀 , 1 (83) − 𝐸 (R)} , (𝑠) 0 󸀠 (𝑠) 𝑀 =𝛾𝑀 , 𝑠 = 2,...,𝑃, 𝐼 𝐼 󸀠 where the centroid adiabaticity parameter 𝛾 belongs to the where 𝑀 are cfi titious masses or unphysical auxiliary param- interval 0<𝛾≤1. eTh se masses arefunctions ofthecompu- eters associated with the nuclei whereas their physical masses tational time. All the other quantities were defined earlier. areidentiefi dby 𝑀 . One should notice that 𝑁×𝑃 ficti- In the next subsection, we address the calculation of (𝑠) 󸀠 (𝑠) tious momenta P =𝑀 R have been formally introduced. 𝐼 𝐼 𝐼 physical observables in the path integral framework. These momenta are required in order to evaluate numerically thepathintegralwithMonteCarlotechniques[30,43].The 9.2. Physical Observables. Anobservable [6,7,45]is ady- reader should notice that the momenta aeff ct neither the namic quantity that may be measured experimentally such as partition function (up to an irrelevant normalisation factor) 𝑠 the average energy or the heat capacity. Numerous observ- nor the physical observables. The ground state energy 𝐸 (R) ables may be obtained directly from the partition func- must be evaluated concurrently with the nuclear energy for tion, for instance, each time-slice. The nuclear coordinates are not linearly independent. (i) the expectation of the energy (average energy): Nevertheless, they may become linearly independent if they 𝜕 lnZ are expressed in terms of their normal modes. eTh normal ⟨𝐸 ⟩=− (84) decomposition [39, 41, 44] is obtained by representing each coordinate in terms of complex Fourier series: (ii) the variance of the energy or energy uc fl tuation: (𝑠) (𝑘) 2𝜋𝑖(𝑠−1)(𝑘−1)/𝑃 𝜕 Z R 󳨀→ ∑a 𝑒 . (79) ⟨(Δ𝐸 )⟩= (85) 𝐼 𝐼 𝑘=1 𝜕𝛽 𝜕𝛽 𝑖𝑗 𝑖𝑗 10 Advances in Chemistry (iii) the heat or thermal capacity: 𝑁 1 (𝑠) (𝑠) (𝑠) 1 (𝑠) (𝑠) 𝑄 𝜂 ̈ =2[𝜇 ∑ ⟨𝜓 |𝜓 ⟩− ]−𝑄 𝜂 ̇ 𝜂 ̇ , 𝑒 1 𝑖 𝑖 𝑒 1 2 2𝛽 𝐶 = ⟨(Δ𝐸 )⟩ (86) 𝑘 𝑇 𝑙 (𝑠) 𝑙−1 (𝑠) ̈ ̇ 𝑄 𝜂 =[𝑄 (𝜂 ) − ] 𝑒 𝑙 𝑒 𝑙−1 𝑙 (𝑠) (𝑠) (iv) the entropy: −𝑄 𝜂 ̇ 𝜂 ̇ (1 − 𝛿 ), 2𝐿 𝑒 𝑙 𝑙+1 𝑙 = 2,...,𝐿; 𝑠 = 1,...,𝑃 (87) 𝑆= (𝑘 𝑇 lnZ) 𝜕𝑇 (90) whereas the nuclear equations of motion transform into (v) the Helmholtz free energy: KS (𝑠) (𝑠) [𝜓 , R ] (𝑠) 𝐴=−𝑘 𝑇 lnZ (88) 𝑀 R =− 𝐼 𝐼 (𝑠) 𝜕R 2 (𝑠) (𝑠+1) (𝑠−1) among others. In addition, if a small perturbation is applied to −𝑀 𝜔 (2R − R − R ) 𝑃 𝐼 𝐼 𝐼 a molecular system, the expectation of the energy associated 󸀠 (𝑠) (𝑠) with this perturbation is ̇ −𝑀 𝜉 R , 𝐼 1 𝐼 (91) 𝐸=𝐸 +𝜆𝑊⇒ 󳨐 1 (𝑠) 󸀠 (𝑠) 1 (𝑠) (𝑠) 0 ̈ ̇ ̇ 𝑄 𝜉 =[∑𝑀 (R ) − ]−𝑄 𝜉 𝜉 , 𝑛 1 𝐼 𝐼 𝑛 1 2 (89) 𝐼 1 𝜕 ⟨𝑊 ⟩=− lnZ, 𝛽 2 1 𝑘 (𝑠) 𝑘−1 (𝑠) 𝑘 (𝑠) (𝑠) ̈ ̇ ̇ ̇ 𝑄 𝜉 =[𝑄 (𝜉 ) − ]−𝑄 𝜉 𝜉 (1 − 𝛿 ), 𝑛 𝑘 𝑛 𝑘−1 𝑛 𝑘 𝑘+1 where 𝜆 is a small running parameter. 𝑘 = 2,...,𝐾 ; 𝑠 = 1,...,𝑃. In the next subsection, we present an isothermal formu- lation of path integrals. eTh quantities appearing in these equations are all den fi ed in Section 8.2. One should notice that number of equations is 9.3. Car–Parrinello Path Integral Molecular Dynamics and quite large due to the evolution parameter swhich is absent Massive Thermostating. Isothermal processes are essential for from the standard Car–Parrinello formulation as reported in two reasons. eTh rfi st one, which is rather evident, is that the Section 9.3. simulation of realistic physical conditions often involves their In the next section, we address some important imple- simulation [34]. eTh second reason must be found in a rather mentational issues. subtle shortcoming of the formalism. eTh massive thermostating of the path integral [34] is 10. Computational Implementation a “sine qua non” condition in order to obtain physically realistic results. Indeed, the harmonic potential leads to In this section, we present some important implementational inefficient and nonergodic dynamics when microcanonical considerations. In particular, we approximate the electronic trajectories are used to generate ensemble averages. interaction with an effective pseudopotential. In addition, In contrast, the thermostats generate ergodic, canonical we demonstrate that the computational complexity may be averages at the expense of introducing sets of auxiliary chain reduced if the orbitals are expressed in terms of a suitable variables which add to the complexity of the calculation, but functional basis. While ab initio molecular dynamics is which are nevertheless required for its precision and physical restricted to small molecules, because of its computational trustworthiness. complexity, the method may be extended to larger systems if As we saw in Section 8.2, it is not possible to generate an a hybrid approach is followed. The latter is introduced in the isothermal process simply by adding extra terms to the La- last subsection. grangian. Rather, the equations of motion, which are obtained from the Euler–Lagrange equations, must be modified 10.1. Pseudopotential. For many calculations, the complete accordingly. eTh refore, friction terms must be added recur- knowledge of the electronic interaction is not essential for sively to the various degrees of freedom. As a result, the electronic equations of motion become the required precision. Consequently, for the sake of com- putational efficiency, the exact electronic potential may be KS (𝑠) (𝑠) approximated by means of an eeff ctive potential, called the [𝜓 , R ] 󵄨 󵄨 (𝑠) (𝑠) (𝑠) 󵄨 󵄨 󵄨 ̈ 󵄨 pseudopotential [46–48]. Moreover, the relativistic effects 𝜇 𝜓 ⟩= − + ∑Λ 𝜓 ⟩ 󵄨 ∗ 󵄨 𝑖 𝑖 󵄨 󵄨 (𝑠) 𝛿(𝜓 ) 𝑗 associated with the core electrons may be implicitly incorpo- rated into the pseudopotential, without recourse to explicit (𝑠) (𝑠) ̇ 󵄨 −𝜇 𝜂 𝜓 ⟩ and intricate approaches. 1 𝑖 𝑖𝑗 𝛿𝐸 𝑘𝐾 𝜕𝜆 𝜕𝐸 Advances in Chemistry 11 A commonly employed pseudopotential, the dual-space projected on a basis, it is entirely determined by its projection Gaussian pseudopotential [46, 47], is composed of two parts, coefficients. us, Th the resulting representation is parametric. a local potential and a nonlocal potential: As a result, the determination of the projection coefficients is equivalent to the determination of the orbitals: the closer the PP 󸀠 𝐿 NL 󸀠 𝑉 (r,r )=𝑉 (𝑟 )+𝑉 (r,r ). (92) basis functions are to the real solution, the more efficient the calculation is. This is due to the fact that fewer coefficients The local potential, which described the local interaction, is are required to adequately represent the underlying orbital. defined as An even more realistic representation may be obtained if the basis functions are centred upon their respective nuclei [13]: 𝑍 𝑟 1 𝑟 𝐿 ion 𝑉 (𝑟 ) Š − erf ( )+ exp (− ( ) ) 𝑟 √ 2 𝑟 2𝑟 𝑓 (r)= ∑𝑐 𝜙 (r− R ), 𝑖 𝑖𝜇 𝜇 𝐼 (99) (93) 𝜇,𝐼 2 4 6 𝑟 𝑟 𝑟 ×(𝐶 +𝐶 ( ) +𝐶 ( ) +𝐶 ( ) ), 1 2 3 4 where 𝜙 is an atomic basis function and R is the location of 𝑟 𝑟 𝑟 𝜇 𝐼 𝐿 𝐿 𝐿 a given nucleus. The orbitals may also be represented by plane waves [14]: where erf is the error function while 𝑟 , 𝐶 , 𝐶 , 𝐶 ,and 𝐶 𝐿 1 2 3 4 are adjustable empirical parameters. eTh nonlocal potential PW 𝑓 (r)= exp [𝑖k ⋅ r], (100) NL 󸀠 k V (r,r ) Ω ℓ (94) where Ω is the volume of the cell associated with the underly- 𝑚 ℓ ℓ ℓ 𝑚 ∗ 󸀠 󸀠 = ∑∑ ∑ 𝑌 (𝜃, 𝜑) 𝑝 (𝑟 )𝑉 𝑝 (𝑟 )(𝑌 ) (𝜃 ,𝜑 ) ing 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. that describes the nonlocal interaction is defined in terms of Otherwise, a large number of plane waves are required for the complex spherical harmonics functions 𝑌 (𝜃, 𝜑) and of an accurate reconstruction of the orbitals. Consequently, the the Gaussian projectors: computational burden associated with the parametric rep- 2 resentation becomes rapidly prohibitive. Nevertheless, if the −(1/2)(𝑟/𝑟 ) ℓ ℓ+2(𝑖−1) orbitals are smooth functions, one may take advantage of the 𝑝 (𝑟 )= 2𝑟 × . (95) ℓ+4𝑖−1/2 √ Fourier transform in order to reduce the complexity involved 𝑟 Γ (ℓ+ (4𝑖 − 1 )/2) in the calculations of the derivatives. For instance, the Lapla- In these equations, ℓ is the azimuthal quantum number and 𝑚 cian of a Fourier transformed function is obtained by multi- is the magnetic quantum number while Γ is the well-known plying this function by the square of the module of the wave gamma function. vector: In the next subsection, we explore the advantages of projecting theorbitalsonafunctional basis. Δ 󳨐⇒ 󳨐 𝑘 , (101) −1 10.2. Orbitals and Basis Functions. Any continuous function, 2 𝑘 󳨐 󳨐 󳨐⇒ Δ, such as an orbital, may be represented as a linear combination of basis functions: whereF is the Fourier transform. If the Fourier transform is 󵄨 calculated with the Fast Fourier Transform algorithm, the 𝜓 ⟩= ∑𝑐 𝑓 ⟩, 𝑖 𝑗 󵄨 (96) complexity of such a calculation, for a 𝑁×𝑁×𝑁 discrete, grid, reduces to where 𝑐 are the projection coefficients and |𝑓 ⟩ are the 󵄨 󵄨 󵄨 󵄨 󵄨 󵄨 𝑂 (𝑁 ) 󳨐⇒ 𝑂 ((𝑁 log 𝑁 ) ) , (102) basis functions. Such decomposition may be either physically 󵄨 󵄨 󵄨 FT 󵄨 FFT motivated, computationally motivated, or both. For instance, ¨ which explains why plane wave basis and the Fourier if physical solutions of the Schrodinger equation are known, transforms are prevalent in ab initio molecular dynamics it is possible to project the orbitals on these solutions [49– simulations. The computational efficiency may be further 51]. The most common Schr odin ¨ ger basis functions are the improved if the Fourier transformations are evaluated with Slater-type basis functions [49]: graphics processing units (GPU), as their architecture makes 𝑆 𝑆 𝑚 𝑦 𝑚 𝑥 𝑧 them particularly suited for these calculations, especially 𝑓 (r)=𝑁 𝑟 𝑟 𝑟 exp [−𝜁 r ] (97) ‖ ‖ m m 𝑥 𝑧 m regarding speed [12]. and the Gaussian-type basis functions [10, 50]: The orbital functions are not always smooth as a result of the strength of the electronic interaction. In this particular 𝐺 𝐺 𝑚 𝑦 𝑚 2 𝑥 𝑧 𝑓 (r)=𝑁 𝑟 𝑟 𝑟 exp [−𝛼 ‖r‖ ]. (98) 𝑦 case,the planewavebasis constitutesaninappropriate choice m m 𝑥 𝑧 m as a very large number of basis elements are required in 𝑆 𝐺 where 𝑁 and 𝑁 are normalisation constants while m order to describe the quasi-discontinuities. Nevertheless, one m m represents the magnetic quantum numbers. Once an orbital is would be inclined to retain the low computational complexity 𝑖𝑗 𝑖𝑗 12 Advances in Chemistry associated with plane waves and the Fourier transform while where 𝑠 and 𝑑 are the wavelet coefficients. Naturally, 𝑘 𝑘 being able to efficiently describe the rougher parts of the higher resolutions may be achieved; the maximum resolution orbitals. This may be achieved with a wavelet basis and the is determined by the resolution of the computational grid. wavelet transforms [20, 48, 52–54]. As opposed to the plane As for the Fourier transform, the wavelet transform admits wave functions, which span the whole space, the wavelets are Fast Wavelet Transform implementation which has the same spatially localised. er Th efore, they are particularly suited for complexity as its Fourier counterpart. the description of discontinuities and fast-varying functions. In the next subsection, we briefly present a hybrid ap- In addition, the wavelet basis provides a multiresolution rep- proachwhichinvolvesbothabinitiomoleculardynamicsand resentation of the orbitals which means that the calculations molecular dynamics. may be performed efficiently at the required resolution level. The wavelet functions are filter banks [55]. In one dimension, 10.3. Hybrid Quantum Mechanics–Molecular Dynamics Meth- they involve two functions: ods. Because of its inherent complexity, the applicability of ab initio molecular dynamics is essentially limited to small (i) the scaling function: molecular domains. Molecular dynamics, which we recently reviewed in [30], is suitable for larger molecular system. 𝜙 (𝑥 )= 2 ∑ ℎ 𝜙 (2𝑥−𝑖 ) (103) A drawback is that it is not adapted for the simulation of 𝑖=1−𝑚 chemical complexes, since it is not based on first principles which is a high-pass filter responsible for the multiresolution as its quantum mechanics counterpart. Rather, molecular dy- aspect of the transform; namics relies on empirical potentials [56] and classical mechanics which allow for the simulation of much larger mo- (ii) the mother wavelet: lecular complexes. Therefore, in order to simulate larger systems, hybrid 𝜅 (𝑥 )= 2 ∑ 𝑔 𝜅 (2𝑥 − 𝑖 ) (104) approaches [12, 57–59] must be followed. In such a method, a 𝑖=1−𝑚 small region of interest is simulated with ab initio molecular whichisaband-passfilterresponsibleforthebasis localisa- dynamics, while the rest of the molecular system is approxi- tion. eTh coefficients of the scaling and the mother wavelet mated with molecular mechanics: function are related by QM/MM QM MM QM-MM L = L + L + L . (109) CP 𝑔 =−1 ℎ . (105) 𝑖 −𝑖+1 An extra Lagrangian term is required in order to ensure a In three dimensions, at the lowest resolution, the scaling proper coupling between the quantum degrees of freedom function is given by and the classical degrees of freedom. This Lagrangian consists 𝑥 𝑦 𝑧 of a bounded and an unbounded part. eTh bounded part con- 𝜙 (r)=𝜙( −𝑖)𝜙( −𝑗)𝜙( −𝑘) (106) sists of the stretching, bending, and torsional terms which are 𝑎 𝑎 𝑎 characterised by their distances, angles, and dihedral angles, while the mother wavelets become respectively. eTh unbounded part contains the electrostatic 𝑥 𝑧 𝜅 (r)=𝜅( −𝑖)𝜙( −𝑗)𝜙( −𝑘), 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. 𝜅 (r)=𝜙( −𝑖)𝜅( −𝑗)𝜙( −𝑘), 𝑎 𝑎 𝑎 11. Conclusions 𝑥 𝑦 𝑧 𝜅 (r)=𝜅( −𝑖)𝜅( −𝑗)𝜙( −𝑘), 𝑎 𝑎 𝑎 In this paper, we have presented a comprehensive, but yet 𝑥 𝑧 4 concise, review of ab initio molecular dynamics from a com- 𝜅 (r)=𝜙( −𝑖)𝜙( −𝑗)𝜅( −𝑘), (107) 𝑎 𝑎 𝑎 putational perspective and from rfi st principles. Although it is always hazardous to speculate about the future, we 𝑥 𝑦 𝑧 𝜅 (r)=𝜅( −𝑖)𝜙( −𝑗)𝜅( −𝑘), foresee two important breakthroughs. Fourier transforms, 𝑎 𝑎 𝑎 which constitute an essential part of many molecular simula- 𝑥 𝑦 𝑧 tions, may be evaluated with high performance on graphical 𝜅 r =𝜙( −𝑖)𝜅( −𝑗)𝜅( −𝑘), ( ) processing units or GPU [54]. eTh same remark applies to 𝑎 𝑎 𝑎 the wavelet transform whose role is essential in describing 𝑥 𝑧 𝜅 (r)=𝜅( −𝑖)𝜅( −𝑗)𝜅( −𝑘), discontinuous orbitals [52]. This opens the door to high per- 𝑎 𝑎 𝑎 formance simulations, while tempering the limitations asso- where 𝑎 is the resolution of the underlying computational ciated with computational complexity. grid. As a result, any orbital function may be approximated For decades, one of the main bottlenecks of molecular by a truncated expansion: simulations has been the calculation of the density func- tionals. As in numerous fields, machine learning constitutes ] ] a promising avenue for the fast and efficient evaluation of 𝜓 (r)= ∑𝑠 𝜙 (r)+ ∑ ∑𝑑 𝜅 (r), (108) 𝑘 𝑘 𝑘 𝑘 these functionals without recourse to explicit calculations 𝑖,𝑗,𝑘 𝑖,𝑗,𝑘 ]=1 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 Advances in Chemistry 13 [17, 24, 25, 60]. Here, the density functional is simply learned [15] K. Kitaura and K. Morokuma, “A new energy decomposition scheme for molecular interactions within the Hartree-Fock from existing examples with machine learning techniques. approximation,” International Journal of Quantum Chemistry, This paves the way for the simulation for larger and more vol. 10, no. 2, pp. 325–340, 1976. complex molecular systems. We plan to review machine [16] H. S. Yu, X. He, and D. G. Truhlar, “MN15-L: A New Local Ex- learning-based approaches in a future publication. change-Correlation Functional for Kohn-Sham Density Func- tional eTh ory with Broad Accuracy for Atoms, Molecules, and Conflicts of Interest Solids,”JournalofChemicalTheoryandComputation ,vol.12,no. 3,pp.1280–1293,2016. eTh authors declare that there are no conflicts of interest re- [17] F.Brockherde,L.Vogt,L.Li,M.E.Tuckerman,K.Burke,andK.- garding the publication of this paper. R. Muller ¨ , “By-passing the Kohn-Sham equations with machine learning,” Nature Communications,vol.8,no. 1, articleno. 872, References [18] A. P. Selvadurai, Partial Dier ff ential Equations in Mechanics 2 , [1] P.Carloni,U.Rothlisberger, andM.Parrinello,“er Th oleand Springer Berlin Heidelberg, Berlin, Heidelberg, 2000. perspective of ab initio molecular dynamics in the study of [19] E. Engel and R. M. Dreizle, DensityFunctionalTheory:An biological systems,” Accounts of Chemical Research,vol.35, no. Advanced Course, Springer, London, UK, 2011. 6, pp. 455–464, 2002. [20] T. D. Engeness and T. A. Arias, “Multiresolution analysis for [2] L. Monticelli and E. Salonen, Eds., Biomolecular Simulations– efficient, high precision all-electron density-functional calcu- Methods and Protocols, Humana Press, Springer, London, UK, lations,” Physical Review B: Condensed Matter and Materials Physics,vol.65, no.16, pp.1–10,2002. [3] T.D.Kuhne ¨ , Ab-Initio Molecular Dynamics, 2013, https://arxiv [21] A. J. Cohen, P. Mori-Sanc ´ hez, and W. Yang, “Challenges for den- .org/abs/1201.5945. sity functional theory,” Chemical Reviews,vol.112,no.1,pp.289– [4] M. E. Tuckerman, “Ab initio molecular dynamics: Basic con- 320, 2012. cepts, current trends and novel applications,” Journal of Physics: [22] H. Ou-Yang and M. Levy, “Theorem for functional derivatives Condensed Matter, vol. 14, no. 50, pp. R1297–R1355, 2002. in density-functional theory,” Physical Review A: Atomic, Molec- [5] J.Broeckhoveand L.Lathouwers,Eds., Time-Dependent Quan- ular and Optical Physics,vol.44, no.1,pp. 54–58, 1991. tum Molecular Dynamics, NATO Advanced Research Workshop [23] P. Mori-Sanc ´ hez, A. J. Cohen, and W. Yang, “Self-interaction- on Time Dependent Quantum Molecular Dynamics: Theory and free exchange-correlation functional for thermochemistry and Experiment,Springer, Snowbird,USA,1992. kinetics,” The Journal of Chemical Physics ,vol.124,no.9,Article [6] F. Jensen, Introduction to Computational Chemistry,Wiley, ID 091102, 2006. Chichester, Uk, 2017. [24] L. Li,T.E.Baker,S.R.White,andK. Burke,“Puredensityfunc- [7] M.E.Tuckerman, Tuckerman, Statistical Mechanics: eor Th y and tional for strong correlation and the thermodynamic limit from Molecular Simulation, Oxford University Press, Oxford, UK, machine learning,” Physical Review B: Condensed Matter and Materials Physics,vol.94, no.24, ArticleID245129, 2016. [8] R. P. Steele, “Multiple-Timestep ab Initio Molecular Dynamics [25] F. Pereira,K. Xiao,D.A.R.S.Latino, C. Wu,Q.Zhang,andJ. Using an Atomic Basis Set Partitioning,” The Journal of Physical Aires-De-Sousa, “Machine Learning Methods to Predict Den- Chemistry A, vol. 119, no. 50, pp. 12119–12130, 2015. sity Functional eTh ory B3LYP Energies of HOMO and LUMO [9] N.Luehr,T.E.Markland,andT.J.Mart´ınez, “Multiple time Orbitals,” Journal of Chemical Information and Modeling,vol.57, step integrators in ab initio molecular dynamics,” The Journal no. 1, pp. 11–21, 2017. of Chemical Physics,vol.140,no.8, ArticleID084116,2014. [26] M. E. Tuckerman and M. Parrinello, “Integrating the Car-Par- [10] H.B. Schlegel,Li. X.,J.M.Millam,G.A.Voth,G.E. Scuseria, rinello equations. I. Basic integration techniques,” The Journal and M. J. Frisch, “Ab initio molecular dynamics: Propagating of Chemical Physics,vol.101,no.2, pp.1302–1315,1994. the density matrix with Gaussian orbitals. III. Comparison with [27] J. Lach, J. Goclon, and P. Rodziewicz, “Structural flexibility of BornOppenheimer dynamics,” Journal of Chemical Physics,vol. the sulfur mustard molecule at finite temperature from Car- 117,no.19,pp.9758–9763,2001. Parrinello molecular dynamics simulations,” Journal of Haz- [11] T. D. Kuhne ¨ , M. Krack, F. R. Mohamed, and M. Parrinello, “Effi- ardous Materials,vol.306,pp.269–277, 2016. cient and accurate car-parrinello-like approach to born-oppen- [28] G. Bussi, T. Zykova-Timan, and M. Parrinello, “Isothermal- heimer molecular dynamics,” Physical Review Letters, vol. 98, isobaric molecular dynamics using stochastic velocity rescal- no.6,ArticleID 066401,2007. ing,” The Journal of Chemical Physics ,vol.130,no.7, ArticleID [12] X. Andrade, A. Castro, D. Zueco et al., “Modified ehrenfest for- 074101, 2009. malism for efficient large-scale ab initio molecular dynamics,” [29] J. Lahnsteiner, G. Kresse, A. Kumar, D. D. Sarma, C. Franchini, Journal of Chemical eor Th y and Computation ,vol.5,no.4,pp. and M. Bokdam, “Room-temperature dynamic correlation be- 728–742, 2009. tween methylammonium molecules in lead-iodine based per- [13] M. Vacher, D. Mendive-Tapia, M. J. Bearpark, and M. A. Robb, ovskites: An ab initio molecular dynamics perspective,” Physical “eTh second-order Ehrenfest method: A practical CASSCF Review B: Condensed Matter and Materials Physics,vol.94, no. approach to coupled electron-nuclear dynamics,” Theoretical 21, Article ID 214114, 2016. Chemistry Accounts,vol.133,no. 7,pp.1–12, 2014. [30] E. Paquet and H. L. Viktor, “Molecular dynamics, monte carlo [14] F.Ding,J.J.Goings,H. Liu,D.B. Lingerfelt,andX. Li,“Ab initio simulations, and langevin dynamics: A computational review,” two-component Ehrenfest dynamics,” The Journal of Chemical BioMed Research International,vol.2015, ArticleID183918, Physics, vol. 143, no. 11, Article ID 114105, 2015. 2015. 14 Advances in Chemistry [31] M. E. Tuckerman, Y. Liu, G. Ciccotti, and G. J. Martyna, “Non- [48] L. Genovese, A. Neelov, S. Goedecker et al., “Daubechies wave- Hamiltonian molecular dynamics: Generalizing Hamiltonian lets as a basis set for density functional pseudopotential calcu- phase space principles to non-Hamiltonian systems,” The Jour- lations,” The Journal of Chemical Physics ,vol.129,no. 1, Article nal of Chemical Physics, vol. 115, no. 4, pp. 1678–1702, 2001. ID 014109, 2008. [49] M. A. Watson, N. C. Handy, and A. J. Cohen, “Density func- [32] M. Ceriotti, G. Bussi, and M. Parrinello, “Colored-noise ther- mostats ala ` Carte,” Journal of Chemical eTh ory and Computa- tional calculations, using Slater basis sets, with exact exchange,” The Journal of Chemical Physics ,vol.119,no. 13,pp. 6475–6481, tion,vol.6,no.4, pp.1170–1180,2010. [33] M. Ceriotti, G. Bussi, and M. Parrinello, “Langevin equation [50] X. Andrade and A. Aspuru-Guzik, “Real-space density func- with colored noise for constant-temperature molecular dynam- tional theory on graphical processing units: Computational ap- ics simulations,” Physical Review Letters,vol.102,no. 2,Article proach and comparison to Gaussian basis set methods,” Journal ID 020601, 2009. of Chemical Theory and Computation ,vol.9,no. 10,pp. 4360– [34] M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Mano- 4373, 2013. lopoulos, “Efficient stochastic thermostatting of path integral [51] S. V. Levchenko, X. Ren, J. Wieferink et al., “Hybrid functionals molecular dynamics,” eTh JournalofChemicalPhysics ,vol.133, for large periodic systems in an all-electron, numeric atom- no. 12, Article ID 124104, 2010. centered basis framework,” Computer Physics Communications, [35] R. Kubo, “eTh fluctuation-dissipation theorem,” Reports on Pro- vol. 192, pp. 60–69, 2015. gress in Physics,vol.29, no.1,articleno.306,pp. 255–284, 1966. [52] T. Deutsch and L. Genovese, “Wavelets for electronic structure [36] G. J. Martyna, A. Hughes, and M. E. Tuckerman, “Molecular dy- calculations,” Collection SFN,vol.12,pp.33–76,2011. namics algorithms for path integrals at constant pressure,” The [53] A. H. Prociuk and S. S. Iyengar, “A multiwavelet treatment of the Journal of Chemical Physics, vol. 110, no. 7, pp. 3275–3290, 1999. quantum subsystem in quantum wavepacket ab initio molecular [37] T. Spura, C. John, S. Habershon, and T. D. Kuhne ¨ , “Nuclear dynamics through an hierarchical partitioning of momentum quantum eeff cts in liquid water from path-integral simulations space,” JournalofChemicalTheoryandComputation ,vol.10,no. using an ab initio force-matching approach,” Molecular Physics, 8, pp. 2950–2963, 2014. vol.113,no.8, pp.808–822,2015. [54] L. Genovese, B. Videau, D. Caliste, J.-F. Meh ´ aut, S. Goedecker, [38] O. Marsalek and T. E. Markland, “Ab initio molecular dynamics andT.Deutsch,“Wavelet-BasedDensity FunctionalTheoryon with nuclear quantum eeff cts at classical cost: Ring polymer MassivelyParallelHybridArchitectures,” Electronic Structure contraction for density functional theory,” The Journal of Chem- Calculations on Graphics Processing Units: From Quantum ical Physics,vol.144,no. 5,ArticleID054112, 2016. Chemistry to Condensed Matter Physics, pp. 115–134, 2016. [39] M. E. Tuckerman, D. Marx, M. L. Klein et al., “Efficient and gen- [55] R. C. Gonzalez and R. E. Woods, Digital Image Processing,Pear- eral algorithms for path integral Car-Parrinello molecular dy- son Prentice Hall, Upper Saddle River, NJ, USA, 2008. namics,” The Journal of Chemical Physics ,vol.104,no.14,pp. [56] K. Vanommeslaeghe, E. Hatcher, and C. Acharya, “CHARMM 5579–5588, 1996. general force field: a force field for drug-like molecules compati- [40] P.Minary,G.J.Martyna,and M. E.Tuckerman, “Algorithms ble with the CHARMM all-atom additive biological force fields,” and novel applications based on the isokinetic ensemble. I. Bio- Journal of Computational Chemistry,vol.31,no.4,pp.671–690, physical andpathintegralmolecular dynamics,” The Journal of Chemical Physics,vol.118,no. 6,pp.2510–2526,2003. [57] J. R. Perilla, B. C. Goh, C. K. Cassidy et al., “Molecular dy- [41] M. Ceriotti, J. More, and D. E. Manolopoulos, “I-PI: A Python namics simulations of large macromolecular complexes,” Cur- interface for ab initio path integral molecular dynamics simu- rent Opinion in Structural Biology,vol.31, pp.64–74,2015. lations,” Computer Physics Communications,vol.185,no.3,pp. [58] N. Sahu and S. R. Gadre, “Molecular tailoring approach: A route 1019–1026, 2014. for ab initio treatment of large clusters,” Chemical Reviews,vol. [42] H. Y. Geng, “Accelerating ab initio path integral molecular dy- 47, no. 9, pp. 2739–2747, 2014. namics with multilevel sampling of potential surface,” Journal of [59] J. Dreyer,B.Giuseppe, E. Ippoliti,V.Genna,M.DeVivo, and Computational Physics, vol. 283, pp. 299–311, 2015. P. Carloni, “First principles methods in biology: continuum [43] K. A. Fichthorn and W. H. Weinberg, “eo Th retical foundations models to hybrid ab initio quantum mechanics/molecular of dynamical Monte Carlo simulations,” The Journal of Chemical mechanics,” in Simulating Enzyme Reactivity: Computational Physics,vol.95,no.2,pp. 1090–1096,1991. Methods in Enzyme Catalysis,I.Tun˜o´nandV.Moliner,Eds., chapter 9, Royal Society of Chemistry, London, UK, 2016. [44] D. Marx and M. Parrinello, “Ab initio path integral molecular dynamics: Basic ideas,” The Journal of Chemical Physics ,vol.104, [60] V. Botu and R. Ramprasad, “Adaptive machine learning frame- no. 11, pp. 4077–4082, 1995. work to accelerate ab initio molecular dynamics,” International JournalofQuantumChemistry,vol.115,no.16,pp.1074–1083, [45] O. Marsalek,P.-Y. Chen,R.Dupuisetal.,“Efficientcalculation of free energy differences associated with isotopic substitution using path-integral molecular dynamics,” Journal of Chemical Theory and Computation ,vol.10,no.4,pp. 1440–1453,2014. [46] S. Goedecker and M. Teter, “Separable dual-space Gaussian pseudopotentials,” Physical Review B: Condensed Matter and Materials Physics,vol.54,no.3,pp. 1703–1710,1996. [47] C. Hartwigsen, S. Goedecker, and J. Hutter, “Relativistic sep- arable dual-space Gaussian pseudopotentials from H to Rn,” Physical Review B: Condensed Matter and Materials Physics,vol. 58,no. 7,pp.3641–3662,1998. Journal of Nanomaterials Journal of International Journal of International Journal of The Scientific Analytical Methods Journal of Photoenergy in Chemistry World Journal Applied Chemistry Hindawi Hindawi Hindawi Publishing Corporation Hindawi Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 http://www www.hindawi.com .hindawi.com V Volume 2018 olume 2013 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Advances in International Journal of Physical Chemistry Medicinal Chemistry Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Submit your manuscripts at www.hindawi.com Journal of Bioinorganic Chemistry and Applications Materials Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Advances in Journal of International Journal of International Journal of BioMed Tribology Chemistry Research International Spectroscopy Electrochemistry Hindawi Hindawi Hindawi Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Enzyme Journal of Journal of International Journal of Biochemistry Analytical Chemistry Spectroscopy Nanotechnology Research Research International Hindawi Hindawi Hindawi Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Nanomaterials

Journal

Advances in ChemistryHindawi Publishing Corporation

Published: Apr 29, 2018

References