Machine Learning Applied to a Variable Charge Atomistic Model for

Subscriber access provided by The University of Melbourne Libraries

Article

Machine Learning a Variable Charge Atomistic Model for Cu/Hf Binary Alloy Oxide Heterostructures Kiran Sasikumar, Henry Chan, Badri Narayanan, and Subramanian K.R.S. Sankaranarayanan Chem. Mater., Just Accepted Manuscript • DOI: 10.1021/acs.chemmater.8b03969 • Publication Date (Web): 10 Apr 2019 Downloaded from http://pubs.acs.org on April 11, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chemistry of Materials

Machine Learning a Variable Charge Atomistic Model for Cu/Hf Binary Alloy Oxide Heterostructures Kiran Sasikumar,1,, Henry Chan1, Badri Narayanan1,3 and Subramanian K. R. S. Sankaranarayanan1,2*   1

Center for Nanoscale Materials, Argonne National Laboratory, Argonne IL, USA 2 Computation

3 Department

Institute, University of Chicago, Chicago IL, USA

of Mechanical Engineering, University of Louisville, Louisville KY, USA 

Equal Contributions

Abstract Alloy oxide heterostructures are essential for a wide variety of energy applications, including anticorrosion coatings, bifunctional catalysis, energy storage, as well as emerging platforms for multi-level non-volatile memory and neuromorphic devices. The key role played by oxide composition, density, and stoichiometry in governing electro-chemical reactions, interface structure, and ionic transport in alloy oxides remains largely unknown. This is primarily owing to the lack of a variable charge interatomic potential models that can adequately describe the various atomic/ionic interactions in alloy oxides within a unified framework. Here, we introduce a charge transfer ionic potential (CTIP) model for Cu/Hf/O alloy system and demonstrate its ability to accurately capture the complex potential energy surface owing to dynamically varying interactions, including metallic (Cu/Hf), ionic (Cu/Hf/Ox), mixed environment (interfaces). We leverage supervised machine learning methods powered by genetic algorithms coupled with local Simplex optimization to efficiently navigate through a high-dimensional parameter space, and identify an optimum set of independent set of CTIP parameters. We train our model against an extensive first-principles based data set that includes lattice constants, cohesive energies, equations of state, and elastic constants of various experimentally observed polymorphs of hafnia, copper oxide, hafnium as well as copper. Our machine learned CTIP model captures the structure, elastic properties, thermodynamics and energetic ordering of various polymorphs. It also accurately predicts the surface properties of both oxides and their metals. To demonstrate the suitability of our CTIP model for investigating dynamic processes, we employed it to identify the atomistic scale mechanisms associated with the initial nanoscale oxidation and surface oxide growth on Cu doped hafnia surfaces. Our machine learned CTIP model can be used to probe

ACS Paragon Plus Environment

Chemistry of Materials 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the dynamic response of Cu/Hf/O alloy interfaces subjected to external stimuli (e.g, electric field, pressure, temperature, strain, etc.), and a variety of atomistic phenomena including the dynamics of switching in emerging neuromorphic platforms.

1. Introduction Fundamental understanding of the initial stages of alloy oxidation and subsequent nanoscale oxide growth is critical to enable groundbreaking advances in a wide range of technologies, including microelectronics, catalysis, wear-protection, anti-corrosion coatings, and neuromorphic computing.1-5 Alloying of metals is a common way of improving the properties (e.g. mechanical, chemical or physical) of materials over their single metal counterpart6. For example, alloying copper with hafnium (~1-2%) leads to exceptional improvements in strength (~120%) as compared to pure copper, especially at elevated temperatures. In most applications, alloys encounter oxidizing conditions during operation; both the oxidation process and the morphology of the formed oxide layer (composition, structure, and density) have a profound impact on the mechanical, electronic, and chemical properties of the overall heterostructure.6-7 Alloy oxidation is typically quite complex, and involves preferential oxidation of components leading to single or multiphase oxide mixtures that may be completely or partially miscible.8-10. Many past experiments have studied the formation mechanism and characterized the microstructure of such oxides. 1113

However, studies on initial alloy oxidation and nanoscale oxide growth remain scarce. 14,15-19 While oxidation and oxide growth on the surface of pure metals has been extensively studied, the

same on alloy surfaces remains less explored.20 An accurate description of the various complex interatomic interactions is required to simulate metal alloy oxidation and the growing metal-oxide heterointerfaces.21-23 In alloys, the oxide growth is complicated due to phenomena such as surface segregation and micromixing of the constituents.19 These can significantly influence the morphology and composition of the growing nanoscale oxides on the alloy surfaces. In such cases, it is highly likely that the oxidation as well as oxide growth mechanisms would be very different from pure or elemental metals.24 Large-scale classical


Page 2 of 38

Page 3 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


molecular simulations with semi-empirical potentials are clearly needed to study initial nanoscale oxide growth on alloys such as Cu-Hf. The treatment of alloy oxidation within the classical MD framework can be quite challenging.25 The success of such metal/oxide atomistic models, especially in the case of alloys, hinges on the ability of the underlying potential model to describe accurately the interactions between atoms in both the metallic and the mixed oxide regions.26 Additionally, the interatomic potential should accurately model the structure, thermodynamic, elastic, chemical, and surface properties of metals (Cu, Hf) as well as those in the various oxide (CuOx and HfOx) phases. To the best of our knowledge, there are no interatomic potential models, especially a reactive model to treat Cu/Hf/O alloy system within a single framework. Here, we introduce a machine learnt charge transfer ionic potential model for alloy metal/oxide system (Cu/Hf/O). Our machine learned model successfully captures the structure, thermodynamics and surface properties of the constituent metals i.e. Cu and Hf, their various oxide polymorphs i.e. hafnia (monoclinic, cubic, and tetragonal), copper oxide (cubic Cu2O) and their interfaces. We choose the potential formalism developed by Zhou and Wadley26 i.e. the charge transfer ionic potential (CTIP); this model dynamically captures the environment-dependent charges on the atoms. 19, 25-26 The determination of forcefield parameters for complex functional forms, for instance, CTIP is non-trivial.27 Given the large parameter space, one requires a systematic approach that can efficiently sample the available parameter space. In some of our recent work, we have demonstrated the success of evolutionary global optimization methods such as genetic algorithms (GA) for training empirical force fields.27,28-31 In this work, we use GA to perform an initial efficient sampling of the parameter landscape. Subsequently, the GA is combined with local minimization methods (Simplex) i.e. we input the populations of parameter sets obtained from GA) to obtain the final optimized CTIP alloy model. We have organized our paper as follows: the variable charge (CTIP) formalism, the training data set, and ML parameterization strategies are described in Section 2. The optimized CTIP parameters obtained using the ML algorithm, their performance and subsequent validation of the predictions of the optimized parameter set i.e. comparison of the various structural and energetic properties against data not



employed in the training set in described in Section 3. In Section 4, we employ our machine learnt CTIP parameter set to investigate the initial stages of Cu-Hf alloy oxidation for several different compositions as well as structure and dynamical evolution of their hetero-interfaces. Finally, Section 5 summarizes the key findings and conclusions of our work.

2. Computational Methods i.

Charge transfer ionic potential (CTIP) for Cu/Hf/O system

To describe the various interatomic interactions in Cu/Hf/O system, we use the modified charge transfer potential model developed by Zhou et al.25-26 CTIP comprises of two parts: ionic interactions in the oxide region are modeled using an electrostatic term (Ees), and the metallic interactions i.e. non-electrostatic contributions modeled using embedded atom method.32 𝐸𝑡𝑜𝑡 = 𝐸𝑒𝑠 + 𝐸𝑚 The CTIP potential formalism overcomes the limitations associated with fixed charge models. It allows for dynamic charge transfer i.e. charges are dependent on the local environment. Such a formalism allows modeling of metal/oxide heterostructures as well as dynamical phenomena such as oxidation and oxide growth. A brief description of the CTIP formalism and the various potential parameters is provided in the supporting information and a more detailed description can be found in the works of Zhou et al.25-26 In Tables 1-5, we list the various optimized EAM parameters obtained by using the machine learning strategy described below. ii. Training data set for parameterization of alloy/oxides For the development of accurate and transferable force fields, we require training and test datasets that sufficiently sample the potential energy landscape. To train our CTIP model for the Cu-Hf-O, we employ the following datasets for training and cross-validation tests. The training and test sets are derived from a combination of experimental (where available) and DFT calculated datasets.


Page 4 of 38

Page 5 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For the Hf metal, we train against the lattice constants, internal coordinates and cohesive energy of HCP, FCC and BCC polymorphs. The relative energetic differences between FCC, BCC and HCP is also included as part of the training. Furthermore, we also include the elastic constants of HCP (most stable polymorph) in the training dataset. The test set includes bulk moduli of the 3 polymorphs as well as the surface energies and the defect formation energies for the most stable phase i.e. HCP. Note that the structural properties for all polymorphs are from the DFT-GGA optimized POSCAR structure files33, cohesive energy and elastic constants of HCP are from experiments [34], and the energy differences between the three polymorphs are from DFT-GGA33. Bulk moduli for HCP and FCC, surface energies for HCP are from DFT-GGA33, and defect energies are from DFT-PBE [34]. The BCC bulk modulus is from Antolin et. al.35. It is based on fast free-energy calculations for unstable high-temperature phases and is consistent with experimental bulk modulus in Trampenau et. al.36.For the Cu metal, we retain the original EAM parameter set optimized by Zhou et al.

25-26

and use the structure, elastic constants, surface energies and the defect

formation energies as test sets. In the case of HfO2, our training set included structure (lattice constants, internal coordinates) and energy of the three different polymorphs i.e. monoclinic, tetragonal and cubic. We also include the relative energy between monoclinic, tetragonal and cubic phases in the training data to preserve the energetic ordering. The elastic constants of cubic was used in training whereas the bulk modulus of cubic is used as a test set. The bulk moduli of monoclinic and tetragonal is also included as part of the training data. We also use the average Hf charge for the 3 polymorphs and the oxygen and hafnium vacancy formation energies in monoclinic as part of the test set. Note that the structural, energy and elastic properties for all polymorphs are calculated using DFT-GGA and obtained from Ref. [33]. In the case of copper oxide i.e. Cu2O, our training data included the structure, energy and elastic constants of the cubic polymorph. The average Cu charge in Cu2O, the bulk modulus and the defect formation energies were included in the test set. Note that the structural and elastic properties for Cu2O are from DFT-GGA33 whereas the cohesive energy is obtained from experiments [37-38].



Once we generate the training data, the next step is to define the objective function, which is the weighted sum of squares of errors of the CTIP predictions. Within our ML workflow, our objective function  i is given as follows: Δ𝑖 =

∑𝑤 (𝑉 𝑗

𝐶𝑇𝐼𝑃 𝑗

2

) ― 𝑉𝐷𝐹𝑇 𝑗

𝑗 DFT

where 𝑉𝐶𝑇𝐼𝑃 and Vj 𝑗

are the CTIP predicted, and DFT predicted properties for any jth structure in the

training set whereas w j represents the corresponding weighting factor. The weight selection can be based either on the difference in magnitude of observables, the different categories of observables (examples include energies, lattice and elastic constants etc.), the number of quantities within a certain observable category and the order in relative importance of the observable for a given problem. In this case, we choose the weights of the various observables (lattice constants, elastic, internal coordinates, etc.) to ensure that the total weighted sum of the squared error for each property type is approximately equal. The weights chosen for the various properties are 2500/Å for the lattice constants, 5000/(eV/atom) for the cohesive energy, 10000/Å for the internal coordinates, 1000/(eV/atom) for the relative cohesive energy and 1/GPa for the elastic constants. iii. Evolutionary strategy for CTIP parameter optimization The training of binary metal oxides under the CTIP formalism involves optimization of ~ 97 parameters. Amongst these, it is worth noting that 24 EAM parameters for Cu-Cu interactions were previously optimized by Zhou et al. 25-26. The Cu-Cu interactions in this work are described with these parameters, and are listed in Table 2. We train the 73 remaining parameters by following an evolutionary approach. Note that these 73 parameters reduce to a further 58 when the appropriate continuity constraints for the splines are applied, and only one of the oxygen Qeq parameters are re-optimized. For the parameterization of the interatomic potential, we calculate the objective function (Δ). Δ is the error between the training and the fitted set. The key to a robust parameterization is the extensive


Page 6 of 38

Page 7 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


sampling of the objective function landscape during the minimization of the objective. In the case of a 58dimensional parameter set, the Δ landscape may be complex and often has multiple minima. In such cases, local optimization can prove to be insufficient here since the optimized parameters obtained vary significantly based on the choice of the initial guess. For efficiently sampling such multi-parameter space, we employ an evolutionary approach that utilizes genetic algorithm (GA) for global parameter optimization. We have shown that this approach can successfully parameterize force field parameterization for a diverse materials class and potential functional forms.27-28, 30-31 Here, our search begins with a random population of Np = 50 parameter sets. The derived (child) parameter sets are obtained using appropriate genetic operations i.e. crossover and mutation.39 We compute Δ, the weighted sum of errors in predicted properties (obtained using MD open source code LAMMPS40) as compared against the training set to assess the fitness of a particular parameter set (individual). Δ, is calculated for each individual and subsequently, Np individuals with the lowest Δ are selected for the subsequent generations. After the GA run has converged, we perform additional local Simplex optimization; we start with 15 promising GA optimized sets that have the lowest errors. We next perform validation on test sets to decide on the final optimized parameter set (see Fig. 1 for a typical schematic of the approach outlined above). The optimal CTIP parameter obtained using this workflow are given in Tables 1-5.



(b) (a)

Δ

Par am et e r

r2 met e Par a

1

(c)

Page 8 of 38

Training:  Structure • Lattice constants • Internal coordinates  Energy  Elastic constants of stable polymorph  Multiple polymorphs

First-principles DFT training and test datasets

Random N p population of m-dimensional parameter space

Test: Bulk moduli Surface energies Defect energies Oxygen vacancy formation energy  Dynamics and structure at 300 K    

Perform genetic operations

Crossover, mutation

Evaluate objective function Δi for parameter set i

No HCP Hf P63/mmc

FCC Hf Fm3#m

BCC Hf Im3#m

FCC Cu Fm3#m

Optimized parameters Yes

Good set? Retain N p parameter sets with lowest objective, Δi

Run test set validations SIMPLEX or local optimization Monoclinic HfO 2 Tetragonal HfO 2 P21/c P42/mnm

Cubic HfO 2 Fm3#m

Yes

Converged?

Cubic Cu2O Pn3#m

Figure 1: Schematic showing the machine learning workflow employed to train the CTIP potential model. (a) Sample objective space landscape for a sample 2-dimensional problem with multiple local minima. The optimized parameters are sensitive to the initial guesses and the objective can often settle into an appropriate local minimum. The multidimensional parameter space can be efficiently sampled with an evolutionary or global search algorithm. (b) A two-stage evolutionary search is utilized here to train the variable charge ionic potential for Hf/Cu/O systems. Structures, energetic ordering of multiple polymorphs, their cohesive energies as well as elastic constants of the stable polymorph for each material system are part of the training set computed from first-principles. (c) The different polymorphs of Hf, Cu, HfO2 and Cu2O used in the training and test sets.

iv. Oxidation Simulation Set-up We use the optimized CTIP model to study nanoscale oxidation of Cu-Hf alloys of various composition by employing an approach similar to that used by Sankaranarayanan et al. for treating oxidation of Ni-Al system19. The 30 Jul 2016 version of LAMMPS40 MD package is used to perform the oxidation


No

Page 9 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simulations. For each alloy composition, we study the oxidation of a slab of approximately 5.75 × 5.53 nm2 whose surface normal are oriented along the crystallographic (0001) direction. A vacuum of 5.67 nm is employed (on each side) along the surface normal (Fig. 2). Periodic boundary conditions are applied along both orthogonal directions of the surface plane. The short-range Coulomb interactions have a sufficiently long cutoff (12 Å). The long-range Ewald summations are used to calculate the Coulomb interaction along the periodic directions. The slabs are equilibrated to the starting temperature of 300 K and annealed in the 0-300 K range with 50 K steps. Thermal equilibration is performed using a Nose-Hoover thermostat41. At each temperature, we perform 5ps of MD using NVT ensemble. We ignore dynamic charge transfer between the alloys during first round of equilibration. Note that the charges for pure metals are zero and hence this approximation is valid. Next, thermalization at the desired temperature (300-600 K) is performed for 100 ps. This includes charge transfer dynamics using Qeq. A Nose-Hoover barostat is used for box relaxation in the x- and y- directions 42. O2 molecules are randomly introduced in the vacuum slab to initiate oxidation of the alloy substrates. We choose initial velocities from a Maxwell-Boltzmann distribution corresponding to the desired temperature (300-600 K). In addition, the z-direction simulation box boundaries have reflecting boundary conditions imposed. During the oxidation simulation, we maintain a constant gas pressure by inserting a new oxygen molecule only after the previous one has formed bonds with the metal. We track the growing oxide-gas front by averaging the positions of the outermost metal ions during the growth process. We insert a new oxygen molecule only after the previous one enters into the oxide/gas interface i.e. it’s z position is less than that of the oxide/gas front. We use the velocity Verlet scheme for integration of equations of motion; a 0.5 fs timestep is used and 200 ps of simulations are performed to study the nanoscale oxidation. Charge relaxation is performed using Qeq at each MD step.


Lz Lz = 5.67 nm HCP Hf (0001) surface Surface area 5.75 x 5.53 nm2

Surface area 5.69 x 5.48 nm2

Surface area 5.63 x 5.43 nm2

Vacuum

(a) 0%Cu

Lz

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 38

Oxygen insertion


(b) 5%Cu

(c) 10%Cu

Figure 2: Schematic depicting the oxidation and oxide growth simulation set-up. Hf substrate with (a) 0%, (b) 5%, and (c) 10% Cu doping, and the surrounding vacuum slab. The Hf surface exposed to gas phase is (0001). 3. Results and Discussion i.

Training of CTIP model for Cu/Hf/O system

Using the supervised machine learning framework depicted in Figure 1, the EAM + QEq (CTIP) parameters are trained for a Hf-Cu-O ternary system. Training set includes structural, thermodynamic, and elastic properties of a) HCP, BCC and FCC Hf, b) FCC Cu, c) the three hafnia polymorphs i.e. monoclinic, tetragonal and cubic, and d) cubic Cu2O. Note that we are fitting 58 model parameters to 189 data points including the 144 internal coordinates of all unit cell atoms in all the polymorphs. Each GA generation involves evolution of 58-parameter space 50-populations. We simulate for a total of 100 GA generations, which is sufficient for sampling the multi-dimensional parameter space. We take the 15 best sets obtained at the end and fine-tune the range chosen for the parameters during the next optimization rounds. We continue this till we obtained a set of parameters that achieve the target properties within a desired range


Page 11 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


i.e. have the lowest errors between the training and test set. During training, the parameters during the

evolutionary process are fed to perform MD simulations and checks are made to ensure energy conservation (and structural stability) with timestep of 0.5 fs. Note that the process is automated and the workflow automatically checks and rejects parameter sets which lead to “dynamically unstable” sets. The EAM + QEq optimal set of parameters obtained from the ML workflow are provided in Tables 1-5. Table 1. CTIP parameters for Cu/Hf/O system. All the oxygen Qeq parameters are kept fixed and all the charge parameters for Cu and Hf are optimized. Element

qmin

qmax

χ (eV)

J (eV)

ξ(Å-1)

Z(e)

O

-2.00

0.00

2.000000

13.657378

2.143957

0.000000

Hf

0.00

4.00

-3.981184

9.956430

1.035324

2.159623

Cu

0.00

4.00

-3.226831

11.401177

0.779276

1.007949

Table 2. EAM parameter set to model interaction of metal-metal. Cu parameters from Refs. [25-26]. Metal

re(Å)

fe

ρe

ρs

Hf

3.387566 2.66353

Cu

2.556162 1.554485 21.175871 21.175871 8.127620 4.334731 0.396620 F0(eV)

F1(eV)

F2(eV)

Hf

-3.412770

-0.080145

Cu

-2.19

0

B(eV)

κ

β

A (eV)

26.829200 26.829200 9.028026 4.703744 0.382736

Metal

Metal

α

F3-(eV)

η(eV)

0.912658 -1.387682

-1.387682

0.474682 -3.416801

0.561830 -2.100595

-2.100595

0.310490 -2.186568

λ

F3+(eV)

Fn0 (eV)


Fn1 (eV)

Fe(eV)

Fn2(eV)

Fn3 (eV)


Page 12 of 38

Hf

0.821277

0.296473

0.984717

-3.375530

-0.380469

1.110565 -1.884495

Cu

0.548085

0.308782

0.756515

-2.170269

-0.263788

1.088878 -0.817603

Table 3. EAM optimized parameter set for pair potentials Pair

re(Å)

α

β

A (eV)

B (eV)

κ

λ

O-Hf

2.102757

9.912055

2.114373

1.174261

1.688205

0.216589

0.561354

O-Cu

1.839583

8.601937

2.136696

0.980096

0.607583

0.301091

0.080479

O-O

3.275047

5.410041

2.557467

0.624804

0.755903

0.336468

0.710976

Table 4. EAM optimal parameters for electron density function of oxygen fe

γ

Ν

1.383394

1.568833

1.013103

Table 5. EAM parameters optimized for embedding energy spline function of oxygen i

F0,i(eV)

F1,i(eV)

0

-1.385712

1 2

F2,i(eV)

F3,i(eV)

ρe,i

ρmin,i

ρmax,i

-2.753538 -1.349940

0.017886

58.452061 0

-1.797200

-2.247050 3.308021

0.000000

68.767131 58.452061 75.008171

-2.093336

-1.164633 5.082761

0.000000

81.249210 75.

58.452061

∞

008171

ii. Performance of the Newly Developed CTIP parameters We compare the CTIP predicted properties against the DFT training and test sets for our newly developed EAM+QEq parameters and assess the performance of ML trained CTIP model. The properties of a) HCP, BCC and FCC Hf, b) FCC Cu, c) the three main hafnia polymorphs i.e. monoclinic, tetragonal and cubic, and d) cubic Cu2O are compared against experiments and ab-initio derived data in Tables 6-11.


Page 13 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Note that most of our training and test data are derived from DFT. Only the cohesive energy of HCP Hf, elastic constants of HCP Hf and the Cu2O cohesive energy are from experiments. Typical error estimates for DFT properties are as follows: Lattice parameter: 1%43, Cohesive Energy: 2%43, Elastic constants: 15%43, and Surface Energy: 10%44. Quantifying typical uncertainties for experimental data is not straightforward and can only be done on a case-by-case basis. The reported experimental cohesive energy (eV/atom), for example, for HCP Hf are as follows: -6.9945, -6.4446, -6.3147 resulting in standard deviation of ~5%. Similarly, available experimental data for the elastic constants of hafnium have a standard deviation of ~4%45, 48-49. For pure hafnium, the lattice parameters obtained using the optimized CTIP are within 3% of the experimentally reported values for HCP, FCC and BCC polymorphs (Table 6). The cohesive energy for hcp is 6.93 eV/atom which agrees well with experimental value of 6.99 eV/atom. The energetic ordering between HCP, FCC, and BCC phases predicted by our model is consistent with that from DFT. Our predicted elastic constants are within 10% of the experimental and DFT reported values (Table 6). The surface energies computed for various surface also agree well with those computed using DFT-GGA. The defect formation energies (vacancy and octahedral interstitial) also agree well with experiments and DFT predictions (Table 6). The CTIP parameters for metallic Cu were retained from the original model by Zhou et al.25-26 and the calculated lattice constants, cohesive energies, elastic constants, surface energy and bulk modulus agree very well with experiments as well as those computed from DFT. Table 7 shows the comparison of CTIP predicted properties for pure Cu with those obtained from other popular interatomic potentials and DFT.

Table 6. Performance of the EAM+QEq (CTIP) interatomic potential parameters developed in this work. We compare properties of hafnium metal derived from CTIP with experiments and DFT calculations. The polymorph is HCP unless specified. Properties Dataset Expt. [34] CTIP DFT-GGA33 DFT-PBE [34] a (Å)

Training

3.1950

3.1874


3.2050

3.1935


Page 14 of 38

b (Å)

Training

3.1950

3.1874

3.2051

3.1935

c (Å)

Training

5.0542

5.1864

5.0570

5.0533

Α

Training

90

90

90

90

Β

Training

90

90

90

90

Γ

Training

120

120

120

120

a (fcc) (Å)

Training

-

4.5010

4.4815

-

a (bcc) (Å)

Training

-

3.5786

3.5416

-

Ec (eV/atom)

Training

-6.99

-6.93

-9.95

-9.96

E (fcc-hcp)

Training

-

0.011

0.072

0.060

Training

-

0.054

0.183

0.155

C11 (GPa)

Training

190

188

184

-

C12 (GPa)

Training

75

87

70

-

C13 (GPa)

Training

66

58

68

-

C33 (GPa)

Training

204

227

194

-

C44 (GPa)

Training

60

48

52

-

C66 (GPa)

Training

58

48

57

-

B (GPa)

Test

110

113

108

115

B (fcc) (GPa)

Test

-

111

101

-

B (bcc) (GPa)

Test

112 [50,51]

106

-

-

γ(0001) (mJ/m2)

Test

-

1922

1710

1133

γ(1100) (mJ/m2)

Test

-

2325

1870

2003

γ(2110) (mJ/m2)

Test

-

2122

1840

3044

(eV/atom) E (bcc-hcp) (eV/atom)


Page 15 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Vacancy

Test

-

3.08

-

2.55

Test

-

6.25

-

6.23

(eV/defect) Interstitial (eV/defect)

Table 7. Properties of FCC copper metal – comparing experiments, DFT and various MD force fields. The CTIP prediction is equivalent to EAM by Zhou et. al. [25] since the Cu-Cu parameters are not refit. Properties Dataset Expt. [37] CTIP EAM – EAM - Mishin COMB2011 DFT Sheng et.

[37, 53]

[37]

[37, 54]

al. [52] a (Å)

Training

3.62 [38]

3.62

3.62

3.62

3.62

3.64

Ec (eV/atom)

Training

-3.54 [55]

-3.54

-3.54

-3.54

-3.54

-3.50

C11 (GPa)

Training

170 [56]

170

175

173

171

173

C12 (GPa)

Training

123 [56]

121

124

123

123

123

C44 (GPa)

Training

76 [56]

76

79

76

48

80

B (GPa)

Test

139 [56]

136

141

140

139

140

γ(100)

Test

1780 [53]

1582

1504

1345

1478

1478

Test

1780

1751

1607

1475

1519

1609

Test

1780

1421

1387

1239

1218

1294

Test

1.28 [57]

1.24

0.99

1.27

1.17

-

(mJ/m2) γ(110) (mJ/m2) γ(111) (mJ/m2) Vacancy (eV/defect)



Octahedral Interstitial

Test

2.8 - 4.2

3.27

3.17

3.06

Page 16 of 38

4.93

-

[58]

(eV/defect)

The CTIP predictions for monoclinic, tetragonal and cubic hafnia properties are given in Tables 8, 9 and 10, respectively. The lattice constants for monoclinic hafnia agree well with experiments and DFT predictions. The cohesive energies and energetic ordering are in excellent agreement with DFT (Fig. 3). The CTIP predicted near equilibrium equation of state (EOS) agrees well with the DFT predictions as shown in Fig. 3 for the various hafnia polymorphs and Cu2O. The CTIP equilibrium charges on hafnium are in good agreement with Bader charges. The oxygen vacancy formation for metal rich environment predicted by CTIP is ~ 3 eV. While this value is lower than that predicted by DFT (~9 eV), it is in excellent agreement with experimental predictions of Guha and co-workers.59 Guha and Narayanan have shown that oxygen vacancy is a dominant intrinsic point defect in nanoscale hafnium oxide dielectric films59. This is consistent with the CTIP predictions since the vacancy formation energies for oxygen are much lower than that for hafnium vacancy (~19 eV). So, the CTIP potential should be able to adequately represent the diffusional characteristics observed in the experiments. In general, the defect formation energy of monoclinic hafnia is in good qualitative agreement with DFT predictions. In the case of tetragonal hafnia, we find that the CTIP lattice constants, cohesive energies, energy ordering and cohesive energies are in excellent agreement with DFT predictions. Likewise, the lattice constants, cohesive energy and elastic constants of CTIP are in good agreement with DFT predictions. The CTIP equilibrium charge on Hf is also in good agreement with the computed Bader charge60. It is worth noting that the bulk modulus of a crystal

can be expressed as a linear combination of elastic constants Cij. In this work, we do not include the bulk modulus in the training set and hence is not a concern for over-fitting. Even if it were, it would be equivalent to increasing the weight for the elastic constant contribution in the objective and result in marginally better elastic constants at the expense of other training set data. In this work, we employed DFT-derived Cij values of Hf and HfO2 phases for training a CTIP for Cu/Hf/O ACS Paragon Plus Environment

Page 17 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


system. Our newly developed CTIP reproduces the Cij values (within ~20 %) and the bulk moduli (within a few GPa) remarkably well. Note also that when we simultaneously train against several properties of the various polymorphs of hafnia, we find that some properties do tend to deviate more significantly than the others. One of the issues might be related to the selection of weights and a multi-objective approach that independently targets the properties of each polymorph may help overcome such issues. Future studies will aim at employing multiobjective GA to train CTIP forcefields. Overall, the CTIP is able to reliably capture the structure, elastic properties and energetics of the various hafnia polymorphs.

Figure 3: Comparison of the predicted CTIP equation of state (EOS) near equilibrium with DFT predictions from Refs. [34] and [37] for HfO2 and Cu2O, respectively. The solid lines correspond to the Murnaghan fit which is used to evaluate the bulk moduli in Tables 6-11. (a) DFT predictions for three polymorphs of HfO2. (b) CTIP predictions for HfO2. (c) CTIP vs. DFT for cubic Cu2O. Here, the energies are as a function of the crystal volumes and are scaled by its equilibrium value computed at the corresponding level of theory.

Table 8. Comparison of monoclinic HfO2 properties predicted by our EAM+QEq (CTIP) interatomic potential parameters with experiments, other MD force fields and ab-initio calculations. Monoclinic (P21/c) Properties

Dataset

Expt. [34]

CTIP

COMB [34]


DFT-GGA33

DFT [34]


Page 18 of 38

a (Å)

Training

5.12 [61-63]

4.99

5.13

5.14

5.12 [64-65]

b (Å)

Training

5.17 [61-63]

5.28

5.21

5.20

5.20 [64-65]

c (Å)

Training

5.30 [61-63]

5.26

5.11

5.33

5.28 [64-65 ]

α

Training

90 [61-63]

90

90

90

90[64-65]

β

Training

99.2 [61-63]

100.3

98.8

99.7

99.7[64-65]

γ

Training

90 [61-63]

90

90

90

90[64-65]

Ec (eV/u.f.)

Training

-

-30.63

-30.89

-30.54

-30.56

E (P21/c -

Training

-

-0.30

-0.13

-0.27

-0.24

B (GPa)

Training

-

391

235

209

251

qHf (e)

Test

-

3.41

3.48

-

3.60 [65]

O vacancy

Test

~3 [59]

3.57

11.95

-

9.34 [66]

Test

-

5.04

13.44

-

9.36

Fm3m) (eV/u.f.)

defectformation energy (eV) – Metal rich O vacancy defectformation energy (eV) – Oxide rich


Page 19 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Hf vacancy

Test

-

19.12

23.53

-

16.9

Test

-

16.20

6.49

-

5.7

defectformation energy (eV) – Metal rich Hf vacancy defectformation energy (eV) – Oxide rich

Table 9. Comparison of tetragonal HfO2 properties predicted by our EAM+QEq (CTIP) potential parameters with DFT calculations. Tetragonal (P42/mnm) Properties

Dataset

CTIP

DFT [33-34]

a (Å)

Training

4.59

4.86

c (Å)

Training

3.19

3.22

Ec (eV/u.f.)

Training

-30.50

-30.47

E (P21/c -

Test

-0.17

-0.192

B (GPa)

Training

257

221

qHf (e)

Test

4.34

3.33 [65]

Fm3m) (eV/u.f.)



Page 20 of 38

Table 10. Comparison of cubic HfO2 properties predicted by our EAM+QEq (CTIP) interatomic parameters versus experiments, other MD force fields and DFT predictions. Cubic (Fm𝟑m) Properties

Dataset

Expt. [34]

CTIP

COMB [34]

DFT-GGA

DFT [34]

[33] a (Å)

Training

5.08 [61-63]

5.19

5.02

5.07

5.05 [64]

Ec (eV/u.f.)

Training

-

-30.33

-30.76

-30.27

-30.32

C11 (GPa)

Training

-

427

561

559

578 [67]

C12 (GPa)

Training

-

71

161

93

121 [67]

C44 (GPa)

Training

-

63

111

69

83 [67]

B (GPa)

Test

-

191

295

248

201-280 [65]

qHf (e)

Test

-

3.15

3.44

-

3.49 [65]

Comparison of the CTIP properties for cubic Cu2O with experiments and those computed using DFT is shown in Table 11. The lattice constant and cohesive energies of cubic Cu2O is in excellent agreement with experiments and DFT-PBE. The elastic constants and the bulk modulus predicted using CTIP are in good agreement with both experiments and DFT. Likewise, the Bulk modulus and the CTIP equilibrium charges on Cu as well as vacancy formation energies (Cu and O) are in good agreement with DFT predictions.

Table 11. Comparison of cubic Cu2O predicted by our EAM+QEq (CTIP) potential parameters with experiments, other MD force fields and DFT calculations. Properties Dataset Expt. [37] CTIP COMB2010 COMB2011 DFT-GGA [33] DFT [68, 37]

[37]

[37, 69-71 68]

a (Å)

Training

4.27 [72]

4.26

4.23-4.25

4.27

4.29

4.27-4.31

Ec (eV/u.f.)

Training

-11.3 [38]

-11.62

-11.9

-11.4

-14.4

-11.4

C11 (GPa)

Training

123 [73-74]

113

105

122

124

123-126


Page 21 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


C12 (GPa)

Training

108 [73-74]

96

89

105

105

106-107

C44 (GPa)

Training

12 [73-74]

28

71

57

8

12-15

B (GPa)

Test

112 [73-74]

103

94

111

112

112-113

qCu (e)

Test

-

0.68

0.79

-

-

0.55

O vacancy

Test

-

1.44

5.68

0.95

-

1.5575

Test

0.4576

1.30

1.34

2.49

-

0.28-1.1775

formation (eV) Cu vacancy formation (eV)

iii. Discussion The evolutionary optimization leads to identification of several local minima and we need to be careful to weed out the sets that appear to be unphysical. This is done in an automated manner by carefully monitoring the evolution of the various populations of the parameter sets and avoiding parameter sets that can lead to instability in an actual dynamical MD run. For instance, while fitting the EAM or CTIP parameters, the electronic density embedding function splines are often hard to fit and need particular attention. Treating all spline parameters as independent may result in the optimization procedure sampling discontinuous splines, such as shown in Fig. S1b. This unnecessarily slows down the convergence. An ideal shape for the metal spline is depicted in Fig. S1a, which plots the splines parameterized by Zhou et al.77 For an efficient search of the parameter space in our global optimization workflow, it is important to impose spline constraints and not treat all parameters as independent. The constraints are chosen so as to avoid discontinuities (see supporting information Fig. S1b). The constraints are chosen to ensure splines are continuous in the 0th, 1st and 2nd order derivatives. Nevertheless, some of the optimal parameters found have splines that have a local maximum as shown in Fig. S1c. Such a parameter set is seen to display



dynamic structural instability when an MD run is performed at 300K. Hence, it is important to test for structural stability at 300K while selecting the final parameter set. Fig. S1d shows the final spline parameters. Coupled with testing for the structural stability in a molecular dynamics simulation and looking for unphysical natures of the interatomic potential, it is possible to weed out incorrect local minima such as set 2 (see Table S1 and supporting information for details). The complex objective function landscape for CTIP, with multiple local minima, has no obvious correlation between parameters (see Fig. S3). The evolutionary algorithm (followed by local optimization) identifies several potential local minima. For a comprehensive sampling of the objective function landscape we start with searching in a broad range for all the parameters. 15 best sets chosen from the optimum of 100 generations (2/15 sets from an intermediate optimization round are shown in Table S2 along with the final published set). Based on the standard deviation of each of the parameters (see Table S2) the ranges to sample for the next round of GA are found. This process is continued till a final search over a narrow range of parameters. A discussion on the comparison on some of the best parameter sets obtained from our optimization procedure is provided in the supplementary information. It is worth noting that the ML workflow is able to come up with a CTIP parameter set that does not require the charge bounds imposed in the original formalism of Zhou et al. For a Hf2Cu alloy test case, we find that the Bader charges for Cu is -1.0841 whereas that for Hf is 0.54205. The CTIP predicts charges of Cu ~-0.34 and Hf~0.17, which is qualitative agreement with DFT calculations. The minimized lattice structures for the Hf2Cu is in excellent agreement with DFT (within 1% error). Note that this alloy was not a part of the training set and the original CTIP potential would destabilize the lattice. Hence, the newly optimized potentials do not require the charge bounds to be imposed which is a significant improvement over the original formalism. Also, note that this work retains existing EAM database for the metals and uses them as a starting point for developing CTIP potentials for the corresponding oxides. In an earlier work, we had developed a CTIP potential for Ta/TaOx system78. The O-O interaction parameters from this earlier work could not be transferred to the Cu/Hf/O system and was found to perform poorly indicating that O-O interaction


Page 22 of 38

Page 23 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


parameters are not readily transferable from one oxide to another. Ideally, it would great to have a universal set of O-O parameters that could be used for all oxides. To be able to do that, one needs to generate a more elaborate training data set that has several multi-component oxides (oxides of Cu, Pt, Pd, Ni, Ta, Hf, Zr, Al etc) and then train a potential against this dataset. Such an approach would generate a universal O-O parameter set that could be used across different oxides. Finally, we note that all the three aspects (training data, model selection and optimization scheme) are critical to develop a transferable and well-formulated potential. There is clearly a need to move from modest and specific fitting data sets to reproduce specific properties that are strictly valid at equilibrium to fitting vast and diverse datasets obtained from first principles and experiments. In particular, to ensure transferability to other oxides and material classes, it is important to have a more elaborate training data set that encompasses features from several different metals (fcc, bcc, hcp) and their oxides. Going forward, it is also preferable parameterization based on quantities that could be measured experimentally (e.g. oxygen vacancy formation energy). MD simulations that employ pre-defined functional forms impose limitations on the extent of physics and chemistry that can be captured by the underlying model. There is a clear need for more flexibility in the functional form to address this challenge. For instance, the work by Gibson et al. is a step in this direction79. Although potential functions with more complex functional form such as ReaxFF and MEAM impart more flexibility compared to simpler forms such as EAM, they will still be limited in their ability to capture the underlying physics and chemistry.80-81 Our work demonstrates that data-driven approaches that take advantage of more elaborate training dataset and advanced numerical optimization can lead to significant improvement in the quality of the potential models. There will, however, be a ceiling limit to the performance that can be achieved using pre-defined functional forms. The existing force-fields that rely on such pre-defined functional forms have limited flexibility. For instance, one can easily transfer parameters from one material class to another (metal to oxide, for example). There is thus a necessity to incorporate more flexibility in the functional form of the potential models that are employed in current MD simulations. One way to circumvent the problem of



utilizing pre-defined functional forms by defining the potential function in terms of a superposition of a linearly-independent basis set of functions:

𝐸𝑡𝑜𝑡𝑎𝑙 =

∑∝ ∑ 𝑓

𝑚(𝑟𝑖𝑗)

𝑚

𝑚

< 𝑖𝑗 >

+

∑𝛽 ∑ 𝑔 (𝑟 𝑛

𝑛

𝑛

< 𝑖𝑗𝑘 >

𝑖𝑗𝑘)

+ ….

where Etotal is the total system energy of

atoms, fm and gn are over-complete sets of basis functions (e.g. Bessel functions, Hermite polynomials, exponentials, trigonometric functions, or combinations thereof), rij and rijk are rotationally-invariant combinations of the atomic coordinates ri, and ’s and ’s are expansion coefficients to be determined by fitting. By defining an over-complete basis set, a sparse (“sparse” in this context means that most of the ’s and ’s are expected to be zero) regression technique rather than Least Squares is needed. We believe that such flexible functional forms combined with more elaborate, well sampled training dataset and advanced optimization schemes are required to develop the next generation of transferable interatomic potentials.

iv. Case study – Oxidation of Cu-Hf alloy systems To exemplify the suitability of our machine learned model, we perform MD of oxidation and nanoscale oxide growth on the surface of pure hafnium and Cu-doped hafnium alloys. We investigate the oxidation behavior at low temperatures (300-600 K) and different concentration (

Machine Learning Applied to a Variable Charge Atomistic Model for

Recommend Documents