Phonon optimized interatomic potential for aluminum

We address the problem of generating a phonon optimized interatomic potential (POP) for aluminum. The POP methodology, which has already been shown to work for semiconductors such as silicon and germanium, uses an evolutionary strategy based on a genetic algorithm (GA) to optimize the free parameters in an empirical interatomic potential (EIP). For aluminum, we used the Vashishta functional form. The training data set was generated ab initio, consisting of forces, energy vs. volume, stresses, and harmonic and cubic force constants obtained from density functional theory (DFT) calculations. Existing potentials for aluminum, such as the embedded atom method (EAM) and charge-optimized many-body (COMB3) potential, show larger errors when the EIP forces are compared with those predicted by DFT, and thus they are not particularly well suited for reproducing phonon properties. Using a comprehensive Vashishta functional form, which involves short and long-ranged interactions, as well as three-body terms, we were ...


I. INTRODUCTION
Aluminum is one of the most abundant and most useful metals in the earth's crust.Micronsized aluminum powder has been routinely used in combustion, propulsion, and energy generation applications, due to its high energy density, low cost of extraction, and relative safety. 1 The past decade, however, has marked a shift toward the use of nano-sized aluminum particles, which have large specific surface area and therefore higher reaction rates. 2,3Unfortunately, the formation of aluminum oxide on the particle surface reduces the rate of heat conduction between the heat release zone and the unburned reactants, which is a key process governing combustion wave propagation. 4or a spherical particle, for instance, the fraction of oxide increases from 6 to 47% as the size of a Corresponding Author Contact Information: Email: ase@gatech.edu,Phone: (404) 894-7514 2158-3226/2017/7(12)/125022/11 7, 125022-1 © Author(s) 2017 the particle decreases from 30 to 3 nm. 5Therefore, the overall burning rates of nano-aluminum powders could be substantially diminished by the oxide surface, through reduced heat conduction rates.
7][8] The thermal conductivity of aluminum oxide, on the other hand, is phonon dominated and is much lower than that of aluminum.It could, hence, become a predominant resistance to heat conduction as the combustion process proceeds towards completion.0][11] Furthermore, not only the two solid phases, but also the interface between them probably plays a crucial role in heat conduction, and the interface may possibly be even more significant than the layers, in terms of resistance.
A series resistance circuit model can be used to estimate the total resistance, R of a nano-aluminum particle with oxide layer, and the contribution of each conductive resistance given as in Eq. 1.
As a rough estimate based on phonon MFP spectra calculated for Al 2 O 3 and the electron MFPs in Al, 12 we assume, k Al = 100 W/m-K; k Al2O3 = 5 W/m-K for a 10 nm aluminum particle with a 3 nm oxidation layer.Lyeo et al. 13 report results for the temperature dependence of interfacial conductance of Al/Al 2 O 3 measured using the time-domain thermoreflectance (TDTR) method.They obtained a room temperature interfacial conductance of ∼200 MW m -2 K -1 .Using this data, the interfacial resistance R int is 5x10 -9 m 2 K/W, and our simplified estimate suggests that the interfacial conductance is likely to dictate ∼ 90% of the total resistance.We therefore expect that the interface will be dominant, the contribution from the oxide layer is < 10%, and the Al comprises < 1%.It is theorized, and there is some evidence that suggests, that electrons also make a negligible contribution to the interfacial resistance, 14 and therefore the phonons are most critical.Classical molecular dynamics (MD) can be used to study phonons and predict both interfacial conductance 15 and phonon thermal conductivity. 16Ab initio methods such as density functional theory (DFT) 17 may also be used, but they are limited by the system size.They also do not decompose the force on atoms in terms of contributions from other atoms; that is, only the total force on atoms results, since the charge density is not easily partitioned into individual atom contributions.Formalisms that provide phonon mode-level information on thermal conductivity 18 and interfacial conductance 19 require force decomposition, so it is not clear how such formalisms could be used in conjunction with DFT.
One promising approach to studying phonon physics is to develop a phonon optimized interatomic potential (POP) that can be used in MD simulations. 20The most widely used potential for Al is the embedded atom method (EAM); 21 the charge optimized many body (COMB3) potential has also recently been used to model Al. 22,23However, these potentials were not specifically designed or optimized to describe phonons, and it is not known how well these potentials will reproduce ab initio MD trajectories and subsequent phonon-related properties.In theory, the complexity in common potential forms should be capable of capturing phonons accurately, but to better quantify the extent to which the existing potentials (EAM and COMB3) can describe phonons, we compare their predicted forces with forces obtained from DFT results for supercells, where the atoms have been randomly displaced from their equilibrium lattice sites.As a quick test, we performed self-consistent field (SCF) calculations on Al (full details given in the methodology section), and the total force errors were compared with those of predictions by EAM 21 and COMB3, 22 using a FCC supercell of Al consisting of 32 atoms.The crystal was first relaxed, then random displacements (maximum = 0.01 Å and RMS = 0.0062 Å) were given to each atom and forces were evaluated.Comparing the forces from the calculations predicted by EAM and COMB3 potentials with the DFT results, the average error in forces for EAM was ∼ 52% and ∼ 48% for COMB3.Accurate force predictions are important for predicting phonon properties, but it is not well understood what value of force error can be tolerated and still yield reasonable phonon properties.Prior works 20,24 on the development of POP for crystalline silicon and germanium suggest ∼10% force error is satisfactory.
In this work, we present a POP for Al based on a training set consisting entirely of DFT calculations.Initial work by Rohskopf et al. 20 showed that the POP methodology can work for semiconductors such as Si and Ge, and suggested 24 that the method is general and can be extended to other materials/classes.Here we apply the POP methodology to a metal (Al) and seek to determine if the POP approach can properly reproduce the phonon contributions to thermal conductivity for Al.Training datasets to sample the phase space were generated from DFT calculations, and consist of five key properties, namely total forces and stresses, energy vs. volume, and harmonic and cubic force constants.As for Si and Ge, a genetic algorithm 25,26 based optimization routine was used to search the space of fitting parameters.The Vashishta 27 functional form was chosen to describe the potential because it consists of short and long-range interactions, and three-body terms to capture phonon interactions.Since the functional form is sufficiently flexible to describe both Al and Al 2 O 3 , in a future study the interfacial conductance between the two phases can be studied in detail.

A. The Vashishta potential functional form
The Vashishta potential has previously been used in MD simulations for thermal transport property evaluations; 28 it was originally developed for silicon carbide 27 and has been extended to Al 2 O 3 . 29he potential is fairly comprehensive, and may be parameterized for describing the Al/Al 2 O 3 interface in future, with the goal of describing each bulk material as well as the interface, all with a single functional form -albeit with different parameters for Al interactions in each phase (i.e. an Al atom in metallic Al will likely be treated as a different species from Al in Al 2 O 3 ).The Vashishta functional form expresses interaction potential energy consisting of two and three-body terms: The two-body term includes steric-size effects, Coulomb interactions, charge-induced dipole, and van der Waals interactions, and can be written as: It is truncated at r c and the shifted potential can be described as: The three-body term is a product of spatial and angular dependent components, In this study, among the 13 independent parameters, 12 parameters H, η, χ, D, ξ, W, r c , B, γ, r 0 , C, and cos θ jik need to be optimized.The charges Z i can be set to zero in the aluminum crystal, since it is charge neutral.

B. Genetic algorithm based optimization
With 12 independent fitting parameters, the corresponding 12 dimensional optimization problem could easily fail if local minima-based optimization strategies are used.Therefore, a genetic algorithm (GA) based optimization technique is used to obtain the parameters.A GA, with its metaheuristic approach, mimics natural selection, providing a more complete and efficient exploration of the multi-dimensional parameter space. 25,30,31The basic approach of a GA is to evolve a potential set of randomly generated candidate solutions (termed "individuals") towards better solutions via multiple evolutionary stages (termed "generations") by minimizing an objective function.The objective function, Z is a measure of the error between the training set values for each configuration and the values for a candidate set of fitting parameters, and is given as: where, z represents generalized errors, and w represents the weights.The subscripts f, e, s, ifc2 and ifc3 represent force, energy, stress, harmonic force constant, and cubic force constants.The errors are given by Rohskopf et al. 20 as: where F DFT ij is the ab initio Hellman-Feynman force on atom i in configuration j and F EIP ij is the force predicted by the potential.The error is summed over all N atoms and M configurations, and is normalized by the number of forces and sum of ab initio forces.The error for energy z e is given by where E DFT j is the ab initio energy and E EIP j is the energy predicted by the fitted potential.Energy is referenced to the energy of the fully relaxed equilibrium configuration.Emphasis is placed on reproduction of the shape of the energy vs. volume curve, and therefore, the energy is described as an energy difference E j = ε j ε coh , where ε j is the total potential energy and ε coh is the total potential energy of the DFT equilibrium configuration with the equilibrium cell dimensions.With this definition for energy, the energy of the relaxed/minimized structure will always be zero.The error for supercell stress z s is then given by, (11)   where σ DFT jk and σ EIP jk are the k th symmetric ab initio and EIP stress tensor components of the j th configuration, respectively.The sum over k stress tensor components considers the six symmetric components σ xx , σ yy , σ zz , σ xy , σ xz , and σ yz .Both of these quantities are important for stability and proper dynamics, which depend strictly on the shape of the energy-volume curve.Fitting of IFCs is accomplished by addition of another error component, z o ifc , where o is the order of the IFC being fit.The IFC error component is given by: where P is the number of IFCs and φ DFT n and φ EIP n are the i th order ab initio and EIP IFCs, respectively.The training set consists of the properties being optimized, which may be obtained from experimental or ab initio results.In the present work, only ab initio DFT results are used to generate the training set, as described in Section III C.After each generation, the individuals with the least error (the lowest objective function value), termed the elite population, are stochastically selected and evolved to the next generation after experiencing some 'mutation' and 'crossover'. 25In a GA, a set of EIP parameters is termed and individual and individuals are represented as binary strings.When represented as a single binary string, one or more bits can be flipped at random to produce a random "mutation".'Crossover' refers to the mixing of multiple individuals to obtain a new 'child' solution.
Here we start with 1200 individuals (that is, 1200 randomly generated parameter sets with 12 parameters each).A group of 1200 individuals is termed a trial, and a trial is allowed to evolve for 250 generations with an elite population of 10%.The elite population is the group of individuals that exhibited the lowest value for the objective function in that generation.The rest of the population is randomly selected and grouped with the elites so that 50% of the total population is mixed with these better performing parameter sets, and they experience crossover and mutation.For a particular generation, this 50% of solutions are termed the "parents".The parents first undergo a crossover operation, where a single random crossover point in the string is chosen for two parent solutions.The binary strings that represent the decimal parameter values of a parent solution are split at the crossover point, and all the strings beyond this point are swapped with corresponding strings from the other parent.This recombination results in a "child" solution that then has a randomly chosen mixture of the two parent binary representations.The next step is to apply the mutation operation, where the binary strings representing decimal values of each "child" solution experience random perturbations from "0" to "1" or vice versa.The mutation operation randomly perturbs existing parameters and it allows the optimization process to explore the full parameter space, reducing the chance of getting stuck in a local minimum.This means that only the top 10% of each generation survive and the remaining 90% of the individuals in the next generation are generated through the aforementioned process.
The objective function value, termed the "fitness", of each parameter set is then evaluated from the error obtained with respect to each fitted quantity in the DFT training set.The properties, such as individual atom forces and super cell stresses, are calculated using the LAMMPS 32 open source code as a shared library with the POP code. 20The respective weights for the objective function components are given in Table I.The best individuals or the elites in each generation, are then used to seed the next generation.In the end a single best solution is not necessarily obtained, but rather many unique solutions with nearly degenerate objective function values are obtained.Here, the best 200 parameter sets were selected for statistical analysis and the resulting parameters cover the ranges given in Table II, which is also reported in the supplementary material.

C. Training data set
Appropriate training datasets that sufficiently sample the relevant phase space were generated from DFT calculations.The DFT calculations were performed employing the local density approximation (LDA) for the exchange-correlation functional, and the normconserving Perdew-Zunger 33 scalar relativistic pseudopotential was used to describe the core electrons.The projector-augmented wave formalism was implemented in QUANTUM ESPRESSO (QE) 34 with a plane wave energy cutoff of 750 eV.The Brillouin zone was sampled using 4x4x4 uniform (Monkhorst-Pack) k-point grids.For the self-consistent field (SCF) calculations, the electron energy convergence threshold was set 1x10 -7 eV, and for the initial structural optimization, the force/atom convergence threshold was set to 1 x10 -4 eV/Å.
To further validate the accuracy of the DFT calculations before proceeding to phonon optimization, we used the open-source ALAMODE code, 35 which is useful for calculating the harmonic and anharmonic force constants from DFT data, and also facilitates the calculation of dispersion relations and the phonon DOS.The phonon DOS was calculated using a 4x4x4 mesh resolution, and dispersion relations were extracted using a fine 150 one-dimensional grid points along each phonon branch.The phonon BTE is solved under a relaxation time approximation (RTA) using a 30x30x30 grid, which produces converged phonon thermal conductivity values that can be compared with our results and with experimental data.
A total of 100 configurations (26 for force constant fitting + 74 for other properties) were used in the training set.Note that the number of configurations is significantly lower than those required for neural network potentials (NNP), Gaussian approximated potentials (GAP), or the spectral neighbor analysis potential (SNAP).In addition, comparing the speeds of these potentials when implemented in LAMMPS, GAP and SNAP are 10 times slower whereas NNP is ∼ 4 times slower than the much simpler Vashishta potential.These are significant computational savings by comparison to more complex potentials, which is one of the key advantages of the POP approach, where the goal is explicitly to study phonons rather than attempting to describe the entire phase/configuration space.

A. Parameter optimization
Table II shows the ranges allowed for the fitting parameters in the Vashishta functional form, where 12 parameters were optimized.The atomic charges in the Coulomb term of the potential were assigned a value of zero, since Al is monatomic.The ranges explored for each parameter are given in the second column and the range of values for each parameter that appear in the best 200 solutions are given in the third column.Initial inputs were chosen over a much broader range, which was manually narrowed down to reduce the search space once several solutions emerged.Column 4 gives the fraction of input parameter range comprised by the solutions.Since there are many solutions, it would be cumbersome to represent all of them in tables as is common practice. 20,22Instead, we have provided LAMMPS input files that refer to specific solutions by a reference index (Vashishta-1, 2, etc.) in the supplementary material.Among these sets of solutions, some may be transferable to describe other properties, such as surfaces, defects, or alloys, although transferability is neither guaranteed, nor an objective of the approach.
The statistics give a sense of how much diversity there is in the solutions.For example, if a parameter has a narrow range, then in all the solutions this parameter was essentially fixed at a nearly unique optimal value.On the other hand, if a parameter varies drastically, the solutions span a wide range and are not degenerate in representation, despite the fact that they all reproduce the forces well and have nearly degenerate objective function values.
Figure 1 shows the error in forces on atoms.The minimum error is ∼32%, and the maximum is ∼45%, excluding a few outliers.While Rohskopf et al. 20 suggest that force error of ∼10% is to be aimed for, it is important to investigate whether the solutions ensure a fair prediction of phonon properties.Therefore, in the following sections, the accuracy of the fitted potential with the lowest force error (i.e., Vashishta-56) is compared against the results of EAM, COMB3, and DFT results.

B. Thermal conductivity accumulation plots
In order to test the accuracy of the fitted potential, it is important to evaluate its efficacy in calculating phonon thermal conductivity.Based on the Wiedemann-Franz law, the thermal conductivity of Al is dominated by electrons (>90%), and it is also quite high (> 200 W/m-K).Since it is dominated by electrons, this puts a large degree of uncertainty on the phonon contribution, if it is to be inferred from experiments.Phonon thermal conductivity can therefore be expected to be less than 20 W/m-K.Here, we calculate the phonon thermal conductivity by solving the phonon Boltzmann transport equation (BTE) under relaxation time approximation (RTA).Solution of BTE was performed using ALAMODE. 35An atomic displacement of 0.01 Å was used for perturbing the atoms from equilibrium positions.Prior researchers 16,36 have also made similar attempts to calculate thermal conductivity of Al, but do not report the magnitude of atomic displacements used.Jain and McGaughey 16 report a thermal conductivity of ∼8 W/m-K at 300 K, whereas Wang et al. 36 predict the phonon thermal conductivity to be ∼4.8W/m-K.Our results, however, indicate a thermal conductivity of ∼ 1 W/m-K at 300 K.Although there is disparity in these three computed values, additional calculations showed that this is likely due to the difference in atomic displacements used in the DFT calculations.Furthermore, because of the large degree of uncertainty associated with the phonon contributions, all the values are arguably consistent with the Wiedemann-Franz Law, since they are all < 20 W/m-K.
The cumulative thermal conductivity (thermal conductivity vs. phonon mean free path) or the accumulation plot is obtained, as shown in Fig. 2 with DFT results.While COMB3 over-predicts the results by about 140% and predictions from the EAM potential are about 40% lower than the DFT results.These potentials are unable to reproduce the thermal conductivity accumulation with respect to the phonon mean free paths, which indicates that they may not be suited for describing phonons.However, the Vashista-56 potential exhibits good agreement in both the total value of thermal conductivity and the accumulation function.

C. Phonon density of states
In addition to the accumulation plots, phonon densities of states (DOS) were also calculated using the different potentials and comparing with the DFT results, as shown in Fig. 3.As can be seen from the figure, phonon DOS predicted by our potential has a better agreement with the DFT results in terms of the proximities of the peaks.The EAM and COMB3 potentials exhibit good qualitative agreement, but the actual peaks in frequencies are off by more than 50% in both height and location.

D. Phonon dispersion relations
Phonon dispersion relations represent the phonon branches.It has been shown 22 that COMB3 and EAM predict dispersion relations well.Figure 4 shows that the Vashishta-56 potential is also quite satisfactory when compared with DFT and experimental 37 results.In theory, an improved functional form such as a Taylor expansion may offer sufficient flexibility to capture all the necessary force constants, particularly for the non-first nearest neighbors.It should also be noted that DFT results can FIG.4. Phonon dispersion relations calculated using Vashishta-56 (this work), EAM, and COMB3 compared with DFT and experimental results. 37ry marginally with the pseudopotential used to describe core electrons, but generally the quantitative and qualitative trends have been captured with high levels of accuracy.

E. Green-Kubo thermal conductivity calculations
Having shown that the Vashishta-56 potential captures the basic phonon properties from the various ALAMODE calculations, it is important to test its applicability in MD simulations.Using Vashishta-56, molecular dynamics simulations were performed using LAMMPS, 32 to test predictions of thermal conductivity (k), using the Green-Kubo relation given by: where k B is the Boltzmann constant and T is temperature.The heat current vector S was derived by Hardy. 38 where V is the volume, p is momentum, m is the atomic mass, U is the potential energy, v is the velocity, and r is the position of atoms represented by i and j.To test for stability, we tried different ensembles and thermostats to make sure the crystal preserved its lattice stability throughout the calculation without any significant deviations in energy or structure.For the thermal conductivity calculation, as a first step, the energy of a FCC crystal with 256 atoms was minimized by the conjugate-gradient (CG) method with a relative tolerance of 1x10 -4 and 1x10 -6 for energy and forces respectively, and the system was then equilibrated in the NPT ensemble for 10 ps.A time-step size of 0.1 fs was used and the Green-Kubo calculations were performed for another 100 ps with a correlation length of 1 ps.The crystal size of 256 atoms and correlation time of ∼ 1 ps were determined to be sufficient after testing for convergence.The temperature dependency of lattice thermal conductivity of the FCC Al crystal was compared with the results from the solution of the phonon Boltzmann transport equation.Results of the calculations are presented in Fig. 5.For each temperature, ∼35 independent results were averaged; the error bars correspond to the standard deviation amongst the independent results.As the figure shows, there is again a good agreement with the BTE results in the temperature range considered.The maximum deviation, ∼10%, is observed at 400 K.For all the other temperatures, the deviation is well within our goal of 10%.This suggests that POP Vashishta-56 can be used to simulate the phonons in Al, with fidelity.

IV. CONCLUSIONS
In this work, phonon-optimized potential (POP) for aluminum has been developed, and this is the first application of this method to a metal.The POP methodology was previously validated for silicon and germanium, and now has been extended to aluminum via the Vashishta functional form, resulting in approximately 200 solutions of interest.The functional form and the parameter sets can be readily used in MD simulations using LAMMPS software.The phonon properties predicted by Vashishta-56 align well with the results from DFT calculations, and the force error is ∼32%, which is somewhat lower than previous EAM and COMB3 potentials.Nonetheless, despite the seemingly large force error, the Vashishta-56 potential is able to accurately reproduce the phonon-phonon interactions, as indicated by the excellent agreement between the thermal conductivity accumulation plots in Fig. 2. When using MD, the thermal conductivity calculated by the Green-Kubo method is within 10% of the BTE-RTA results.In the future, this same approach can be extended to aluminum oxide (Al 2 O 3 ) and to the interface between Al and Al 2 O 3 , which is of interest to the combustion community.
FIG.1.Error (%) in forces on atoms for the various sets (trials) of Vashishta potential parameters.

TABLE I .
Fitting properties and weights.

TABLE II .
Input and output parameter ranges, fraction of search space obtained as output, and the specific values obtained for the Vashishta-56 parameter set.