Assorted Topics [Partly based on Chapters 9 & 10, Sholl & Steckel]

1. Assorted Topics[Partly based on Chapters 9 & 10, Sholl & Steckel] Molecular Dynamics Monte Carlo – Revisited Accuracy of Standard DFT Some beyond-DFT methods Multiscale Modeling Final Thoughts

2. Molecular DynamicsA bird’s eye view In the absence of external forces and damping … Time, t = 0 Time, t = dt v2 dx = v1dt v1 v2 v2dt v1 In the presence of external forces … Time, t = 0 Time, t = dt F1 dx = v1dt+0.5(F1/m)dt2 v1 v1+(F1/m)dt

3. Molecular Dynamics (MD) • External forces could be because of external electric/magnetic fields, but mainly due to interatomic interactions; in the latter case • Fi = dE(R1,R2,…RN)/dRi • E(R1,R2,…RN) may be obtained from DFT (ab initio MD) or empirical potentials (classical MD) • Average kinetic energy  defines temperature • 1.5NkT = Si0.5mivi2 • Time step has to be of the order of 10-15 seconds (fs) so that trajectories during vibrations can be captured

4. Ensembles • Microcanonical Ensemble • E, N, V constant during the simulation (closed, isolated system) • Temperature fluctuates with time • Canonical Ensemble • N, V, T constant during the simulation (system in contact with a heat bath that maintains a constant temperature) • At the end of every few time steps, velocity is rescaled so that the temperature is maintained • Molecular Statics • If the velocity is rescaled to zero at the end of every time stem, we reach the zero temperature ground state (this is the same as geometry optimization)

5. Ab Initio MD • Electrons are treated quantum mechanically, but nuclear motion is treated classically (Newtonion mechanics) • At the present time, these simulations can extend only up to nanoseconds (if that) • If the time step is fs, this implies 106 iterations • Classical MD can be extended to microseconds-milliseconds as energy and force evaluations are cheaper (but also less accurate)

6. Hyperdynamics • Standard MD is impractical to solve long time (> microseconds) diffusion problems, as as it may involve larger than 1010 time steps! • A way to accelerate the dynamics is by making the potential wells at atomic sites “shallow” (proposed by Art Voter); this artificially slows down the atomic vibrations (hence timesteps can be large) by several orders of magnitude while preserving the basics physics of diffusion

7. Standard Implementations • Euler method • xi(t + dt) = xi(t) + vi(t) dt • vi(t + dt) = vi(t) + (Fi/mi) dt • Correct up to first order, and hence requires very small dt • Verlet method • Taylor expand: xi(t + dt) = xi(t) + vi(t) dt + (1/2)(Fi/mi) dt2 + (1/6)(d3xi(t)/dt3) dt3 + O(dt4) • Note: xi(t - dt) = xi(t) - vi(t) dt + (1/2)(Fi/mi) dt2 - (1/6)(d3xi(t)/dt3) dt3 + O(dt4) • Add: xi(t + dt) = 2xi(t) - xi(t - dt) + (Fi/mi) dt2 • Correct up to 3rd order! • To determine the position at the next time step, we need: current position, previous position and current force (what do we do for the 2nd step of the simulation?) • Velocity still given by Euler formula • No dependence of position evolution on velocity  good for constant energy MD simulations • If we want to perform constant T simulations, then we need to rescale velocities, and hence the position and velocity evolution has to be coupled • How do we do this?

8. MD Example I: Liquids and Amorphous Systems: InPLewis, De Vita & Car, PRB 57 1594 (1998) • The melt & quench approach Each plateau: “equilibration”

9. Radial Distribution Function: InP • g(r) = r(r)/(4pr2ravg) • P-P bonds found in the amorphous phase (not found in crystalline phase)

10. Example II: Equilibrium geometries of ultrasmall nanocrystals • 13-atom Pt cluster (particularly stable “magic size” cluster, and occurs in catalysts) • High T MD simulations, followed by “quenching” (i.e., standard geometry optimization) of specific low-energy configurations • Several disordered clusters which we couldn’t have guessed based on symmetry, are as stable as symmetric geometries; MD is thus useful in the identification low-energy structures

11. Monte Carlo (MC) RevisitedThe case of diffusivity • MC methods are stochastic, and involves ensemble averages of properties (in contrast to MD simulations which involve time averages of properties once equilibrium is reached) • Let us consider diffusion (in 1-d): For the case of vacancies, we will attempt to relate the macroscopic diffusivity and measured activation barriers for long range diffusion (Q) to atomic-level site-to-site hopping barriers (E) of vacancies • D = D0exp[-Q/kT] • <x2> = Dt (i.e., D is related to an ensemble average) • Probability of hopping, p = exp[-E/kT] • D, Q, < x2> are macroscopic quantities, while p and E are atomic-level properties • E may be calculated using DFT for a variety of situations, followed by MC simulations of large scale diffusion vs. microstructure

12. t = 0 t = 1 t = 2 t =3 <x2> = 0 <x2> = 1 <x2> = 2 <x2> = 3 MC Ensemble Averaging • 1-d lattice with a vacancy • We consider an ensemble of 8 lattices • At high T, p = 1 • It is clear that <x2> is linearly related to t  at T = , D = 1 • This simple scheme can be used to determine the macroscopic diffusivity if the microscopic information is available

13. Density Functional TheoryFormidable problem  Manageable problemTrue! But what to expect? Density Functional Theory (DFT) [W. Kohn, Chemistry Nobel Prize, 1999] 1-electron wave function (function of 3 variables!) The “average” potential seen by electron i 1-electron energy (band structure energy) Energy can be obtained from r(r), or from i and ei (i labels electrons) DFT provides a practical prescription for computing the total ground state energy When combined with statistical thermodynamics, DFT provides free energies too

14. Predictions of geometry • Structural details predicted typically to within 1% of experiments for a wide variety of solids and molecules • Results from various sources [collected in Ramprasad, Shi and Tang, in Physics and Chemistry of Nanodielectrics, Dielectric Polymer Nanocomposites (Springer)]

15. Predictions of other basic properties

16. Band offsets

17. The Exchange-Correlation FunctionalPerdew’s Ladder & the problem with choices! • Unfortunately, a step up the ladder does not guarantee greater accuracy! LDA GGA-PW91 GGA-PBE GGA-BP88 GGA-BLYP GGA-RPBE GGA-revPBE GGA-OLYP GGA-OPBE B3LYP B1LYP KMLYP O3LYP X3LYP LDA+X GGA+X GGA-PBE0 HSC LDA+U GGA+U LDA+G0W0 GGA+G0W0 LDA+GW0 GGA+GW0 LDA+GW GGA+GW QMC HF HF-CI MPn CCSD

18. Beyond-DFT • (Not so) quick fixes • Self-interaction correction (VH≠ 0 for H2+ in DFT) • DFT + X (X: exact exchange; or is it?) • DFT + U (inspired by the previous 2 methods) • Beyond the 1-electron approximation • GW (G: Green function; W: screened Coulomb interaction) • Replacing W by the bare Coulomb interaction results in the HF exchange potential • Quantum Monte Carlo (QMC) • Direct solution of the many-electron Schrodinger equation • Multi-dimensional integrals evaluated using MC methods • HF + Configuration Interaction (CI) & related strategies • Largely used by chemists and has been applied only to molecules so far • The “gold standard” • Order-N methods • Scaling: • Standard DFT: N2-N3 • HF: N4 • GW, QMC: N4-N5 • CI: N6 • Linear scaling (or Order-N) methods take advantage of localization Physics Today, April 2008

19. Modeling vs. Scales Time Engineering Design Years • If more than one “box” is involved in a computation  multi-scale modeling Hours Finite element Analysis (Continuum/classical) Minutes Seconds Mesoscale Modeling (Semi-classical) Microseconds Molecular Mechanics Nanoseconds Picoseconds Quantum Mechanics Femtoseconds mm nm mm m Å Distance e.g., density functional theory (DFT)

20. Multi-scale modeling:Two perspectives • Implicit multi-scale approach • Perform computations at a higher level, with parameters and constants occurring in the equations determined at a more accurate “lower” level of theory (e.g., solving diffusion equations with diffusion coefficients determined using atomic level computations) • Not new! • Explicit multi-scale approach • Simultaneously perform different levels of computation for the same system (e.g., atomic level computation in a small important region plus finite element type computation further away) • Very recent, and very difficult, approach

21. Implicit Multi-scale modeling:Historical examples • At any level of theory, the parameters found in equations conceptualize events occurring at smaller scales not directly treated by the theory • Elasticity theory: sij = Cijklekl (elasticity theory assumes the existence of the elastic constants Cijkl which contain the sum total of relevant atomic level information) • Electromagnetism: although classical electromagnetic theory does not explicitly consider electrons, permeability (m), permittivity (e) and conductivity (s) are input material properties, whose existence is determined by electronic structure • Kinetic theory of gases: PV = nRT • Phase transformations, statistical thermodynamics • Diffusion coefficients, kinetic rate constants (kMC) • Etc., etc.

22. Explicit Multiscale ApproachExample: Crack propagation in silicon • “Concurrent coupling of length scales: Methodology and applications”, Broughton et al, Phys. Rev. B60, 2391 (1999) • The “gold standard” in the area of multi-scale computation • Silicon was chosen in this study due to • Its industrial relevance • The availability of robust & well-tested empirical potentials (for molecular dynamics simulations) • The accuracy of the tight-binding method for Si (which is computationally more efficient than fully first principles methods) • Same approach can be used for other materials (metals, organics, ceramics, etc.) • Combined tight-binding, molecular dynamics & finite element method • Uses parallel processing (especially for the tight-binding part)

23. The multi-scale framework • Tight-binding at the crack tip • Semi-empirical potentials (Stillinger-Webber) + molecular dynamics in the neighborhood of the crack tip (to model plastic behavior) • Finite element method elsewhere (elasticity theory)

24. “Hand-shaking”Handling the boundary regions • Within each region, a Si atom behaves like a Si atom • But, at the boundary, we need to have a silicon-hydrogen hybrid (silogen) • Within the MD region, we have a collection of Si atoms • At the boundary between the MD and the finite element (FE) regions, the FE mesh is equivalent to the the Si unit cell, but away from the boundary, the mesh sizes can be large

25. Results • Can reliably predict the stress fields in materials • Provides a framework for the modeling of the physical or mechanical breakdown of materials • All length/time scales treated simultaneously (example of explicit multi-scale modeling) • This may not be the most desired approach • Perhaps can break the problem into pieces that can be handled separately

26. Final thoughts … • Several types of computational methods at various length & time scales • Need to use judgment in determining the right type of method or model to study a given problem • No computational method should be used as a “black box” • Challenges • Geometric complexity; Chemical complexity; Disorder; Multiple phenomena; Models of reality • DFT • Typical scaling: 1980s: few-several atoms; 1990s: few-several 10s of atoms; 2000s: few-several 100s of atoms; future: ??? • Future directions: Materials Design using existing methods; ab initio thermodynamics; beyond-DFT methods • Regarding DFT predictions … • DFT is in general complementary to experiments • Can provide insights, trends and guidance • Sometimes, computation may be the only option • … Hypothesis non fingo … [Isaac Newton] Science is primarily an experimental investigation, and so experiments are the ultimate authority in determining the validity of computational theory! • “The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws lead to equations much too complicated to be soluble.” [Paul Dirac]