Ab Initio Molecular Dynamics

Ab Initio Molecular Dynamics

(Parte 1 de 7)

John von Neumann Institute for Computing

Ab Initio Molecular Dynamics: Theory and Implementation

Dominik Marx and Jurg Hutter published in

Modern Methods and Algorithms of Quantum Chemistry, Proceedings, Second Edition, J. Grotendorst (Ed.), John von Neumann Institute for Computing, Julich, NIC Series, Vol. 3, ISBN 3-0-05834-6, p. 329-477, 200.

c© 2000 by John von Neumann Institute for Computing Permission to make digital or hard copies of portions of this work for personal or classroom use is granted provided that the copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise requires prior specific permission by the publisher mentioned above.


Lehrstuhl fur Theoretische Chemie, Ruhr–Universitat Bochum

Universitatsstrasse 150, 44780 Bochum

Germany E-mail: dominik.marx@theochem.ruhr-uni-bochum.de

Organisch–chemisches Institut, Universitat Zurich

Winterthurerstrasse 190, 8057 Zurich

Switzerland E-mail: hutter@oci.unizh.ch

The rapidly growing field of ab initio molecular dynamics is reviewed in the spirit of a series of lectures given at the Winterschool 2000 at the John von Neumann Institute for Computing, Julich. Several such molecular dynamics schemes are compared which arise from following various approximations to the fully coupled Schrodinger equation for electrons and nuclei. Special focus is given to the Car– Parrinello method with discussion of both strengths and weaknesses in addition to its range of applicability. To shed light upon why the Car–Parrinello approach works several alternate perspectives of the underlying ideas are presented. The implementation of ab initio molecular dynamics within the framework of plane wave–pseudopotential density functional theory is given in detail, including diagonalization and minimization techniques as required for the Born–Oppenheimer variant. Efficient algorithms for the most important computational kernel routines are presented. The adaptation of these routines to distributed memory parallel computers is discussed using the implementation within the computer code CPMD as an example. Several advanced techniques from the field of molecular dynamics, (constant temperature dynamics, constant pressure dynamics) and electronic structure theory (free energy functional, excited states) are introduced. The combination of the path integral method with ab initio molecular dynamics is presented in detail, showing its limitations and possible extensions. Finally, a wide range of applications from materials science to biochemistry is listed, which shows the enormous potential of ab initio molecular dynamics for both explaining and predicting properties of molecules and materials on an atomic scale.

1 Setting the Stage: Why Ab Initio Molecular Dynamics ?

Classical molecular dynamics using “predefined potentials”, either based on empirical data or on independent electronic structure calculations, is well established as a powerful tool to investigate many–body condensed matter systems. The broadness, diversity, and level of sophistication of this technique is documented in several monographs as well as proceedings of conferences and scientific schools 12,135,270,217,69,59,177. At the very heart of any molecular dynamics scheme is the question of how to describe – that is in practice how to approximate – the interatomic interactions. The traditional route followed in molecular dynamics is to determine these potentials in advance. Typically, the full interaction is broken up into two–body, three–body and many–body contributions, long–range and short– range terms etc., which have to be represented by suitable functional forms, see

Sect. 2 of Ref. 253 for a detailed account. After decades of intense research, very elaborate interaction models including the non–trivial aspect to represent them analytically were devised 253,539,584.

Despite overwhelming success – which will however not be praised in this review – the need to devise a “fixed model potential” implies serious drawbacks, see the introduction sections of several earlier reviews 513,472 for a more complete digression on these aspects. Among the most delicate ones are systems where (i) many different atom or molecule types give rise to a myriad of different interatomic interactions that have to be parameterized and / or (i) the electronic structure and thus the bonding pattern changes qualitatively in the course of the simulation. These systems can be called “chemically complex”.

The reign of traditional molecular dynamics and electronic structure methods was greatly extended by the family of techniques that is called here “ab initio molecular dynamics”. Other names that are currently in use are for instance Car– Parrinello, Hellmann–Feynman, first principles, quantum chemical, on–the–fly, direct, potential–free, quantum, etc. molecular dynamics. The basic idea underlying every ab initio molecular dynamics method is to compute the forces acting on the nuclei from electronic structure calculations that are performed “on–the–fly” as the molecular dynamics trajectory is generated. In this way, the electronic variables are not integrated out beforehand, but are considered as active degrees of freedom. This implies that, given a suitable approximate solution of the many–electron problem, also “chemically complex” systems can be handled by molecular dynamics. But this also implies that the approximation is shifted from the level of selecting the model potential to the level of selecting a particular approximation for solving the Schrodinger equation.

Applications of ab initio molecular dynamics are particularly widespread in materials science and chemistry, where the aforementioned difficulties (i) and (i) are particularly severe. A collection of problems that were already tackled by ab initio molecular dynamics including the pertinent references can be found in Sect. 5. The power of this novel technique lead to an explosion of the activity in this field in terms of the number of published papers. The locus can be located in the late–eighties, see the squares in Fig. 1 that can be interpreted as a measure of the activity in the area of ab initio molecular dynamics. As a matter of fact the time evolution of the number of citations of a particular paper, the one by Car and Parrinello from 1985 entitled “Unified Approach for Molecular Dynamics and Density–Functional Theory” 108, parallels the trend in the entire field, see the circles in Fig. 1. Thus, the resonance that the Car and Parrinello paper evoked and the popularity of the entire field go hand in hand in the last decade. Incidentally, the 1985 paper by Car and Parrinello is the last one included in the section “Trends and Prospects” in the reprint collection of “key papers” from the field of atomistic computer simulations 135. That the entire field of ab initio molecular dynamics has grown mature is also evidenced by a separate PACS classification number (71.15.Pd “Electronic Structure: Molecular dynamics calculations (Car–Parrinello) and other numerical simulations”) that was introduced in 1996 into the Physics and Astronomy Classification Scheme 486. Despite its obvious advantages, it is evident that a price has to be payed for



Figure 1. Publication and citation analysis. Squares: number of publications which appeared up to the year n that contain the keyword “ab initio molecular dynamics” (or synonyma such as “first principles MD”, “Car–Parrinello simulations” etc.) in title, abstract or keyword list. Circles: number of publications which appeared up to the year n that cite the 1985 paper by Car and Parrinello 108 (including misspellings of the bibliographic reference). Self–citations and self–papers are excluded, i.e. citations of Ref. 108 in their own papers and papers coauthored by R. Car and / or M. Parrinello are not considered in the respective statistics. The analysis is based on the CAPLUS (“Chemical Abstracts Plus”), INSPEC (“Physics Abstracts”), and SCI (“Science Citation Index”) data bases at STN International. Updated statistics from Ref. 405.

putting molecular dynamics on ab initio grounds: the correlation lengths and relaxation times that are accessible are much smaller than what is affordable via standard molecular dynamics. Another appealing feature of standard molecular dynamics is less evident, namely the “experimental aspect of playing with the potential”. Thus, tracing back the properties of a given system to a simple physical picture or mechanism is much harder in ab initio molecular dynamics. The bright side is that new phenomena, which were not forseen before starting the simulation, can simply happen if necessary. This gives ab initio molecular dynamics a truly predictive power.

Ab initio molecular dynamics can also be viewed from another corner, namely from the field of classical trajectory calculations 649,541. In this approach, which has its origin in gas phase molecular dynamics, a global potential energy surface is constructed in a first step either empirically or based on electronic structure calculations. In a second step, the dynamical evolution of the nuclei is generated by using classical mechanics, quantum mechanics or semi / quasiclassical approx- imations of various sorts. In the case of using classical mechanics to describe the dynamics – the focus of the present overview – the limiting step for large systems is the first one, why so? There are 3N − 6 internal degrees of freedom that span the global potential energy surface of an unconstrained N–body system. Using for simplicity 10 discretization points per coordinate implies that of the order of 103N−6 electronic structure calculations are needed in order to map such a global potential energy surface. Thus, the computational workload for the first step grows roughly like ∼ 10N with increasing system size. This is what might be called the “dimensionality bottleneck” of calculations that rely on global potential energy surfaces, see for instance the discussion on p. 420 in Ref. 254.

What is needed in ab initio molecular dynamics instead? Suppose that a useful trajectory consists of about 10M molecular dynamics steps, i.e. 10M electronic structure calculations are needed to generate one trajectory. Furthermore, it is assumed that 10n independent trajectories are necessary in order to average over different initial conditions so that 10M+n ab initio molecular dynamics steps are required in total. Finally, it is assumed that each single–point electronic structure calculation needed to devise the global potential energy surface and one ab initio molecular dynamics time step requires roughly the same amount of cpu time. Based on this truly simplistic order of magnitude estimate, the advantage of ab initio molecular dynamics vs. calculations relying on the computation of a global potential energy surface amounts to about 103N−6−M−n. The crucial point is that for a given statistical accuracy (that is for M and n fixed and independent on N) and for a given electronic structure method, the computational advantage of “on–the–fly” approaches grows like ∼ 10N with system size.

Of course, considerable progress has been achieved in trajectory calculations by carefully selecting the discretization points and reducing their number, choosing sophisticated representations and internal coordinates, exploiting symmetry etc. but basically the scaling ∼ 10N with the number of nuclei remains a problem. Other strategies consist for instance in reducing the number of active degrees of freedom by constraining certain internal coordinates, representing less important ones by a (harmonic) bath or friction, or building up the global potential energy surface in terms of few–body fragments. All these approaches, however, invoke approximations beyond the ones of the electronic structure method itself. Finally, it is evident that the computational advantage of the “on–the–fly” approaches diminish as more and more trajectories are needed for a given (small) system. For instance extensive averaging over many different initial conditions is required in order to calculate quantitatively scattering or reactive cross sections. Summarizing this discussion, it can be concluded that ab initio molecular dynamics is the method of choice to investigate large and “chemically complex” systems.

Quite a few review articles dealing with ab initio molecular dynamics appeared in the nineties 513,223,472,457,224,158,643,234,463,538,405 and the interested reader is referred to them for various complementary viewpoints. In the present overview article, emphasis is put on both broadness of the approaches and depth of the presentation. Concerning the broadness, the discussion starts from the Schrodinger equation. Classical, Ehrenfest, Born–Oppenheimer, and Car–Parrinello molecular dynamics are “derived” from the time–dependent mean–field approach that is ob- tained after separating the nuclear and electronic degrees of freedom. The most extensive discussion is related to the features of the basic Car–Parrinello approach but all three ab initio approaches to molecular dynamics are contrasted and partly compared. The important issue of how to obtain the correct forces in these schemes is discussed in some depth. The most popular electronic structure theories implemented within ab initio molecular dynamics, density functional theory in the first place but also the Hartree–Fock approach, are sketched. Some attention is also given to another important ingredient in ab initio molecular dynamics, the choice of the basis set.

Concerning the depth, the focus of the present discussion is clearly the implementation of both the basic Car–Parrinello and Born–Oppenheimer molecular dynamics schemes in the CPMD package 142. The electronic structure approach in CPMD is Hohenberg–Kohn–Sham density functional theory within a plane wave / pseudopotential implementation and the Generalized Gradient Approximation. The formulae for energies, forces, stress, pseudopotentials, boundary conditions, optimization procedures, parallelization etc. are given for this particular choice to solve the electronic structure problem. One should, however, keep in mind that a variety of other powerful ab initio molecular dynamics codes are available (for instance CASTEP 116, CP-PAW 143, fhi98md 189, NWChem 446, VASP 663) which are partly based on very similar techniques. The classic Car–Parrinello approach 108 is then extended to other ensembles than the microcanonical one, other electronic states than the ground state, and to a fully quantum–mechanical representation of the nuclei. Finally, the wealth of problems that can be addressed using ab initio molecular dynamics is briefly sketched at the end, which also serves implicitly as the “Summary and Conclusions” section.

2 Basic Techniques: Theory

2.1 Deriving Classical Molecular Dynamics

The starting point of the following discussion is non–relativistic quantum mechanics as formalized via the time–dependent Schrodinger equation in its position representation in conjunction with the standard Hamiltonian

I,i for the electronic {ri} and nuclear {RI} degrees of freedom. The more convenient atomic units (a.u.) will be introduced at a later stage for reasons that will soon become clear. Thus, only the bare electron–electron, electron–nuclear, and nuclear– nuclear Coulomb interactions are taken into account.

The goal of this section is to derive classical molecular dynamics 12,270,217 starting from Schrodinger’s wave equation and following the elegant route of Tully 650,651. To this end, the nuclear and electronic contributions to the total wavefunction Φ({ri},{RI};t), which depends on both the nuclear and electronic coordinates, have to be separated. The simplest possible form is a product ansatz

where the nuclear and electronic wavefunctions are separately normalized to unity at every instant of time, i.e. 〈χ;t|χ;t〉 = 1 and 〈Ψ;t|Ψ;t〉 = 1, respectively. In addition, a convenient phase factor

was introduced at this stage such that the final equations will look nice; ∫ drdR

refers to the integration over all i = 1,and I = 1,... variables {ri} and {RI},

respectively. It is mentioned in passing that this approximation is called a one– determinant or single–configuration ansatz for the total wavefunction, which at the end must lead to a mean–field description of the coupled dynamics. Note also that this product ansatz (excluding the phase factor) differs from the Born–Oppenheimer ansatz 340,350 for separating the fast and slow variables

even in its one–determinant limit, where only a single electronic state k (evaluated for the nuclear configuration {RI}) is included in the expansion. Inserting the separation ansatz Eq. (3) into Eqs. (1)–(2) yields (after multiplying from the left by 〈Ψ| and 〈χ| and imposing energy conservation d〈H〉/dt ≡ 0) the following relations

This set of coupled equations defines the basis of the time–dependent self–consistent field (TDSCF) method introduced as early as 1930 by Dirac 162, see also Ref. 158. Both electrons and nuclei move quantum–mechanically in time–dependent effective potentials (or self–consistently obtained average fields) obtained from appropriate averages (quantum mechanical expectation values 〈...〉) over the other class of degrees of freedom (by using the nuclear and electronic wavefunctions, respectively). Thus, the single–determinant ansatz Eq. (3) produces, as already anticipated, a mean–field description of the coupled nuclear–electronic quantum dynamics. This is the price to pay for the simplest possible separation of electronic and nuclear variables.

The next step in the derivation of classical molecular dynamics is the task to approximate the nuclei as classical point particles. How can this be achieved in the framework of the TDSCF approach, given one quantum–mechanical wave equation describing all nuclei? A well–known route to extract classical mechanics from quantum mechanics in general starts with rewriting the corresponding wavefunction in terms of an amplitude factor A and a phase S which are both considered to be real and A〉0 in this polar representation, see for instance Refs. 163,425,535. After transforming the nuclear wavefunction in Eq. (7) accordingly and after separating the real and imaginary parts, the TDSCF equation for the nuclei

is (exactly) re–expressed in terms of the new variables A and S. This so–called “quantum fluid dynamical representation” Eqs. (9)–(10) can actually be used to solve the time–dependent Schrodinger equation 160. The relation for A, Eq. (10), can be rewritten as a continuity equation 163,425,535 with the help of the identification of the nuclear density |χ|2 ≡ A2 as directly obtained from the definition

Eq. (8). This continuity equation is independent of and ensures locally the con- servation of the particle probability |χ|2 associated to the nuclei in the presence of a flux.

More important for the present purpose is a more detailed discussion of the relation for S, Eq. (9). This equation contains one term that depends on , a contribution that vanishes if the classical limit

is taken as → 0; an expansion in terms of would lead to a hierarchy of semi- classical methods 425,259. The resulting equation is now isomorphic to equations of motion in the Hamilton–Jacobi formulation 244,540

of classical mechanics with the classical Hamilton function defined in terms of (generalized) coordinates {RI} and their conjugate momenta {PI}. With the help of the connecting transformation dPI

can be read off. Thus, the nuclei move according to classical mechanics in an effective potential V Ee due to the electrons. This potential is a function of only the nuclear positions at time t as a result of averaging He over the electronic degrees of freedom, i.e. computing its quantum expectation value 〈Ψ|He|Ψ〉, while keeping the nuclear positions fixed at their instantaneous values {RI(t)}. However, the nuclear wavefunction still occurs in the TDSCF equation for the electronic degrees of freedom and has to be replaced by the positions of the nuclei for consistency. In this case the classical reduction can be achieved simply by replacing

classical nuclei as given by Eq. (15). This yields e.g. for the position operator∫ the required expectation value. This classical limit leads to a time–dependent wave equation for the electrons

which evolve self–consistently as the classical nuclei are propagated via Eq. (15).

Note that now He and thus Ψ depend parametrically on the classical nuclear posi- tions {RI(t)} at time t through Vn−e({ri},{RI(t)}). This means that feedback between the classical and quantum degrees of freedom is incorporated in both directions (at variance with the “classical path” or Mott non–SCF approach to dynamics 650,651).

The approach relying on solving Eq. (15) together with Eq. (18) is sometimes called “Ehrenfest molecular dynamics” in honor of Ehrenfest who was the first to address the question a of how Newtonian classical dynamics can be derived from Schrodinger’s wave equation 174. In the present case this leads to a hybrid or mixed approach because only the nuclei are forced to behave like classical particles, whereas the electrons are still treated as quantum objects.

Although the TDSCF approach underlying Ehrenfest molecular dynamics clearly is a mean–field theory, transitions between electronic states are included aThe opening statement of Ehrenfest’s famous 1927 paper 174 reads: “Es ist wunschenswert, die folgende Frage moglichst elementar beantworten zu konnen: Welcher Ruckblick ergibt sich vom Standpunkt der Quantenmechanik auf die Newtonschen Grundgleichungen der klassischen Mechanik?” in this scheme. This can be made evident by expanding the electronic wavefunction Ψ (as opposed to the total wavefunction Φ according to Eq. (5)) in terms of many electronic states or determinants Ψk

(Parte 1 de 7)