New potentials for multiscale simulations of liquid metals

We present the technique and results for finding norm-conserving pseudopotentials and EAM potentials that can be used to recover atomic and electronic structure of liquid iron at and above the melting point. Pseudopotentials were found by minimizing the energy differences of our results with all-electron reference methods; EAM potentials — by the modified hybridization method proposed earlier by Belashchenko. We show that these potentials are at least as accurate in describing liquid iron as the established potentials in the field.


Introduction
In the last decade multiscale computer simulations have transitioned from the category of curiosities to the toolbox of practicing material scientists. Increasing computational resources and no need for setting of expensive experiments make such simulations an essential step in the development of new materials. Multiscale modelling is typically understood [1] as the modelling from the smallest scale (i.e., electrons and atoms) to the full system level (e.g., autos). Every step in this multiscale hierarchy requires its own simulation method -from quantum mechanics to continuum modelling -which depends on a set of parameters obtainable from previous steps. Each method, in its turn, has its own assumptions and limitations, that can greatly influence the obtained results.
The first two steps in the multiscale modelling hierarchy usually comprise ab initio electronic simulations, done with methods based on the density functional theory (DFT) [2], and atomic scale simulations with classical molecular dynamics. Ab initio simulation methods are based on the solution of Schrödinger equation for a system of particles. The DFT methods suggest searching for the solution of the equation in the form of one-electron wavefunction, which can be expanded in terms of some basis -typically either atomic orbitals or plane waves are used. If the electron density is close to the atomic nucleus, large gradients occur in it, that can be described properly only by enlarging the expansion basis substantially. For such electrons pseudopotentials are used, that replace the Coulomb repulsion usually present in Schödinger equation. The choice of pseudopotentials and basis sets are two main sources of errors within the DFT framework. Classical molecular dynamics, in its turn, relies heavily on the choice of interatomic potentials, which are usually constructed largely basing on the author's intuition, as the choice of the fitting parameters is very wide. We tried to tackle the problem sequentially: first, we estimated the error of ab initio calculations of liquid iron; then, we tried to minimize it; and finally on the basis of ab initio data we built the series of EAM potentials and chose the most adequate one. In this paper we present the findings.

Computational details
A prerequisite for the error analysis is the repeatability of results -different software implementations of the same DFT formalism must give identical predictions for identical systems. Recently Lejaeghere et al [3] developed a technique to consistently determine the Δ-factorthe single number that describes the difference in the predictions of different methods. The technique was applied to a variety of DFT methods, however, the majority of them use plane waves as a basis. We used the technique to find Δ-factor for SIESTA software [4], that uses atomic orbitals, and wrote the PseudoGenerator code [5] to minimize Δ-factor in comparison with de-facto standard Wien2K and VASP software packages.
Later, atomic and electronic structure of crystalline and liquid iron was found by ab initio molecular dynamics (AIMD). AIMD simulations were made with SIESTA software package for supercells of 200 atoms of iron at 4 temperatures from the range of liquid phase existence: 1833K, 1873K, 1923K and 2023K. Siesta implements Car-Parrinello molecular dynamics in the basis of numerical atomic orbitals, with core electrons represented by nonlocal norm-conserving pseudopotentials [6]. Standard DZP basis was employed, with GGA as the approximation of exchange and correlation. As local magnetic moments on iron atoms exist also above Curie temperature, we employed spin polarization for all AIMD runs. Previously we performed the study of influence of local magnetic moments on the structure of liquid iron [7] and found that it is necessary to properly account for magnetic effects to get the structure RDF close to experimental. The same conclusion was made by Ganesh and Widom in [8].
Based on ab initio atomic structure, we constructed a set of EAM potentials using a modification of the method developed by Belashchenko in [9] for recovering an EAM potential from an RDF. The details of algorithm will be published elsewhere; here we will give only an outline.
To reconstruct EAM potentials, first we needed to establish a set of hybrid effective pair potentials (EPPs). It is known that one can approximately represent an EAM potential as the effective pair potential corresponding to the average electron density in a system under fixed thermodynamic conditions: Here φ EAM (r) is the effective pair representation of EAM potential, and ρ is the mean electron density in the system at the given conditions. The key point of the Belashchenko method is the choice of the pair part recovered from RDF as φ(r). He proved, that if some pair potentials reconstruct the RDF for the liquid, then the linear combination of the found potentials also reconstructs the target RDF: given that φ i (r) here are the basis effective pair potentials, and λ i are arbitrary real numbers. We employed a method by Schommers [10], which is used to generate pair potentials from a known RDF, to obtain three linearly independent pair potentials. The potentials are exploited later to get a set of hybrid pair potentials from Eq. (2), which differ by the choice of λ i . Any set of λ i coefficients from the known range (range bounds are determined by the thermodynamical stability of the system) does not distort significantly the model RDF curve in NVT ensemble. λ i can be selected to fit several physical properties of the modelled system. Namely, Figure 1. Δ-factor for SIESTA software package compared to all-electron calculations. The coloured symbols distinguish the elements from different periods of the Periodic Table (see legend for correspondence between symbols and period numbers); atomic radii are also shown (empty star symbols). we use one of these coefficients to optimize hybrid EPPs to experimental atomic density in NPT ensemble. With condition (3), we have one free parameter λ left. In this paper, we used this parameter to fit the potential to a number of physical properties of the system, both experimentally measured and found by ab initio simulations. With such a pair part, EAM potential also reproduces target RDF, provided the derivative of embedding function in (1) is equal to or sufficiently close to zero. Embedding function was computed from the universal equation of state for metals by Rose et al [11]. Effective electron density ψ(r) was chosen in the form of decaying exponent.

Norm-conserving pseudopotential for iron
Employing the technique described in [3], we calculated Δ-factors for more than 50 different elements with SIESTA software package, using default pseudopotentials (taken from [12]); the results are shown in Fig. 1 along with atomic radii. It can be seen, that the difference between SIESTA and all-electron calculation results grows with the period number. Apart from this, for each period Δ-factor is larger for the elements with the smallest atomic radius, i.e. situated in the middle of the period. Two notable outliers are iron and platinum.
It can be seen from the table, that the change of pseudopotential improves both the lattice constant and the magnetic moment of iron; the resulting magnetic moment appearing even closer to the experimental data than that from Wien2K and VASP calculations. To check if the pseudopotential can be applied to the calculations of structure in the liquid state, we have simulated the structure of liquid iron, which is shown in Fig. 2 (along with experimental data and RDFs made from classical MD calculations). It can be seen, that the structure recovered with the new pseudopotential is at least as close to the X-ray diffraction data as the old one.

EAM potential for iron
To check the influence of EAM potential on the structure of the liquid we have calculated RDFs of liquid iron for potentials corresponding to several λ parameters ranging from -2 to 3. The results are presented in Fig. 2. We can see, that the hybridization principle standsthe differences in RDFs are very subtle, that is, different potentials resulted in the same liquid structure.
To determine how model properties depend on a hybrid EAM potential choice, we calculated several iron properties both in liquid state and during melting, versus hybridization parameter λ. We examined the following liquid iron characteristics at T = 1833K: selfdiffusion coefficient D, dynamic viscosity coefficient η and volume thermal-expansion coefficient β. Melting temperature T melt , enthalpy ΔH melt and changing the volume at melting ΔV melt have been studied during melting process. Table 2 presents all mentioned characteristics, which were calculated for the potentials with λ = 1 (the best fit to properties), compared to the potentials from λ range boundaries (λ = −2 and λ = 3). The reference data are also presented for comparison. As a reference we used experimental and ab initio data; besides, we compared data obtained with our potentials to that from EAM potentials for iron reported by other authors. We used several EAM potentials: potentials No.2 and No.5 by Mendelev et al [17] (below referenced to as M2 and M5), potential by Belashchenko [9] (referenced to as B) and potential reported by Chamati et al [18] (referenced to as C).
When comparing the results for different values of λ, given in Tab. 2, we can see that every physical property is dependent on λ; however, different properties depend on the parameter to a different extent. The most sensitive property appears to be the volume change at melting ΔV melt , which changes by 237% across the λ interval; the least sensitive property is T melt , with 9% change. We optimized λ against the set of properties using least squares, and found λ = 1 to be the parameter that leads to the most optimal potential. Comparing the results for λ = 1 with reference data, we can find that the relative deviation from the target values for D and η is comparatively small (D -5%, η -9%). The value of β, on the contrary, deviates greatly from the experimental value, with relative error of 59% with EAM and 29% with the pair potential. As for the characteristics at melting point, the calculated values of T melt and ΔH melt are very close to reference for EAM with λ = 1.0; relative errors are just 0.2% and 1.2% respectively. ΔV melt was recovered with sufficiently larger error of 34%; however, it is the best result compared to other authors.

Conclusions
In the paper we have found the Δ-factor for SIESTA software package and compared it with other packages implementing DFT. We showed, that the Δ-factor depends on the choice of pseudopotential, and suggested a procedure to fit pseudopotential parameters taking the data on Δ-factor into account. Using proposed procedure, we modified the pseudopotential and showed that the new pseudopotential performs better in describing structure and properties of both crystalline and liquid iron. After that, taking structural data from ab initio calculations as reference, we built a set of EAM potentials that differ by the only parameter λ and each of which accurately reproduces the RDF of liquid iron from ab initio calculations; most of the physical characteristics of liquid iron, as well as at melting, appeared to have strong dependence on the parameter. We chose one of the potentials (corresponding to λ = 1) as the most adequate based on the RDF recovery error and found it to recover the whole set of properties at least as accurate as the majority of existing EAM potentials for iron.