1Reconsidering Car–Parrinello molecular dynamics using direct propagation of molecular orbitals developed upon Gaussian type atomic orbitals
Type: Goal |
Advantage: None |
Novelty: None |
ConceptID: Goa1
1
Reconsidering Car–Parrinello molecular dynamics using direct propagation of molecular orbitals developed upon Gaussian type atomic orbitals
Type: Goal |
Advantage: None |
Novelty: None |
ConceptID: Goa1
2
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.
Type: Goal |
Advantage: None |
Novelty: None |
ConceptID: Goa1
3
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.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod1
4
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.
Type: Conclusion |
Advantage: None |
Novelty: None |
ConceptID: Con1
5
First tests calculations are presented by considering simple systems (water and monothiooxalic acid molecules and a (H2O)20 water cluster).
Type: Experiment |
Advantage: None |
Novelty: None |
ConceptID: Exp1
6
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.
Type: Conclusion |
Advantage: None |
Novelty: None |
ConceptID: Con2
Introduction
7
What is a chemical reaction from a theoretical point of view?
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac1
8
One answer to this apparently simple question may be the knowledge of the associated potential energy surface (PES) and the temperature of the reaction.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac1
9
The accuracy of the calculated PES is the target of most theoretical methods.
Type: Motivation |
Advantage: None |
Novelty: None |
ConceptID: Mot1
10
However, exploration of PES is often based on the characterisation of singular points, such as minima and transition states.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac2
11
Although some of these structures remain difficult to localise, in most cases, the search for these peculiar points relies on chemical intuition.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac2
12
The characteristic values of chemical reactions (kinetic constant, Gibbs energy, etc.) deduced from the knowledge of these particular points are sometimes approximated.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac3
13
For instance, estimation of entropy contributions depends on the harmonic approximation.1
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac4
14
Another way to explore the PES is the molecular dynamics (MD) where the evolution of the system during time is simulated.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac5
15
The standard technique was initially based on the use of classical interatomic potentials parametrised to experimental or ab initio data.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac5
16
Recently, the advances in ab initio molecular dynamics have been pronounced.2,3
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac6
17
Rather than using empirical interaction potentials, the electronic energy and its derivatives, the forces on nuclei, are computed by quantum mechanics.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac6
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac7
19
In the BO approach the propagation involves the nuclear degrees of freedom which are treated classically.
Type: Method |
Advantage: None |
Novelty: Old |
ConceptID: Met1
20
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.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac8
21
According to this scheme, the minimisation energy procedure must be fulfilled at each step in order to obtain accurate electronic energy and gradients.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac8
22
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
Type: Method |
Advantage: None |
Novelty: Old |
ConceptID: Met2
23
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.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod2
24
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.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod2
25
Thus, the CPU time is reduced and larger molecular systems might be treated as compared to the BO approach.
Type: Method |
Advantage: Yes |
Novelty: Old |
ConceptID: Met2
26
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.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod2
27
These algorithms were developed to study and predict equilibrium and non-equilibrium properties of condensed-matter systems.7
Type: Method |
Advantage: None |
Novelty: Old |
ConceptID: Met2
28
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.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod3
29
On the other hand, the restriction to DFT methods is sometimes not judicious for chemical investigations.
Type: Method |
Advantage: No |
Novelty: Old |
ConceptID: Met2
30
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.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac9
31
It is well accepted that theoretical investigations of chemical reactivity need the use of hybrid functionals in the DFT approach.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac9
32
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.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod4
33
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.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac10
34
This type of basis set has numerous advantages.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac10
35
The number of basis functions needed to describe molecular orbitals with a good accuracy is far smaller than with plane wave functions.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod5
36
Although the necessary integrals are more difficult to calculate, particularly two-electron integrals, efficient algorithms are available to perform this task.
Type: Method |
Advantage: None |
Novelty: Old |
ConceptID: Met2
37
Using only atom centred basis functions within the CP scheme has been initially proposed by Field9 in the context of simulated annealing calculations.
Type: Method |
Advantage: None |
Novelty: Old |
ConceptID: Met2
Type: Method |
Advantage: None |
Novelty: Old |
ConceptID: Met2
39
Nevertheless, this scheme has shown some difficulty with energy conservation.
Type: Method |
Advantage: No |
Novelty: Old |
ConceptID: Met2
40
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
Type: Method |
Advantage: No |
Novelty: Old |
ConceptID: Met2
41
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
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod6
42
This method has proven its ability in treating complex chemical systems with good accuracy.
Type: Method |
Advantage: None |
Novelty: Old |
ConceptID: Met1
43
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
44
Instead of using a density matrix scheme, we use the original CP Lagrangian.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod7
45
Consequently, the dynamics variables chosen to represent the electronic degree of freedom in this approach are the molecular orbitals ψi.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod8
46
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.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod9
47
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
48
Accordingly, our approach corresponds neither to strict CP nor to the BO method.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac11
49
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.
Type: Method |
Advantage: Yes |
Novelty: New |
ConceptID: Met3
50
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
51
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.
Type: Method |
Advantage: Yes |
Novelty: New |
ConceptID: Met3
52
Moreover in this approach no discontinuity is observed in the total energy.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs1
53
For MD simulation the accuracy of the gradient appears to be an important issue.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac12
54
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.
Type: Method |
Advantage: Yes |
Novelty: New |
ConceptID: Met3
55
Thus, one can more easily consider the accurate treatment of chemical problems with a smaller computational cost.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac13
56
In Section II, we outline the theoretical framework for the algorithm and detail some particular aspects linked to the method.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac14
57
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.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac15
58
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.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac16
59
In Section IV, concluding remarks and future perspectives of development are presented.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac17
General method
60
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}.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod10
61
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.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod10
62
(Note that μi does not have the unit of mass but the unit of energy times a squared time for reasons of dimensionality).
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod10
63
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:
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod11
64
The equations of propagation may be integrated using the Verlet algorithm19 that takes the following form:
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod12
65
The values of Λij are fixed by imposing the orthogonality constraints of the molecular orbitals.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod12
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac18
67
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).
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod13
68
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).
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod13
69
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 = S′S−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)〉.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod14
70
Then, the orbitals are corrected by adding the force due to orthonormality constraint:where Xij = (δt)2/2μiΛij.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod15
71
The unknown Lagrange multipliers are calculated by imposing the constraint condition eqn. (2.2) on the molecular orbitals at time (t +δ t) expressed using eqn. (2.10).
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod16
72
The matrix X should then satisfy the following equation:XX† + XB + B†X† = I − Awhere Aij = 〈i(t + δt)|j(t + δt)〉 and Bij = 〈ψi(t)|j(t + δt)〉.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod17
73
This equation can be solved iteratively using:
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod18
74
with the initial guess:Eqn.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod19
75
(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
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod20
76
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).
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod21
77
One may note that this change of basis set corresponds to the first one that will occur in the next step of the propagation.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod22
78
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).
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod22
79
This set of new molecular orbitals and nuclear positions are then used to compute the new electronic potential energy ε(t + δt) and nuclear gradients.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod23
80
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.
Type: Experiment |
Advantage: None |
Novelty: None |
ConceptID: Exp2
81
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
82
This procedure slightly improves the wave function and allows us to obtain more accurate gradients.
Type: Method |
Advantage: Yes |
Novelty: New |
ConceptID: Met3
83
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.
Type: Method |
Advantage: Yes |
Novelty: New |
ConceptID: Met3
84
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.
Type: Method |
Advantage: Yes |
Novelty: New |
ConceptID: Met3
85
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.
Type: Method |
Advantage: Yes |
Novelty: New |
ConceptID: Met3
86
It is also important to note that fictitious masses μi are different according to molecular orbitals.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
87
It is well established that core orbitals would be less affected than valence orbitals during dynamics simulation.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac19
88
Then, a similar mass weighting scheme as Iyengar et al.,17 is used.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac20
89
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
90
At each step of the simulation, the fictitious mass of each molecular orbital is verified according to the precedent scheme.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
91
In order to start one simulation, initial velocity vectors have to be applied on nuclei.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
92
Consequently, a given total kinetic energy has been initially distributed equivalently among all degrees of freedom, removing global translation and rotation of the system.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
93
The values of initial total kinetic energy will be precised in each subsection.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac21
94
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
95
Thus, two converged SCF calculations (BO like) are performed at time t = 0 and (t – δt) in order to initialise the simulation.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
96
The classical propagation of the wave function starts at time t = 2δt.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
Computational tests and results
97
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.
Type: Experiment |
Advantage: None |
Novelty: None |
ConceptID: Exp1
98
As stated in the introduction, the objective of these calculations is only to test our code on different dynamical problems.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac16
99
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.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac16
100
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).
Type: Goal |
Advantage: None |
Novelty: None |
ConceptID: Goa2
101
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met4
Water molecule
102
For the first system, an isolated water molecule, calculations have been performed using the Hartree–Fock (HF) method.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod24
103
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.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod24
104
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.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac22
105
All BO trajectories presented for this system have been realized with a simple Verlet integrator.
Type: Method |
Advantage: None |
Novelty: Old |
ConceptID: Met5
106
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.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs2
107
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met4
108
The chosen value for the fictitious mass parameter is usual for CP simulations.7,16
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac23
109
The energy difference oscillates regularly around the zero value with a maximum deviation of 1.6 × 10−4 au.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs2
110
This indicates that the quality of the PES explored by the CP dynamics is very close to the one explored by the BO dynamics.
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res1
111
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.
Type: Goal |
Advantage: None |
Novelty: None |
ConceptID: Goa3
112
For this purpose, δt varies from 0.05 fs to 0.25 fs with a total simulation time of 100 fs.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met4
113
The maximum deviation of the total energy, compared to the initial total energy, and its average are reported in Table 1.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs3
114
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.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs3
115
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.
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res2
116
An increase of the average error of the total energy is also observed with a comparable variation.
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res3
117
However, the range of CP dynamics’ error still remains in very good limits (≤10−4 au) and is very close to BO dynamics’.
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res3
118
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.
Type: Conclusion |
Advantage: None |
Novelty: None |
ConceptID: Con3
119
The CPU time ratio between the BO and CP dynamics for different time steps is reported in Table 1.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs4
120
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.
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res4
121
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
Type: Method |
Advantage: No |
Novelty: Old |
ConceptID: Met2
122
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.
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res5
123
Then, one should used a time step of 0.5 fs for BO dynamics in order to obtain a similar CPU cost.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac24
124
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.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac24
125
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
126
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met4
127
Contrary to the CP dynamics, the restricted BO method shows a non-conservation of the total energy, as one can expect.
Type: Method |
Advantage: None |
Novelty: Old |
ConceptID: Met5
128
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.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod25
129
This guess vector corresponds to the variational solution adapted to the molecular structure expressed at time (t).
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod25
130
Consequently it is not fully adapted to the one expressed at time (t + δt).
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod25
131
In the CP approach, the molecular orbitals are propagated together with the nuclei.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod26
132
One can then consider that the molecular orbitals guess vectors are closer to the variational solution than in the BO approach.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac25
133
Consequently, a small number of iterations are sufficient to reach a quasi-convergence of the electronic energy and the nuclear gradients are accurately calculated.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac26
The monothiooxalic acid
134
The second molecule used to test our algorithm is the monothiooxalic acid for which simple and double intra molecular proton transfer reaction can occur.
Type: Experiment |
Advantage: None |
Novelty: None |
ConceptID: Exp3
135
This system has been chosen in order to test the ability of our program to reproduce more complicated dynamical behaviour (especially bond breaking).
Type: Goal |
Advantage: None |
Novelty: None |
ConceptID: Goa4
136
This system has already led to theoretical investigations:25 the PES has been explored and stationary points have been mostly localised.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac27
137
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.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac28
138
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.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac29
139
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
140
More precise dynamical investigations of this molecular system will be considered in future works.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac30
141
All electronic structure calculations have been performed with the hybrid DFT method, using the B3LYP exchange correlation functional.28,29
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac31
142
Core electrons of carbon, oxygen and sulfur were replaced by pseudo-potentials with their associated double-ζ basis set.24
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met4
143
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met4
144
The fictitious mass parameter was = 0.1 amu bohr2 and a 6-311++G(d,p) basis set was employed for hydrogen.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod27
145
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).
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met4
146
This configuration is less stable by 0.1 au than the most stable conformation.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs5
147
The initial total kinetic energy was 5.0 × 10−3 a.u.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs6
148
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.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac32
149
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.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs7
150
In this more complicated system the degree of accuracy is lower than in the case of a water molecule.
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res6
151
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.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs8
152
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.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs9
153
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.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs10
154
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.
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res7
155
Moreover, as shown in Fig. 4, the difference of potential energy between BO and CP increases with time.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs11
156
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.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac33
157
As a consequence, explored areas of the potential energy surface are not strictly identical.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs12
158
The increase of difference between BO and CP kinetic nuclear energies is less pronounced than on electronic energies.
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res8
159
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 Å.
Type: Goal |
Advantage: None |
Novelty: None |
ConceptID: Goa5
160
This fact confirms that the explored path by the CP dynamics is close to the BO one, but not strictly the same.
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res9
161
As a consequence, differences between BO and CP, in terms of energies and nuclear structures may increase with time (see Figs. 4–6).
Type: Conclusion |
Advantage: None |
Novelty: None |
ConceptID: Con4
162
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.
Type: Goal |
Advantage: None |
Novelty: None |
ConceptID: Goa6
163
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.
Type: Goal |
Advantage: None |
Novelty: None |
ConceptID: Goa6
164
Results are compared with BO trajectories with the same propagation characteristics, except the integrator.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met4
165
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.
Type: Experiment |
Advantage: None |
Novelty: None |
ConceptID: Exp4
166
Differences between CP and BO dynamics for different basis sets are shown in Tables 2 and 3.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs13
167
Each component of the discrepancies on the total energy are in the same order of magnitude.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs13
168
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).
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs13
169
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.
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res10
170
With respect to the BO simulations, maximum deviations are three times larger, but they still remain in acceptable limits.
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res10
171
This demonstrates that the CP scheme suits different curvatures of PES with a significant degree of accuracy.
Type: Conclusion |
Advantage: None |
Novelty: None |
ConceptID: Con5
172
As stated in the previous paragraph, the total energy is well conserved in the CP approach.
Type: Conclusion |
Advantage: None |
Novelty: None |
ConceptID: Con6
173
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.
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res11
A (H2O)20 water cluster
174
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).
Type: Experiment |
Advantage: None |
Novelty: None |
ConceptID: Exp5
175
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.
Type: Model |
Advantage: None |
Novelty: None |
ConceptID: Mod28
176
The initial total kinetic energy was 0.17 au.
Type: Observation |
Advantage: None |
Novelty: None |
ConceptID: Obs14
177
The simulation has been run for 500 fs using a 0.25 fs time step (2000 steps).
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met4
178
This trajectory can be regarded as a significant test case to evaluate the performance of the algorithm.
Type: Goal |
Advantage: None |
Novelty: None |
ConceptID: Goa7
179
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.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac34
180
In this example, one should notice that the comparison with BO trajectory is more expensive, the BO simulation leading to prohibitive computational time.
Type: Method |
Advantage: No |
Novelty: Old |
ConceptID: Met1
181
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).
Type: Result |
Advantage: None |
Novelty: None |
ConceptID: Res12
182
This final computational test makes us confident on the robustness and the applicability of our algorithm.
Type: Conclusion |
Advantage: None |
Novelty: None |
ConceptID: Con7
Conclusion
183
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
184
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.
Type: Method |
Advantage: None |
Novelty: New |
ConceptID: Met3
185
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.
Type: Conclusion |
Advantage: None |
Novelty: None |
ConceptID: Con7
186
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.
Type: Conclusion |
Advantage: None |
Novelty: None |
ConceptID: Con2
187
This glimpses dynamical studies of bigger chemical systems within an acceptable computational time.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac35
188
With more sophisticated tools, the propagation scheme may be more efficient, as it has already implemented in all other molecular dynamics suites of programs.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac36
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac36
190
These aspects are currently under development.
Type: Background |
Advantage: None |
Novelty: None |
ConceptID: Bac37