The Nuts and Bolts of First-Principles Simulation

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

2. Overview of Lecture • Why bother? • What can it tell you? • How does it work? • Damped MD in CASTEP • BFGS in CASTEP • Future directions • Conclusions Lecture 15: Structural Calculations and Pressure

3. Why bother? • Want to find ground state of system • Need to minimise energy of electronic structure at fixed ionic positions and then optimise the ionic positions and/or the unit cell shape and size (particularly if external pressure applied) • Theoretical minimum depends on choice of pseudopotential, plane-wave cut-off energy, choice of XC functional, etc. • Often not exactly the same as experiment! Fully converged calculation should get agreement to better than 1% Lecture 15: Structural Calculations and Pressure

4. What Can It Tell You? • Equilibrium bond lengths and angles • Equilibrium cell parameters • Discriminating between competing structures • Elastic constants • Surface reconstructions • Pressure-driven phase transitions • Starting point for many advanced investigations … Lecture 15: Structural Calculations and Pressure

5. How Does It Work? • Electrons adjust instantly to position of ions so have a multi-dimensional potential energy surface and want to find global minimum • Treat as an optimisation problem • Simplest approach is steepest descents • More physical approach is damped MD • More sophisticated approaches are conjugate gradients or BFGS • Could use simulated annealing to avoid getting stuck in local minima Lecture 15: Structural Calculations and Pressure

6. Steepest Descents (I) • Simplest minimiser • Step downhill in direction of local steepest gradient using trial step length • Line minimisation to find the optimal step length • Repeat to convergence Lecture 15: Structural Calculations and Pressure

7. Steepest Descents (II) Traversing a long, narrow valley Enlargement of a single step showing the line minimisation in action – the step continues until the local energy starts to rise again whereupon a new direction is selected which is orthogonal to the previous one Lecture 15: Structural Calculations and Pressure

8. Steepest Descents (III) • Advantages • Easy to implement • Very robust – hard to get it confused! • Reliable – will always get to the minima eventually • Disadvantages • Often very slow to converge • Can get stuck in local minima Lecture 15: Structural Calculations and Pressure

9. Damped MD (I) • Improved Steepest Descents • Move ions using velocities as well as forces  second-order equation of motion  more efficient than steepest descents • Need to add damping term ‘-gv’ to forces • Initially v=0 • Use physical insight to get ‘optimal damping’ factor and to adjust the time step s.t. get rapid convergence Lecture 15: Structural Calculations and Pressure

10. Damped MD (II) x Over-damped Under-damped Critically damped t Lecture 15: Structural Calculations and Pressure

11. Damped MD (III) • Advantages • Easy to code, robust, much more efficient than either steepest descents or simulated annealing • Can do with Car-Parrinello or conjugate gradient minimisation of electrons. • If not CP then can accelerate the ab initio part of the calculation using wavefunction extrapolation • Disadvantages • Convergence rate depends on damping factor • Can get stuck in local minima Lecture 15: Structural Calculations and Pressure

12. Conjugate Gradients (I) • Improved Steepest Descents • Moves ions according to the gradient with a line minimisation to find the best step length • But not just the locally steepest gradient – constructs a direction which is conjugate to all previous directions so it does not undo earlier minimizations at later times  big improvement in rate of convergence • See previous talks on electronic minimisation strategies for details Lecture 15: Structural Calculations and Pressure

13. Conjugate Gradients (II) Traversing a long, narrow valley The initial search direction is given by steepest descents. Subsequent search directions are constructed to be orthogonal to the previous and all prior search directions. Lecture 15: Structural Calculations and Pressure

14. Conjugate Gradients (III) • Advantages • Rapid rate of convergence – in a quadratic energy landscape, each iteration should converge one degree of freedom • Low storage requirements • Disadvantages • More complex to code and cannot use with CP • No knowledge of the Hessian explicitly generated • Can get stuck in local minima Lecture 15: Structural Calculations and Pressure

15. BFGS (I) • Basic idea: • Energy surface around a minima is quadratic in small displacements and so is totally determined by the Hessian matrix A (the matrix of second derivatives of the energy): • so if we knew A then could move from any nearby point to the minimum in 1 step! Lecture 15: Structural Calculations and Pressure

16. BFGS (II) • The Problem • We do not know Aa priori • Therefore we build up a progressively improving approximation to A (actually H=A-1) as we move the ions according to the BFGS (Broyden-Fletcher-Goldfarb-Shanno) algorithm. • Also known as a quasi-Newton or variable metric method. • Positions updated according to: Lecture 15: Structural Calculations and Pressure

17. BFGS (III) • In a perfectly quadratic system then l=1 is the optimal step length. In a real system need a line minimization to find the optimal l: F.DX start best l 0 1 trial Lecture 15: Structural Calculations and Pressure

18. BFGS (IV) • Davidon-Fletcher-Powell (DFP) variant • Older method • Mathematically equivalent to BFGS • Less tolerant of round-off error or inexact line minimisation • Direct Hessian Method • Calculates A rather than H • Uses Cholesky decomposition to keep similar speed but should guarantee that A remains non-singular. Very rare problem in practice. Lecture 15: Structural Calculations and Pressure

19. BFGS (V) • DFP: Hessian updated as • BFGS: Hessian updated as Lecture 15: Structural Calculations and Pressure

20. BFGS (VI) • Advantages • Convergence rate similar (or better) than CG • Extra physical information generated from Hessian • Disadvantages • Complex to code and cannot use with CP • Storage of Hessian ~ (number d.o.f.)^2 which is prohibitive for electronic problem but OK for ionic • Can get stuck in local minima Lecture 15: Structural Calculations and Pressure

21. Simulated Annealing (I) • Stochastic method • Monte-Carlo style exploration of energy landscape • Always accept steps that go downhill in energy • Sometimes accept uphill steps depending on Boltzman factor combining energy change with current temperature of system • Progressively reduce temperature and iterate to ground state Lecture 15: Structural Calculations and Pressure

22. Simulated Annealing (II) U(x) start stop x Lecture 15: Structural Calculations and Pressure

23. Simulated Annealing (III) • Advantages • Very robust and reliable • Reasonably immune to local minima • Disadvantages • Exceedingly slow to converge • Need to be careful in cooling rate so as not to quench and trap in local minimum if too fast, and not to waste too much time if too slow. • Cannot guarantee will find exact global minimum Lecture 15: Structural Calculations and Pressure

24. Lecture 15: Structural Calculations and Pressure

25. Damped MD in CASTEP (I) • Need to calculate damping g that corresponds to critical damping • For simple harmonic oscillator, this can be found from the characteristic frequency w0 • For a system with a dominant natural mode then this can be found and used to set g • As convergence is approached, the mode with largest w will be converged first whereupon g can be re-evaluated to accelerate the convergence of the slower modes • In the same way, the time step of the MD can be linked to the dominant mode w and therefore progressively increased as convergence is approached. Lecture 15: Structural Calculations and Pressure

26. Damped MD in CASTEP (II) • CASTEP has several different algorithms for calculating the optimal g and time step • Whilst it is to be expected that BFGS will be superior, experience has shown that the old CASTEP BFGS can sometimes be beaten • Although DMD will require more iterations than BFGS, each one will be cheaper because of wavefunction extrapolation • Sometimes it is more stable and succeeds in finding a minimum where old CASTEP BFGS fails • A useful fall-back but not so necessary with new CASTEP and the improved BFGS algorithm? Need more experience to tell! Lecture 15: Structural Calculations and Pressure

27. Lecture 15: Structural Calculations and Pressure

28. BFGS in old CASTEP (I) • Simultaneous optimisation of ions and positions  seeks to minimise the enthalpy H=E+pW • Works in space of cell parameters (a,b,c,a,b,g) and ionic positions with optional external pressure • Convergence criteria • Simultaneous convergence in enthalpy, RMS force, RMS stress component, RMS displacement of ions but has ‘escape route’ if everything except stress is converged Lecture 15: Structural Calculations and Pressure

29. BFGS in old CASTEP (II) • Problems • Uses DFP not BFGS update • Uses heuristic for step length rather than rigorous line minimisation • Starts with initial H = I (unit matrix) • Periodically resets H to I to reduce error accumulation whether strictly necessary or not • Does not guarantee that a step will go down in enthalpy Lecture 15: Structural Calculations and Pressure

30. BFGS in new CASTEP (I) • Also seeks to minimise the enthalpy H=E+pW • Works in space of fractional ionic positions and strains with optional external pressure • Unified approach allows better strain/coordinate coupling • Does rigorous line minimisation to find optimal step length • Only recalculate ground-state at the new structure if sufficiently different to the trial structure • Also bisection search if the line minimisation fails • Also quadratic step if necessary • Builds up H using BFGS update • Only reset H if run out of search directions Lecture 15: Structural Calculations and Pressure

31. BFGS in new CASTEP (II) • Starts with more complex initial guess at H • Block diagonal using input/default guess at bulk modulus (B) for cell part and input / default guess at the characteristic frequency (w0) for the ionic part • Guaranteed to preserve symmetry if so desired • Analyse H to get updated estimates for B and w0 • Only allows ‘uphill’ steps if all else has failed (including resetting H) • Stringent convergence criteria • Simultaneous convergence in enthalpy, max. modulus of force on any ion, max. component of stress tensor, max. displacement of any ion • Over a ‘window’ of successive iterations Lecture 15: Structural Calculations and Pressure

32. Lecture 15: Structural Calculations and Pressure

33. Lecture 15: Structural Calculations and Pressure

34. Lecture 15: Structural Calculations and Pressure

35. Lecture 15: Structural Calculations and Pressure

36. Future Directions • Transition State Search • Find saddle point structures not just stable ones • Redundant Internal Coordinates • Transform to normal coordinates and then only optimise the non-redundant degrees of freedom • Recent theoretical breakthrough in application to crystals/extended systems • Need to see how it compares to the new BFGS • Variable cell damped MD Lecture 15: Structural Calculations and Pressure

37. Conclusions • Two independent algorithms have been implemented within new CASTEP for structural optimisation • Damped MD • Can be very efficient depending on g • Only works with fixed cell parameters (so far!) • BFGS • Can do either/or/both ionic positions and cell parameter optimisation • Can re-use calculated B and/or w0 from earlier runs to accelerate convergence in later runs Lecture 15: Structural Calculations and Pressure