Performance of the standard exchange-correlation functionals in predicting melting properties fully from first principles Application to Al and magnetic Ni

We apply the efﬁcient two-optimized references thermodynamic integration using Langevin dynamics method [Phys. Rev. B 96 , 224202 (2017)] to calculate highly accurate melting properties of Al and magnetic Ni from ﬁrst principles. For Ni we carefully investigate the impact of magnetism on the liquid and solid free energies including longitudinal spin ﬂuctuations and the reverse inﬂuence of atomic vibrations on magnetic properties. We show that magnetic ﬂuctuations are effectively canceling out for both phases and are thus not altering the predicted melting temperature. For both elements, the generalized gradient approximation (GGA) and the local-density approximation (LDA) are used for the exchange-correlation functional revealing a reliable ab initio conﬁdence interval capturing the respective experimental melting point, enthalpy of fusion, and entropy of fusion.

(Received 17 January 2020; revised manuscript received 24 March 2020; accepted 3 April 2020; published 30 April 2020) We apply the efficient two-optimized references thermodynamic integration using Langevin dynamics method [Phys. Rev. B 96, 224202 (2017)] to calculate highly accurate melting properties of Al and magnetic Ni from first principles. For Ni we carefully investigate the impact of magnetism on the liquid and solid free energies including longitudinal spin fluctuations and the reverse influence of atomic vibrations on magnetic properties. We show that magnetic fluctuations are effectively canceling out for both phases and are thus not altering the predicted melting temperature. For both elements, the generalized gradient approximation (GGA) and the local-density approximation (LDA) are used for the exchange-correlation functional revealing a reliable ab initio confidence interval capturing the respective experimental melting point, enthalpy of fusion, and entropy of fusion. DOI: 10.1103/PhysRevB.101.144108

I. INTRODUCTION
Melting points, T m , of materials have always been an important topic in materials science. The theoretical determination of T m is by now possible with ab initio molecular dynamics (MD) simulations based on density functional theory (DFT). A specific approach for that purpose is the free-energy approach in which the melting temperature is determined by the crossing point of the Gibbs energies of solid and liquid. After its first application to Si in 1995 [1], a number of works have followed utilizing this approach to various materials [2][3][4][5]. However, even though the ab initio based free-energy approach provides in principle DFT accuracy, it faces two critical challenges-the high computational cost and the reliance on a good reference system for the involved thermodynamic integration-which have hindered a more general application.
An important strategy to improve the ab initio based freeenergy approach is to design a good reference system to speed up the thermodynamic integration during the solid and liquid free-energy calculations. Alfè and co-workers [3][4][5] made developments in this direction. Specifically, they applied an inverse power potential as an intermediate reference between the Lennard-Jones liquid and the ab initio system. Using a similar idea, Duff et al. [6] developed the twostage up-sampled thermodynamic integration using Langevin dynamics (TU-TILD) method for efficiently calculating the anharmonic free energy of the solid phase. The critical feature of TU-TILD is that the reference potential-for which an embedded atom method (EAM) parametrization was used in the original work-is fitted to a set of DFT MD energies so that it reproduces the ab initio system very closely at finite temperatures. Extending this idea, we have recently proposed the two-optimized references thermodynamic integration using Langevin dynamics (TOR-TILD) methodology [7] for melting property calculations, in which a second optimized reference potential is fitted to DFT MD energies of the liquid phase. Additionally, TOR-TILD provides a special strategy for calculating the free-energy surface on the reference potential level by utilizing the interface method. Meanwhile, the idea of designing a better reference system to speed up the thermodynamic integration has been extended to the machine learning potentials [8] and artificial neural network potentials [9]. In the present work we apply the TOR-TILD method utilizing two bespoke EAM potentials to calculate the melting properties of aluminum and nickel.
The melting point of Al has been the subject of several theoretical studies, but the results show significant scatter depending on the employed method as shown in Table I. In the early days interatomic potentials were generally used to calculate the melting point. Moriarty et al. [10] predicted a melting point of 1050 K with the generalized pseudopotential theory. Mei and Davenport [11] applied an embedded atom method (EAM) potential fitted to structural properties and obtained a value of 800 K. Straub et al. [12] employed a TABLE I. Theoretical studies on the melting point of Al utilizing the free-energy approach in comparison to experiment; GPT = generalized pseudopotential theory; EAM = embedded atom method; OF-DFT = orbital-free density functional theory; LDA = local-density approximation; GGA = generalized gradient approximation; PW91 = Perdew-Wang functional [14]. See text for further details. two-body potential fitted to various ab initio results and predicted a melting point of 955 K. It is apparent that the predicted melting point strongly depends on the chosen potential. In 1998, de Wijs et al. [2] calculated the melting point of Al fully from ab initio, utilizing the local-density approximation (LDA) for the exchange-correlation functional. They [2] predicted a melting temperature of 890 K, i.e., about 40 K below the experimental value of 933 K. Jesson and Madden [13] applied a special orbital-free version of DFT based on LDA input and predicted a rather low melting point of 615 K. Two years later, Vočadlo and Alfè [5] used again Kohn-Sham based DFT (as de Wijs et al. [2]) but in the generalized gradient approximation (GGA). A melting point of 786 K was obtained, which was subsequently corrected to 912 K by employing experimental lattice constant information [5]. From this collection of previous results, it becomes apparent that not only do the T m results from interatomic potentials show substantial scatter, but the direct ab initio predictions do as well. In the present work we therefore calculate the melting point (and other melting properties) of Al with the highly accurate TOR-TILD method utilizing both the LDA and GGA for the exchange-correlation functional. We compare our results to those of de Wijs et al. [2] and Vočadlo and Alfè [5] in Sec. III C.
TOR-TILD is not limited to calculating the melting properties of nonmagnetic phases, such as Al (present work) or Cu [7]. It can be extended to magnetic materials such as Ni by the inclusion of terms related to magnetic excitations. As regards Ni, ab initio studies on its melting properties are rare, in particular at low pressure. Pozzo and Alfè [17] utilized the coexistence approach (also known as the interface method [18]) to calculate the melting curve of Ni as a function of pressure (up to 100 GPa) with GGA. They obtained a zero pressure melting temperature of 1637 K [17], which agrees reasonably well with the experimental value of 1728 K [19]. Spin polarization was not included in Ref. [17] by anticipating that nonmagnetic calculations could provide a reasonable description at sufficiently high temperatures. However, magnetism is known to significantly contribute to the free energy of the solid phases of Ni, although cancellation is observed when considering phase stabilities [20][21][22][23][24]. The impact of magnetism on the liquid phase and on the melting properties remains to be investigated.
In our present work, we therefore apply the TOR-TILD methodology to calculate the melting properties of magnetic Ni with both LDA and GGA. In all MD simulations spin polarization of the Ni atoms and the coupling to electrons and atomic vibrations is included. As Ni becomes paramagnetic above its Curie temperature of 628 K a recently developed longitudinal spin fluctuations theory [25] is employed to compute the effective local moments of Ni at high temperatures in both the solid and liquid structure to clarify the magnetic contribution to the free energy. The combined impact of vibrations and magnetic fluctuations on the electronic free energy is also investigated.

II. METHODOLOGY
Within the free-energy approach the Gibbs energies of solid and liquid need to be calculated. As fixed volume conditions are well amenable for ab initio calculations, we start with a calculation of the Helmholtz free energy F (V, T ) in an NV T ensemble as a function of volume V and temperature T . Once F (V, T ) is available for the relevant V and T range, the Gibbs energy G(P, T ) can be obtained by a Legendre transformation, where the pressure P is obtained from the free-energy surface as P = −(∂F/∂V ) T . The melting point T m is determined by the condition of equality of the Gibbs energies of solid, G solid (P, T ), and liquid, G liquid (P, T ). To get T m fully from ab initio, we therefore need to calculate F solid (V, T ) and F liquid (V, T ) from first principles and apply Eq. (1) to both. Any other thermodynamic equilibrium quantity, such as the enthalpy or entropy, can be derived by proper thermodynamic relations [26]. Within the free-energy Born-Oppenheimer approximation [27,28] the total free energy of a solid can be decomposed as where E (V ) is the total electronic energy at T =0 K, F el (V, T ) the electronic contribution for a static lattice, F qh (V, T ) the quasiharmonic contribution, F ah (V, T ) the anharmonic contribution, and F el-vib (V, T ) the adiabatic electron-phonon coupling. The last term F mag-el,vib (V, T ) applies to magnetic materials and can be decomposed as where F mag is the magnetic contribution for the static lattice, F mag-el its coupling to electrons, and F mag-vib its coupling to atomic vibrations. In the present work, for the solid phases, we resort to the results from Ref. [20] in which the free-energy surfaces of fcc Al and fcc Ni were obtained with high accuracy using the TU-TILD method [6]. These free-energy surfaces contain all the relevant contributions for both elements according to the above prescription. Noteworthy, the impact of the atomic vibrations on the electronic free energy, F el-vib (V, T ) in Eq. (2), was fully taken into account in Ref. [20] and shown to be of importance for Ni. Further, the magnetic contribution for solid Ni, F mag-el,vib (V, T ), was intensively analyzed in Ref. [20].
In the here considered temperature window from 1500 to 2000 K, Ni shows paramagnetism since its Curie temperature is much lower, 628 K. Therefore, we assume that magnetic short-range order can be neglected and that a hightemperature approximation is applicable. The magnetic free energy, is in this case dominated by the magnetic entropy S mag . To include the contributions due to transverse fluctuations of the local magnetic moments m, we utilize an often employed mean-field approximation (MFA) for S mag as where k B is the Boltzmann constant. Note that Eq. (5) has been originally derived for an idealized magnetic system with localized m's, i.e., with a fixed magnitude. Longitudinal spin fluctuations (LSFs) in Ni and their explicit magnetic entropy contribution are challenging to model explicitly as discussed in Ref. [21]. In order to approximate the impact of LSFs we have thus applied the following procedure. We have combined Eq. (5) with magnetic moments m derived from our MD calculations. Note that although Eq. (5) derives from a localized system, it allows us to study, at least qualitatively, the impact of LSFs. A conceptual difficulty in combination with supercell methods as employed in the present work is, however, that approaches utilizing only directionally constrained local magnetic moments cannot be applied to Ni because a standard supercell calculation with disordered magnetic moments will converge into a nonmagnetic solution. To circumvent this difficulty, the local magnetic moments were obtained from spin polarized DFT MD simulations in the ferromagnetic state for different volumes and temperatures with the electronic temperature adjusted to the temperature of the atomic vibrations. In this way the impact of atomic motion on the site-averaged local magnetic moments, i.e., m MD , is included as well as the impact of the thermal expansion, such as the coupling term F mag-vib in Eq. (3). We have also studied the impact of explicit magnetic disorder, its interplay with atomic vibrations and electronic density of states (DOS), by employing a recently developed LSF theory [25] in combination with the disordered-local moment (DLM) approach. A detailed analysis will be given in Sec. III D.
For the liquid free energy, a decomposition as performed in Eq. (2) for the solid free energy is not possible because a static reference lattice is missing. Instead, one needs to fully rely on thermodynamic integration from an appropriate reference system. According to the TOR-TILD prescription, two optimized reference potentials are utilized in order to speed up the liquid free-energy calculations. Thus, the liquid free energy can be written as Here, F liquid ref1 is the free energy of the first reference system fitted to the ab initio solid phase (labeled "ref1"), F liquid ref1→ref2 is the free-energy difference between "ref1" and the second reference system fitted to the ab initio liquid phase (labeled "ref2"), and F liquid ref2→DFT is the free-energy difference between "ref2" and the ab initio system. For the calculation of the first term, F liquid ref1 , TOR-TILD provides a specially designed procedure by connecting the liquid free energy to the solid free energy at the melting point of "ref1" [7]. The last two terms in Eq. (6) are obtained by thermodynamic integration. As for F liquid ref2→DFT , the number of expensive DFT calculations is reduced to a minimum because "ref2" reproduces the ab initio liquid phase accurately. Importantly, the determination of the first two terms in Eq. (6) involves only efficient empirical potential calculations and is therefore extremely efficient, even in large supercells to overcome finite-size effects [7]. For the "ref1" potentials we utilize here the optimized potentials developed for the solid anharmonic free-energy calculations in Ref. [20]. The "ref2" potentials have been newly fitted for the present work. The fitting information for all employed reference potentials will be provided in Sec. III A.
Within our approach electronic excitations and their adiabatic coupling to atomic motion are fully included for the liquid phase at the stage of calculating F liquid ref2→DFT . For liquid Ni, spin polarization of the Ni atoms and the coupling of magnetism to electrons and atomic vibrations are included during the calculation of F liquid ref2→DFT by the same approach as used for the corresponding free energy term F mag-el,vib of the solid. In particular, the high-temperature approximation is applied for the treatment of the magnetic contribution of liquid Ni. The usage of ferromagnetic calculations in both solid and liquid Ni and its impact on the melting point will be discussed in Sec. III D.
Note that for the calculation of F ah and F liquid ref2→DFT , the upsampling technique [29] is applied to efficiently coarse grain the configuration space and thereby reduce the number of the very expensive, highly converged DFT calculations.

A. Computational details
For the DFT calculations we used the projector-augmented wave (PAW) method [30] as implemented in the VASP software package [31][32][33][34]. PAW potentials containing three valence electrons for Al and ten valence electrons for Ni were used. LDA and GGA were employed as exchange-correlation functionals, with the Perdew-Burke-Ernzerhof (PBE) [35] parametrization for GGA. The sets of explicitly DFT computed volume and temperature points for the solid and liquid free-energy surfaces are given in Table II for Al and in  Table III for Ni. For details on the fitting procedure of the 144108-3 corresponding DFT solid free-energy surface and the respective convergence parameters we refer to our study in Ref. [20]. The DFT liquid calculations, i.e., calculations for F liquid ref2→DFT in Eq. (6), were performed in a 3×3×3 supercell with 108 atoms for both Al and Ni. The explicitly computed DFT points were used as input to fit polynomials up to third order to obtain an analytical description of the free-energy surface as a function of volume and temperature. As mentioned in the methodology part, the upsampling technique [29] was applied to speed up the thermodynamic integration calculations. For Al the plane-wave cutoff and k-point mesh (Monkhorst-Pack [36]) were set to 250 eV and 3×3×3 for the low converged calculations, and 400 eV and 8×8×8 for the high converged calculations. For Ni the plane-wave cutoff and k-point mesh (Monkhorst-Pack [36]) were set to 300 eV and 2×2×2 for the low converged calculations, and 450 eV and 5×5×5 for the high converged calculations. The electronic structure calculations for modeling the paramagnetic state of Ni were performed employing the exactmuffin-tin-orbital (EMTO) method [37,38] in the Lyngby version of the code [39] within LDA. A 4×4×4 k-point mesh was used for the 108-atom supercells. The basis functions in the calculations were expanded up to l max = 2 (s, p, and d orbitals). To improve the convergence with partial waves, higher tails (up to l = 3) were included in the charge density calculations. For the integration the Green's functions were calculated in combination with the Fermi smearing at a high temperature of 1750 K. Longitudinal spin fluctuations (LSFs) were modelled combining the disordered local-moment approach and the recently proposed LSF theory [25].
For the reference potential calculations, i.e., calculations for F liquid ref1 and F liquid ref1→ref2 in Eq. (6), we used the LAMMPS software package [40]. The coexistence approach [18] was employed to calculate the melting point of the "ref1" potential utilizing a tetragonal 10×10×20 supercell with 8000 atoms resulting in T m ref 1 = 824 K for Al (experiment: 933 K) and T m ref 1 = 1472 K for Ni (experiment: 1728 K). The other reference potential calculations were performed in a cubic 10×10×10 supercell with 4000 atoms. For fitting the liquid free-energy surface of our potential, i.e., F liquid ref1 (V, T ), we used the same volumes as for the DFT calculations (see Tables II  and III), but at a denser temperature mesh (steps of 5 K).
Both "ref1" and "ref2" potentials were fitted with MEAMfit code [41]. The volumes and temperatures used for the fitting are summarized in Tables IV and V. The resulting potential parameters are provided in Tables VI and VII. We used in particular three pair terms for each of the density contributions and pair potentials. For details on the meaning of the potential parameters we refer to our previous work [6].
For the MD simulations we used a time step of 5 fs and the Langevin thermostat with a friction parameter of 0.01 fs −1 to control the temperature.

B. Computational efficiency
The efficiency of an optimized reference potential can be quantified by the standard deviation of the difference to DFT energies during the MD simulations as follows: where E is the energy difference between the reference potential and the ab initio liquid calculated with low converged DFT parameters. The smaller σ is, the less of the expensive DFT MD calculations are required to statistically converge the thermodynamic integration. Importantly, since the standard error σ n is proportional to the square root of the number of (uncorrelated) sampling steps n, σ n = σ / √ n, reducing σ for example by half results in a speed-up factor of 4, i.e., only 1/4 MD steps are required to get the same standard error (for which one typically targets at 1 meV/atom). Therefore,  in the following we evaluate the efficiency of the reference potentials by the square of the standard deviation. We concentrate on the evaluation of the "ref2" potentials for the liquid thermodynamic integration to DFT. We recall our previous results for Cu [7] to enhance the comparison. In Ref. [7] we investigated the efficiency of the "ref2" potential (fitted to DFT MD trajectories of liquid Cu) by comparing it to the "ref1" potential (fitted to DFT MD trajectories of solid Cu) and a well-developed EAM potential available from the literature [42] (labeled "lit"). We found [7] that with our highly optimized "ref2" potential the efficiency of the thermodynamic integration for the liquid phase is drastically improved, especially at the coupling parameter λ = 1. At λ = 1 the standard deviations of "ref2", "ref1", and "lit" were respectively 1.3, 2.2, and 5.5 meV/atom. Thus, our "ref2" potential is more efficient than the "ref1" potential by a factor of ∼3, i.e., (2.2/1.3) 2 , and much faster than "lit" by an order of magnitude, i.e., (5.5/1.3) 2 , in the CPU timing.
We have done similar efficiency tests for liquid Al and Ni with "ref1", "ref2", and "lit" at λ = 1 with LDA ("lit" for Al from Ref. [43] and for Ni from Ref. [42]). The thermodynamic integration calculations were performed at a lattice constant of 4.05 Å and a temperature of 1200 K for Al and 3.54 Å and 2000 K for Ni. The standard deviations of "ref2", "ref1", and "lit" for Al are 1.5, 1.9, and 4.6 meV/atom, respectively. For Ni we have 2.7, 5.9, and 13.8 meV/atom, respectively. The resulting speed-up factors of the "ref2" potentials defined by (σ /σ ref 2 ) 2 are given in Fig. 1 along with the results for Cu. We can see that the largest speed-up is achieved for the "ref2" potential of Ni, with a speed-up factor of almost 30 as compared to the "lit" potential.
Note that the comparison to literature potentials is done here only for the purpose of a quantitative evaluation of the  1. Speed-up factors of the "ref2" potentials in the thermodynamic integration to DFT with respect to well-developed EAM potentials from the literature ("lit"; red bars) and the "ref1" potentials fitted to the DFT solid phase (blue bars). The speed-up factor is defined as (σ /σ ref 2 ) 2 , where σ is the standard deviation of the corresponding potential and σ ref 2 is the standard deviation of the "ref2" potential. efficiency of our methodology. The TOR-TILD method is fully independent of the availability of empirical potentials from the literature. This is important as for other material systems empirical potentials might not be readily available. Table VIII compile  For Al our GGA calculations yield a value of 888 K, which is 45 K lower than the experimental melting point of 933 K [15,16]. Considering the LDA calculations, we find a melting temperature of 972 K which is 39 K higher than the experimental result. In our previous work [7], we found the same phenomenon for Cu in that the experimental melting point is located close to the middle of the predictions from GGA and LDA. The reason can be traced back to the overbinding property of LDA, which generally provides smaller lattice constants, stiffer bulk moduli, and phonon frequencies compared to experiment and GGA [7,51]. The stiffer LDA system is more resistant to melting and exhibits a higher melting temperature. Similarly as found for Cu [7], we find here the same behavior for Al, i.e., that GGA and LDA provide a confidence interval for the experimental temperature.

Figures 2 and 3 and
Based on our results, it is possible to understand why the prediction of Vočadlo and Alfè [5] (GGA-PW91, 786 K) is about 104 K lower than that of de Wijs et al. [2] (LDA, 890 K). However, both of these values are about 100 K lower than our prediction for the corresponding exchangecorrelation functional, and as a consequence the experimental melting point does not fall within them. We have performed additional calculations in order to understand the origin of this discrepancy. These calculations show that the discrepancy 144108-5  VI. Parametrization of the "ref1" potentials (fit to the DFT solid) for Al and Ni. For the meaning of the parameters we refer to Ref. [6].    arises mainly from a different treatment of the inner core electrons. In Refs. [2] and [5] the ultrasoft pseudopotential method was used, which enables more efficient DFT calculations at the expense of accuracy. In the present work, we have instead employed the PAW method which can be generally considered to provide a higher accuracy because smaller radial cutoff radii around the core can be used and because the exact valence wave functions with all nodes in the core region are reconstructed. For Ni our predictions from GGA (1557 K) and LDA (1884 K) again provide a reliable ab initio confidence interval with respect to the experimental melting point 1728 K [19]. We found that the impact of magnetism on the melting temperature is very small, which supports the conjecture in Ref. [17] that nonmagnetic calculations can reasonably well describe free-energy differences between the solid and liquid phases. If we neglect the magnetic contribution, the temperature window predicted by GGA and LDA becomes slightly smaller, 1570 K for GGA and 1860 K for LDA, i.e., 13 K higher and 24 K lower than the respective values including magnetism. Importantly, each of the phases separately, solid and liquid, receives an appreciable contribution from magnetism to the free energy. The reason behind the small impact on the melting point is compensation: Both solid and liquid exhibit similar local magnetic moments with similar temperature dependencies as shown exemplary for GGA in Fig. 4(a). The similar magnetic moments are a direct consequence of the strong impact of atomic vibrations that smoothen the electronic DOS of both phases similarly as shown in Fig. 4(b). Such a smoothening was observed previously for solid phases [52,53]. Our calculations show that atomic vibrations affect the solid and liquid phase similarly. As a consequence, the magnetic free-energy contributions for both solid and liquid are almost identical, with a difference of only 1.5 meV for GGA and 2.5 meV for LDA at the corresponding melting point. As we will discuss in Sec. III D, these small differences are even further diminished TABLE VII. Parametrization of the "ref2" potentials (fit to the DFT liquid) potentials of Al and Ni. For the meaning of the parameters we refer to Ref. [6].     when the impact of LSFs on the electronic free energy is taken into account. Figure 5(a) summarizes the melting point results for Al and Ni and also for Cu from Ref. [7]. The figure nicely highlights that the concept of the ab initio confidence interval applies to all three investigated elements, i.e., the experimental melting points fall in between the GGA and LDA results. Figure 5 Besides the melting temperature, with TOR-TILD we have also access to other important equilibrium melting properties, e.g., the enthalpy of fusion, the entropy of fusion, and the volume change at the melting point, all of which are difficult to measure experimentally. All computed melting properties are visualized in Figs. 2 and 3 as a function of temperature. In Table VIII we quantify the deviation from experimental values at the respective melting points.
A careful inspection of Table VIII reveals that the concept of the ab initio confidence interval also extends to the enthalpy and entropy of fusion for both elements. To stress it again, this means that the experimental values are likely to be found inside of the range predicted by GGA and LDA, although not necessarily in the middle of the confidence interval. The same observation applies to Cu as found previously in Ref. [7].
Matters are more subtle for the volume change V m at the melting point. Considering separately the absolute volumes of solid V m solid and liquid V m liquid for which a sufficient number of experiments are available, we can clearly see that the   Table VIII values computed from distinct experiments for V m solid and V m liquid (numbers in brackets). We see that some values do fall into the GGA-LDA prediction, others not. Thus, due to the limited experimental situation, we can only speculate that the concept of the ab initio confidence interval is likely to apply as well for V m of Al and Ni. Again, a similar statement applies to Cu [7].

D. Special considerations for magnetism in Ni
Recall that we have applied ferromagnetic DFT calculations in the paramagnetic regime for both solid and liquid nickel, including the magnetic entropy arising from transverse magnetic fluctuations according to Eq. (5). In general LSF contributions, i.e., variations in the magnitude of the magnetic moments, which are only partially captured by the ferromagnetic calculations, are important for an accurate description of the high-temperature paramagnetic state of fcc bulk Ni [22]. As discussed above, LSFs caused by the atomic motion are already included in our scheme by utilizing local magnetic moments computed from the MD runs. In principle LSFs could also have a reverse impact on lattice vibrations. It has, however, been shown that in the special case of fcc Ni, the impact of magnetic excitations on lattice vibrations is negligible [23].
Previous studies have also revealed a negligible impact of LSFs in Ni on stacking fault energies [20,21] and vacancy formation energies [21] due to an effective cancellation. So far, however, only free-energy differences between solid phases of Ni have been considered (e.g., fcc vs hcp [20]), whereas a similar estimation for the liquid structure has been lacking. In the spirit of the previous work on stacking fault energies [20] we therefore provide in the following several supporting calculations that these finite-temperature magnetic effects are unlikely to modify the energy difference between solid and liquid phases and thus the predicted melting properties.
The LSFs could impact the electronic free-energy differences (by modifying the electronic DOS) as well as the explicit transverse magnetic contributions (by modifying the effective local magnetic moments). As discussed, our ferro-magnetic calculations reveal very similar magnetic moments for both the solid and liquid phase. This is one of the main reasons why the magnetic free-energy contributions stemming from the transverse spin degree of freedom within the current treatment cancel each other to a large extent providing only a negligible impact on the melting point. To estimate the impact of LSFs on the local magnetic moments in the paramagnetic state and the corresponding impact on the electronic freeenergy contributions, we computed the electronic DOS and effective local moments of Ni at high temperature (1750 K) in both solid and liquid employing the LSF theory from Ref. [25]. For this we selected from our ferromagnetic MD runs three snapshots each for the solid and liquid structure (108 atom cell) at T = 1750 K. These snapshots were employed in subsequent exact muffin-tin orbitals (EMTO) [37,38] calculations in combination with the recently developed LSF theory [25] as also discussed in Ref. [24]. To better compare the impact we performed the calculations for both solid and liquid at the same volume of V ≈ 12.22 Å 3 /atom.
The impact of LSFs on the electronic DOS is exemplary shown for a representative snapshot of the solid and the liquid phase in Fig. 6. It is found that the LSFs further broaden the electronic DOS of the individual structures and also decrease the DOS near the Fermi level. LSFs thus decrease the explicit electronic free energy for both phases. Inspecting the relevant energy range at the considered temperature it is found that, when comparing the electronic free-energy difference with and without LSFs, the impact is slightly larger on the liquid phase as compared to the solid phase. Considering the total electronic free-energy differences it is found that the inclusion of LSFs stabilizes the electronic free energy of the solid phase by about 2 meV/atom as compared to the liquid phase. On the other hand we find that the LSF stabilized local magnetic moments in the liquid phase are slightly larger as compared to the ones in the solid phase. This is also consistent with our ferromagnetic calculations shown in Fig. 6. The LSF-included magnetic contributions are thus stabilizing the liquid phase by about 2 meV/atom which is almost exactly the value we approximated from the ferromagnetic calculations [see Fig. 4(a)]. Given both contributions together, the explicit LSF contributions to the magnetic free energy and the impact on electronic free-energy differences, we find thus that the impact of LSFs on the solid-liquid phase transition is effectively canceling out.

IV. CONCLUSIONS
We have applied the TOR-TILD method to compute the melting points of Al and magnetic Ni with high accuracy utilizing GGA and LDA as exchange-correlation functionals. To this end, we have extended the TOR-TILD method to be applicable to magnetic materials. We have demonstrated the excellent performance of the method. The largest speed-up factor has been achieved for Ni, for which the optimized reference potential for the liquid phase is about a factor of 30 more efficient than the available literature potential.
Based on the computed melting points for Al and Ni, and taking additionally the previous results for Cu into account, we can draw an important conclusion. Melting points computed with the GGA and LDA functional suggest an ab initio based confidence interval for the prediction of experimental melting points. Further investigations for a larger number and variety of examples are, however, required to further generalize these findings. The GGA melting points provide a lower bound and the LDA melting points an upper bound. As explicitly shown, the deviations in the computed melting points can be traced back to the underbinding and overbinding property of GGA and LDA, respectively. The concept of the ab initio confidence interval also applies to the enthalpy and entropy of fusion. For the volume change at the melting point the situation could not be finally settled due to limited experimental information.
For Ni, for which a full magnetic treatment is at present out of reach, we have in detail analyzed the magnetic contribution at different levels of approximation. The main result is that, even though each of the phases, solid and liquid, is affected rather strongly by magnetism, there is a cancellation of the magnetic impact on the melting properties which derive from free-energy differences. A cancellation was observed previously for solid phases of Ni. Here we could show that the cancellation also applies to differences between solid and liquid phases.
The physical reason for the cancellation is a substantial broadening of the electronic densities of states of all phases due to similarly strong thermal vibrations. The resulting shapes of the densities of states as a function of energy become very similar. Including the impact of longitudinal spin fluctuations further enhances this effect and renders the electronic density of states of fcc and liquid Ni almost indistinguishable.
Finally, it is worth noting that, although TOR-TILD provides computationally efficient access to high accuracy melting properties, its practical application requires at present sophisticated user knowledge and suffers from many technicalities. Therefore, the human work load is a serious bottleneck for a general application of TOR-TILD. To overcome this challenge, we are implementing a fully automatized TOR-TILD simulation protocol into pyiron [54,55], which is designed to automatize complex computational tasks in material science. This development will be published elsewhere.