Real-World Predictions from Ab Initio Molecular Dynamics Simulations - Multiscale Molecular Methods in Applied Chemistry - page 120

optimization, i.e., the wavefunction derivatives are neglected, energy and forces are
also no longer fully consistent. This would constitute, for energies sufficiently close
to the Born-Oppenheimer energy, no real problem if these errors did not depend on
the optimization method and initial guess wavefunction. Initial wavefunctions
selected in a non-time reversible way, e.g., by taking the last optimized wavefunction,
will transfer this property with the force error to the nuclear dynamics. The nuclear
dynamics is therefore no longer time reversible, despite a seemingly time reversible
integration algorithm, and very poor stability of integration results. To cure this
problem a time reversible propagation of the initial wavefunction has to be chosen.
A time reversible BO molecular dynamics scheme based on the propagation
of one-particle density matrices has been proposed by Niklasson [ 8 ]. The equation
of motion for the density matrix is

P

2

¼ o ð D P Þ: (30)

Here, D is the self-consistent, optimized density matrix. Using a time-reversible
Verlet scheme we get an explicit integration of the form

t 2

2

P ð t þ d t Þ¼ 2 P ð t Þ P ð t d t Þþd o ½ D ð t Þ P ð t Þ: (31)

t 2

2

k ¼ d o ¼

2 as the original form proposed, using stability analysis the
largest possible value of

With

k k

can be determined. A larger value of is desirable, as for

d t a larger k

a given time step corresponds to a stiffer harmonic potential keeping

the propagated density matrix P closer to the optimized density matrix D . The
largest value of k that is consistent with stability is the one that guarantees that the
distance D ( t ) P ( t ) does not diverge. For the Verlet family of integrators this

optimal value is in fact 2.

The propagation of wavefunctions expanded in atom centered basis functions
needs special care. It is best to use an extrapolated contra-covariant density matrix
PS as a projector on to the occupied subspace

k ¼

X

K 1

C p C T ð t þ d t Þ¼ B k C ð t m d t Þ ð t m d t Þ S ð t m d t Þ C ð t d t Þ; (32)

m ¼ 0

where S ( t ) is the overlap matrix and C ( t ) are the orbital expansion coefficients at
time t . An approximate time reversible predictor-corrector method proposed by
Kolafa [ 17 ], always stable predictor-corrector (ASPC), originally proposed for
classical polarizable force fields, can be used. The extrapolation coefficients for
this method are

2 K þ 2

K m

m

:

B k ¼ð 1 Þ ð m þ 1 Þ (33)

2 K
K