Materials Process Design and Control Laboratory

1. Using multi-body energy expansions from ab-initio calculations for computation of alloy phase structures V. Sundararaghavan and Prof. Nicholas Zabaras Materials Process Design and Control Laboratory Sibley School of Mechanical and Aerospace Engineering188 Frank H. T. Rhodes Hall Cornell University Ithaca, NY 14853-3801 Email: zabaras@cornell.edu URL: http://mpdc.mae.cornell.edu Materials Process Design and Control Laboratory

2. PREDICTION OF STABLE STRUCTURES Ca Cu hP6 oP12 Stable Pt clusters (Doye and Wales, New J. Chem., 1998) Stable configurations of adsorbed species • Computational techniques • -Exhaustive or heuristic search aided by DFT calculations • Cluster expansion

3. Comparison with CE Cluster expansion • Only configurational degrees of freedom • Relaxed calculation required but only a few calculations required • Periodic lattices, Explores superstructures of parent lattice Multi-body expansion • Configurational and positional degrees of freedom • Relaxed DFT calculations are not required • Periodicity is not required • Requires a large number of cluster energy evaluations • Convergence issues Materials Process Design and Control Laboratory

4. Hybrid cluster expansions • Allow positional degrees of freedom in cluster expansions • For periodic lattices Cluster expansion for the fixed lattice Pair potentials for local relaxations Geng, Sluiter et al, Phys Rev B 2006 Materials Process Design and Control Laboratory

5. Multi-body expansion ∑ ∑ ∑ = + + + … Position and species Total energy Symmetric function JW Martin - Journal of Physics C, 1975, Empirical potentials (3 body): Murrell-Mottram (Mol. Phys 1990) Materials Process Design and Control Laboratory

6. Multi-body expansion Example of calculation of multibody potentials E1(X1) = V(1)(X1) E1(X2) = V(1)(X2) E2(X1,X2) = V(2)(X1,X2) + V(1)(X1) + V(1)(X2) Evaluate (ab-initio) energy of several two atom structures to arrive at a functional form of E2(X1,X2) Inversion of potentials V(2)(X1,X2) = E2(X1,X2) - (E1(X1)+ E1(X2)) = Increment in energy due to pair interactions Drautz, Fahnle, Sanchez, J Phys: Condensed matter, 2004 Materials Process Design and Control Laboratory

7. Multi-body expansion Inversion of potentials EL is found from ab-initio energy database, L << M Calculation of energies Drautz, Fahnle, Sanchez, J Phys: Condensed matter, 2004 Materials Process Design and Control Laboratory

8. Fitting energy surfaces • To calculate the energy of a 3 body structure (E3), we need to identify E2 and E1, values. • Two body energy E2(X1,X2) is the energy of an isolated cluster of 2 atoms at positions X1 and X2. • The database may not contain this energy since the energy values have only been obtained for atoms at locations (xi,yi) that are different from (X1,X2) • We use interpolation methods for retrieving energy at (X1,X2) from the database of energies at (xi,yi) . For example, we can use a polynomial interpolation of the form: Interpolation allows us to compute a large number of energies from a well-sampled database Materials Process Design and Control Laboratory

9. Smolyak algorithm Extensively used in statistical mechanics Provides a way to construct interpolation functions based on minimal number of points Uni-variate interpolation Multi-variate interpolation Smolyak interpolation Accuracy the same as tensor product Within logarithmic constant Increasing the order of interpolation increases the number of points sampled

10. Results in multiple orders of magnitude reduction in the number of points to sample Smolyak algorithm: reduction in points For 2D interpolation using Chebyshev nodes Left: Full tensor product interpolation uses 256 points Right: Sparse grid collocation used 45 points to generate interpolant with comparable accuracy For multi-atom systems, sample all combinations of atoms (eg. E(A-A-A), E(A-A-B), E(A-B-B),E(B-B-B) and construct interpolants.

11. 4 5 4 a b b 2 3 1 2 3 1 5 a CLUSTER REPRESENTATION Specification of clusters of various order by position variables • Convex hull technique to represent all atoms in the positive z-direction • Use independent coordinates to represent the cluster geometry A point in 6 dimensional space Materials Process Design and Control Laboratory

12. Computations were performed in parallel on a 64 node quad-processor LINUX cluster CLUSTER ENERGY COMPUTATIONS • Executables • Cluster coordinates • Energy interpolation • Batch input for PWSCF • Read energies from PWSCF • Energy calculation • Plane-wave electronic density functional program ‘quantum espresso’ (http://www.pwscf.org) calculations are used to compute energies given the atomic coordinates and lattice parameters. • These calculations employ LDA and use ultra-soft pseudopotentials. • Single k-point calculations were used for isolated clusters, the cell size was selected so that the effect of periodic neighbors are negligible. • For multi-component systems, a constant energy cutoff equal to cutoff for the hardest atomic potential (e.g. B in B-Fe-Y-Zr) is used. • MP smearing (ismear=1, sigma=0.2) is used for the metallic systems. Materials Process Design and Control Laboratory

13. LINKING THE MULTIBODY EXPANSION TO OTHER SOFTWARE Multi Body Expansion (MBE) The multibody expansion software written in C++ Two parts: potential generation & energy computation Energy computation part is the Hamiltonian Molecular dynamics- LAMMPS Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) is a classical molecular dynamics (MD) code developed by S. Plimpton et. al (Sandia national lab) Directly linked energy computation part in LAMMPS with MBE Useful for molecular dynamics and energy minimization Monte Carlo- MCCCS Towhee Monte Carlo for Complex Chemical Systems (MCCCS) developed by M. G. Martin, J. I. Siepmann et. al. Available at http://towhee.sourceforge.net/ Fortran based code. Linked Towhee and MBE using a library Performs a variety of calculations in all ensembles

14. ENERGY SURFACES FOR ISOLATED CLUSTERS Y X Platinum isolated cluster energies computed using multi-body potentials a = (1.5*X+0.5)*7.5 Bohr b = (1.5*Y+0.5)*7.5 Bohr b a Materials Process Design and Control Laboratory

15. COMPUTATION OF CLUSTER ENERGIES The complete potential surface for a 3 Pt cluster. Figure (a) shows computed Platinum three-atom cluster energies, while (b) shows extension of energies using pair potential terms beyond the cutoff. Distance between atoms 1-3 (Bohr) Distance between atoms 1-3 (Bohr) Distance between atoms 1-2 (Bohr) Distance between atoms 1-2 (Bohr) (b) (a) Materials Process Design and Control Laboratory

16. Convergence results for different energy functions Energy of the system scales as n1/2 where n is the number of atoms Order of interactions necessary for full convergence: 4 Energy of the system scales as n2 where n is the number of atoms Order of interactions necessary for full convergence: 2 5 body term contribution =0 3 body term contribution =0 Materials Process Design and Control Laboratory

17. Convergence results for different energy functions Using pair potentials: Lennard Jones for Helium atoms Order of interactions necessary for full convergence is 2 as expected (since it is a pair potential) 3 body term contribution =0 Materials Process Design and Control Laboratory

18. Oscillations in MB energy for complex energy functionals EAM potentials: Platinum system Energies (En) calculated from an n-body expansion correct energy • Energies oscillate around the true energy • -Approach: Low pass filtering (convolution operation) that cuts off high frequency oscillations. • -Compute the energy at the minima using self consistent field calculation Materials Process Design and Control Laboratory

19. Computation of MBE energy filters + + .. + Weighted MBE Is the total energy correlated with structural energies of clusters ? Materials Process Design and Control Laboratory

20. Weighted MB energy a1 a2 True energy a3 Materials Process Design and Control Laboratory

21. Extrapolatory tests on weighted MBE 16 atom Au-Cu FCC cluster AuCu3 4 unit cell, 4 at/cell True energy MBE 4th order Weighted MBE energies once built for a small set of configurations provide accurate energy fit for various different inter-atomic distances within that configuration.

22. Selection of order of expansion Test various MBE orders in extrapolatory modes Weighted 4th order MBE Weighted 2nd order MBE True energies Cohesive energy (Ryd) Cohesive energy (Ryd) True energies Weighted 3rd order MBE Weighted MBE expansion coefficients are fitted using 12 atom cluster energies and the results are presented for a 16 atom cluster. Energies may differ but the weighted MBE captures the energy minima within 4th order expansion. Cohesive energy (Ryd) True energies Materials Process Design and Control Laboratory

23. Number of isolated cluster calculations Depth of interpolation 0 Weighted MBE 4th order -1 4 120 Actual energy -2 Cohesive energy (Ryd) 4 560 -3 -4 + 4 1820 -5 -6 + 6 6.5 7 7.5 8 8.5 9 Lattice parameter (Bohr) Energy minima Platinum clusters 16 atom FCC cluster • Coefficients obtained using an 12 atom cluster energies at different lattice parameters Materials Process Design and Control Laboratory

24. 1 Number of isolated cluster calculations Depth of interpolation 0 Weighted MBE 4th order -1 4 276 -2 -3 Cohesive energy (Ryd) -4 4 2024 -5 -6 + -7 4 10626 -8 -9 6 6.5 7 7.5 8 8.5 9 Lattice parameter (Bohr) + Energy minima Platinum clusters 24 atom FCC cluster Actual energy Materials Process Design and Control Laboratory

25. Number of isolated cluster calculations Depth of interpolation -1 4 276 -2 -3 -4 4 2024 Cohesive energy (Ryd) -5 -6 + 4 10626 -7 -8 + -9 6 6.5 7 7.5 8 8.5 9 Lattice parameter (bohr) Platinum clusters A random 24 atom configuration Weighted MBE 4th order Actual energy Materials Process Design and Control Laboratory

26. Stable phase structures of Au-Cu alloy Super-cell approach For computing stable structures of periodic lattices, a 4x4x4 supercell (216 atoms) is used as an approximation. Weighted MBE is several orders of magnitude faster than a relaxed DFT calculation. Useful for amorphous structures Small cluster calculations are used to compute the weights in the weighted MBE expansion FCC structures are considered here for Au-Cu. Materials Process Design and Control Laboratory

27. Stable phase structures of Au-Cu alloy AuCu3 cell relaxation Au3Cucell relaxation a = 6.62 bohr = 3.50 A a = 7.3 bohr = 3.86 A 3x3x3 supercell AuCu3 lattice parameter: 3.76 A Au3Cu lattice parameter: 4.04 A Materials Process Design and Control Laboratory

28. Stable phase structures of Au-Cu alloy AuCu3 cell relaxation Au3Cucell relaxation a = 6.71 bohr = 3.55 A a = 7.4 bohr = 3.92 A 4x4x4 supercell AuCu3 lattice parameter: 3.76 A Au3Cu lattice parameter: 4.04 A Materials Process Design and Control Laboratory

29. APPLICATION TO SURFACE PHENOMENA Minimum energy surface of h on Pt(111) Plot of minimum energy in z direction for the primitive cell Highly anharmonic potential energy surface FCC->HCP (55 mev), FCC->TOP (160 mev) H confined to FCC-HCP-FCC valleys FCC site G.Kallen,G.Wahnstrom, Quantum treatment of H on a Pt(111) surface, Phys Rev B, 65 (2001) (Baskar and Zabaras, 2007)

30. Conclusions • MB expansion provides atom position dependent potentials that are used to identify stable structures. • Ab-initio database of cluster energies are created and interpolation for various cluster positions are generated using efficient sparse grid interpolation algorithms. • Weighted MBE is fast and captures the energy minima within a small order of expansion. • Technique is applicable to study stability of amorphous systems, molecules and clusters. Materials Process Design and Control Laboratory