280 Quantum molecular dynamics
the smooth evolution of the orbitals! However, if the rotation is always close to the
unit transformation,
U= 1 +h^2 A (9.53)whereAis a Hermitian transformation of order one, varying smoothly with time,
the Verlet algorithm will still work: apart from the motion governed by the equation
of motion, the algorithm might induce some extra forces which cause the orbitals to
rotate smoothly in Hilbert space, and this latter motion can be dealt with perfectly
by the Verlet algorithm. It is difficult to see whether orthogonalisation algorithms
satisfy these requirements and it is therefore easiest to construct the algorithm such
that it is equivalent to the unambiguous time evolution resulting from the Verlet
algorithm (without extra rotation) to a precision of at least orderh^4 , which is the
overall precision of the Verlet algorithm.
A method which is based on the Verlet algorithm and which solves the (^) klin
(9.50)by the orthogonality requirements is the iterative algorithm called SHAKE
by Ryckaertet al.[5], which was mentioned inSection 8.6.3. This algorithm was
used in the original work of Car and Parrinello [11]. It is straightforward and does
not introduce rotations of the set of orbitals. Moreover, it orthogonalises the states to
arbitrary precision (depending on the number of iterations performed). For details
we refer to the cited literature.
Most other methods first predict the form of the (orthonormal)ψkat the next time
step with some precision and then perform an additional orthonormalisation of these
predicted orbitals by constructing orthonormal linear combinations of them. The
idea behind this is that if the prediction is accurate, only a few orthonormalisation
iterations are needed. As the Verlet algorithm prescribes an orthonormalisation by
mixing in theψk(t)through the Lagrange multipliers (seeEq. (9.50)), and not the
ψk(t+h), such a final re-orthonormalisation can only be justified if the changes
involved are of orderh^4 which is the overall accuracy of the Verlet algorithm for
a single step. Therefore these algorithms must first predict the new values to order
O(h^4 )and the re-orthonormalisation should yield the new states lying close to the
old states. Note that after each step orthonormality is of then satisfied to machine
precision whereas the error in the integration algorithm is of orderh^4.
Let us now consider one such algorithm in detail. Over a time steph, the orbitals
ψkshift by an amount of orderhowing to the fact that they have a velocity. The
acceleration, which is caused by the force due to the derivative of the total energy
and by the force due to the Lagrange multiplier terms, then gives an additional shift
of orderh^2. The Hamiltonian force term occurring in the Car–Parrinello equations
of motion is given by−Hψk. First we neglect this force term for simplicity and
find Lagrange multipliers which guarantee orthonormality in the absence of forces.