Computational Modelling of Materials

1. Computational Modelling of Materials Recent Advances in Contemporary Atomistic Simulation Or: Understanding the physical and chemical properties of materials from an understanding of the underlying atomic processes http://secamlocal.ex.ac.uk/people/staff/ashm201/Atomistic/ Lectures Introduction to computational modelling and statistics 1 Potential models 2 Density Functional (quantum) 1 3 Density Functional 2 4

2. Introduction The increased power of computers have allowed a rapid advance in the use of simulation techniques for modelling the properties of materials. Why do it? • Interpret of experiment • Extrapolate experimental data • Empirical Search • Prediction of New Effects But how? The answer depends on the length and time-scale

3. Atomistic Simulation: Choices Which Technique? • Energy Minimisation • Molecular Dynamics • Monte Carlo • Genetic Algorithms How do you calculate the forces? • Interatomic Potentials • Quantum Mechanics What Conditions? • Select Ensemble • Select Periodic Boundary Conditions

4. Simulation of Forces All the atomistic simulation techniques require that the total interaction energy is evaluated and are more efficient if the forces between every atom is evaluated. • Interatomic Potentials (Force fields) • Parameterised equations describing forces - fast • Empirical Derivation • Non-empirical Derivation • Quantum Mechanics • direct solution of the Schrodinger Equation – slow+reliable? • Semi-empirical • Density functional approach • Molecular Orbital approach

5. Simulation Techniques: • Energy Minimisation • Calculate Lowest Energy Structure • Gives structural, mechanical and dielectric properties • Molecular Dynamics • Calculates the effect of Temperature • Gives dynamics e.g. diffusivity • Monte Carlo • Calculates a range of structures • Gives the thermally averaged properties • Genetic Algorithms • Calculates a range of structures • Efficient search for global minimum

6. Atomistic simulation - Dynamics Summary • Molecular Dynamics can provide reliable models • Effect of Temperature • Time evolution of system • Highly suited to liquids and molecular systems • Calculate dynamical properties, e.g. diffusivity • PROVIDED • Reliable potential models • Molecular Dynamics • Robust and reliable for solids and their surfaces • BUT • Takes a long time to search configurational space • Does not easily allow atoms to pass over large energy barriers • Can use constrained methods – but usually need to know where the atom/molecule needs to go. } Monte Carlo + Genetic Algorithms

7. Monte Carlo • In the widest sense of the term, Monte Carlo (MC) simulations mean any simulation (not even necessarily a computer simulation) which utilizes random numbers in the simulation algorithm. • The term “Monte Carlo” comes from the famous casinos in Monte Carlo. • Another closely related term is stochastic simulations, which means the same thing as Monte Carlo simulations.

8. Monte Carlo • Metropolis MC • A simulation algorithm, central to which is the formula which determines whether a process should happen or not. Originally used for simulating atom systems in an NVT thermodynamic ensemble, but nowadays generalized to many other problems. • Simulated annealing • The Metropolis MC idea generalized to optimization, i.e. finding minima or maxima in a system. This can be used in a very wide range of problems, many of which have nothing to do with materials. • Thermodynamic MC • MC when used to determine thermodynamic properties, usually of atomic systems.

9. Monte Carlo • Kinetic MC, KMC • MC used to simulate activated processes, i.e. processes which occur with an exponential probability • e−Ea/kT • Quantum Monte Carlo, QMC • A sophisticated electronic structure calculation method. • Diffusional Monte Carlo (stochastic projector technique, which solves the imaginary time-dependent Schroedinger equation). In theory DMC is exact!

10. Metropolis Monte Carlo THE JOURNAL OF CHEMICAL PHYSICS VOLUME 21, NUMBER 6 JUNE, 1953 • The approach is to calculate energy Ei then • randomly move an atom or molecule to give a new energy, Ej • Then decide whether to accept or reject move • Can easily extract Thermodynamic properties within NVT -Canonical Ensemble

11. Selection U Local Minima • Metropolis Monte Carlo • If the new energy is lower (i.e. a more stable structure) then accept the move • If the new energy is higher (less stable) then • generate a random number between 0 and 1 • calculate: Pij = exp(-(Ej -Ei)/kT) • only accept the move if, Pij is higher than the random number. • This enables the system to focus on the important configurations Global Minimum q

12. Example of Use: TiO2 • Particularly powerful when used with energy minimisation • Prediction of crystal structure without prior knowledge of atom positions • Freeman etal J.Materials Chem, 1993, 3, 531 • used Monte Carlo to select a number of likely structures • followed by energy minimisation of each candidate to locate the precise atom positions • Successfully found all the phases of TiO2

13. Example of Use:Template Design Lewis, et al, Nature , 382, 604.

14. Predicted New Template for Levyne • ZEBEDDE suggests 1,2-dimethylcyclohexane as a template for LEV • Using 2-methylcyclohexylamine, a LEV structured CoAlPO (DAF-4) is formed • Barratt et al, Chem Commun,1996, 2001

15. Computer Designed Template • Bi-cyclohexane motif • Amine derivative • 4-piperidino piperidine Co-AlPO4 Preparations • 170oC, 4hours • Chabazitic structure • NO competing phase

16. Problems with Monte Carlo • The major problem is that computer resources • A lot of configurations need to be sampled to obtain reasonable statistics • A lot of configurations need to be sampled to ensure that you have found the global minimum • Hence need to keep rejection rate down • Has no ‘memory’ of good solutions

17. Problem:Structure of Clusters and Nuclei • Clusters span a wide range of particle sizes – from molecular (well separated, quantized states) to micro-crystalline (quasi-continuous states). • How do properties change as they grow ? • Clusters constitute new materials (nanoparticles) which may have properties that are distinct from those of discrete molecules or bulk matter. • New chemistry ?

18. Nucleation of Zinc SulphideS.H. Gomez, E. Spano, C.R.A. Catlow • Generating Nuclei via molecular dynamics • Start with individual atoms are monitor how and they assemble. • ZnS • In the bulk both ions 4-fold coordinated • But get 3-fold coordinated clusters. (ZnS)25 (ZnS)12

19. Comparison of Stability • Although still small – show continued stability of ‘bubble’ structures • (ZnS)47 Bulk like cluster (300 kJ/mol less stable) Bubble-like Bulk-like CHEM COMMUN (7): 864-865 APR 7 2004 + J AM CHEM SOC 127 (8): 2580-2590 MAR 2 2005

20. Alternative Approach: Genetic Algorithms for Cluster Geometry Optimisation • GA procedure is for optimising a function, structure or process which depends on a large number of variables. • Developed by computer scientists in the 1970’s. • Based on principals of natural evolution. • Works through a combination of mating, mutation and “natural selection”. Roy L. Johnston, University of Birmingham DALTON T (22): 4193-4207 2003

21. gene B D A C A allele chromosome GA Definitions • Chromosome – a string of variables (genes) corresponding to a trial solution. • Allele – the value of a particular gene (i.e. variable).

22. GA Approach • Take a Population – the set of trial solutions. • Measure of the quality of each member of the population - Fitness (usually by calculating the total interaction energy) • Proceed with mating - the overall process of selecting strings (parents) and exchanging their genes to produce new strings (offspring).

23. Roulette Wheel Selection: parents are chosen with a probability proportional to their fitness: Selection Process

24. parents + Single Point Crossover offspring + Generating new structures • Crossover – the process of exchanging genes between chromosomes. • Some offspring will be fitter than their parents. • Due to crossover the GA effectively explores the parameter space in parallel.

25. Single Point Mutation Possible Problem • It is possible to get stagnation – where certain structures can appear to be ‘frozen-in’. • Overcome by introducing new genetic material which ensures population diversity – preventing in-breeding and stagnation. • Mutation – randomly changing certain genes in selected members of the population.

26. Protein folding G.A. Cox, T. V. Mortimer-Jones, R. P. Taylor and R. L. Johnston, Theor. Chem. Acc.112, 163-178 (2004). Crystal structure solution K.D.M. Harris, R.L. Johnston and B.M. Kariuki, Acta Cryst. A54, 632-645 (1998). Spectral deconvolution Conformational analysis I GA 2 Some Other Applications of GAs A variety of GAs have now been written for cluster geometry optimization.

27. Apply “cut and paste”crossover operator One new cluster generated from each mating operation. Perform energy minimisation using BFGS algorithm The Birmingham Cluster GARoy L. Johnston • Mutation achieved by randomly moving a fraction (N/3) of atoms. Mutation probability:  Pmute = 0.1 • The mutation operator acts on the offspring.

28. Ionic MgO Clusters Rigid Ion Model • First term – long-range electrostatic Coulomb energy. • Second term – short-range repulsive Born-Mayer potential, which reflects the short range repulsive energy due to overlap of the ions. PHYS CHEM CHEM PHYS 3 (22): 5024-5034 2001

29. Formal charges q =  2

30. Formal charges q =  1

31. Variation of Structure with Magnitude of Formal Ion Charge q (MgO)8 (MgO)9 (MgO)12

32. Conclusions – GA • The GA is an efficient technique for searching for global minima –a variety of potentials (LJ, Morse, Ionic, MM, Gupta, TB, EAM …) have been studied. • As with Monte Carlo the chief problem is the time taken to investigate the different possible structures • When particles become much bigger, e.g. beyond 10nm, most efficient is Molecular Dynamics • Care needed in generating structures

33. Cl Cl- O = C = O H-F K+ + + Cl - - + - + - C Cl Cl - - - - ++ ++ Electrostatic Forces (Multipolar Forces) • Most molecules have an uneven distribution of charge, e.g. ions quadrupole octopole dipolar This leads to electrostatic (Coulomb) forces between the molecules. If we approximate the charge distribution as a collection of discrete charges qi, where qi are charges in molecule 1 and qj are those in molecule 2

34. Potential Models (Force Fields) • Potential models rely on Born-Oppenheimer, ignore electronic motions and calculate the energy of a system as a function of nuclear positions only • • Potential modelsrely on: • – Relatively “simple” expressions that capture the essentials of the interatomic and intermolecular interactions. Such as stretching of bonds, the opening and closing of angles, rotations about bonds, etc. • – Transferability: the ability to apply a given form for a potential model to many materials by tweaking parameters (e.g. MgO vs CeO2) taken from Dr. S. C. Glotzer’s lectures on Computational Nanoscience of Soft Materials, University of Michigan http://www.engin.umich.edu/dept/cheme/people/glotzertch.html

35. a c b d e f Composite Pair Potentials for Small Molecules • For small molecules (e.g. Ar, N2, CO2) many neglect molecular flexibility and treat the molecule as rigid. • Commonly used models include: - Lennard-Jones (12,6) e.g. CO2 LJ + Coulomb taken from Prof. K. Gubbins lectures on Computer simulation , NC State Univ http://chumba.che.ncsu.edu/

36. Flexible molecules • Total pair energy breaks into a sum of terms Intramolecular only • UvdW- van der Waals • Uel - electrostatic • Upol - polarization • Ustr - stretch • Ubend - bend • Utors - torsion • Ucross - cross Mixed terms See Leach 2nd ed., ch. 4; also, Gubbins and Quirke, pp. 25-27, 28-33

37. A Typical Force Field taken from Dr. S. C. Glotzer’s lectures on Computational Nanoscience of Soft Materials, University of Michigan http://www.engin.umich.edu/dept/cheme/people/glotzertch.html

38. A (More Complicated) Force Field Analytic expression for the CFF 95 force field

39. Some Commonly Used Models • • There are many different Potentials in the literature, particularly for organics. • In most cases, they are developed to treat a particular class of systems. • Some commonly used FFs are: (in blue: original systems studied; in red, • some useful references and/or websites) • - MM2, MM3 and MM4 (N. L. Allinger et al.) • → small organic molecules • → http://europa.chem.uga.edu/index.html • - MMFF (Merck Molecular Force Field, proposed by T. A. Halgren) • → biomolecules • → T.A. Halgren, J. Comput. Chem. 17, 490 (1996) • - AMBER (Assisted Model Building with Energy Refinement, by P. A. • Kollman et al.) • → biomolecules • → http://www.amber.ucsf.edu/amber/amber.html • - CVFF (A. Hagler -> Biosym -> MSI -> Accelrys) -> COMPASS • → biomolecules -> more general • → Dauber-Osguthorpe& Hagler

40. Some Commonly Used Models - OPLS (Optimized Potentials for Liquid Simulation, W. L. Jorgensen et al) → organic liquids →W. Damm, A. Frontera, J. Tirado-Rives, W.L. Jorgensen, J. Comput. Chem. 18, 1955 (1997); http://zarbi.chem.yale.edu/ - CHARMM (Chemistry at HARvard Macromolecular Mechanics, by M. Karplus and coworkers) → biomolecules → http://www.charmm.org/ - ECEPP (Empirical Conformational Energy Program for Peptides, by H. A. Scheraga et al.) → biomolecules → http://www.tc.cornell.edu/Research/Biomed/CompBiologyTools/eceppak/ http://www.chem.cornell.edu/has5/ - GROMOS (GROningen MOlecular Simulation, by W. F. van Gunsteren and coworkers) → biomolecules → http://www.igc.ethz.ch/gromos/

41. Other Models • • There are also potential models, such as • MOMEC (P. Comba and T. W. Hambley) and • SHAPES (V. S. Allured et al) that were developed for transition metal complexes • There are also models developed with the purpose of treating the full periodic table, such as • UFF (Universal Force Field, by A. K. Rappe et al.), • RFF (Reaction Force Field, by A. K. Rappe et al.) and • DREIDING (by S. L. Mayo et al.)

42. Problems: Unlike-Atom Interactions(non-bonding) • “Mixing rules” give the potential parameters for interactions of atoms that are not the same type • no ambiguity for Coulomb interaction • for effective potentials (e.g., LJ) it is not clear what to do • Lorentz-Berthelot is a widely used choice

43. B B B B A A A + A Problems: Unlike-Atom Interactions(bonding) • Conservation of equilibrium bond distance and energy. On altering for example, charge, adjust short range parameters to maintain distance and energy. • Issue for simple force fields • Bond energy: U = 0.5 k (r AB – r 0AB)2 If new bond is approx the equilibrium bind length then the energy of reaction about 0 energy. • Treatment is a very weak link in quantitative applications of molecular simulation

44. More Potentials for Solids • Even for polar/ionic solids there are a vast array of models, (e.g. see refs by Bush, Catlow, de Leeuw, Dove, Gale, Lewis, Jackson, Parker and Woodley) that are based on the shell model and for models based on three body potentials see refs by (S. Garofalini et al.) • There other models for metals (Finnis and Sinclair) [Phil.Mag. A 50 (1984) 45; for an improvement see Phil. Mag. A 56 (1987) 15]. • Semiconductors [Tersoff, Phys. Rev. Lett. 56 (1986) 632] extended by Brenner [D. W. Brenner, Phys. Rev. B 42 (1990) 9458] for conjugated systems, see further extensions [Stuart et al., J. Chem. Phys. 112 (2000) 6472]and [Che et al., Theor. Chem. Acc. 102 (1999) 346]. where

45. Shell Model Potential For example: polar solids • Electrostatic • despite simple expression (q1q2/r12) it has poor convergence - use methods by Ewald, Parry and Madelung etc. • Short-range • includes repulsion + dispersion A12exp(-r12/p12) - C12/r126 • where A, p and C are needed for each pair of atoms • Electronic polarisability • Via Shell model • specify shell charge and spring constant • Angle dependent forces • For polyanions

46. Ewald Method • Approach for calculating the Coulombic interaction energy • Replace point charges (charge density – delta functions) by Gaussians. • Gives • difference between Gaussians and delta functions • Interacting Gaussians • remove interaction of Gaussian with self q q q q q

47. Shell Model – many body forces • Valence electrons Massless shell • distorted by electric field, size of distortion dependent on strength of spring, i.e. variable polarisability • For quadrupolar distortions see work by P.A. Madden etal • Shell charge remains symmetric U = 0.5 k d2 Y (shell charge) Free ion polarisability a = Y2/k k (spring constant)

48. Partially covalent solids For example: work by S.H. Garofalini Two-Body Term Three-Body Term Tetrahedral B. P. Feuston and S. H. Garofalini, J. Chem. Phys., 89 (1988) 5818 (note error in Table I, where beta headings are mixed) R. G. Newell, B. P. Feuston, and S. H. Garofalini, J. Materials Research, 4 (1989) 434. S. Blonski and S. H. Garofalini, Surf. Sci. 295 (1993) 263.

49. Issues when using Potential Models • The main problem in fitting a general model is to ensure its transferability while using a reasonable number of parameters; in order to be useful the model has to be able to predict correctly properties for compounds that fall outside the set used to fit the parameters • How different models are linked together is still an area of debate – are the results meaningful? • When using a potential model, it is important to know what is being included and how, and what isn’t. • Leach, AR Molecular Modelling: Principles and Applications; 2nd Edn (2001)Pearson Prentice Hall

50. Derivation of parameters • Empirical fitting • to crystal structure, elastic and dielectric constants • problems with • validation (must not use all exptal data) e.g. ir and raman • interatomic separations far from those used in fitting e.g. at high temperatures and pressures • overcome with…. • Non-empirical fitting • to electronic structure calculations • problems with • incomplete description of forces e.g. dispersion • open shell atoms (e.g. transition metals)