So far, total energy calculations have been described for only fixed set of atomic coordinates. A natural consequence of the Car-Parrinello formalism is that electronic and ionic relaxation can take place simultaneously. For this to be performed, the total forces on the ions must be found in order to update the ionic positions.
The force on an ion, I, at position can be obtained from the full derivative of the total energy, E, by
As the electronic wavefunctions change, the force on the ions will be altered, therefore the full derivative has to be expressed in terms of changes in the wavefunction:
it follows that is just . But are electronic eigenstates with eigenvalues . Therefore after some algebraic manipulation the final two terms can be written
which is zero since is the normalisation constant. This shows that when each electronic orbital is an eigenstate of the Hamiltonian, then the partial derivative of the total energy with respect to the ionic positions is the force felt by the ions. This is the Hellmann-Feynman theorem[49, 50] and can be generalised to any order of derivative of the total energy.
In practice, the wavefunctions are only calculated to a given tolerance therefore they are never exact eigenstates of the Kohn-Sham Hamiltonian. As a result the forces calculated using the Hellmann-Feynman theorem incur error. This can be seen when the expression for the force is written formally as
where are all external energies such as the Ewald energy. The final term, called the variational force, vanishes when the wavefunctions are completely converged so that the conditions of the Hellmann-Feynman theorem are satisfied. It is in this term that the errors occurring from unconverged wavefunctions or incomplete basis sets are incurred. The error on the force is first order with respect to error in the wavefunction (as compared to errors in the total energy which are second order). It follows that the Hellmann-Feynman theorem can only be implemented when the wave functions are very close to self consistency. Only then can the ionic equations of motion be integrated and the ionic positions updated.
This leads to the use of the Born-Oppenheimer approximation. This is based on the fact that typical electronic velocities are much greater than that of the ions. It can therefore be assumed that the ions move so slowly relative to the electrons that at any point in time, the electrons will be in their ground state for that instantaneous ionic configuration. A schematic view of such an ionic trajectory in given in Figure 3.5.
Figure 3.5: A schematic view of the potential surface of . A typical ionic relaxation trajectory is shown by the solid line which must lie close to the Born-Oppenheimer surface (dashed line) for the Hellmann-Feynman theorem to be implemented without incurring large errors.
The need for the system be be kept close to the Born-Oppenheimer surface governs the method in which Car-Parrinello molecular dynamics is used in practice. The ground state wavefunction is first calculated self-consistently for a fixed set of ionic coordinates. In this method an initial wavefunction has to be assumed - this is given by initializing the expansion coefficients of the wavefunction by random numbers. Although this initial state is far from the Kohn-Sham eigenstates it does not assume any initial symmetries which, if incorrect, could lead to prohibitively long calculation to find the correct wavefunctions. The total energy is reduced by the conjugate gradients method (ensuring the wavefunctions are orthogonalised at the end of each iteration) until the wavefunctions are sufficiently converged to give correct forces. New ionic positions are then calculated under the influence of these forces either by integration of the equations of motion or by direct minimisation techniques. After each change in the position of the ions, the wavefunctions will no longer be eigenstates of the new ionic positions. In order to find the wavefunctions for this better set of ionic positions it is not necessary to randomise the wavefunctions. The motion of the ions will be sufficiently small that the old wavefunctions will be a good place to initialise the iterative procedure again.
This process of calculating the electronic wavefunctions, using the Hellmann - Feynman theorem to calculate the forces on the ions, moving the ions under the influence of these forces and updating the wavefunction will be a will keep the trajectory of the molecular dynamics simulation on the Born-Oppenheimer surface if a sufficiently small ionic timestep is used in conjunction with well converged wavefunctions.
A summary of the computational procedure is shown in Figure 3.6.
Figure 3.6: This flow chart illustrates the procedure used in ab initio molecular dynamics simulations.
This process forms the method which is used in later chapters to find a 0K local minimum of the ionic positions and the ground state electronic structure of materials. But this method also lends itself to finite temperature ab initio molecular dynamics simulations. Instead of moving the ions in order to reduce the forces acting (and therefore reduce the total energy of the system) the ions are initially given a random velocity, with a total kinetic energy given by per degree of freedom. Temperature is a time averaged quantity as energy interchanges between potential and kinetic forms, therefore it is not physically correct to scale the velocities every timestep to remain at a fixed temperature. The time averaged temperature is maintained by means of a Nosé thermostat. Integration of the ionic equations of motion continue at that temperature/kinetic energy governed by the thermostat allowing the ions to move under the influence of the Hellmann-Feynman forces. The system is returned as close as possible to the Born-Oppenheimer surface before the ions are moved at each timestep.