WIXLIB

AA BriefBrief IntroductionIntroduction toto AbAb InitioInitioMolecularMolecular DynamicsDynamicsMatt ProbertCondensed Matter Dynamics GroupDepartment of Physics,University of York, rgCMD GroupDepartment of Physics

Overview of Talk In this talk I hope to give you someideas as to why you might want to doMD and what it can tell you. I hope to pass on some practical tipsand advice, and answer some of yourquestions, particular w.r.t. CASTEP I shall illustrate with examples wherepossible. Time is short CMD GroupDepartment of Physics

Why MD? Atoms move!– We may be interested in studying time dependentphenomena, such as molecular vibrations, phonons,diffusion, etc.– We may be interested in studying temperaturedependant phenomena, such as free energies,anharmonic effects, etc. Ergodic Hypothesis– One of the key principles behind the usefulness ofMD for statistical mechanics studies– Iff our MD trajectory is “good enough” then a timeaverage over the trajectory is equivalent to anensemble average – hence MD averages are useful.CMD GroupDepartment of Physics

Alternatives Monte Carlo– can do thermal averages– hard to do time dependant things Hybrid MD/MC– bad MD as good MC– generate configurations using poor/cheap/fastMD but then evaluate contribution toensemble average using MCCMD GroupDepartment of Physics

Types of ab initio MD Classical Motion– We use classical mechanics to move the atoms Born-Oppenheimer approximation decouples nucleus andelectrons– But using forces and stresses derived from theelectronic wavefunction– No quantum fluctuations, tunneling, zero pointmotion, etc. Quantum Motion– Can include ZPM etc using ab initio Path Integral MD Damped MD as a geometry optimizer– BFGS ought to be a lot better but not always – seeProbert, J. Comput. Phys. 191, 130 (2003)CMD GroupDepartment of Physics

Choice of Ensemble NVE––––Micro-canonical ensembleConstant Number of atoms, Volume and EnergyCorresponds to Newtonian mechanicsGood for non-equilibrium situations, e.g. watching abond vibrate or doing impact movies NVT– Canonical ensemble – constant Temperature– More physical as it allows energy exchange with aheat bath– Good for simulating thermal equilibrium– Choice of thermostating algorithmsCMD GroupDepartment of Physics

Choice of Ensemble NPH– constant pressure P and enthalpy H– Choice of barostats to handle pressure:– Andersen can allow cell to change size isotropically(liquids) whilst Parrinello-Rahman can allow changesin size and shape (solids)– External pressure can be isotropic (hydrostatic) oranisotropic (shear stress etc). NPT– Most physically relevant as system is now connectedto a piston and a heatbath.– Again, choice of thermostats and barostats µVT - constant chemical potential µCMD GroupDepartment of Physics

How do you do it? NVE Integrate classical equations of motion– discretize time time step– different integration algorithms, e.g. VelocityVerlet:f (t ) 23r (t δt ) r (t ) v(t ).δt ( ).δt O δt2mf (t ) f (t δt )v(t δt ) v(t ) .δt O δt 22m( )– trade-off time step vs. stability vs. accuracy– need accurate forces (high cutoff energiesand good k-point sampling)CMD GroupDepartment of Physics

Other Ensembles Other ensembles can be simulated by usingappropriate equations of motion– Usually derived from an extended Lagrangian (e.g.Nosé-Hoover, Parrinello-Rahman)– Recent developments in Liouvillian formulation havebeen very successful in deriving new symplecticintegration schemes Langevin schemes need to be deriveddifferently as non-Hamiltonian!– Need Focker-Planck & Liouville equation– see Quigley & Probert, J. Chem. Phys. 120, 11432(2004) or my last talk at the FHI Berlin!CMD GroupDepartment of Physics

Simple Example: N2 Naïve Materials Studio approach:– put 2 N atoms in a 5 A box at (0.4,0.5,0.5)and (0.6,0.5,0.5)– Use Gamma point for BZ sampling (it is anisolated molecule after all!)– Use default settings, e.g. “medium” Ecut.– Run NVT dynamics at default T 273 K usingLangevin thermostat with default “Langevintime” of 0.1 ps and default time step of 1.0 fs– What do you see?CMD GroupDepartment of Physics

Simple N2 MovieCMD GroupDepartment of Physics

Temperature ?200001800016000T me (ps)CMD GroupDepartment of Physics0.50.60.7

Constant of Motion ?-540-540.5Eham e (ps)CMD GroupDepartment of Physics0.50.60.7

What is Going On? Why is the temperature not constant if it issupposed to be NVT? The initial conditions were a long wayfrom equilibrium. Doing a simple fixed-cellgeometry optimisation relaxed 2 eV. This excess PE is turned into KE by theMD – hence the huge initial temperaturesbefore the thermostat is able to control it. The 2 eV excess PE shows up in thechange in “constant of motion”CMD GroupDepartment of Physics

What is this “constant of motion”? It certainly does not seem very constant!– It depends on the ensemble but is essentially theclosest thing to the “value of the Hamiltonian” whichshould be a conserved quantity:NVE : EHam Eelectrons KEionsNVT : EHam Eelectrons KEions PE NHC KE NHCNPH : EHam Eelectrons KEions pextV KEcellNPT : EHam Eelectrons KEions pextV KEcell PE NHC KE NHCMay fluctuate on short times but no long-term drift!CMD GroupDepartment of Physics

Ignoring initial T transient1000900800T (K)70060050040030020010000.40.450.50.55time (ps)CMD GroupDepartment of Physics0.60.65

Ignoring initial T transient-542.8-542.82Eham 96-542.980.40.450.50.55time (ps)CMD GroupDepartment of Physics0.60.65

So Better but still some wobble in T – why? T is only strictly defined as a macroscopicquantity – what you are seeing is theinstantaneous KE of a 2-particle system! Hence it is the average T that is importantand should be conserved: T 217 140 K And that will have a stat. mech. finite sizevariation given by δT2 T* 273 129 KCMD GroupDepartment of PhysicsT 3 N ions

CASTEP MD keywordsMost set in the param file: task Molecular Dynamicsmd num iter 10000md delta t 1.0 fsmd ensemble NVE or NVT or NPH or NPTmd temperature 300 Kmd thermostat Langevin or Nose-Hoovermd barostat Andersen-Hoover orParrinello-Rahmanshould be obvious but what about md ion t ormd extrap? What do they do?CMD GroupDepartment of Physics

Nosé-Hoover keywords Nosé-Hoover chains are a standarddeterministic way of thermostating system– Add an extra degree of freedom to the Lagrangian, torepresent heat-bath with coupling depending on theinstantaneous and target temperatures– But is not guaranteed to be ergodic One way to improve this is to add a thermostatto the thermostat etc resulting in a NoséHoover chain– md nhc length 5 sets the length of this chain– md ion t 100 fs sets the characteristic time forthe feedback – for most efficient thermostating youwant to set this time to resonate with dominant periodof your systemCMD GroupDepartment of Physics

Langevin keywords Langevin dynamics are an alternative andstochastic way of thermostating system– Implements a heat bath via FluctuationDissipation theorem– md ion t 100 fs sets the characteristictime for the feedback - set this to be longerthan the dominant period of your system– Typically 5*md ion t is sufficient to lose alltrace of initial conditions and be in equilibrium– Guaranteed to be ergodic if run long enoughCMD GroupDepartment of Physics

Barostat keywords What about the barostat? How is thatcontrolled? In all MD schemes, the barostat isimplemented by giving something afictitious “mass”– Andersen-Hoover uses log(V/V0) whilstParrinello-Rahman uses the cell h-matrix In both cases, this “mass” is set bymd cell t which sets the time scale forrelaxations of the cell motion. Should beslow CMD GroupDepartment of Physics

Extrapolation ExplainedBackground With ab initio MD, the forces and stresses arederived from the wavefunction ϕ– Hence need a converged ϕ at each time step With CPMD, this is achieved by integrating thewavefunction and the ionic positions together CASTEP uses BOMD and hence must reminimise ϕ each time, which is costly Wavefunction extrapolation is a useful speedup:– instead of using ϕ(t) as the initial guess at the newϕ(t δt) we extrapolate forwards in time using the MDintegrator as a frameworkCMD GroupDepartment of Physics

Extrapolation keywords BUT we do not know the acceleration of ϕ– Approximate it using known change in ϕ over previous time steps If we use the current value of ϕ and that at the previoustime step, then we have a 1st order extrapolation scheme:md extrap first Using the pre-previous time as well leads to a 2nd orderscheme: md extrap second We can also switch between 1st and 2nd on alternate stepsas a compromise: md extrap mixed The extrapolation can be done using coefficients fitted tothe instantaneous behaviour of the ionic MD(md extrap fit true) or using constant coefficients(md extrap fit false)CMD GroupDepartment of Physics

Go-faster Stripes CASTEP uses convergence window todetermine SCF convergence– default is for elec convergence win 3SCF iterations to be within elec energy tol(default 10-5 eV/atom) With well-behaved MD this can be over-kill– The extrapolation saves many SCF cycles– Hence can use md elec energy tol andmd elec convergence win to slackentolerances if all is well.CMD GroupDepartment of Physics

Doing N2 “properly” 1) Do a proper convergence test for cut-offenergy at fixed k-sampling400 eV 2) Check for finite size interactions5x5x5 A, 0.01 charge isosurfaceCMD GroupDepartment of urfaceisosurface6x5x5 A,A, 0.0010.01 charge

Doing N2 “properly” Now do geometry optimisation:δE 0.1 meV, final freq. est. 2387.5 cm-1(this is automatic from BFGS analysis)τ 1/(100.c.ν) 15 fsec so δt 1 fsec OK? Can change units of CASTEP input/output– e.g. energy unit kcal/mol– e.g. frequency unit THz, etc Now do NVE run – best for testing qualityof MD – using default T 273 K:CMD GroupDepartment of Physics

Doing N2 “properly”-543.41091Eham 1096-543.4109700.020.040.06time (ps)CMD GroupDepartment of Physics0.080.1

Doing N2 “properly”500450400T (K)35030025020015010050000.020.040.060.08time (ps)CMD GroupDepartment of Physics0.10.120.14

Now what? Problem is in the velocity initialisation:– assigning a temperature means a randomvelocity to each degree of freedom– this leads to motion in arbitrary directions– hence molecule rotates– hence falls foul of small size in y & zSolution is to usea 7x7x7 box orcontrol the initialvelocitiesCMD GroupDepartment of Physics

Default Nosé-Hoover in 7A3 box T150-300 299 256 K Eham -543.40 0.02 eV160014001200T (K)1000800600400200000.050.10.15time (ps)CMD GroupDepartment of Physics0.20.250.3

Velocity Control If doing NVE or NPH then can set T 0 K– But not if doing NVT or NPT!– So any initial velocity comes from the initialstrain w.r.t. equilibrium, or by user input Can set up any condition by editing the.cell file, e.g.%block IONIC VELOCITIESang/ps12.7 12.7 12.70.0 0.0 0.00.0 0.0 0.00.0 0.0 0.0 etc %endblock IONIC VELOCITIESCMD GroupDepartment of PhysicsNB 3*12.7 Ang/ps speed of sound in siliconHence can simulate highvelocity shock, nonequilibrium MD, etc

Non-Equilbrium MD Movie generated in FHI using PyMol and MovieMaker Bottom-most atom given initial velocity, others at rest CMD GroupDepartment of Physics

NPT Statistical Mechanics Phosphorus (II) iodide soft molecular crystal with triclinic cell Ecut-off 300 eV, 3x2x2 k-points T 250 K, P 50 MPa Highly flexible cell - βT 5.4 0.1 GPaCMD GroupDepartment of Physics

Path Integral MD Hydrogen defect in silicon– Important defect with strong coupling ofquantum ZPM to surrounding silicon latticemd use pathint truemd num beads 16md pathint staging truenum farms 16 Movie generated from .md file usingPovRay to render each timestepCMD GroupDepartment of Physics

Path Integral MDCMD GroupDepartment of Physics