Computational Methods for Ab Initio Molecular Dynamics

Paquet, Eric; Viktor, Herna L.

doi:https://doi.org/10.1155/2018/9839641

Advances in Chemistry

On this page

Abstract Introduction Conclusions Conflicts of Interest References Copyright Related Articles

Review Article | Open Access

Volume 2018 | Article ID 9839641 | https://doi.org/10.1155/2018/9839641

Computational Methods for Ab Initio Molecular Dynamics

Eric Paquet¹and Herna L. Viktor²

Academic Editor: Yusuf Atalay

Received23 Feb 2018

Accepted20 Mar 2018

Published29 Apr 2018

Abstract

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

1. Introduction

Ab initio molecular dynamics (AIMD) is an irreplaceable technique for the realistic simulation of complex molecular systems and processes associated with biological organisms [1, 2] such as monoclonal antibodies as illustrated in Figure 1. With ab initio molecular dynamics, it is possible to predict in silico phenomena for which an in vivo experiment is either too difficult or expensive, or even currently deemed infeasible [3, 4]. Ab initio molecular dynamics essentially differs from molecular dynamics (MD) in two ways. Firstly, AIMD is based on the quantum Schrödinger equation while its classical counterpart relies on Newton’s equation. Secondly, MD relies on semiempirical effective potentials which approximate quantum effects, while AIMD is based on the real physical potentials.

In this paper, we present a review of ab initio molecular dynamics from a computational perspective and from first principles. Our main aim is to create a document which is self-contained, coherent, and concise but still as complete as possible. As the theoretical details are essential for real-life implementation, we have provided the equations for the relevant physical aspect and approximations presented. Most important, we expose how these equations and concepts are related to one another.

The paper is organised as follows. In Section 2, quantum mechanics is reviewed from a molecular dynamics perspective. Two important approximations are introduced, namely, the adiabatic and the Born–Oppenheimer approximations. Subsequently, in Section 3, the Ehrenfest molecular dynamics is presented, which allows for a semiclassical treatment of the nuclei. This opens the door, in Section 4, to the Born–Oppenheimer molecular dynamics formulation. This is followed, in Section 5, by the Hartree–Fock formulation which approximates the atomic antisymmetric wave function by a single determinant of the orbitals. Sections 6 and 7 introduce the Kohn–Sham formulation whose primary objective is to replace the interacting electrons by a fictitious, but equivalent, system of noninteracting particles. Electrons are reintroduced as dynamical quantities, in Section 8, through the Car–Parrinello formulation. Isothermal and isobaric processes are addressed in this section. In Section 9, a path integral formulation of molecular dynamics is presented, which makes allowance for a quantum treatment of the nuclear degrees of freedom. Finally, in Section 10, some important computational issues are addressed such as a simplification of the electronic interaction with the pseudopotential, the representation of orbitals in terms of a functional basis, the use of the Fourier and wavelet transform in order to reduce the computational complexity, and the simulation of larger systems with hybrid molecular dynamics. This is followed by a conclusion in Section 11.

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

We review some important notions of quantum mechanics from a molecular perspective. In quantum mechanics [5–7], the Hamiltonian [5–7] is the operator corresponding to the total energy of the molecular system associated with the electrons and the nuclei. The total Hamiltonian is the sum of the kinetic energies of all the particles plus the potential energy of the constituent particles; in occurrence the electrons and the nucleiwhere is the Planck constant and is the mass of a given nucleus. The electronic and nuclear Cartesian coordinates for a given particle are defined, respectively, asThe ensemble of all electronic coordinates is represented by while the ensemble of all nuclear coordinates is denoted by . In addition, the gradient and the Laplacian operators for a given electron and nucleus are defined, respectively, asThe gradient operator, a vector, is related to the momentum while the Laplacian operator, a scalar quantity, is associated with the kinetic energy [5–7]. The electronic Hamiltonian is defined as the sum of the kinetic energy of the electrons, the potential energy associated with each pair of electrons, the potential energy associated with each pair of electron-nucleus, and the potential energy concomitant to each pair of nuclei:where is the electronic mass, is the electronic charge, is the atomic number (the number of protons in a given nucleus), and is the vacuum permittivity. The quantum molecular system is characterised by a wave function whose evolution is determined by the time-dependent Schrödinger equation [5–7]:where . The Schrödinger equation admits a stationary or time-independent solution:where is the energy associated with the electronic wave function and is a set of quantum numbers that labelled the Eigenstates or wave functions as well as the Eigenvalues or energies associated with the stationary Schrödinger equation. The total wave function may be expended in terms of time-dependent nuclear wave functions and stationary electronic wave functions :It should be noticed that this expansion is exact and does not involve any approximation. The electronic wave functions are solution of the stationary Schrödinger equation, an Eigen equation which involves the electronic Hamiltonian. If one substitutes this expansion in the time-dependent Schrödinger equation, one obtains a system of equations for the evolution of the nuclear wave functions:As may be seen, the nuclear wave functions are coupled to the electronic wave functions. The strength of the coupling is determined by the coupling coefficients which form a matrix:This system of equations is too complex to be solved directly. Consequently, in the next subsection, we consider two important approximations of the Schrödinger equation, namely, the adiabatic and the Born–Oppenheimer approximations, which aim to reduce such a complexity. These approximations, as well as those that later follow, reduce substantially the duration of the calculations allowing for larger molecular systems to be simulated and longer time-scales to be explored [8, 9].

2.1. Adiabatic and Born–Oppenheimer Approximations

The above-mentioned system of differential equations is too complex to be solved directly. Some approximations must be performed in order to reduce the computational complexity, while maintaining the predictive accuracy of the calculations. The first approximation, which is called the adiabatic approximation [5–7], assumes that the electronic wave functions adapt quasi-instantaneously to a variation of the nuclear configuration. This approximation is justified by the fact that the electrons are much lighter than the nucleus. Consequently, such an approximation is valid, unless the motion of the electron becomes relativistic, which may happen if the electrons are too close to the nucleus. Mathematically, this approximation implies thatAs a result, the coupling matrix becomes diagonal:where is the well-known Kronecker delta.

The second approximation, which is called the Born–Oppenheimer approximation [6, 10, 11], assumes that the electronic and the nuclear motions are separable as a result of the difference between nuclear and electronic masses. As a result, the expansion reduces to a single Eigen state:In addition, as the energy associated with the wave function is usually much larger than the corresponding coupling constant:the time-dependent nuclear Schrödinger further reduces toFurthermore, it is assumed that the electronic wave function is in its ground or nonexcited state:This occurs if the thermal energy is lower than the energy difference between the ground state and the first excited state, that is, if the temperature is low.

In the next section, we introduce another approximation, in which the motion of heavy nuclei is described by a semiclassical equation.

3. Ehrenfest Molecular Dynamics

Another possible approximation is to assume that the nuclear motion is semiclassical as it is the case when the nuclei are relatively heavy. This implies that, instead of being determined by the Schrödinger equation, the average nuclear motion is determined by Newton’s equation. The right form for this equation is given by the Ehrenfest theorem [6, 12–14] which states that the potential, which governs the classical motion of the nucleons, is equal to the expectation or average, in the quantum sense, of the electronic Hamiltonian with respect to the electronic wave function:where the bra–ket notation is to be understood asfor any operator and where . This equation may be further simplified if the adiabatic and the Born–Oppenheimer approximations are enforced:In the specific context of Ehrenfest molecular dynamics, the electrons follow a time-dependent Schrödinger equation:The fact that this equation is time-dependent ensures that the orthogonality of the orbitals is maintained at all times. The motivation is to be found in the fact that the electronic Hamiltonian is a Hermitian operator [6]. As stated earlier, this is not the case in the Born–Oppenheimer approximation in which the electronic wave function is governed by the stationary Schrödinger equation which does not maintain the orthonormality over time.

In the next section, we apply the adiabatic and Born–Oppenheimer approximation in the context of ab initio molecular dynamics.

4. Born–Oppenheimer Molecular Dynamics

In Born–Oppenheimer molecular dynamics [6, 10, 11], it is assumed that the adiabatic and the Born–Oppenheimer approximations are valid and that the nuclei follow a semiclassical Newton equation whose potential is determined by the Ehrenfest theorem. It is further assumed that the electronic wave function is in its ground state (lowest energy). Since the stationary Schrödinger equation is used for the electronic degrees of freedom, the orthogonality condition must be enforced with a Lagrange multiplier as this condition is not preserved by the stationary Schrödinger equation over time. The orthogonality of the orbitals is a physical requirement which must be enforced at all time in order to be able to compute real observable quantities. The Born–Oppenheimer molecular dynamics is characterised by a Lagrangian [6] which is defined as the difference between the kinetic and the potential energy:The first term of the Lagrangian corresponds to the kinetic energy of the nuclei, the second term corresponds to the nuclear potential as obtained from the Ehrenfest theorem, and the last term is a Lagrange function which ensures that the orbitals remain orthonormal at all time. The Lagrange multipliers are denoted by . Since the electrons follow a Fermi–Dirac statistics [6], they obey the Pauli Exclusion Principle (two electrons cannot be in the same quantum state) which means that the electronic wave function must be an antisymmetric function of its orbitals, namely, the wave functions associated with the individual electrons.

The equations of motion associated with this Lagrangian are obtained from the Euler–Lagrange equations [6]. There is one system of Euler–Lagrange equations for the nuclear degrees of freedom:and one system of Euler–Lagrange equations for their electronic counterpart:From the Euler–Lagrange equations, one obtainsfor the nuclear equations of motion, andfor the electronic equations of motion. One should notice the presence of the Lagrange multiplier in the electronic equations of motion which ensures that the orbital remains orthonormal at all time.

In the next section, we introduce another approximation in order to further reduce the complexity of the electronic Hamiltonian.

5. Hartree–Fock Molecular Dynamics

From now on, we shall adopt the atomic unit system:in order to alleviate the notation. In Hartree–Fock molecular dynamics [15], the antisymmetric electronic wave function is approximated by a single determinant of the electronic orbitals :Two operators are associated with the electronic interaction. The first one, the Coulomb operator, corresponds to the electrostatic interaction between the orbitals:The second operator, the exchange operator [6, 15], is associated with the exchange energy, a quantum mechanical effect that occurs as a result of the Pauli Exclusion Principle:From the Coulomb and the exchange operators, an electronic Hamiltonian may be inferred:where the external potential is defined asThis Hamiltonian determines, in turn, the electronic Lagrangian:A Lagrange multiplier term is added to the Lagrangian in order to ensure that the orbitals remain orthonormal at all time: a quantum mechanical requirement [6]. From the Euler–Lagrange equations, one obtains the equations of motion for the electrons:which differ from (24) by the Hamiltonian.

In the next subsection, we seek to replace the interacting electrons by a fictitious but equivalent system of noninteracting particles in order to further reduce the computational complexity.

6. Kohn–Sham Molecular Dynamics

The Kohn–Sham formulation [16, 17] aims to replace the interacting electrons by a fictitious system of noninteracting particles that generate the same electronic density as the physical system of interacting particles. The electronic density is defined aswhere is the occupation number of a given orbital.

In this formulation of AIMD, the Ehrenfest term is replaced by the Kohn–Sham energy:The latter is define aswhere the first term is the total electronic kinetic energy associated with the electrons, the second term is the electrostatics interaction energy between the electronic density and the external potential, and the third term is the self-electrostatic interaction energy associated with the electronic density. The latter involves the interaction of the electronic density with it self-created electrostatic potential. This potential, called the Hartree potential, is obtained by solving the Poisson equation [18]:The last term is the celebrated exchange-correlation energy or density functional [16, 19–21] which takes into account the residual electronic interaction, that is, the self-interaction of the electrons. Unfortunately, the density functional has no closed-form solution but many approximations are known. Some of these approximations are considered in the next section.

From the exchange-correlation energy, one may define the exchange-correlation potential:which is the functional derivative [22] of the exchange-correlation energy with respect to the electronic density. In addition, one may define the Kohn–Sham potential:As for the other approaches, a Lagrangian is defined in which the orthonormality conditions are enforced with Lagrange multipliers. The electronic equations of motion, which are determined by the Euler–Lagrange equations arewhere the fictitious one-particle Hamiltonian is given byA canonical form may be obtained if a unitary transformation is applied to the previous equation:In the next section, we introduce some density functionals in order to replace the interacting electrons by an equivalent but yet simpler system of noninteracting particles.

7. Exchange-Correlation Energy

The detailed analysis of density functionals, also known as exchange-correlation energies [23], is beyond the scope of this review. We refer the interested reader to [21] for an exhaustive analysis. In this section, we briefly introduce a few common density functionals. As mentioned earlier, the aim of the density functional is to express the electronic interaction in terms of the sole electronic density.

In the simplest case, the exchange-correlation energy is a functional of the electronic density alone for which the most important representative is the local density approximation:One may take into consideration, in addition to the electronic density, the gradient of the electronic density, in which case the approach is called the generalised gradient approximation:where the dimensionless reduced gradient is defined asAmong these approximations is the B88 approximation:where refers to the spin and the Perdew, Burke, and Ernzerhof (PBE) approximation:The exchange energy, which is associated with the Pauli Exclusion Principle, may be characterised by the Hartree–Fock energy: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:where , , and are empirical parameters. Obviously, these parameters are application dependent.

Recently, it has been proposed to evaluate the functional density with machine learning techniques. The functional density is learned by examples instead of directly solving the Kohn–Sham equations [17, 24, 25]. As a result, substantially less time is required to complete the calculations allowing for larger system to be simulated and longer time-scales to be explored.

In the next section, we further improve the precision of the calculations by reintroducing the dynamic electronic degrees of freedom which, until now, have been absent from the equations of motion.

8. Car–Parrinello Molecular Dynamics

8.1. Equations of Motion

The Kohn–Sham dynamics, as formulated in the previous sections, does not take the dynamics of the electrons into account despite the fact that it is present: an unattractive unphysical feature. Indeed, one obtains from the Lagrangian and the Euler–Lagrange equationswhich means the orbitals are not dynamical fields of the molecular system. Nevertheless, the Lagrangian may be modified in order to also include the dynamic nature of electronic degrees of freedom through the introduction of a fictitious electronic kinetic energy term. The extended Lagrangian, which is called the Car–Parrinello Lagrangian [11, 26, 27]differs from the original Kohn–Sham Lagrangian by the second term which associates a fictitious kinetic energy to the electronic orbitals. The parameter acts as a fictitious electronic mass or inertia. As in the previous sections, the nuclear equations of motionand the electronic equations of motionmay be inferred from the Euler–Lagrange equations. These equations involve a second order time derivative of the orbitals which means that the latter are now proper dynamical quantities.

In the next two subsections, we consider ab initio molecular dynamics at constant temperature and at constant pressure.

8.2. Massive Thermostating and Isothermal Processes

In the previous sections, we have implicitly assumed that the molecular system under consideration was isolated. Of course, this is not compatible with most experimental conditions [27–29] in which the system is either kept at constant temperature (isothermal) in a heat bath or at constant pressure (isobaric) as a consequence, for instance, of the atmospheric pressure. In this subsection and in the next, we explain how to address these important issues.

As we have recently explained in a computational review on molecular dynamics [30], an isothermal process cannot be obtained by adding an extra term to the Lagrangian. Rather, the equations of motions must be modified directly as a result of the non-Hamiltonian nature of the isothermal process [31]. In order to obtain an isothermal molecular system, a friction term must be recursively added to each degree of freedom [28, 29]. Therefore, the electronic equations of motion must be modified as follows:while the nuclear equations of motion must take a similar form:where , is the temperature of the heat bath, is the Boltzmann constant, the are frictional electronic degrees of freedom, the are the friction coefficients associated with the electronic orbitals, the are frictional nuclear degrees of freedom, the are the friction coefficients associated with the nuclei, and denotes the number of dynamical degrees of freedom to which the nuclear thermostat chain is coupled. The friction terms are proportional to the velocity of the corresponding degrees of freedom as it is customary. Not only must friction terms be added to the physical degrees of freedom, but each nonphysical friction term, in turn, must be thermalised by another nonphysical friction term until an isothermal state is properly achieved. The terms and may be considered as dynamical friction coefficients for the physical degrees of freedom.

In the next section, we present an alternative approach based on the generalised Langevin equation.

8.3. Generalised Langevin Equation and Isothermal Processes

Massive thermostating is not the sole approach to simulate an isothermal process. Indeed, the latter may be achieved by means of the generalised Langevin equation [32–34]. In this subsection, we restrict ourselves to one generalised coordinate in order not to clutter the notation. The generalisation to more than one coordinate is immediate. The generalised Langevin equation, which is a differential stochastic equation, may be written aswhere is the drift matrix: is the diffusion matrix, is a generalised coordinate associated with an atom or a nucleus (position), is the corresponding generalised momentum, and is a set of hidden nonphysical momenta. The structure of the diffusion matrix is similar to the one of the corresponding drift matrix. The random process is characterised by an uncorrelated Gaussian noise:The hidden momenta, associated with the generalised Langevin equation, may be marginalised because the evolution of the variables , in the free particle limit, is described by a linear Markovian stochastic differential equation. As a result, the generalised Langevin equation becomes equivalent to a Langevin equation with memory kernel and noise correlation [32–34]:The memory kernel, which is a dissipative term associated with friction, is given by the expressionwhereas the noise correlation, which characterised the fluctuations associated with the random noise, is provided bywhere the matrixes and are defined asrespectively. The canonical, isothermal ensemble is obtained by applying the fluctuation-dissipation theorem [35]. The fluctuation-dissipation theorem states that the Fourier transforms of the noise correlation and the memory kernel must be related byThe fluctuation-dissipation theorem implies in turn thatTherefore, the fluctuation-dissipation theorem fixes the diffusion matrix once the drift matrix is determined. As a result, an isothermal process follows immediately.

In the next subsection, we address isobaric molecular processes, which are very common as most experiments are performed at atmospheric pressure.

8.4. Isobaric Processes

As opposed to the isothermal process, the isobaric process is a Hamiltonian process which means that it may be obtained by adding extra terms to the Lagrangian [28, 36].

In order to model an isobaric molecular process, the volume occupied by the molecular system is divided into congruent parallelepipeds called Bravais cells. These cells are characterised by three oriented vectors which correspond to the three primitive edges spanning their volume. These vectors are concatenated in order to form a Bravais matrix:The Bravais cell is characterised by its volume and by a local metric:In order to introduce the volume as a dynamic quantity into the Lagrangian, the nuclear and the electronic coordinates are reformulated in terms of the Bravais matrix and of the scaled coordinates:The Bravais matrix and the scale coordinates constitute distinct degrees of freedom. In order to obtain an isobaric system, the Car–Parrinello Lagrangian must be supplemented with three additional terms which appear on the third line of the Lagrangian:The first term represents the kinetic energy associated with the scale coordinates. One should notice the presence of the metric or quadratic form in the inner product. The second term represents the kinetic energy associated with the Bravais matrix. The matrix norm is the square of the Frobenius norm while is a fictitious mass or inertia attributed to the Bravais cells. The last term is associated with the pressure of the barostat (ambient pressure). The equations of motion are obtained, as usual, from the Euler–Lagrange equations. In particular, the Euler–Lagrange equations for the Bravais cells take the form:In the next section, we present a path integral formulation of ab initio molecular dynamics which allows a quantum formulation of both the electronic and nuclear degrees of freedom.

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

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

9.1. Path Integral Formulation

The quantum path integral formulation [36, 39–42] is based on the partition function of the quantum statistical canonical ensemble which is defined as the trace of the exponential of the Hamiltonian operator:The partition function describes the statistical properties of the molecular system. Since the operators associated with the electrons and the nuclei do not commute, the exponential of the Hamiltonian operator must be evaluated with the Trotter factorization:If the completeness relationwhich is an identity operator, is recursively inserted into the Trotter formula and if the adiabatic approximation is performed:the partition function becomesThe quantitiesare the amplitude and the angular frequency associated with the quantum harmonic oscillators. The latter appear as a consequence of the path integral formulation. The computational time is a discrete evolution parameter associated with the evolution of the molecular system. As such, it represents a particular time-slice. The path integral associated with the partition function is a weighted sum over all possible nuclear trajectories. The nuclear configuration at a particular time-slice is provided by the ensemble of all the individual nuclear configurations at this specific instant. The weighting function corresponds to the exponential factor which consists of two terms: the first one is related to the harmonic potential energy associated with the nuclei while the second is the energy associated with the electrons.

As in the previous sections, the Born–Oppenheimer approximation is legitimate if the thermal energy is much smaller than the difference between the electronic ground state and the excited states:It follows that the partition function becomesFrom this partition function, an extended Lagrangian may be defined:where are fictitious masses or unphysical auxiliary parameters associated with the nuclei whereas their physical masses are identified by . One should notice that fictitious momenta have been formally introduced. These momenta are required in order to evaluate numerically the path integral with Monte Carlo techniques [30, 43]. The reader should notice that the momenta affect neither the partition function (up to an irrelevant normalisation factor) nor the physical observables. The ground state energy must be evaluated concurrently with the nuclear energy for each time-slice.

The nuclear coordinates are not linearly independent. Nevertheless, they may become linearly independent if they are expressed in terms of their normal modes. The normal decomposition [39, 41, 44] is obtained by representing each coordinate in terms of complex Fourier series:The complex Fourier coefficients are further expanded in terms of their real and imaginary parts:As for the electronic structure, it may be obtained from one of the previously introduced approaches. For instance, for a Car–Parrinello technique formulated in terms of a Kohn–Sham Hamiltonian, the path integral Lagrangian in normal coordinates becomeswhere the normal mode masses are defined asto which the fictitious normal mode masses are closely related:where the centroid adiabaticity parameter belongs to the interval These masses are functions of the computational time. All the other quantities were defined earlier.

In the next subsection, we address the calculation of physical observables in the path integral framework.

9.2. Physical Observables

An observable [6, 7, 45] is a dynamic quantity that may be measured experimentally such as the average energy or the heat capacity. Numerous observables may be obtained directly from the partition function, for instance,(i)the expectation of the energy (average energy):(ii)the variance of the energy or energy fluctuation:(iii)the heat or thermal capacity:(iv)the entropy:(v)the Helmholtz free energy:

among others. In addition, if a small perturbation is applied to a molecular system, the expectation of the energy associated with this perturbation iswhere is a small running parameter.

In the next subsection, we present an isothermal formulation of path integrals.

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

Isothermal processes are essential for two reasons. The first one, which is rather evident, is that the simulation of realistic physical conditions often involves their simulation [34]. The second reason must be found in a rather subtle shortcoming of the formalism.

The massive thermostating of the path integral [34] is a “sine qua non” condition in order to obtain physically realistic results. Indeed, the harmonic potential leads to inefficient and nonergodic dynamics when microcanonical trajectories are used to generate ensemble averages.

In contrast, the thermostats generate ergodic, canonical averages at the expense of introducing sets of auxiliary chain variables which add to the complexity of the calculation, but which are nevertheless required for its precision and physical trustworthiness.

As we saw in Section 8.2, it is not possible to generate an isothermal process simply by adding extra terms to the Lagrangian. Rather, the equations of motion, which are obtained from the Euler–Lagrange equations, must be modified accordingly. Therefore, friction terms must be added recursively to the various degrees of freedom. As a result, the electronic equations of motion becomewhereas the nuclear equations of motion transform intoThe quantities appearing in these equations are all defined in Section 8.2. One should notice that number of equations is quite large due to the evolution parameter which is absent from the standard Car–Parrinello formulation as reported in Section 9.3.

In the next section, we address some important implementational issues.

10. Computational Implementation

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

10.1. Pseudopotential

For many calculations, the complete knowledge of the electronic interaction is not essential for the required precision. Consequently, for the sake of computational efficiency, the exact electronic potential may be approximated by means of an effective potential, called the pseudopotential [46–48]. Moreover, the relativistic effects associated with the core electrons may be implicitly incorporated into the pseudopotential, without recourse to explicit and intricate approaches.

A commonly employed pseudopotential, the dual-space Gaussian pseudopotential [46, 47], is composed of two parts, a local potential and a nonlocal potential:The local potential, which described the local interaction, is defined aswhere is the error function while , , , , and are adjustable empirical parameters. The nonlocal potentialthat describes the nonlocal interaction is defined in terms of the complex spherical harmonics functions and of the Gaussian projectors:In these equations, is the azimuthal quantum number and is the magnetic quantum number while is the well-known gamma function.

In the next subsection, we explore the advantages of projecting the orbitals on a functional basis.

10.2. Orbitals and Basis Functions

Any continuous function, such as an orbital, may be represented as a linear combination of basis functions:where are the projection coefficients and are the basis functions. Such decomposition may be either physically motivated, computationally motivated, or both. For instance, if physical solutions of the Schrödinger equation are known, it is possible to project the orbitals on these solutions [49–51]. The most common Schrödinger basis functions are the Slater-type basis functions [49]:and the Gaussian-type basis functions [10, 50]:where and are normalisation constants while represents the magnetic quantum numbers. Once an orbital is projected on a basis, it is entirely determined by its projection coefficients. Thus, the resulting representation is parametric. As a result, the determination of the projection coefficients is equivalent to the determination of the orbitals: the closer the basis functions are to the real solution, the more efficient the calculation is. This is due to the fact that fewer coefficients are required to adequately represent the underlying orbital. An even more realistic representation may be obtained if the basis functions are centred upon their respective nuclei [13]:where is an atomic basis function and is the location of a given nucleus.

The orbitals may also be represented by plane waves [14]:where is the volume of the cell associated with the underlying discrete grid and is the wave vector associated with the plane wave. From a physical point of view, plane waves form an appropriate basis when the orbitals are smooth functions. Otherwise, a large number of plane waves are required for an accurate reconstruction of the orbitals. Consequently, the computational burden associated with the parametric representation becomes rapidly prohibitive. Nevertheless, if the orbitals are smooth functions, one may take advantage of the Fourier transform in order to reduce the complexity involved in the calculations of the derivatives. For instance, the Laplacian of a Fourier transformed function is obtained by multiplying this function by the square of the module of the wave vector: where is the Fourier transform. If the Fourier transform is calculated with the Fast Fourier Transform algorithm, the complexity of such a calculation, for a discrete, grid, reduces towhich explains why plane wave basis and the Fourier transforms are prevalent in ab initio molecular dynamics simulations. The computational efficiency may be further improved if the Fourier transformations are evaluated with graphics processing units (GPU), as their architecture makes them particularly suited for these calculations, especially regarding speed [12].

The orbital functions are not always smooth as a result of the strength of the electronic interaction. In this particular case, the plane wave basis constitutes an inappropriate choice as a very large number of basis elements are required in order to describe the quasi-discontinuities. Nevertheless, one would be inclined to retain the low computational complexity associated with plane waves and the Fourier transform while being able to efficiently describe the rougher parts of the orbitals. This may be achieved with a wavelet basis and the wavelet transforms [20, 48, 52–54]. As opposed to the plane wave functions, which span the whole space, the wavelets are spatially localised. Therefore, they are particularly suited for the description of discontinuities and fast-varying functions. In addition, the wavelet basis provides a multiresolution representation of the orbitals which means that the calculations may be performed efficiently at the required resolution level. The wavelet functions are filter banks [55]. In one dimension, they involve two functions: (i) the scaling function:which is a high-pass filter responsible for the multiresolution aspect of the transform; (ii) the mother wavelet:which is a band-pass filter responsible for the basis localisation. The coefficients of the scaling and the mother wavelet function are related byIn three dimensions, at the lowest resolution, the scaling function is given bywhile the mother wavelets becomewhere is the resolution of the underlying computational grid. As a result, any orbital function may be approximated by a truncated expansion:where and are the wavelet coefficients. Naturally, higher resolutions may be achieved; the maximum resolution is determined by the resolution of the computational grid. As for the Fourier transform, the wavelet transform admits Fast Wavelet Transform implementation which has the same complexity as its Fourier counterpart.

In the next subsection, we briefly present a hybrid approach which involves both ab initio molecular dynamics and molecular dynamics.

10.3. Hybrid Quantum Mechanics–Molecular Dynamics Methods

Because of its inherent complexity, the applicability of ab initio molecular dynamics is essentially limited to small molecular domains. Molecular dynamics, which we recently reviewed in [30], is suitable for larger molecular system. A drawback is that it is not adapted for the simulation of chemical complexes, since it is not based on first principles as its quantum mechanics counterpart. Rather, molecular dynamics relies on empirical potentials [56] and classical mechanics which allow for the simulation of much larger molecular complexes.

Therefore, in order to simulate larger systems, hybrid approaches [12, 57–59] must be followed. In such a method, a small region of interest is simulated with ab initio molecular dynamics, while the rest of the molecular system is approximated with molecular mechanics:An extra Lagrangian term is required in order to ensure a proper coupling between the quantum degrees of freedom and the classical degrees of freedom. This Lagrangian consists of a bounded and an unbounded part. The bounded part consists of the stretching, bending, and torsional terms which are characterised by their distances, angles, and dihedral angles, respectively. The unbounded part contains the electrostatic interaction between the molecular mechanics atoms and the quantum density as well as the steric interaction which may be approximated, for instance, by a van der Waals potential.

11. Conclusions

In this paper, we have presented a comprehensive, but yet concise, review of ab initio molecular dynamics from a computational perspective and from first principles. Although it is always hazardous to speculate about the future, we foresee two important breakthroughs. Fourier transforms, which constitute an essential part of many molecular simulations, may be evaluated with high performance on graphical processing units or GPU [54]. The same remark applies to the wavelet transform whose role is essential in describing discontinuous orbitals [52]. This opens the door to high performance simulations, while tempering the limitations associated with computational complexity.

For decades, one of the main bottlenecks of molecular simulations has been the calculation of the density functionals. As in numerous fields, machine learning constitutes a promising avenue for the fast and efficient evaluation of these functionals without recourse to explicit calculations [17, 24, 25, 60]. Here, the density functional is simply learned from existing examples with machine learning techniques. This paves the way for the simulation for larger and more complex molecular systems. We plan to review machine learning-based approaches in a future publication.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

References

P. Carloni, U. Rothlisberger, and M. Parrinello, “The role and perspective of ab initio molecular dynamics in the study of biological systems,” Accounts of Chemical Research, vol. 35, no. 6, pp. 455–464, 2002.
View at: Publisher Site | Google Scholar
L. Monticelli and E. Salonen, Eds., Biomolecular Simulations–Methods and Protocols, Humana Press, Springer, London, UK, 2013.
T. D. Kühne, Ab-Initio Molecular Dynamics, 2013, https://arxiv.org/abs/1201.5945.
M. E. Tuckerman, “Ab initio molecular dynamics: Basic concepts, current trends and novel applications,” Journal of Physics: Condensed Matter, vol. 14, no. 50, pp. R1297–R1355, 2002.
View at: Publisher Site | Google Scholar
J. Broeckhove and L. Lathouwers, Eds., Time-Dependent Quantum Molecular Dynamics, NATO Advanced Research Workshop on Time Dependent Quantum Molecular Dynamics: Theory and Experiment, Springer, Snowbird, USA, 1992.
F. Jensen, Introduction to Computational Chemistry, Wiley, Chichester, Uk, 2017.
M. E. Tuckerman, Tuckerman, Statistical Mechanics: Theory and Molecular Simulation, Oxford University Press, Oxford, UK, 2013.
R. P. Steele, “Multiple-Timestep ab Initio Molecular Dynamics Using an Atomic Basis Set Partitioning,” The Journal of Physical Chemistry A, vol. 119, no. 50, pp. 12119–12130, 2015.
View at: Publisher Site | Google Scholar
N. Luehr, T. E. Markland, and T. J. Martínez, “Multiple time step integrators in ab initio molecular dynamics,” The Journal of Chemical Physics, vol. 140, no. 8, Article ID 084116, 2014.
View at: Publisher Site | Google Scholar
H. B. Schlegel, Li. X., J. M. Millam, G. A. Voth, G. E. Scuseria, and M. J. Frisch, “Ab initio molecular dynamics: Propagating the density matrix with Gaussian orbitals. III. Comparison with BornOppenheimer dynamics,” Journal of Chemical Physics, vol. 117, no. 19, pp. 9758–9763, 2001.
View at: Google Scholar
T. D. Kühne, M. Krack, F. R. Mohamed, and M. Parrinello, “Efficient and accurate car-parrinello-like approach to born-oppenheimer molecular dynamics,” Physical Review Letters, vol. 98, no. 6, Article ID 066401, 2007.
View at: Publisher Site | Google Scholar
X. Andrade, A. Castro, D. Zueco et al., “Modified ehrenfest formalism for efficient large-scale ab initio molecular dynamics,” Journal of Chemical Theory and Computation, vol. 5, no. 4, pp. 728–742, 2009.
View at: Publisher Site | Google Scholar
M. Vacher, D. Mendive-Tapia, M. J. Bearpark, and M. A. Robb, “The second-order Ehrenfest method: A practical CASSCF approach to coupled electron-nuclear dynamics,” Theoretical Chemistry Accounts, vol. 133, no. 7, pp. 1–12, 2014.
View at: Publisher Site | Google Scholar
F. Ding, J. J. Goings, H. Liu, D. B. Lingerfelt, and X. Li, “Ab initio two-component Ehrenfest dynamics,” The Journal of Chemical Physics, vol. 143, no. 11, Article ID 114105, 2015.
View at: Publisher Site | Google Scholar
K. Kitaura and K. Morokuma, “A new energy decomposition scheme for molecular interactions within the Hartree‐Fock approximation,” International Journal of Quantum Chemistry, vol. 10, no. 2, pp. 325–340, 1976.
View at: Publisher Site | Google Scholar
H. S. Yu, X. He, and D. G. Truhlar, “MN15-L: A New Local Exchange-Correlation Functional for Kohn-Sham Density Functional Theory with Broad Accuracy for Atoms, Molecules, and Solids,” Journal of Chemical Theory and Computation, vol. 12, no. 3, pp. 1280–1293, 2016.
View at: Publisher Site | Google Scholar
F. Brockherde, L. Vogt, L. Li, M. E. Tuckerman, K. Burke, and K.-R. Müller, “By-passing the Kohn-Sham equations with machine learning,” Nature Communications, vol. 8, no. 1, article no. 872, 2017.
View at: Publisher Site | Google Scholar
A. P. Selvadurai, Partial Differential Equations in Mechanics 2, Springer Berlin Heidelberg, Berlin, Heidelberg, 2000.
View at: Publisher Site
E. Engel and R. M. Dreizle, Density Functional Theory: An Advanced Course, Springer, London, UK, 2011.
T. D. Engeness and T. A. Arias, “Multiresolution analysis for efficient, high precision all-electron density-functional calculations,” Physical Review B: Condensed Matter and Materials Physics, vol. 65, no. 16, pp. 1–10, 2002.
View at: Publisher Site | Google Scholar
A. J. Cohen, P. Mori-Sánchez, and W. Yang, “Challenges for density functional theory,” Chemical Reviews, vol. 112, no. 1, pp. 289–320, 2012.
View at: Publisher Site | Google Scholar
H. Ou-Yang and M. Levy, “Theorem for functional derivatives in density-functional theory,” Physical Review A: Atomic, Molecular and Optical Physics, vol. 44, no. 1, pp. 54–58, 1991.
View at: Publisher Site | Google Scholar
P. Mori-Sánchez, A. J. Cohen, and W. Yang, “Self-interaction-free exchange-correlation functional for thermochemistry and kinetics,” The Journal of Chemical Physics, vol. 124, no. 9, Article ID 091102, 2006.
View at: Publisher Site | Google Scholar
L. Li, T. E. Baker, S. R. White, and K. Burke, “Pure density functional for strong correlation and the thermodynamic limit from machine learning,” Physical Review B: Condensed Matter and Materials Physics, vol. 94, no. 24, Article ID 245129, 2016.
View at: Publisher Site | Google Scholar
F. Pereira, K. Xiao, D. A. R. S. Latino, C. Wu, Q. Zhang, and J. Aires-De-Sousa, “Machine Learning Methods to Predict Density Functional Theory B3LYP Energies of HOMO and LUMO Orbitals,” Journal of Chemical Information and Modeling, vol. 57, no. 1, pp. 11–21, 2017.
View at: Publisher Site | Google Scholar
M. E. Tuckerman and M. Parrinello, “Integrating the Car-Parrinello equations. I. Basic integration techniques,” The Journal of Chemical Physics, vol. 101, no. 2, pp. 1302–1315, 1994.
View at: Publisher Site | Google Scholar
J. Lach, J. Goclon, and P. Rodziewicz, “Structural flexibility of the sulfur mustard molecule at finite temperature from Car-Parrinello molecular dynamics simulations,” Journal of Hazardous Materials, vol. 306, pp. 269–277, 2016.
View at: Publisher Site | Google Scholar
G. Bussi, T. Zykova-Timan, and M. Parrinello, “Isothermal-isobaric molecular dynamics using stochastic velocity rescaling,” The Journal of Chemical Physics, vol. 130, no. 7, Article ID 074101, 2009.
View at: Publisher Site | Google Scholar
J. Lahnsteiner, G. Kresse, A. Kumar, D. D. Sarma, C. Franchini, and M. Bokdam, “Room-temperature dynamic correlation between methylammonium molecules in lead-iodine based perovskites: An ab initio molecular dynamics perspective,” Physical Review B: Condensed Matter and Materials Physics, vol. 94, no. 21, Article ID 214114, 2016.
View at: Publisher Site | Google Scholar
E. Paquet and H. L. Viktor, “Molecular dynamics, monte carlo simulations, and langevin dynamics: A computational review,” BioMed Research International, vol. 2015, Article ID 183918, 2015.
View at: Publisher Site | Google Scholar
M. E. Tuckerman, Y. Liu, G. Ciccotti, and G. J. Martyna, “Non-Hamiltonian molecular dynamics: Generalizing Hamiltonian phase space principles to non-Hamiltonian systems,” The Journal of Chemical Physics, vol. 115, no. 4, pp. 1678–1702, 2001.
View at: Publisher Site | Google Scholar
M. Ceriotti, G. Bussi, and M. Parrinello, “Colored-noise thermostats àla Carte,” Journal of Chemical Theory and Computation, vol. 6, no. 4, pp. 1170–1180, 2010.
View at: Publisher Site | Google Scholar
M. Ceriotti, G. Bussi, and M. Parrinello, “Langevin equation with colored noise for constant-temperature molecular dynamics simulations,” Physical Review Letters, vol. 102, no. 2, Article ID 020601, 2009.
View at: Publisher Site | Google Scholar
M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Manolopoulos, “Efficient stochastic thermostatting of path integral molecular dynamics,” The Journal of Chemical Physics, vol. 133, no. 12, Article ID 124104, 2010.
View at: Publisher Site | Google Scholar
R. Kubo, “The fluctuation-dissipation theorem,” Reports on Progress in Physics, vol. 29, no. 1, article no. 306, pp. 255–284, 1966.
View at: Publisher Site | Google Scholar
G. J. Martyna, A. Hughes, and M. E. Tuckerman, “Molecular dynamics algorithms for path integrals at constant pressure,” The Journal of Chemical Physics, vol. 110, no. 7, pp. 3275–3290, 1999.
View at: Publisher Site | Google Scholar
T. Spura, C. John, S. Habershon, and T. D. Kühne, “Nuclear quantum effects in liquid water from path-integral simulations using an ab initio force-matching approach,” Molecular Physics, vol. 113, no. 8, pp. 808–822, 2015.
View at: Publisher Site | Google Scholar
O. Marsalek and T. E. Markland, “Ab initio molecular dynamics with nuclear quantum effects at classical cost: Ring polymer contraction for density functional theory,” The Journal of Chemical Physics, vol. 144, no. 5, Article ID 054112, 2016.
View at: Publisher Site | Google Scholar
M. E. Tuckerman, D. Marx, M. L. Klein et al., “Efficient and general algorithms for path integral Car-Parrinello molecular dynamics,” The Journal of Chemical Physics, vol. 104, no. 14, pp. 5579–5588, 1996.
View at: Google Scholar
P. Minary, G. J. Martyna, and M. E. Tuckerman, “Algorithms and novel applications based on the isokinetic ensemble. I. Biophysical and path integral molecular dynamics,” The Journal of Chemical Physics, vol. 118, no. 6, pp. 2510–2526, 2003.
View at: Publisher Site | Google Scholar
M. Ceriotti, J. More, and D. E. Manolopoulos, “I-PI: A Python interface for ab initio path integral molecular dynamics simulations,” Computer Physics Communications, vol. 185, no. 3, pp. 1019–1026, 2014.
View at: Publisher Site | Google Scholar
H. Y. Geng, “Accelerating ab initio path integral molecular dynamics with multilevel sampling of potential surface,” Journal of Computational Physics, vol. 283, pp. 299–311, 2015.
View at: Publisher Site | Google Scholar
K. A. Fichthorn and W. H. Weinberg, “Theoretical foundations of dynamical Monte Carlo simulations,” The Journal of Chemical Physics, vol. 95, no. 2, pp. 1090–1096, 1991.
View at: Publisher Site | Google Scholar
D. Marx and M. Parrinello, “Ab initio path integral molecular dynamics: Basic ideas,” The Journal of Chemical Physics, vol. 104, no. 11, pp. 4077–4082, 1995.
View at: Publisher Site | Google Scholar
O. Marsalek, P.-Y. Chen, R. Dupuis et al., “Efficient calculation 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.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
C. Hartwigsen, S. Goedecker, and J. Hutter, “Relativistic separable dual-space Gaussian pseudopotentials from H to Rn,” Physical Review B: Condensed Matter and Materials Physics, vol. 58, no. 7, pp. 3641–3662, 1998.
View at: Publisher Site | Google Scholar
L. Genovese, A. Neelov, S. Goedecker et al., “Daubechies wavelets as a basis set for density functional pseudopotential calculations,” The Journal of Chemical Physics, vol. 129, no. 1, Article ID 014109, 2008.
View at: Publisher Site | Google Scholar
M. A. Watson, N. C. Handy, and A. J. Cohen, “Density functional calculations, using Slater basis sets, with exact exchange,” The Journal of Chemical Physics, vol. 119, no. 13, pp. 6475–6481, 2003.
View at: Publisher Site | Google Scholar
X. Andrade and A. Aspuru-Guzik, “Real-space density functional theory on graphical processing units: Computational approach and comparison to Gaussian basis set methods,” Journal of Chemical Theory and Computation, vol. 9, no. 10, pp. 4360–4373, 2013.
View at: Publisher Site | Google Scholar
S. V. Levchenko, X. Ren, J. Wieferink et al., “Hybrid functionals for large periodic systems in an all-electron, numeric atom-centered basis framework,” Computer Physics Communications, vol. 192, pp. 60–69, 2015.
View at: Publisher Site | Google Scholar
T. Deutsch and L. Genovese, “Wavelets for electronic structure calculations,” Collection SFN, vol. 12, pp. 33–76, 2011.
View at: Google Scholar
A. H. Prociuk and S. S. Iyengar, “A multiwavelet treatment of the quantum subsystem in quantum wavepacket ab initio molecular dynamics through an hierarchical partitioning of momentum space,” Journal of Chemical Theory and Computation, vol. 10, no. 8, pp. 2950–2963, 2014.
View at: Publisher Site | Google Scholar
L. Genovese, B. Videau, D. Caliste, J.-F. Méhaut, S. Goedecker, and T. Deutsch, “Wavelet-Based Density Functional Theory on Massively Parallel Hybrid Architectures,” Electronic Structure Calculations on Graphics Processing Units: From Quantum Chemistry to Condensed Matter Physics, pp. 115–134, 2016.
View at: Publisher Site | Google Scholar
R. C. Gonzalez and R. E. Woods, Digital Image Processing, Pearson Prentice Hall, Upper Saddle River, NJ, USA, 2008.
K. Vanommeslaeghe, E. Hatcher, and C. Acharya, “CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields,” Journal of Computational Chemistry, vol. 31, no. 4, pp. 671–690, 2010.
View at: Publisher Site | Google Scholar
J. R. Perilla, B. C. Goh, C. K. Cassidy et al., “Molecular dynamics simulations of large macromolecular complexes,” Current Opinion in Structural Biology, vol. 31, pp. 64–74, 2015.
View at: Publisher Site | Google Scholar
N. Sahu and S. R. Gadre, “Molecular tailoring approach: A route for ab initio treatment of large clusters,” Chemical Reviews, vol. 47, no. 9, pp. 2739–2747, 2014.
View at: Publisher Site | Google Scholar
J. Dreyer, B. Giuseppe, E. Ippoliti, V. Genna, M. De Vivo, and P. Carloni, “First principles methods in biology: continuum models to hybrid ab initio quantum mechanics/molecular mechanics,” in Simulating Enzyme Reactivity: Computational Methods in Enzyme Catalysis, I. Tuñón and V. Moliner, Eds., chapter 9, Royal Society of Chemistry, London, UK, 2016.
View at: Google Scholar
V. Botu and R. Ramprasad, “Adaptive machine learning framework to accelerate ab initio molecular dynamics,” International Journal of Quantum Chemistry, vol. 115, no. 16, pp. 1074–1083, 2015.
View at: Publisher Site | Google Scholar

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.

PDF Download Citation

Download other formats

Order printed copies

Views

15046

Downloads

844

Citations