Atomistic Simulations of Pure Tin Based on a New Modified Embedded-Atom Method Interatomic Potential

A new interatomic potential for the pure tin (Sn) system is developed on the basis of the second-nearest-neighbor modified embedded-atom-method formalism. The potential parameters were optimized based on the force-matching method utilizing the density functional theory (DFT) database of energies and forces of atomic configurations under various conditions. The developed potential significantly improves the reproducibility of many fundamental physical properties compared to previously reported modified embedded-atom method (MEAM) potentials, especially properties of the β phase that is stable at the ambient condition. Subsequent free energy calculations based on the quasiharmonic approximation and molecular-dynamics simulations verify that the developed potential can be successfully applied to study the allotropic phase transformation between α and β phases and diffusion phenomena of pure tin.


Introduction
The materials systems with tin (Sn) have received great attention owing to their importance in many modern technologies. In the electronics industry, soldering is the most important technique of connecting the substrate and electronic devices. Among various kinds of solders, Sn-Pb alloys are the most popular and traditional solders. Recently, Sn-based lead-free solders such as Sn-Ag-Cu alloys have replaced the traditional solders [1]. In the electrochemical industry, Sn-Li alloys are notable materials due to their applications to electrodes [2]. Moreover, liquid Sn-Li alloys are currently being considered as candidates for plasmafacing materials in tokamak fusion reactors [3,4].
For the understanding of detailed phenomena in Sn and its alloys, atomic scale simulations are effective to complement experiments providing detailed insights into the atomic scale processes. Among those simulations, density functional theory (DFT) calculation is a widely used method with many advantages such as a straightforward selectivity of alloy components and the accuracy. However, the DFT calculation is computationally demanding and cannot be applied in all applications. If the material phenomena are related to a comparably large system which cannot be covered by the DFT calculation, large-scale atomistic simulations such as molecular dynamics (MD) and Monte Carlo (MC) are highly desirable. Because the predictive accuracy of such atomistic simulations is critically dependent on the accuracy of interatomic potentials, the availability of reliable interatomic potentials of target systems is of crucial importance. In particular, a relevant description of pure Sn should take precedence because it can be a basis for the development of interatomic potentials for concerned Sn-based multi-component alloy systems.
Until now, the development of an accurate interatomic potential for pure Sn has been a difficult task because of the complexity of the system due to the presence of allotropic phase transformation. In the periodic table, Sn is located at the borderline between covalent and metallic bonding elements. Above 286 K, Sn crystallizes into a body-centered tetragonal crystal structure (β-Sn) having a characteristic of metallic bonding [5]. Below 286 K, Sn crystallizes into a diamond cubic structure (α-Sn) having a characteristic of covalent bonding [5]. So far, several interatomic potentials for pure Sn have been developed in previous works. There are available potentials based on pair potential models [6][7][8], but the practical utilization of these potentials is greatly limited because these models ignore many-body bonding characteristics of pure Sn. There is an available many-body interatomic potential [9] based on embedded-atom method (EAM) model [10], but this potential cannot predict the β-Sn phase as a stable phase at the ambient condition because target properties of this potential are confined to phases at a high-pressure. Alternatively, there are available many-body interatomic potentials based on the modified embedded-atom method (MEAM) model [11] by Ravelo and Baskes [5], by Vella et al. [12], and by Etesami et al. [13]. It seems reasonable to apply the MEAM model to the Sn system since this model was originally devised to well describe various types of atomic bonds, including the metallic and covalent bonds, in a single formalism. However, these MEAM potentials were developed mostly focusing on properties of liquid phase and showed deficiencies in reproducing physical properties of solid phases, especially for properties of the β-Sn phase.
In the present study, we provide a new interatomic potential based on the second nearest-neighbor modified embedded-atom method (2NN MEAM) model [14][15][16]. The 2NN MEAM is a model that overcomes several shortcomings of the original MEAM model by partially considering 2NN interactions of atoms. In addition, the present study considered a relevant method of the parameter optimization to improve the reliability of the developed potential. In previous works for MEAM potentials [5,12,13], potential parameters were optimized focusing on specific physical properties obtained by experimental and density-functional theory information. In contrast, the present optimization was performed based on the force-matching method proposed by Ercolessi and Adams [17]. This method considers not a physical property itself, but forces and energies related to various atomic configurations including configurations at finite temperatures expected from the DFT calculation. It has been confirmed that this method can greatly improve the accuracy and transferability of developed interatomic potentials [18][19][20].
The remainder of the present article is organized as follows. Section 2 describes detailed processes of the DFT calculations and optimization of potential parameters. In Section 3, the accuracy and transferability of the developed potential is presented with suitable examples. The conclusion is finally drawn in Section 4.

Methods
The optimization of the present interatomic potential was performed based on the force-matching method which considers both the structural energy and the forces acting on each atom as a fitting target. The overall procedure follows previous works [20,21] in which 2NN MEAM potentials for pure Ni, Ti and Li were similarly determined based on the same method in a systematic manner. First, a database of atomic energies and forces of various atomic configurations was established based on the DFT calculation. Then, an optimum set of potential parameters was determined by minimizing errors between the expectation by potential parameters and the DFT database.

Construction of a DFT Database
A series of DFT calculations were performed using Vienna Ab initio Simulation Package (VASP) code [22][23][24] based on the projector augmented wave (PAW) method [25]. For the exchange-correlation functional, the Perdew-Burke-Ernzerhof generalized-gradient approximation (GGA) [26] was used. A plane-wave kinetic-energy cutoff of 400 eV and the Methfessel-Paxton smearing method with a width of 0.1 eV were used. All calculations were performed with a density of the k-point mesh equivalent to a 21 × 21 × 21 mesh for the face-centered cubic (fcc) primitive unit cell, and the corresponding similar density of the k-point mesh were employed for other unit cells and supercells. In the present DFT calculation, atomic configurations resulting from various conditions were considered for the construction of the DFT database to ensure the sufficient transferability of the developed interatomic potential to various possible applications. Detailed information on cell structures and corresponding k-point meshes is summarized in Table 1. Table 1. Atomic configurations considered in the present density functional theory (DFT) calculations. In the Stability column, "stable" indicates that the structure is reported in an equilibrium phase diagram. Other phases are labeled as "hypothetical". The Strain column indicates the strain applied to the supercells, where H, O, and M denote hydrostatic, orthorhombic, and monoclinic strains, respectively.

N Sn atoms N vacancies
The equilibrium lattice constants and bulk modulus at 0 K were calculated by employing the Birch-Murnaghan equation of state [27,28]. For the calculation of properties involved with defects, positions of each atom were relaxed at a constant volume and cell shape. The vacancy migration energy was calculated with a suitable saddle-point configuration utilizing the nudged elastic band (NEB) method [29,30]. For the calculation of the surface energy, rectangular cells with a stacking of 12 Å to 13 Å thick slab and a vacuum region of 10 Å were employed without the relaxation of cell dimensions parallel to free surfaces. The convergence criteria for energies and forces of all defect calculations were 10 −6 eV/atom and 10 −2 eV/Å, respectively. Phonon calculations were performed using the "Phonopy" code [31,32] based on the direct force constant approach [33]. Supercells of 64 atoms were used for the phonon calculation for α and β phases with the convergence criteria for energy and forces of 10 −8 eV and 10 −4 eV/Å, respectively. To calculate the cohesive energy, a reference energy with a single atom was obtained by considering the spin-polarized calculation.
To obtain a sufficient number of effective force data at finite temperatures, two-step DFT calculations were conducted. First, ab initio MD simulations [22] were performed for a total of 2000 steps with a timestep of 1.5 fs using relatively low convergence parameters such as a single k-point and a default value of the cutoff energy. To determine accurate energies and forces of each configuration, well-converged calculations with a higher cutoff energy (400 eV) and denser k-point mesh were then followed using randomly extracted configurations from ab initio MD simulations.

Optimization of Potential Parameters
For the development of a unary potential based on the 2NN MEAM formalism, an optimization of 14 independent potential parameters is required. Four parameters [the cohesive energy (E c ), the equilibrium nearest-neighbor distance (r e ), the bulk modulus (B) of the reference structure and the adjustable parameter d] are involved with to the universal equation of state. Seven parameters [the decay lengths (β (0) , β (1) , β (2) , and β (3) ) and the weighting factors (t (1) , t (2) , and t (3) )] are involved with the electron density. The parameter A is required for the embedding function and two parameters [C min and C max ] are required for the many-body screening. Detailed explanations on these potential parameters can be found in literature [14][15][16].
The target configurations in the DFT database consist of configurations of stable phases (α, β, and liquid) and hypothetical phases (fcc, bcc, and hcp) at various temperatures. In particular, the inclusion of configurations at finite temperatures is indispensable for the force-matching because these configurations can provide effective atomic force data for the optimization. A total of 63 configurations resulting from various temperatures and strain conditions as well as various defect configurations were considered for the fitting. Detailed information on atomic configurations used for the fitting process is listed in Table 1.
The optimization was performed by comparing energies of target configurations and forces on each atom expected by a candidate set of potential parameters and those expected by the DFT calculation. The optimization started with specifying a reference structure of the potential parameters, a radial cutoff distance, and fitting weights of each configuration. An optimizer based on the genetic algorithm then adjusted the candidate set of parameters to reduce the sum of energy errors of each configuration and force errors of each atom. If the derived set of potential parameters does not satisfactorily reproduce overall physical properties, another optimization was performed with a different set of the reference structure, radial cutoff value, and fitting weights of configurations. The trial reference structures were diamond cubic, fcc and bcc structures, and the trial radial cutoff values were 4.5, 5.0, 5.5 and 6.0 Å. During the optimization process, we realized that the consideration of similar weighting between configurations of α and β phases results in a general worsening of the reproducibility of various physical properties. Considering the importance of the β phase stable at ambient condition, the final optimization was performed using increased fitting weights for atomic configurations of the β phase. Figure 1 shows results of the final optimization represented by statistical correlations between the calculated energies and forces from the developed potential and the DFT calculation. Each correlation was obtained using atomic configurations expected by the DFT calculation without further atomic relaxations. The resultant root mean square (RMS) errors of energies and forces are 0.044 eV/atom and 0.097 eV/Å, respectively. These values are significantly higher than the results of previous optimization for pure Li [21]. It seems that this is because the handling of phases with different types of atomic bonding (covalent and metallic) simultaneously is significantly more difficult than the handling of phases with single type of atomic bonding (metallic) based on a single potential formalism. parameters [ min and max ] are required for the many-body screening. Detailed explanations on these potential parameters can be found in literature [14][15][16]. The target configurations in the DFT database consist of configurations of stable phases (α, β, and liquid) and hypothetical phases (fcc, bcc, and hcp) at various temperatures. In particular, the inclusion of configurations at finite temperatures is indispensable for the force-matching because these configurations can provide effective atomic force data for the optimization. A total of 63 configurations resulting from various temperatures and strain conditions as well as various defect configurations were considered for the fitting. Detailed information on atomic configurations used for the fitting process is listed in Table 1.
The optimization was performed by comparing energies of target configurations and forces on each atom expected by a candidate set of potential parameters and those expected by the DFT calculation. The optimization started with specifying a reference structure of the potential parameters, a radial cutoff distance, and fitting weights of each configuration. An optimizer based on the genetic algorithm then adjusted the candidate set of parameters to reduce the sum of energy errors of each configuration and force errors of each atom. If the derived set of potential parameters does not satisfactorily reproduce overall physical properties, another optimization was performed with a different set of the reference structure, radial cutoff value, and fitting weights of configurations. The trial reference structures were diamond cubic, fcc and bcc structures, and the trial radial cutoff values were 4.5, 5.0, 5.5 and 6.0 Å . During the optimization process, we realized that the consideration of similar weighting between configurations of α and β phases results in a general worsening of the reproducibility of various physical properties. Considering the importance of the β phase stable at ambient condition, the final optimization was performed using increased fitting weights for atomic configurations of the β phase. Figure 1 shows results of the final optimization represented by statistical correlations between the calculated energies and forces from the developed potential and the DFT calculation. Each correlation was obtained using atomic configurations expected by the DFT calculation without further atomic relaxations. The resultant root mean square (RMS) errors of energies and forces are 0.044 eV/atom and 0.097 eV/Å , respectively. These values are significantly higher than the results of previous optimization for pure Li [21]. It seems that this is because the handling of phases with different types of atomic bonding (covalent and metallic) simultaneously is significantly more difficult than the handling of phases with single type of atomic bonding (metallic) based on a single potential formalism. We finally confirmed that a reference structure of fcc provides the optimum result in reproducing various physical properties of pure Sn. A cutoff value of 5.0 Å was confirmed to be large We finally confirmed that a reference structure of fcc provides the optimum result in reproducing various physical properties of pure Sn. A cutoff value of 5.0 Å was confirmed to be large enough to reproduce various physical properties with an acceptable computational efficiency. In the subsequent section, all calculations based on the new MEAM potential were performed using this radial cutoff distance. Table 2 lists the final set of potential parameters for pure Sn. Table 2. Optimized 2NN MEAM potential parameter set for the pure Sn system. The following properties are dimensionful: the cohesive energy E c (eV/atom), the equilibrium nearest-neighbor distance r e (Å), and the bulk modulus B (10 12 dyne/cm 2 ). The reference structure is fcc Sn.

Results and Discussion
In this section, results of atomistic simulations performed by the developed potential are presented. The results are compared to corresponding experiments and DFT calculations to evaluate the accuracy and transferability of the developed potential. In addition, the performance of the developed potential was compared with previous MEAM potentials by Ravelo and Baskes [5], by Vella et al. [12], and by Etesami et al. [13]. We used the potential of Ravelo and Baskes [5] in a form adopted by Vella et al. [12]. For the clear comparisons between the present potential and the previous potentials, all physical properties were recalculated in the same manner, whether or not some of properties were already reported in previous studies [5,12,13]. All atomistic simulations were performed based on the LAMMPS code [34]. If not specially designated as MD simulations, obtained properties represent results of molecular statics simulations at 0 K using atomic configurations of at least 4000 atoms. All MD runs were performed starting from initial configurations optimized by corresponding molecular statics simulations at 0 K using a timestep of 0.002 ps. Nosé-Hoover thermostat and barostat [35,36] were used for controlling temperature and pressure, respectively.
In Table 3, calculated bulk, elastic and defect properties of pure Sn at 0 K using the present potential are compared with experimental data, DFT calculations and results using previous MEAM potentials. The experimental cohesive energy, which is also in good agreement with the DFT expectation, is well reproduced by the present potential. The ground state structure of pure Sn expected by the DFT calculation is the α structure. The present potential well reproduces such tendency and the structural energy differences between various stable and hypothetical solid phases as well. Among structural energy differences, the difference between the low temperature stable α phase and the high temperature stable β phase is important for the expectation of the allotropic transformation. The present potential accurately reproduces this trend while previous potentials by Vella et al. [12] and by Etesami et al. [13] show significantly higher energy differences. The structural parameters such as lattice constants of the α phase and lattice constants and c/a ratio of the β phase are also well reproduced by the present potential in closer agreement with experimental data.
In the elastic properties, all MEAM potentials generally present a difficulty in reproducing the C 44 of the β phase. Especially, previous potentials by Ravelo and Baskes [5] and by Etesami et al. [13] indicate the deviation more than one order of magnitude. The present potential significantly improves the reproducibility of the C 44 and other elastic details as well. For example, it was reported that the β phase lattice is stiffer in the c-direction compared to aor b-direction (C 33 > C 11 ). The present potential correctly reproduces this trend as well as the absolute values of each elastic constant.
The defect properties of the β phase such as the properties related to the vacancy and the free surface are also examined and listed in Table 3. In the present study, the activation energy of vacancy diffusion (Q vac ) is defined as a sum of the vacancy formation energy (E vac f ) and the vacancy migration energy (E vac m ) considering the reported substitutional diffusion mechanism of the β phase [37]. Based on previous experimental study [38] and the present DFT calculation, the self-diffusivity of the β phase is anisotropic along crystal directions, and the activation energy of vacancy diffusion along the c-direction ( c) is higher than that along the aor b-direction (⊥c). The present potential successfully reproduces this trend. In the calculation of the vacancy migration energy and resultant activation energy of vacancy diffusion, we realized that the calculation along the aor b-direction by previous potentials indicates unphysical results such as the negative migration energy and the instability of the saddle point position. These values are thus not presented in Table 3. The surface energies of the β phase are also correctly reproduced only by the present potential while previous potentials show a significant overestimation. Table 3. Calculated bulk, elastic and defect properties of pure Sn using the present 2NN MEAM potential, in comparison with experimental data, DFT data, and the calculation results using previous MEAM potentials by Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13]. The following quantities are listed: the cohesive energy E c (eV/atom), the lattice constant a (Å), the bulk modulus B and the elastic constants C 11 , C 12 , C 13 , C 33 , C 44 and C 66 (10 12 dyne/cm 2 ), the structural energy differences ∆E (eV/atom), the vacancy formation energy E vac   [39]. b Ref. [40]. c Ref. [41]. d Ref. [42]. e Ref. [38]. f Present DFT calculation. g Present calculation using the reported potential parameters. h The potential parameters adopted by Vella et al. [12] are used.
We further examined the transferability of the developed potential by comparing a group of properties at finite temperatures. Table 4 lists the thermal properties of the β phase (the thermal expansion coefficient and the specific heat), in comparison with experimental data. These properties were obtained based on the MD simulations using an isobaric-isothermal (NPT) ensemble at the target temperature and zero pressure. For the thermal expansion coefficient of the β phase, the potential by Ravelo and Baskes [5] shows a better agreement with the experiment than other potentials. Instead, other potentials show a better reproducibility than the potential by Ravelo and Baskes [5] in the case of the specific heat of the β phase. Table 4. Calculated thermal properties of pure Sn using the present 2NN MEAM potential, in comparison with experimental data and the calculation results using previous MEAM potentials by Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13]. The listed quantities correspond to the thermal expansion coefficient ε (10 −6 /K), the heat capacity at constant pressure C P (J/mol K), the melting temperature T m (K), the enthalpy of melting ∆H m (kJ/mol), and the volume change upon melting ∆V m /V solid (%).

Property
Exp  [41]. b Ref. [43]. c Ref. [44]. d Present calculation using the reported potential parameters. e The potential parameters adopted by Vella et al. [12] are used. Figure 2 shows the temperature dependence of the atomic volume of pure Sn. Initially, the β phase with 16,224 atoms was equilibrated at 5 K, and the temperature was then gradually increased 1000 K with a heating rate of 1.0 K/ps using the NPT ensemble at zero pressure. Each potential indicates a discontinuous jump in the volume at a certain temperature. This jump indicates the occurrence of the melting (β → liquid). Because the heating simulation was performed without any heterogeneous nucleation sites, the melting occurs at a temperature much higher than the equilibrium melting temperature. Therefore, this temperature can be regarded as an overheated melting point. The equilibrium melting temperature at zero pressure was further calculated using the interface method [45,46] employing a simulation cell consisting of liquid and solid phases and listed in Table 4. The table also lists the enthalpy change and the volume change upon melting which were obtained at the calculated equilibrium melting temperature. For the previous potential by Ravelo and Baskes [5], it was reported that the liquid part of the interface simulation is crystallized into a structure different from the β phase [12], and thus the equilibrium melting point and resultant enthalpy and volume changes are not presented. The developed potential shows discrepancies in the properties associated with the melting while the potential by Etesami et al. [13] indicates the best agreement in the melting point with the experiment. This can be interpreted by the fact that the potential [13] was developed focusing mostly on the melting phenomenon without consideration to the reproducibility of many other important properties.
As listed in Table 4, the present potential is somewhat deficient in describing properties related to the melting compared to previous potentials. The present potential underestimates the enthalpy of melting and overestimates the volume change upon melting compared to experimental values. Considering that the present potential satisfactory describes physical properties of β phase, it is necessary to further investigate the properties associated with the structure of the liquid phase. Figure 3 shows the calculated radial distribution function of liquid structures at different temperatures. Although the present potential shows a generally acceptable quality to reproduce the height and position of the first peak, it shows a deficiency in reproducing characteristics of other peaks. We attribute these deficiencies to the use of increased fitting weights of the β phase compared to other phases during the fitting process.  Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13]. The heating simulations were started with the β phase at 5 K. The discrete jumps represent the occurrence of the melting.  [7]) and the calculation results using previous MEAM potentials by Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13].
We further examined the reproducibility of phonon properties and resultant properties at finite temperature based on the harmonic approximation (HA) and quasiharmonic (QHA) approximation. Figure 4 shows the phonon density of states (DOS) of α and β phases expected by the present and previous potentials compared to the DFT expectation. As shown in Figure 4a, the phonon DOS of the α phase is closely reproduced by the present potential. In contrast, previous potentials indicate a significant overestimation of the phonon DOS at high frequencies compared to the DFT expectation. A similar trend is also presented for the phonon DOS of the β phase as shown in Figure 4b. Moreover, previous potentials cannot reproduce the stability of the β phase under perturbative forces. As shown in the magnified figure of the phonon DOS of the β phase (Figure 4c), previous potentials indicate dynamical instability (imaginary phonon frequencies) in contrast to expectations by the present potential and DFT calculations. Even though this problem of previous potentials is expected to interfere with various possible applications of previous potentials, the present potential is free from this problem.  Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13]. The heating simulations were started with the β phase at 5 K. The discrete jumps represent the occurrence of the melting.  Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13]. The heating simulations were started with the β phase at 5 K. The discrete jumps represent the occurrence of the melting.
We further examined the reproducibility of phonon properties and resultant properties at finite temperature based on the harmonic approximation (HA) and quasiharmonic (QHA) approximation. Figure 4 shows the phonon density of states (DOS) of α and β phases expected by the present and previous potentials compared to the DFT expectation. As shown in Figure 4a, the phonon DOS of the α phase is closely reproduced by the present potential. In contrast, previous potentials indicate a significant overestimation of the phonon DOS at high frequencies compared to the DFT expectation. A similar trend is also presented for the phonon DOS of the β phase as shown in Figure 4b. Moreover, previous potentials cannot reproduce the stability of the β phase under perturbative forces. As shown in the magnified figure of the phonon DOS of the β phase (Figure 4c), previous potentials indicate dynamical instability (imaginary phonon frequencies) in contrast to expectations by the present potential and DFT calculations. Even though this problem of previous potentials is expected to interfere with various possible applications of previous potentials, the present potential is free from this problem.  [7]) and the calculation results using previous MEAM potentials by Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13].
We further examined the reproducibility of phonon properties and resultant properties at finite temperature based on the harmonic approximation (HA) and quasiharmonic (QHA) approximation. Figure 4 shows the phonon density of states (DOS) of α and β phases expected by the present and previous potentials compared to the DFT expectation. As shown in Figure 4a, the phonon DOS of the α phase is closely reproduced by the present potential. In contrast, previous potentials indicate a significant overestimation of the phonon DOS at high frequencies compared to the DFT expectation. A similar trend is also presented for the phonon DOS of the β phase as shown in Figure 4b. Moreover, previous potentials cannot reproduce the stability of the β phase under perturbative forces. As shown in the magnified figure of the phonon DOS of the β phase (Figure 4c), previous potentials indicate dynamical instability (imaginary phonon frequencies) in contrast to expectations by the present potential and DFT calculations. Even though this problem of previous potentials is expected to interfere with various possible applications of previous potentials, the present potential is free from this problem.  Therefore, only the present potential can be further utilized to expect the physical properties at finite temperatures based on the QHA. One simple example is calculations of the thermal expansion coefficient and the heat capacity as presented in Figure 5. The present potential exhibits an acceptable quality comparable to the DFT calculation when reproducing such properties. Another important example is the expectation of the allotropic phase transformation between α and β phases. Because the present potential can accurately reproduce the phonon properties without suffering from the imaginary phonon modes, a vibrational contribution on the free energy can be calculated based on the QHA. The Gibbs free energy at the given temperature T and pressure P can be obtained from the Helmholtz free energy ( , ) at the given temperature and volume V through the transformation, where the right-hand side of this equation represents a value when the minimum value is found in square brackets by varying the volume. In the present study, the ( , ) was computed by a sum of the DFT total energy ( ) at the given volume, and the vibrational contribution to the free energy ( , ) at the given temperature and volume.
( , ) = ( ) + ( , ) A detailed explanation on the determination of the free energy based on the phonon calculation is given in Ref. [47]. Therefore, only the present potential can be further utilized to expect the physical properties at finite temperatures based on the QHA. One simple example is calculations of the thermal expansion coefficient and the heat capacity as presented in Figure 5. The present potential exhibits an acceptable quality comparable to the DFT calculation when reproducing such properties. Another important example is the expectation of the allotropic phase transformation between α and β phases. Because the present potential can accurately reproduce the phonon properties without suffering from the imaginary phonon modes, a vibrational contribution on the free energy can be calculated based on the QHA. The Gibbs free energy at the given temperature T and pressure P can be obtained from the Helmholtz free energy F(T, V) at the given temperature and volume V through the transformation, where the right-hand side of this equation represents a value when the minimum value is found in square brackets by varying the volume. In the present study, the F(T, V) was computed by a sum of the DFT total energy E DFT (V) at the given volume, and the vibrational contribution to the free energy F vib (T, V) at the given temperature and volume.
A detailed explanation on the determination of the free energy based on the phonon calculation is given in Ref. [47]. compared with those using the present DFT calculation and experimental data (Refs. [43,48]). Figure 6 shows the calculated temperature dependence of Gibbs free energies of α and β phases when the reference state is set to the β phase at 0 K. As stated, both DFT and the present MEAM potential expect that the α phase is a ground state of pure Sn. The Gibbs free energy of each phase decreases with increasing temperature due to the increasing contribution of the vibrational entropy. Because the vibrational contribution is more significant in the β phase than in the α phase, the β phase starts to become more stable than the α phase at a certain temperature. We define this temperature as the phase transformation temperature between α and β phases. The calculated transformation temperature by the present MEAM potential (460 ± 5 K) agrees well with that by the present DFT calculation (470 ± 5 K) while these values are higher than experimental data (286 K) [5]. Although this discrepancy seems to be caused by the limitation of the QHA, the present potential at least reproduces the DFT expectation.
As a final confirmation of the transferability of the developed potential, the self-diffusion of liquid Sn was examined using the MD simulation. In the atomic scale, the MD simulation provides a time dependency of the mean square displacement (MSD), and this is related to the bulk diffusivity based on the following equation (Einstein relation), where D is the diffusivity, < 2 ( ) > is the MSD of atoms, and t is the time.  The results using the present 2NN MEAM potential are compared with those using the present DFT calculation and experimental data (Refs. [43,48]). Figure 6 shows the calculated temperature dependence of Gibbs free energies of α and β phases when the reference state is set to the β phase at 0 K. As stated, both DFT and the present MEAM potential expect that the α phase is a ground state of pure Sn. The Gibbs free energy of each phase decreases with increasing temperature due to the increasing contribution of the vibrational entropy. Because the vibrational contribution is more significant in the β phase than in the α phase, the β phase starts to become more stable than the α phase at a certain temperature. We define this temperature as the phase transformation temperature between α and β phases. The calculated transformation temperature by the present MEAM potential (460 ± 5 K) agrees well with that by the present DFT calculation (470 ± 5 K) while these values are higher than experimental data (286 K) [5]. Although this discrepancy seems to be caused by the limitation of the QHA, the present potential at least reproduces the DFT expectation. compared with those using the present DFT calculation and experimental data (Refs. [43,48]). Figure 6 shows the calculated temperature dependence of Gibbs free energies of α and β phases when the reference state is set to the β phase at 0 K. As stated, both DFT and the present MEAM potential expect that the α phase is a ground state of pure Sn. The Gibbs free energy of each phase decreases with increasing temperature due to the increasing contribution of the vibrational entropy. Because the vibrational contribution is more significant in the β phase than in the α phase, the β phase starts to become more stable than the α phase at a certain temperature. We define this temperature as the phase transformation temperature between α and β phases. The calculated transformation temperature by the present MEAM potential (460 ± 5 K) agrees well with that by the present DFT calculation (470 ± 5 K) while these values are higher than experimental data (286 K) [5]. Although this discrepancy seems to be caused by the limitation of the QHA, the present potential at least reproduces the DFT expectation.
As a final confirmation of the transferability of the developed potential, the self-diffusion of liquid Sn was examined using the MD simulation. In the atomic scale, the MD simulation provides a time dependency of the mean square displacement (MSD), and this is related to the bulk diffusivity based on the following equation (Einstein relation), where D is the diffusivity, < 2 ( ) > is the MSD of atoms, and t is the time.  As a final confirmation of the transferability of the developed potential, the self-diffusion of liquid Sn was examined using the MD simulation. In the atomic scale, the MD simulation provides a time dependency of the mean square displacement (MSD), and this is related to the bulk diffusivity based on the following equation (Einstein relation), where D is the diffusivity, R 2 (t) is the MSD of atoms, and t is the time.
The self-diffusivity of liquid Sn was calculated using NPT ensemble MD runs starting with the β phase of 16,224 atoms. First, the initial structure was melted by maintaining the cell at a high temperature (2000 K) for a sufficient time (40 ps). The temperature was then changed to each target temperature within a range of 600 and 2000 K, and the relaxation was performed for 40 ps. The simulation cell was further maintained at the target temperature for a total simulation time of 2000 ps, and the time evolution of the MSD was counted at every 4 ps. We confirmed that the MSD of atoms accumulated through these conditions shows a clear linear relationship with the time. Figure 7 shows resultant self-diffusivities of the liquid phase at various temperatures (600-2000 K) compared with experimental data. Despite a deviation between experimental results, the present potential successfully reproduces the general trend of an experimental result by Bruson et al. [49]. The calculated self-diffusivities using the present potential are similar to those using the previous potential by Vella et al. [12] even though the values are slightly underestimated at low temperatures. Compared to the previous potentials by Ravelo and Baskes [5] and by Etesami et al. [13], the present potential indicates significantly improved results, especially at high temperatures. It is interesting to note that the self-diffusivities of the liquid phase are sufficiently reproduced despite the discrepancy of properties related to the liquid structure. The result seems to emphasize the effectiveness of the force-matching method to obtain a robust potential that demonstrates a good transferability to kinetic properties such as the migration energy and the attempt frequency of diffusing atoms. However, it should be noted that the transferability of the present potential to kinetic properties of low-coordinated atoms on surfaces (surface diffusion) is not guaranteed because the present fitting was mostly performed by using atomic configurations in the bulk state. present DFT calculation are illustrated. Each energy value was calculated with respect to the reference state of the β phase at 0 K.
The self-diffusivity of liquid Sn was calculated using NPT ensemble MD runs starting with the β phase of 16,224 atoms. First, the initial structure was melted by maintaining the cell at a high temperature (2000 K) for a sufficient time (40 ps). The temperature was then changed to each target temperature within a range of 600 and 2000 K, and the relaxation was performed for 40 ps. The simulation cell was further maintained at the target temperature for a total simulation time of 2000 ps, and the time evolution of the MSD was counted at every 4 ps. We confirmed that the MSD of atoms accumulated through these conditions shows a clear linear relationship with the time. Figure  7 shows resultant self-diffusivities of the liquid phase at various temperatures (600-2000 K) compared with experimental data. Despite a deviation between experimental results, the present potential successfully reproduces the general trend of an experimental result by Bruson et al. [49]. The calculated self-diffusivities using the present potential are similar to those using the previous potential by Vella et al. [12] even though the values are slightly underestimated at low temperatures. Compared to the previous potentials by Ravelo and Baskes [5] and by Etesami et al. [13], the present potential indicates significantly improved results, especially at high temperatures. It is interesting to note that the self-diffusivities of the liquid phase are sufficiently reproduced despite the discrepancy of properties related to the liquid structure. The result seems to emphasize the effectiveness of the force-matching method to obtain a robust potential that demonstrates a good transferability to kinetic properties such as the migration energy and the attempt frequency of diffusing atoms. However, it should be noted that the transferability of the present potential to kinetic properties of lowcoordinated atoms on surfaces (surface diffusion) is not guaranteed because the present fitting was mostly performed by using atomic configurations in the bulk state. Figure 7. Calculated self-diffusivity of liquid Sn using the present 2NN MEAM potential, in comparison with experimental data (Refs. [49][50][51][52]) and the calculation results using previous MEAM potentials by Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13].
To summarize, we have shown that the present potential generally reproduces various physical properties of pure Sn, even though there is a difference in the performance for each property. The developed potential performs very well in describing structural, elastic and thermal properties of the β phase, phonon properties of solid phases, and diffusion properties of solid and liquid phases. It is less suited to describing the melting behavior and the liquid structure compared to previous MEAM potentials [5,12,13], which were developed to target specific physical properties of the liquid phase. This fact should be kept in mind in future applications of the present potential. The LAMMPS files implementing the developed potential can be obtained from an online repository [53].  [49][50][51][52]) and the calculation results using previous MEAM potentials by Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13].
To summarize, we have shown that the present potential generally reproduces various physical properties of pure Sn, even though there is a difference in the performance for each property. The developed potential performs very well in describing structural, elastic and thermal properties of the β phase, phonon properties of solid phases, and diffusion properties of solid and liquid phases. It is less suited to describing the melting behavior and the liquid structure compared to previous MEAM potentials [5,12,13], which were developed to target specific physical properties of the liquid phase. This fact should be kept in mind in future applications of the present potential. The LAMMPS files implementing the developed potential can be obtained from an online repository [53].

Conclusions
We provide a new robust interatomic potential for pure Sn on the basis of the 2NN MEAM model. The potential is developed based on the force-matching method utilizing the DFT database of energies and forces of atomic configurations under various conditions. The developed potential significantly improves the reproducibility of various physical properties compared to previously reported MEAM potentials, especially for properties of the β phase. The present potential exhibits superior transferability to the phonon properties and the properties related to diffusion phenomena. As a possible example, the allotropic phase transformation between α and β phases is investigated based on the QHA. The results indicate that the present potential can be successfully used for the expectation of the transformation temperature at least at the level of the DFT expectation. The present potential also successfully reproduces the self-diffusivity of liquid Sn. The developed potential can be a suitable basis for implementing atomistic simulations of technically important Sn-based alloy systems.