The Nuts and Bolts of First-Principles Simulation

1. The Nuts and Bolts of First-Principles Simulation Lecture 18: First Look at Molecular Dynamics Durham, 6th-13th December 2001 CASTEP Developers’ Groupwith support from the ESF k Network

2. Overview of Lecture • Why bother? • What can you it tell you? • How does it work? • Practical tips • Future directions • Conclusions Lecture 18: First look at MD

3. Why Bother? • Atoms move! • Time dependant phenomena • Ionic vibrations (phonons, IR spectra, etc) • Diffusion, transport, etc. • Temperature dependant phenomena • Equilibrium thermodynamic properties • Catalysis and reactions • Free energies • Temperature driven phase transitions, melting, etc Lecture 18: First look at MD

4. Radiation damage in zircon T=300 K T=600 K Lecture 18: First look at MD

5. Na+ diffusion in quartz Lecture 18: First look at MD

6. What Can It Tell You? • Ensemble Averages • Temperature, pressure, density, configuration energy, enthalpy, structural correlations, time correlations, elastic properties, etc. • Correlation Functions • Time dependent, e.g. velocity auto-correlation function Cvv(t) • Spatially dependent, e.g. radial distribution function g(r) • Fluctuations • Energy fluctuations  Cv, enthalpy fluctuations  Cp, etc. • Distribution Functions • E.g. velocity distribution function, energy distribution function, etc. Lecture 18: First look at MD

7. Velocity Auto-Correlation Function Cvv 1 Gas 0 t Solid Liquid Lecture 18: First look at MD

8. Radial Distribution Function g(r) Solid Liquid 1 Gas r/a0 0 1 2 3 4 Lecture 18: First look at MD

9. How Does It Work? • Classical dynamics of ions using ab initio forces derived from the electronic structure • Integrate classical equation of motion • Discretise time  time step • Different integration algorithms • Trade-off time step • long-term stability vs. short-time accuracy • Ergodic Hypothesis • MD trajectory samples phase space • time average = ensemble average Lecture 18: First look at MD

10. Integration Algorithms (I) • Euler • Simplest method but unstable to error propagation • Runge-Kutta • Excellent stability but too many force evaluations and not symplectic (time reversible) • Predictor-Corrector • Old CASTEP – not symplectic  unsuitable for MD • Verlet • Position Verlet – not explicit velocities so using thermostats is not straightforward • Velocity Verlet – current and new CASTEP Lecture 18: First look at MD

11. Integration Algorithms (II) • Multiple time / lengths scale algorithms • Recent theoretical development • Excellent results in special cases but hard to apply in general purpose code • Car-Parrinello • Combines electron and ion MD • Time step dominated by electrons not ions • Cannot handle metals • Iterative ab initio methods such as CASTEP require more effort to minimize the electrons but compensate by taking larger time steps based upon ions – even better with constraints … Lecture 18: First look at MD

12. Time Step • Should reflect physics not algorithm • e.g. smallest phonon period/10 • Effects the conservation properties of system and long-time stability • Typically ~ femto-sec for ab initio calculations • Limitation on time scale of observations – total run-length ~pico-sec routine, nano-sec exceptional • Use of constraints to increase time step • Freeze motions that are not of interest Lecture 18: First look at MD

13. Types of MD • Micro-canonical = constant NVE • Simplest MD - purely Newtonian dynamics • Canonical = constant NVT • Closer to experiment but need to add a thermostat • Isobaric-Isothermal = constant NPT • Closest to experimental conditions but need to add a barostat as well • Grand Canonical = constant mVT • Cannot do with ab initio MD but has been used with MC Lecture 18: First look at MD

14. Micro-Canonical Ensemble • Suitable for investigating time dependent phenomena • E.g. simple way to sample a single normal modes/ vibrational frequency of complex systems • Set temperature=0 • Tweak relevant bond • Watch the system evolve Lecture 18: First look at MD

15. Lecture 18: First look at MD

16. Canonical Ensemble • Ensemble of choice for investigating finite temperature phenomena • E.g. diffusion • Set appropriate temperature, let system evolve and monitor MSD • E.g. vibrational spectra • Set appropriate temperature, let system evolve and calculate Fourier transform of velocity auto-correlation function Lecture 18: First look at MD

17. Choice of Thermostat • Velocity rescaling • Simple, but breaks the smooth evolution of the system and without theoretical foundation • Not used in CASTEP • Nosé-Hoover • Couples system to external heat bath using an auxiliary variable • Deterministic evolution but not always ergodic • Langevin • Based on fluctuation-dissipation theorem and coupling to an external heat bath • Stochastic evolution but always ergodic Lecture 18: First look at MD

18. Nosé-Hoover Thermostat (I) • Extended Lagrangian and Hamiltonian • Modified equations of motion Lecture 18: First look at MD

19. Nosé-Hoover Thermostat (II) • Need to specify the thermostat ‘mass’ Q • Choose Q so as to cause thermostat-system coupling frequency to resonate with characteristic frequency of system – tricky! • New CASTEP – input coupling frequency instead and code then estimates appropriate Q Lecture 18: First look at MD

20. Lecture 18: First look at MD

21. Langevin Thermostat (I) • Modified equation of motion • Fluctuation • Has proper statistical properties, e.g. thermal fluctuations of system obey Lecture 18: First look at MD

22. Langevin Thermostat (II) • Time-scale of thermal fluctuations depends on the Langevin damping time tL • Need to choose s.t. tL is greater than the characteristic period tc of the system s.t. short-time dynamics is accurately reproduced • “Rule of 10s” • Choose time step s.t. tcdt*10 • Choose Langevin damping time s.t. tL tc*10 • Choose run length s.t. truntL*10 Lecture 18: First look at MD

23. Lecture 18: First look at MD

24. Influence of Electronic Minimizer • All-Bands • Variational minimization  accurate forces • Problems with metals • Density Mixing • Non-variational minimization  need higher accuracy Y to get same accuracy forces and need to correct forces  less accurate MD • OK with metals • Ensemble DFT • Variational minimization  accurate forces • Great with metals Lecture 18: First look at MD

25. Wavefunction Extrapolation (I) • Advantages • Generate better guess for Yat new ionic configuration • Less work for electronic minimizer  faster • Assumes can extrapolation Y forwards in time in similar manner to ionic positions • Can either do first or second order extrapolation • Can either used fixed values for (a,b) or those which minimize difference between MD and extrapolated coordinates Lecture 18: First look at MD

26. Wavefunction Extrapolation (II) x Y- Y0 t Y+ Lecture 18: First look at MD

27. Wavefunction Extrapolation (III) • All bands / Ensemble DFT • Extrapolate Y only • Density-Mixing • Must extrapolate Y and  independently else Residual = 0 and not ground state! • Decompose  into atomic and non-atomic contributions • Move the atomic charges onto the new ionic coordinate • Extrapolate the non-atomic part only Lecture 18: First look at MD

28. Lecture 18: First look at MD

29. Practical Tips (I) • Equilibration • Sensitivity of system to initial conditions • Depends on quantity of interest • E.g. if after equilibrium average, then must allow system to evolve to equilibrium before start data collection • Auto-correlation functions give useful information on the “memory” of the system to the quantity of interest. Lecture 18: First look at MD

30. production equilibration Lecture 18: First look at MD

31. Practical Tips (II) • Sampling • It is very easy to over-sample the data and consequently under-estimate the variance • Successive configurations are highly correlated - not independent data points • Need to determine optimal sampling frequency of the quantity of interest ‘A’ and either save data at appropriate intervals or adjust error bars  need to analyse variance in blocks of size tb Lecture 18: First look at MD

32. Practical Tips (III) Lecture 18: First look at MD

33. Future Directions • Isothermal-Isobaric Ensemble • Variable cell MD • Allow cell size and shape to evolve under internal stress and external pressure • Closest to experimental conditions • Important for generic phase transitions Lecture 18: First look at MD

34. Conclusions • Two ensembles can be simulated in CASTEP using Velocity Verlet integration • Large time step, excellent long term energy conservation and stability • Micro-Canonical (NVE) • Canonical (NVT) • Nosé-Hoover thermostat – deterministic • Langevin thermostat – stochastic • Can be used to study many phenomena - see later talks for example applications! Lecture 18: First look at MD