The Nuts and Bolts of First-Principles Simulation

1. The Nuts and Bolts of First-Principles Simulation Lecture 14: Forces and Stresses Durham, 6th-13th December 2001 CASTEP Developers’ Groupwith support from the ESF k Network

2. Overview of Lecture • Why bother? • Theoretical background • CASTEP details • Symmetry and User Constraints • Conclusion Lecture 14: Forces and Stresses

3. Why bother? (I) • Structure optimisation • Minimum energy corresponds to zero force • Much more efficient than just using energy alone • Equilibrium bond lengths, angles, etc. • Minimum enthalpy corresponds to zero force and stress • Can therefore minimise enthalpy w.r.t. supercell shape due to internal stress and external pressure • Pressure-driven phase transitions Lecture 14: Forces and Stresses

4. Why bother? (II) • Molecular dynamics • Can do classical dynamics of ions using forces derived from abinitio electronic structure • Copes with unusual geometry, bond-breaking, chemical reactions, catalysis, diffusion, etc • Incorporates effects of finite temperature of ions • Can generate thermodynamic information from ensemble averaging • Time dependent phenomena • Temperature driven phase transitions Lecture 14: Forces and Stresses

5. Theoretical Background • Hellman-Feynman Theorem • basic Quantum Mechanics • Density Functional Theory • how it applies in DFT Lecture 14: Forces and Stresses

6. Hellman-Feynman Theorem (I) • Classically we have the force F at position R is determined from the potential energy as • Quantum mechanically we therefore expect where Lecture 14: Forces and Stresses

7. Hellman-Feynman Theorem (II) • If we write the three unit cell vectors a, b, c as the columns of a matrix h then the effect of an applied strain is to change the shape of the unit cell: • We then have the stress tensor s related to the strain tensor e by: where is the volume of the unit cell. Lecture 14: Forces and Stresses

8. c c b b a a b b g g aa+da a Stress and strain in action sxx c b a NB Much messier if non-orthogonal cell sxy b g g+dg a Lecture 14: Forces and Stresses

9. Hellman-Feynman Theorem (III) • The Hellman-Feynman Theorem states that for any perturbation l we have which obviously includes the case we are interested in. We have assumed that the wavefunction is properly normalised and is an exact eigenstate of H. Lecture 14: Forces and Stresses

10. Hellman-Feynman Theorem (IV) • To evaluate <E> for an unknown wavefunction Y we first expand it in terms of a complete set of fixed basis functions j and then use the Variational Principle to find the set of complex coefficients ci that minimise the energy. • If the basis set is incomplete then we arrive at an upper-bound for the energy. Lecture 14: Forces and Stresses

11. Hellman-Feynman Theorem (V) • If we have an approximate eigenstate Y, for example from using an incomplete basis set, then we must keep all 3 terms in the general expression. • If our basis set depends upon the ionic positions, such as atomic centred Gaussians, then the other derivatives in the general expression will contribute so-called Pulay forces (stresses). • Note that Pulay forces (stresses) will vanish in the limit of a complete basis set, but that this is never realized in practice, or if position independent basis functions, such as plane-waves, are used. Lecture 14: Forces and Stresses

12. Hellman-Feynman Theorem (VI) • If we choose plane-waves as our basis functions, then because these functions are independent of the ionic coordinates, it can easily seen that the general expression for the forces becomes: Lecture 14: Forces and Stresses

13. Hellman-Feynman Theorem (VII) • That is, we can calculate the forces using the same expansion coefficients as we used to variationally minimise the energy, using matrix elements of the ionic derivative of the Hamiltonian. • This makes calculation of the forces relatively cheap once the variational energy minimisation has been completed if we are using a plane-wave basis set. • Similar expressions can also be derived for the stresses. Lecture 14: Forces and Stresses

14. Density Functional Theory (I) • In DFT we have the Kohn-Sham Hamiltonian: • Therefore we only get contributions to the forces from the electron-ion (pseudo)potential and the ion-ion Coulomb interaction (the Ewald sum). • Also contribution from exchange-correlation potential if using non-linear core corrections. • However, for the stresses, we also get a contribution from the kinetic energy and Hartree terms. Lecture 14: Forces and Stresses

15. Density Functional Theory (II) • As we do not have a complete basis, the wavefunction will not be exact even within DFT. • If we use a variational method to minimize the total energy, then we know that the energy and hence the wavefunction will be correct to second order errors. • However, the forces will only be correct to first order  need a larger basis set for accurate forces than for energies. Lecture 14: Forces and Stresses

16. Lecture 14: Forces and Stresses

17. Density Functional Theory (III) • However, if we use a non-variational minimisation technique, such as density mixing, then such statements cannot be made. • We can no longer guarantee that the energy found is an upper-bound on the true ground state energy. • This complicates the application of the Hellman-Feynman theorem. • Consequently, non-variational forces and stresses are less reliable. Lecture 14: Forces and Stresses

18. CASTEP details (I) • Forces and stresses are almost the highest level functionality • Single subroutine call to first derivatives module is all that is required to return the forces for the current model if ground state is already known. Ditto stress. • firstd puts together the different contributions to the force (stress) using other functional modules so the physics is obvious. • These in turn call down to operations on charge densities and potentials. Even at this level the physics is obvious and the details of USPs etc are hidden. Lecture 14: Forces and Stresses

19. CASTEP details (II) • The use of Ultra-Soft Pseudopotentials further complicates things, as there are now additional contributions to both the forces and the stresses from the charge augmentation. • However, the modular design of new CASTEP completely hides this from the higher level programmer. Lecture 14: Forces and Stresses

20. CASTEP details (III) • There is a problem with the application of the Hellman-Feynman theorem with non-variational minimisers • Consequently the CASTEP code contains a first-order correction to the forces derived from density mixing. • However, the corresponding correction to the stresses is not known. • This has implications for structure optimisation and molecular dynamics with density mixing • Therefore the more recent Ensemble DFT approach (which is fully variational) is to be preferred. Lecture 14: Forces and Stresses

21. Lecture 14: Forces and Stresses

22. CASTEP details (IV) • If the unit cell changes shape then the number of plane-waves and the FFT grid will, at some point, change discontinuously. • Consequently, it becomes difficult to compare results at different cell sizes at the same nominal basis set size (cut-off energy) as the effective quality of the basis set is not the same, unless the basis set is fully converged (impossible). • This can be countered by using the Finite Basis Set Correction, which calculates the change in total energy upon changing the basis set size at a fixed cell size, and then uses this to correct the total energy and stress at nearby cell sizes. Lecture 14: Forces and Stresses

23. Lecture 14: Forces and Stresses

24. Symmetry • If the symmetry of the system has been calculated (controlled by keyword in input file) then this can be used to symmetrise the forces and stresses. • This ensures that the forces (stresses) have the same symmetry as the model. • Consequently, the symmetry of the system will be preserved in any structural relaxation or dynamics. Lecture 14: Forces and Stresses

25. User Constraints (I) • Sometimes it is desired to impose additional constraints on any structural relaxation or dynamics. • Currently, CASTEP can apply any arbitrary number of linear constraints on the atomic coordinates, up to the number of degrees of freedom. • E.g. fixing an atom, constraining an atom to move in a line or plane, fixing the relative positions of pairs of atoms, fixing the centre of mass of the system, etc. Lecture 14: Forces and Stresses

26. User Constraints (II) • Both force and stress constraints work in the same way – by Lagrange multipliers in the extended Lagrangian of the system: • So for a given set of constraints S(q) we need to know the derivatives of the constraints and then if we can determine l we have the constraint force. Lecture 14: Forces and Stresses

27. User Constraints (III) • Linear constraints for the ionic motion makes the matrix of constraint derivatives trivial and therefore determining l: Lecture 14: Forces and Stresses

28. User Constraints (IV) • If we have multiple constraints, S(q), R(q), etc. then we may satisfy each constraint simultaneously if the constraint matrices are all mutually orthogonal and so we use Gram-Schmidt on the coefficients. • With non-linear constraints, both the constraints and their derivatives need to be specified, and it is not possible in general to satisfy them all simultaneously. Consequently, iterative procedures such as SHAKE or RATTLE are required. Lecture 14: Forces and Stresses

29. User Constraints (V) • Constraints can also be applied to the unit cell lengths and angles (giving 6 degrees of freedom). • Any length (angle) can be held constant, or tied to one or both of the others. • This is most useful if only a subset of the symmetries of the original unit cell is to be enforced. Lecture 14: Forces and Stresses

30. User Constraints (VI) • Cell constraints are implemented in a slightly different way. The stress tensor is transformed from its normal symmetric representation in Cartesian coordinates into the space of cell lengths and angles. • There is then a 1:1 correspondence between the stress components and the cell degrees of freedom and so the constraints may be trivially applied to the stress and hence to the evolution of the system. Lecture 14: Forces and Stresses

31. Conclusion • Hellman-Feynman Theorem gives us a simple recipe for calculating ab initio forces and stresses • plane-wave basis has big advantage • DFT implementation • Can be combined with symmetry and/or constraints • Major use in structural relaxation and molecular dynamics Lecture 14: Forces and Stresses