Computational Physics

9.4 Orthonormalisation; conjugate gradient and RM-DIIS techniques 285

descent direction:
ζk=− 2

[

Hψk−

∑

l

〈ψk|H|ψl〉ψk

]

. (9.75)

In a method proposed by Stichet al.[14], this approximation is used in the line
minimisations. Gillan[15]has devised a conjugate gradients method in which this
approximation is not made. Finally, Teteret al.[16]have proposed a particularly
efficient method in which only one stateψkis updated at a time. See Problems 9.3
and 9.5 for examples.
The conjugate gradients method is known to converge slowly for this problem.
The reason is that the steepest descent direction, which should be close to the
difference between the current orbitals{ψk}and those that minimise the energy,
{ψ ̃k}, may not have this property. To see what goes wrong, let us expandδψk=
ψk−ψ ̃kin the eigenstatesξpof the Hamiltonian:

δψk=

∑

p

αpkξp. (9.76)

Note that whereas the indiceskandlin the previous discussion run only over the
occupied electronic states, the indexpruns overallthe eigenstates of the variational
Hamiltonian matrix, and their number is equal to the number of states in the vari-
ational basis set used. Working out the steepest descent direction for this state, we
have
ζk=− 2

∑

p

(αpkεp−

∑

l

(^) klαpk)ξp (9.77)
(note that the part involvingψ ̃kon the right hand side vanishes). The right hand
sign will contain important contributions from the high energy levels, and they will
spoil the desired proportionality between steepest descent direction andδψk. This
can be remedied by an extension to the method, calledpre-conditioning[17]. We
shallnottreatthisindetailhere–fordetailsseeRefs.[ 8 , 16 , 17 ].
The conjugate gradients method can be readily applied to energy minimisation
problems. However, if we want to perform a dynamic simulation for the atoms, that
is, without friction acting on the nuclei, a rather subtle problem arises. To see this,
consider a nucleus which is moving in the positivex-direction with an electronic
charge distribution around it, which has converged to the ground state. As the
nucleus moves, the charge cloud starts lagging behind and this discrepancy might
grow larger and larger with time. But this is not possible in the Verlet algorithm as
this symplectic algorithm does not allow for energy drift (seeSection 8.4). In the
Verlet simulation, the electron charge cloud keeps oscillating around the nucleus so
that errors remain bounded, and even cancel out on average. This is an important
advantage of the molecular dynamics approach of Car and Parrinello [ 18 – 20 ]. There
is, however, the possibility that the electron system absorbs energy from the nuclear