Reconsidering Car–Parrinello molecular dynamics using direct propagation of molecular orbitals developed upon Gaussian type atomic orbitals

A reconsideration of Car–Parrinello molecular dynamics using only atom centred basis functions is presented by doing direct propagation of molecular orbitals in conjunction with propagation of nuclei.

The electronic degree of freedom chosen for the propagation are the coefficients of the linear combination of atomic orbitals development where atomic orbitals are expressed upon Gaussian functions.

Considering that the wave function is propagated, we show that only very few iterations are sufficient to calculate nuclear gradients and electronic potential energy with good accuracy, instead of fulfilling the minimisation energy procedure at each step of the dynamics simulation.

First tests calculations are presented by considering simple systems (water and monothiooxalic acid molecules and a (H2O)20 water cluster).

These tests present the good behaviour of the propagation scheme if one considers the conservation of the total energy and the accuracy of the potential energy surface compared with Born–Oppenheimer calculations.


What is a chemical reaction from a theoretical point of view?

One answer to this apparently simple question may be the knowledge of the associated potential energy surface (PES) and the temperature of the reaction.

The accuracy of the calculated PES is the target of most theoretical methods.

However, exploration of PES is often based on the characterisation of singular points, such as minima and transition states.

Although some of these structures remain difficult to localise, in most cases, the search for these peculiar points relies on chemical intuition.

The characteristic values of chemical reactions (kinetic constant, Gibbs energy, etc.) deduced from the knowledge of these particular points are sometimes approximated.

For instance, estimation of entropy contributions depends on the harmonic approximation.1

Another way to explore the PES is the molecular dynamics (MD) where the evolution of the system during time is simulated.

The standard technique was initially based on the use of classical interatomic potentials parametrised to experimental or ab initio data.

Recently, the advances in ab initio molecular dynamics have been pronounced.2,3

Rather than using empirical interaction potentials, the electronic energy and its derivatives, the forces on nuclei, are computed by quantum mechanics.

In this way, Born–Oppenheimer (BO)4–7 molecular dynamics, as well as Car–Parrinello (CP)8 molecular dynamics, fall into this category.

In the BO approach the propagation involves the nuclear degrees of freedom which are treated classically.

For each nuclear configuration reached by the molecular system during the dynamics, the electronic energy and the corresponding nuclear gradients are calculated on the fly and are used to compute the next step of the dynamics.

According to this scheme, the minimisation energy procedure must be fulfilled at each step in order to obtain accurate electronic energy and gradients.

Two decades ago, Car and Parrinello proposed a new approach that involved a new degree of freedom in the dynamics, namely the wave function.8

The associated Lagrangian includes all the usual terms (the electronic and kinetic nuclear energy) plus a term that allows the classical propagation of the wave function.

Consequently, the nuclear degrees of freedom and the wave function evolve according to the Newton law and the minimisation energy procedure can be normally avoided.

Thus, the CPU time is reduced and larger molecular systems might be treated as compared to the BO approach.

Since the appearance of the CP method, most of the algorithms are based on the density functional theory (DFT) approach coupled with a plane wave basis set description of the electronic structure.

These algorithms were developed to study and predict equilibrium and non-equilibrium properties of condensed-matter systems.7

While this approach presents many advantages, the use of plane waves can not be well suited for the study of non-periodic chemical systems and to describe cusps and high electronic inhomogeneities.

On the other hand, the restriction to DFT methods is sometimes not judicious for chemical investigations.

As a matter of fact, one may prefer other quantum chemical methods (MPn, hybrid DFT, etc.) in order to increase the accuracy of the potential energy surface explored during the simulation.

It is well accepted that theoretical investigations of chemical reactivity need the use of hybrid functionals in the DFT approach.

Even if it can be viewed as a small modification of the Fock matrix, hybrid functionals have the advantage of including a large part of the electronic correlation contrary to the HF method.

If these two arguments are taken into account, the use of atom centred basis functions appears to be a good alternative and the Gaussian type basis set is a natural choice for quantum chemists.

This type of basis set has numerous advantages.

The number of basis functions needed to describe molecular orbitals with a good accuracy is far smaller than with plane wave functions.

Although the necessary integrals are more difficult to calculate, particularly two-electron integrals, efficient algorithms are available to perform this task.

Using only atom centred basis functions within the CP scheme has been initially proposed by Field9 in the context of simulated annealing calculations.

First CP dynamics using this type of basis set have been accounted for within the Hartree–Fock (HF)10–12 and generalised valence bond (GVB)13,14 methods.

Nevertheless, this scheme has shown some difficulty with energy conservation.

Moreover, Carter’s group decided to give up the CP molecular dynamics considering that it is less efficient than the BO style dynamics.14,15

Recently a different CP molecular dynamics approach has been proposed where the dynamics variables chosen to represent the electronic degree of freedom are the individual elements of the reduced one-particle density matrix expressed with Gaussian orbitals.16–18

This method has proven its ability in treating complex chemical systems with good accuracy.

We decided to reconsider the Carter et al., approach by combining the advantages of using Gaussian type orbitals for electronic structure calculation and the benefit of coupled propagation of the nuclei and the wave function.

Instead of using a density matrix scheme, we use the original CP Lagrangian.

Consequently, the dynamics variables chosen to represent the electronic degree of freedom in this approach are the molecular orbitals ψi.

By means of a Verlet scheme,19 the nuclear coordinates and the coefficients of the LCAO development of each molecular orbital are jointly propagated in order to reproduce an accurate dynamical behaviour of the molecular system.

It is important to note that the implementation presented here is not a real CP method in the sense that a partial convergence of the minimisation energy procedure is achieved at each step.

Accordingly, our approach corresponds neither to strict CP nor to the BO method.

In the present procedure, we are using the advantage of propagating the wave function by means of a CP type Lagrangian in order to reduce the computational time needed for an accurate description of the molecular dynamics.

Because of the propagation of the wave function between two different nuclear structures, a small fixed number of iterations (between two and three) is sufficient to obtain enough accuracy on gradients evaluation.

This approach is different from converging the wave function at each given propagation step (for example 50 or 100) in the sense that numerical error accumulation does not appear abruptly.

Moreover in this approach no discontinuity is observed in the total energy.

For MD simulation the accuracy of the gradient appears to be an important issue.

Moreover this partial convergence has the supplementary advantage of obtaining a dynamical accuracy that is close to the BO approach where the minimisation energy procedure is fulfilled at each step.

Thus, one can more easily consider the accurate treatment of chemical problems with a smaller computational cost.

In Section II, we outline the theoretical framework for the algorithm and detail some particular aspects linked to the method.

In Section III, first tests are presented considering illustrative simple examples: the dynamics of a single water molecule, the simulation of proton transfer reaction in monothiooxalic acid and the dynamics behaviour of (H2O)20 water cluster.

Here, we do not pretend to give insight into the chemical aspects of the problem but rather to test the efficiency in energy conservation by comparison with BO dynamics.

In Section IV, concluding remarks and future perspectives of development are presented.

General method

The Car–Parrinello method is based on an extended Lagrangian,8 in which the electrons, represented by a set of molecular orbitals {ψi(r⃑)}, follow a fictitious dynamics coupled with the motion of the nuclei, represented by a set of positions {RI}.

In an orthonormal basis of molecular orbitals, this Lagrangian can be written as:where MI are the nuclear masses and μi the fictitious masses associated with the molecular orbitals.

(Note that μi does not have the unit of mass but the unit of energy times a squared time for reasons of dimensionality).

The electronic potential energy of the system for a given nuclear structure is represented by the term ε and the holonomic constraints eqn. (2.2) are satisfied through the set of Lagrange multipliers Λij.∫dr⃑ψi*(r⃑)ψj(r⃑) = δijThe Euler–Lagrange equations derived from the Car–Parrinello Lagrangian are:

The equations of propagation may be integrated using the Verlet algorithm19 that takes the following form:

The values of Λij are fixed by imposing the orthogonality constraints of the molecular orbitals.

The algorithm for the constrained electronic equations of motion is based on the approach of Ryckaert et al.,20 and this algorithm was clearly described by Tuckerman and Parrinello.21

Using the Verlet algorithm, the unconstrained molecular orbitals evolve as a first step according to:and the nuclear coordinates are jointly propagated according to eqn. (2.3).

At this stage of the propagation the unconstrained molecular orbitals at time (t + δt) are developed upon atomic orbitals ϕν expressed in term of nuclear coordinates at time (t) and (t − δt).

Consequently, one must perform an update of the basis set for homogeneity reasons in order to develop the molecular orbitals |i(t + δ t)〉 only with atomic orbitals expressed at time (t):U = SS−1where S is the overlap matrix at time (t): Sμν = 〈ϕμ(t)|ϕν(t)〉 and S′ the overlap matrix between time (t − δt) and (t): Sμν = 〈ϕμ(t − δt)|ϕν(t)〉.

Then, the orbitals are corrected by adding the force due to orthonormality constraint:where Xij = (δt)2/2μiΛij.

The unknown Lagrange multipliers are calculated by imposing the constraint condition eqn. (2.2) on the molecular orbitals at time (tt) expressed using eqn. (2.10).

The matrix X should then satisfy the following equation:XX + XB + BX = IAwhere Aij = 〈i(t + δt)|j(t + δt)〉 and Bij = 〈ψi(t)|j(t + δt)〉.

This equation can be solved iteratively using:

with the initial guess:Eqn.

(2.12) can usually be iterated to a tolerance of 10−8 over the RMS of the matrix X within a few (4–6) iterations.21

Once the orthogonality constraints are satisfied, a new update of the basis set is performed in order to express molecular orbitals at time (t + δt) as a function of atomic orbitals at time (t + δt).

One may note that this change of basis set corresponds to the first one that will occur in the next step of the propagation.

Consequently, the algorithm required in fact only one of this type of procedure and only one matrix inversion (that is the most time consuming routine of the propagation algorithm).

This set of new molecular orbitals and nuclear positions are then used to compute the new electronic potential energy ε(t + δt) and nuclear gradients.

Our propagation code is linked to the Gaussian suite of programs22 that performs energies and forces calculation, including in particular the molecular orbitals orthonormalisation constraints derivatives,23 by means of calculating the derivatives over nuclear coordinates of electronic potential energy without any approximation.

At this level, our implementation deviates from the CP method in a sense that the electronic energy calculation is achieved after a small and fixed number of iterations.

This procedure slightly improves the wave function and allows us to obtain more accurate gradients.

Moreover these few cycles have in fact two advantages: the errors on the nuclear gradients are considerably reduced and the precision of the electronic potential energy is greatly improved leading to results close to the BO method as it will be presented in the next section.

The advantage of the method lies in the fact that the propagated molecular orbitals used as a guess for calculation of the electronic energy are closer to the variational solution than the ones used in the BO method.

Thus, when it needs two or three iterations to approach energy convergence within a tolerance of 10−4–10−5 au, the BO dynamics needs a much larger number of cycles and consequently is more time consuming.

It is also important to note that fictitious masses μi are different according to molecular orbitals.

It is well established that core orbitals would be less affected than valence orbitals during dynamics simulation.

Then, a similar mass weighting scheme as Iyengar et al.,17 is used.

The scheme chosen has the following form:μi =  if εi > −2 auμi =  |εi| if εi < −2 auwhere εi is the orbital energy of the molecular orbital ψi and  is a constant predefined by the user.

At each step of the simulation, the fictitious mass of each molecular orbital is verified according to the precedent scheme.

In order to start one simulation, initial velocity vectors have to be applied on nuclei.

Consequently, a given total kinetic energy has been initially distributed equivalently among all degrees of freedom, removing global translation and rotation of the system.

The values of initial total kinetic energy will be precised in each subsection.

Concerning the wave function, according to the Verlet algorithm, it can be evaluated at time t if the wave functions corresponding to the nuclear conformation at time (t – δt) and (t – 2δt) are known.

Thus, two converged SCF calculations (BO like) are performed at time t = 0 and (t – δt) in order to initialise the simulation.

The classical propagation of the wave function starts at time t = 2δt.

Computational tests and results

In order to illustrate the implementation of the approach discussed above, we have considered the dynamics of simple systems: a single water molecule, the monothiooxalic acid (Fig. 1) and a (H2O)20 water cluster.

As stated in the introduction, the objective of these calculations is only to test our code on different dynamical problems.

We do not intend to explore any chemical or physical problem but rather to show the good behaviour of the chosen propagation scheme especially by comparing BO trajectories with the dynamics performed with this method.

For this purpose, we will focus our attention on the quality of the PES explored by the dynamics and on the conservation of the real total energy (i.e. the kinetic nuclear energy plus the electronic potential energy).

One may note that the fictitious kinetic energies of molecular orbitals is not taken into account for the total energy analysis and all calculations presented in the present paper have been performed without any thermostat.

Water molecule

For the first system, an isolated water molecule, calculations have been performed using the Hartree–Fock (HF) method.

Hydrogen atoms have been represented by a 6-311++G(d,p) basis set and core electrons of oxygen were replaced by a pseudo-potential24 with its associated double-ζ basis set.

This molecular system and the level of theory used to describe its electronic structure has been chosen for its simplicity and for its low computational cost.

All BO trajectories presented for this system have been realized with a simple Verlet integrator.

The BO and CP difference of the potential energy during the propagation is presented in Fig. 2 where the initial conditions are identical in both simulations: the initial total kinetic energy distributed is 1.5 × 10−3 au.

The time step was δt = 0.05 fs, the fictitious mass parameter  = 0.1 amu bohr2 (≈170 au) and the total time 100 fs which represents 2000 steps.

The chosen value for the fictitious mass parameter is usual for CP simulations.7,16

The energy difference oscillates regularly around the zero value with a maximum deviation of 1.6 × 10−4 au.

This indicates that the quality of the PES explored by the CP dynamics is very close to the one explored by the BO dynamics.

In order to gain more insight into the stability of our algorithm, the behaviour of the CP dynamics was compared with the BO dynamics by increasing the propagation time step.

For this purpose, δt varies from 0.05 fs to 0.25 fs with a total simulation time of 100 fs.

The maximum deviation of the total energy, compared to the initial total energy, and its average are reported in Table 1.

For both methods, the maximum deviation increases linearly as a function of the time step δt, ranging from 5.9 × 10−4 to 2.91 × 10−3 au.

The increase of the time step by a factor of five induces an increase of the maximum error by at least an order of magnitude.

An increase of the average error of the total energy is also observed with a comparable variation.

However, the range of CP dynamics’ error still remains in very good limits (≤10−4 au) and is very close to BO dynamics’.

Thus, within the same propagator (a simple Verlet here) the CP dynamics is able to reproduce the BO dynamics with a significant saving of CPU time.

The CPU time ratio between the BO and CP dynamics for different time steps is reported in Table 1.

This indicates that even with longer time steps, propagated molecular orbitals are closer to the variational solution than the default initial guess used in the BO approach.

Moreover, the previous CP approach using only atom centred basis functions have considered that CP style dynamics is less efficient: BO can afford longer time steps because forces are exact.14,15

Results gathered in Table 1 lead us to consider CP more effective than BO at the same level of accuracy: for instance, the CP dynamics with a time step of 0.25 fs is two times computationally less expensive than a BO dynamics with the same time step.

Then, one should used a time step of 0.5 fs for BO dynamics in order to obtain a similar CPU cost.

However, saving of computational time by using a larger time step (i.e. 0.5 fs) within BO dynamics induced larger errors compared with CP dynamics using a time step of 0.25 fs: the total energy average is 317 × 10−6 au and the maximum deviation of the total energy is 577 × 10−5.

One of the main features of our propagation routine is the partial convergence of the minimisation energy procedure during the calculation of the electronic structure at a given geometry.

The propagation scheme has been confronted to a “restricted” BO dynamics where the number of SCF cycles is restricted to the one adopted in the CP approach.

Contrary to the CP dynamics, the restricted BO method shows a non-conservation of the total energy, as one can expect.

In the BO approach, the molecular orbitals guess vectors used to calculate the electronic potential energy and nuclear gradients at time (t + δt) come from the calculation performed with the time (t) nuclear structure.

This guess vector corresponds to the variational solution adapted to the molecular structure expressed at time (t).

Consequently it is not fully adapted to the one expressed at time (t + δt).

In the CP approach, the molecular orbitals are propagated together with the nuclei.

One can then consider that the molecular orbitals guess vectors are closer to the variational solution than in the BO approach.

Consequently, a small number of iterations are sufficient to reach a quasi-convergence of the electronic energy and the nuclear gradients are accurately calculated.

The monothiooxalic acid

The second molecule used to test our algorithm is the monothiooxalic acid for which simple and double intra molecular proton transfer reaction can occur.

This system has been chosen in order to test the ability of our program to reproduce more complicated dynamical behaviour (especially bond breaking).

This system has already led to theoretical investigations:25 the PES has been explored and stationary points have been mostly localised.

Other intra molecular proton transfer reactions have been theoretically studied26,27 and have concluded that dynamical investigations are required to especially discriminate between the concerted or not concerted double proton transfer mechanism.

In many cases, quantum dispersion and tunnelling effects play an important role in proton transfer processes and consequently impose the introduction of nuclear quantum effects in the dynamical treatment.

In our approach, the nuclei are treated classically and we do not pretend to give any answer to the different questions that concern this chemical aspect.

More precise dynamical investigations of this molecular system will be considered in future works.

All electronic structure calculations have been performed with the hybrid DFT method, using the B3LYP exchange correlation functional.28,29

Core electrons of carbon, oxygen and sulfur were replaced by pseudo-potentials with their associated double-ζ basis set.24

The first test corresponds to BO and CP trajectories calculated for a total time of 50 fs with a time step of 0.1 fs.

The fictitious mass parameter was  = 0.1 amu bohr2 and a 6-311++G(d,p) basis set was employed for hydrogen.

For each simulation, the same molecular structure has been used as a starting point: the initial geometry was near a second order saddle point, that corresponds to a double proton transfer (O⋯H⋯O and S⋯H⋯O).

This configuration is less stable by 0.1 au than the most stable conformation.

The initial total kinetic energy was 5.0 × 10−3 a.u.

One should stress that CP dynamics using the simple Verlet propagator are now compared with BO dynamics within a Velocity Verlet scheme that is well known to be more robust and more efficient.

The conservation of total energy along the trajectory, shown in Fig. 3 and in Table 2, reflects the accuracy of the CP method compared to the BO approach.

In this more complicated system the degree of accuracy is lower than in the case of a water molecule.

The maximum deviation of the total energy is 8.6 × 10−3 au and 2.9 × 10−3 au, the average deviation of total energy is 1.7 × 10−4 au and 5.4 × 10−5 au for CP and BO methods respectively.

In order to detail the contribution of each component to the discrepancy on the total energy, CP and BO electronic energy and nuclear kinetic energy differences are reported in Fig. 4 and 5 and in Table 3.

The kinetic nuclear energy difference, 5.5 × 10−3 au, is one order of magnitude larger than for the electronic energy which is 2.5 × 10−4 au.

This indicates that differences on the wave function have a larger effect on nuclear gradients than on the electronic energy during the propagation as one can expect.

Moreover, as shown in Fig. 4, the difference of potential energy between BO and CP increases with time.

Indeed a CP trajectory should not exactly reproduce a BO trajectory: the molecular orbitals are not obtained from converged self-consistent field calculations and then, errors are committed on nuclear gradients.

As a consequence, explored areas of the potential energy surface are not strictly identical.

The increase of difference between BO and CP kinetic nuclear energies is less pronounced than on electronic energies.

In order to gauge the accuracy of the method and to demonstrate that the paths followed by BO and CP dynamics are close enough, the root mean square difference (rmsd) of all atomic distances di,j (see eqn. (3.16)) has been calculated as a function of time and is reported in Fig. 6.One can observe that nuclear structures between BO and CP trajectories move away during the simulation: the rmsd of all lengths increase with time, the maximum and the average are respectively 3.9 × 10−4 Å and 1.9 × 10−4 Å.

This fact confirms that the explored path by the CP dynamics is close to the BO one, but not strictly the same.

As a consequence, differences between BO and CP, in terms of energies and nuclear structures may increase with time (see Figs. 4–6).

In order to obtain more insight into the efficiency of our method, the behaviour of the CP dynamics was tested for different hydrogen basis sets.

Different descriptions of hydrogen induce changes in the curvature of the PES and we wish to test if the method is able to reproduce this effect.

Results are compared with BO trajectories with the same propagation characteristics, except the integrator.

For all trajectories, the total time of simulation is 50 fs with a time step of 0.1 fs, the fictitious mass parameter  was 0.1 amu bohr2 for CP simulations.

Differences between CP and BO dynamics for different basis sets are shown in Tables 2 and 3.

Each component of the discrepancies on the total energy are in the same order of magnitude.

The maximum difference of the kinetic nuclear energy is one order of magnitude larger (≈10−3 au) than the maximum difference of electronic energy (≈10−4 au).

In the same way, the maximum deviations of the total energy are in the same order of magnitude for the different basis sets used for CP trajectories.

With respect to the BO simulations, maximum deviations are three times larger, but they still remain in acceptable limits.

This demonstrates that the CP scheme suits different curvatures of PES with a significant degree of accuracy.

As stated in the previous paragraph, the total energy is well conserved in the CP approach.

For the different basis sets used, the average of the total energy is of the order of 10−4 au for CP trajectories which is one order of magnitude larger than for BO trajectories.

A (H2O)20 water cluster

Finally, we have investigated a larger system compared to the two previous examples by performing one dynamics trajectory on a (H2O)20 water cluster (Fig. 7).

A double-ζ valence type orbital basis set augmented by polarisation functions on all atoms, namely 3-21G(d,p), has been used leading to 380 primitives Gaussian functions.

The initial total kinetic energy was 0.17 au.

The simulation has been run for 500 fs using a 0.25 fs time step (2000 steps).

This trajectory can be regarded as a significant test case to evaluate the performance of the algorithm.

Indeed, regardless of its apparent simplicity, the corresponding potential energy surface is essentially governed by weak interactions (mainnly hydrogen bonding) which is then rather complicated for ab initio dynamics simulation.

In this example, one should notice that the comparison with BO trajectory is more expensive, the BO simulation leading to prohibitive computational time.

If we consider the average total energy deviation with respect to the initial one (3.6 × 10−4 au) and the maximum energy deviation (1.1 × 10−3 au), these results are almost comparable to the one obtained on smaller systems (see above).

This final computational test makes us confident on the robustness and the applicability of our algorithm.


A reconsideration of CP molecular dynamics using only atom centred basis functions is presented by combining the advantages of using Gaussian type orbitals for electronic structure calculation and the benefit of coupled propagation of the nuclei and the wave function.

Instead of using a density matrix approach, the dynamics variables chosen to represent the electronic degrees of freedom are the coefficients of the LCAO development.

Although previous works have concluded to give up the CPMD style dynamics with Gaussian type orbitals, we show in this paper the efficiency of the algorithm.

With respect to BO simulations, the total energy is conserved with a good degree of accuracy and the BO trajectories have been almost reproduced by the CP scheme with significant computational time profits.

This glimpses dynamical studies of bigger chemical systems within an acceptable computational time.

With more sophisticated tools, the propagation scheme may be more efficient, as it has already implemented in all other molecular dynamics suites of programs.

The time step may be increased with the use of symplectic reversible integrators,30 constant temperature molecular dynamics may be envisaged with Nosé–Hoover chains of thermostats,21 as well as the determination of free energy pathways using constrained molecular dynamics.31

These aspects are currently under development.