Efﬁcient approach to compute melting properties fully from ab initio with application to Cu

Applying thermodynamic integration within an ab initio -based free-energy approach is a state-of-the-art method to calculate melting points of materials. However, the high computational cost and the reliance on a good reference system for calculating the liquid free energy have so far hindered a general application. To overcome these challenges, we propose the two-optimized references thermodynamic integration using Langevin dynamics (TOR-TILD) method in this work by extending the two-stage upsampled thermodynamic integration using Langevin dynamics (TU-TILD) method, which has been originally developed to obtain anharmonic free energies of solids, to the calculation of liquid free energies. The core idea of TOR-TILD is to ﬁt two empirical potentials to the energies from density functional theory based molecular dynamics runs for the solid and the liquid phase and to use these potentials as reference systems for thermodynamic integration. Because the empirical potentials closely reproduce the ab initio system in the relevant part of the phase space the convergence of the thermodynamic integration is very rapid. Therefore, the proposed approach improves signiﬁcantly the computational efﬁciency while preserving the required accuracy. As a test case, we apply TOR-TILD to fcc Cu computing not only the melting point but various other melting properties, such as the entropy and enthalpy of fusion and the volume change upon melting. The generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional and the local-density approximation (LDA) are used. Using both functionals gives a reliable ab initio conﬁdence interval for the melting point, the enthalpy of fusion, and entropy of fusion. DOI:


I. INTRODUCTION
Theoretical predictions of the melting point T m of materials have a long history. Since 1910 two empirical rules have been available to determine T m , the Lindemann criterion [1] and the Born criterion [2]. The Lindemann criterion relates melting to a vibrational instability, i.e., melting occurs when the root mean square displacement of the atoms reaches a critical fraction of the nearest-neighbor distance. The Born criterion relates melting to an elastic instability, i.e., melting happens when the shear modulus vanishes. More recently, it became understood that these empirical rules hold true only for a restricted class of systems [3]. For example, in aluminum the elastic shear constants remain finite at T m in contradiction to the Born prediction [4,5]. An alternative to the empirical rules are simulations with classical model potentials for explicitly determining T m . Such model potentials are either fitted to experimental data or to ab initio calculations. The advantage of using classical model potentials is that they can be efficiently applied to large systems and for long simulation times. However, the main disadvantage is that these classical potentials are generally lacking transferability so that the accuracy of the predictions is questionable.
In the past twenty years, with increasing computing power ab initio molecular dynamics (MD) simulations based on density functional theory (DFT) have been used as a general simulation tool for determining T m . In 1995 Sugino and Car [6] for the first time calculated T m of Si from first principles without relying on any adjustable parameters. They used the local-density approximation (LDA) to explicitly calculate the Gibbs energies of the solid and liquid by thermodynamic integration. T m was determined by the crossing point of the Gibbs energies (cf. Fig. 1). Somewhat later, Alfè called this method the free-energy approach (FE), applying it to Fe and Al [7,8]. Since the work of Sugino and Car [6], FE became a standard method for calculating T m and a number of works based on it have followed as indicated in Table I.
While it provides in principle DFT accuracy, FE suffers from two main problems which hindered a more general application. The first one is the high computational cost of the free energy calculations. The free energies of solid and liquid need to be calculated with high precision because they cross at a shallow angle (Fig. 1). The difference in the slopes is the entropy change upon melting, i.e., entropy of fusion, which critically affects the precision of the T m calculations. For example, the entropy of fusion of fcc Cu is about 1.2 k B , which means that an error of 1 meV/atom in either the liquid or solid Gibbs energy can cause an error of ∼10 K in T m . The second problem is that the computational efficiency strongly relies on the reference system used for the thermodynamic integration. FE only works well if one is able to find a reference system for which the free energy can be easily calculated and which is as close as possible to the ab initio system so that the computational effort needed to calculate the free energy difference between the two systems can be reduced to a minimum. As calculations of solid free energies are nowadays mature (Sec. II A), the problem of FE is how to efficiently and accurately calculate the liquid free energy. Up to now, different reference systems have been used for the liquid free energy calculations. Car and Sugina [6] used a Stillinger-Weber potential as the reference system for liquid Si. de Wijs et al. [9] applied the Lennard-Jones fluid as the reference system for liquid Al. Alfè et al. [7,[10][11][12] used inverse-power potentials for liquid Fe and Al. However, these reference systems are limited to special materials. A general approach which is easy and efficient to apply for obtaining a proper reference system is missing.
Because of the problems with FE, researchers started to resort to other methods to determine T m from ab initio MD. The two most commonly used methods are the coexistence approach by Alfè (aka interface method) [7,8] and the Z method [21]. The coexistence approach is based on the simultaneous simulation of solid and liquid in coexistence. Originally the coexistence approach was carried out by using classical model potentials [30]. Alfè has proposed an on-top ab initio correction to obtain the same level of accuracy as with FE. The idea is to start with coexistence simulations for a reference sys-tem, e.g., embedded-atom method (EAM) potential, and then to correct the resulting free energy by a perturbation approach employing DFT calculations [7,13,14,16]. The crucial requirement is that the model potential is reasonably close to the ab initio system, i.e., the accuracy of the final results is influenced by the choice of the model potential. It is worth noting that for FE the reference system only influences the computational efficiency, while the final results are independent of the choice of the reference system (provided that they are statistically converged). To avoid the dependence of the coexistence approach on the model potential, direct DFT calculations with the coexistence approach have been performed for a few cases [15][16][17][18][19]. However, because of the extreme computational effort these calculations were based on very low convergence parameters, e.g., only the point was taken into account to sample the Brillouin zone, or on small system sizes [20]. For general applications a full quantum mechanical simulation of solid and liquid in coexistence seems still out of reach.
The Z method is based on a one-phase approach. It was developed by Belonoshko et al. first by using empirical potentials [31] and later ab initio calculations [21]. The idea is to simulate the solid in the NV E ensemble corresponding to a temperature of the superheated limit. If the simulation is long enough, the solid will melt spontaneously and the temperature will drop back to T m . The corresponding pressure-temperature relation follows approximately the shape of the letter Z, explaining the name of the method. Even though the Z method has been successfully applied to a few systems [22][23][24][25], its strong dependence on simulation size and length are still a matter of debate [27,[32][33][34]. Table I shows that in the last decade there was a strong focus on the coexistence approach and the Z method. More recently, TABLE I. Representative ab initio studies of melting points; FE=the free-energy approach; Coexist+Perturb=the coexistence approach+DFT correction; DFT-Coexist=direct DFT simulation of coexistence; DFT-CoexistSmall=coexistence using small DFT systems; 2PT-MF=two-phase thermodynamic model using a memory function formalism; LDA=the local-density approximation; GGA=the generalized gradient approximation; PW91=the Perdew-Wang GGA functional [28]; PBE=the Perdew-Burke-Ernzerhof GGA functional [29]. Overall, our survey reveals that-after an initial period of active research-FE studies became quite rare, which can be traced back to the aforementioned difficulties. Meanwhile there was, however, some interesting progress in computing anharmonic free energies in solids. An example is the upsampled thermodynamic integration using Langevin dynamics (UP-TILD) method [37] which provides an efficient hierarchical scheme to coarse grain the configuration space during thermodynamic integration from ≈10 7 to ≈10 2 configurations retaining ab initio accuracy for the anharmonic free energy of solids. With this technique the calculation of the anharmonic free energy for solid Al became orders of magnitude more efficient. In 2015 an improved version of UP-TILD, the two-stage upsampled thermodynamic integration using Langevin dynamics (TU-TILD), has been proposed by Duff et al. [38] and applied to calculate the anharmonic free energy of solid ZrC. Within TU-TILD the thermodynamic integration is performed in two stages. The first stage switches the system from the quasiharmonic reference state to a highly optimized empirical interatomic potential. In the second stage the empirical potential is switched to the full DFT system, keeping the upsampling step of the UP-TILD method. The key feature of TU-TILD is that the empirical potential is fitted only to the most relevant part of the phase space (e.g., fcc phase in a certain temperature window) while transferability to other parts of the phase space (e.g., other phases) is not required. With this development DFT accuracy is combined with high efficiency. For the case of ZrC, a speed-up factor of 50 compared to the UP-TILD method was achieved. So far UP-TILD and TU-TILD have been only applied for computing free energies of solids.
Our aim in this paper is to develop a general, accurate and efficient methodology within an ab initio MD framework that will enable the calculation of various melting properties (melting point, enthalpy, and entropy of fusion) for a wide range of elements as required for example in phase diagram databases. The basic idea is to extend the TU-TILD method to liquid free energy calculations combining the accuracy of FE with the efficiency of TU-TILD. Within our methodology a flexible and efficient approach for generating good reference systems close to the ab initio system will be used for both the solid and the liquid, and we thus refer to the new method as the two-optimized references thermodynamic integration using Langevin dynamics (TOR-TILD) method. As an example, the melting point and other melting properties of fcc Cu are calculated with TOR-TILD.

A. Free energy approach
The free energy approach is based on explicitly calculating the Helmholtz free energy F (V ,T ) in an NV T ensemble of both solid and liquid as a function of volume V and temperature T . The reason to start with the Helmholtz free energy as the thermodynamic potential is that the corresponding fixed volume conditions are well amenable to first principles calculations. 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 ) ( Fig. 1). 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 [39].
Within the free energy Born-Oppenheimer approximation [40,41], we can decompose the total free energy of a nonmagnetic solid as: where E 0K (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, F el-vib (V ,T ) the adiabatic electron-phonon coupling, and F vac (V ,T ) the vacancy contribution. The sum of F qh (V ,T ) and F ah (V ,T ) constitutes the full vibrational free energy of the solid. The first three terms of Eq. (2) can be calculated straightforwardly and are understood well [42]. The calculation of the anharmonic free energy is more involved requiring sophisticated techniques [37,38,[43][44][45][46][47][48][49][50].
Here we rely on the TU-TILD method [38] as introduced above. The coupling between electrons and phonons, including anharmonicity, can be captured within the free energy Born-Oppenheimer approximation in an adiabatic fashion (see Refs. [51,52] for details). The vacancy term can be calculated according to our previous developments [37,53], also including anharmonicity and electron-phonon coupling. Note however that the contribution of vacancies to the total free energy is rather small (≈−0.1 meV/atom at the melting point [37,53]). A decomposition as performed in Eq. (2) for the solid free energy is not possible for the case of the liquid, because a static reference lattice is missing. Instead, one needs to fully rely on thermodynamic integration from an appropriate reference system: is the free energy of the reference system and F liquid ref→DFT (V ,T ) the free energy difference between the 224202-3 reference system and the ab initio system, obtained by thermodynamic integration. Alfè [11,12] used for example an inverse power potential as a reference system and calculated F liquid ref of Fe and Al by thermodynamic integration to the Lennard-Jones liquid, for which the free energy has been tabulated [54]. Our approach is different as we will use an optimized empirical potential as a reference system. Further, we will calculate F liquid ref by connecting it to the solid free energy using a specifically designed procedure (Sec. II B). Electronic excitations and their adiabatic coupling to atomic motion within the liquid are included in our approach at the stage of calculating F liquid ref→DFT (V ,T ).

B. TOR-TILD: The improved free energy approach
Our goal is to efficiently and accurately calculate the free energy surface of the liquid, F liquid DFT (V ,T ). We assume that we have already computed an accurate DFT based free energy surface of the solid, F solid DFT (V ,T ), using the standard TU-TILD method [38]. This means that we have an optimized classical reference potential available (EAM or modified EAM), fitted to MD trajectories of the DFT solid. We also have the free energy surface of this potential for the solid phase available, which we call F solid ref1 (V ,T ). The main extension to TOR-TILD is that we have to fit a second optimized reference potential, this time however to MD trajectories of the DFT liquid. Thus, we rely on two optimized references (hence TOR), one optimized for the solid (labelled "ref1") and one optimized for the liquid ("ref2"). Based on these preliminaries, TOR-TILD works as sketched in Fig. 2 and described in the following.
Step 1: We use the coexistence approach [30] to locate T m ref1 , the melting point of the "ref1" potential at constant pressure. Since the reference potential is a classical potential, the calculation of T m ref1 can be performed very efficiently and accurately. In particular we use sufficiently large supercells (several thousands of atoms) and long timescales (several tens of ps) to get a well converged T m ref1 . Since at the (constant pressure) melting point the Gibbs energies of solid and liquid are equal, we can also easily relate the free energies of the solid, F solid ref1 (V m,sol ref1 ,T m ref1 ), and liquid, F Step 2: . We now perform a thermodynamic integration (λ integration) from the "ref1" to the "ref2" potential within the liquid phase at V m,liq ref1 and T m ref1 to obtain: where E liquid ref2 and E liquid ref1 are the potential energies corresponding to the two reference potentials, λ is a coupling constant, and . λ,T denotes the thermodynamic average at a temperature T and coupling λ. The calculations in this step rely only on classical potentials and can be therefore performed with high accuracy and precision. Details on the different steps are given in the main text. Note that for visualization purposes the curvature along the volume axis has been made somewhat smoother in the shown free energy surfaces.
Step 3: From Step 2 we know one point on the liquid free energy surface of the "ref2" potential, i.e., F liquid ref2 (V m,liq ref1 ,T m ref1 ). Starting from this point, we integrate the pressure, P (V ,T ), of the liquid of the "ref2" potential to lower and higher volumes (V integration) using [55] F liquid ref2 and the internal energy, U (T ,V ), of the liquid of the "ref2" potential to lower and higher temperatures (T integration) using [55] where k B is the Boltzmann constant, to get F liquid ref2 (V ,T ) at any required volume and temperature. The calculations in this step also rely only on classical reference potentials and can be therefore performed efficiently and with high accuracy.
Step 4: From F liquid ref2 (V ,T ) we perform a thermodynamic integration to the DFT liquid phase where E liquid DFT-low represents the potential energy of the DFT system. For the DFT calculations we use "low" convergence parameters (k points, cutoff) that are still accurate enough to describe the phase space distribution well. This strategy is in the spirit of the original UP-TILD method [37]. Since the "ref2" potential has been optimized to describe the DFT liquid phase (using specifically low parameters), the thermodynamic integration in Eq. (7) is very efficient.
Step 5: Finally, we correct for the error introduced by the low converged parameters by upsampling to high converged parameters using a perturbative approach: where with the sum running over N uncorrelated snapshots extracted from the DFT-low MD and with E liquid DFT,i and E liquid DFT-low,i corresponding to the DFT potential energies of the ith snapshot calculated using high and low converged parameters, respectively. The upsampling strategy is the same as used in the original UP-TILD and TU-TILD methods [37,38]. From the condition G liquid DFT (T m DFT ) = G solid DFT (T m DFT ) the melting point of the DFT system T m DFT can be obtained.

C. TOR-TILD: Discussion and details
An important advantage of our method is that the calculations requiring only the reference potentials, i.e., steps 1, 2, and 3, can be performed with high accuracy and precision, i.e., in large supercells (several thousand of atoms) and over long time scales (several tens of ps). This allows us to accurately average over all degrees of freedom of the liquid, at least on the level of the reference potential.
For the determination of the melting point of the reference potential using the coexistence approach a special procedure turned out to be necessary. Running a single coexistence calculation resulted in an appreciable error in the melting point of the reference potential which propagated into the final DFT melting temperature. Our strategy to remove this error is as follows: We start with many (in the range of a hundred) coexistence calculations in parallel, with the only difference lying in a different initial random seed. The resulting statistical distribution of melting points closely resembles a Gaussian function and the mean can be used as a precise prediction of the melting point of the reference potential. The corresponding standard error is well below 1 K and can be neglected in the final DFT melting point.
The optimization of the reference potentials is a critical point of the TOR-TILD method. We use the recently developed MEAMfit code for that purpose [56]. This code takes as input MD trajectories and corresponding energies and/or forces to fit classical EAM or reference free modified EAM potentials. Our fitting strategy is to include trajectories of a few hundred MD steps at a few volumes and temperatures in the relevant range and to use the DFT energies as target quantities. The used supercells are in the range of a few tens to a hundred atoms. The resulting datasets are rather large consisting of several thousand geometries and energies, but they can be nevertheless efficiently processed by the MEAMfit code [56].

A. Computational details
For the DFT calculations we used the projector-augmented wave method [57] as implemented in the VASP software package [58][59][60][61] with the PAW Cu potential containing 11 valence electrons. LDA and GGA were employed as exchangecorrelation functionals, with the Perdew-Burke-Ernzerhof (PBE) [29] parametrization for GGA. The set of explicitly DFT computed volume and temperature points for the solid and liquid free energy surfaces is given in Table II. Details on the fitting of the DFT solid free energy surface and the respective convergence parameters will be discussed elsewhere [62]. The vacancy contribution was taken from our study in Ref. [53]. The DFT liquid calculations were performed in a 3 × 3 × 3 supercell with 108 atoms (see Sec. III B for finite size effects). The explicitly DFT computed points were fitted with second order polynomials in both volume and temperature. The plane wave cutoff and k point mesh (Monkhorst-Pack [63]) were set to 300 eV and 2 × 2 × 2 for the low converged calculations. For the high converged calculations, a cutoff of 450 eV and a k point mesh of 6 × 6 × 6 were used. The volumes and temperatures used for fitting the "ref1" and "ref2" potentials are summarized in Table III. We found that to obtain a "ref1" potential which is not only close to the DFT calculations but also can be applied to the coexistence approach we need configurations from several temperatures for fitting. For "ref2" which only serves as an efficient reference for the DFT liquid the highest relevant temperature is sufficient for fitting. The resulting potential parameters are provided in Table IV. We used in particular three pair terms for each of the density contributions and pair potentials. For details on the potential parameters we refer to our previous work [38].
For the reference potential calculations we used the LAMMPS software package [64]. The coexistence approach (Step 1) was performed in a tetragonal 10 × 10 × 20 supercell with 8000 atoms (cf. Fig. 3) resulting in T m ref1 = 1274 K with the "ref1" potential fitted to LDA energies and 1144 K when fitting to PBE energies (experiment: 1358 K). The other reference potential calculations were performed in a cubic 10 × 10 × 10 supercell with 4000 atoms (see Sec. III B for finite size effects). For the MD simulations we used a time step of 5 fs and a Langevin thermostat with a friction parameter of 0.01 fs −1 to control the temperature. For fitting the liquid free energy surface from our potential we used the same volumes as for the DFT calculations (see Table II), but a denser temperature mesh (steps of 5 K).

B. Finite size effects
An important advantage of the TOR-TILD method is that the impact of finite size effects can be easily tested on the optimized potential level. The red line in Fig. 4 shows corresponding results for liquid Cu at 1400 K and one can observe that already a supercell of 6 × 6 × 6 repetitions (in terms of the fcc cubic cell) with 864 atoms is sufficient to obtain a convergence of well below 1 meV/atom. We have also tested the convergence of the full DFT liquid free energy up to a supercell of 4 × 4 × 4 with 256 atoms. The black line in Fig. 4 shows nicely that for the DFT accessible supercells the full DFT free energy closely follows the convergence behavior of the potential (red line). This is an important finding as it means that the difference between the optimized potential and DFT converges rapidly with supercell size, as additionally emphasized by the blue line in the inset of Fig. 4.

C. Computational efficiency
The key point of our method is the highly optimized reference potentials that are fitted to the most relevant part of the DFT phase space; "most relevant" with respect to the thermodynamic integration from the reference potential to the DFT system. A consequence of such a fitting procedure is that our potentials have low transferability to other parts of the phase space, e.g., different structural phases or structural defects that are purposefully not included into the fitting  database. Our fitting strategy stands therefore in clear contrast to the usual construction strategy of empirical potentials where good transferability is an important criterion. Yet, our reference potentials do not need to be transferable, they serve only as an efficient bridge to DFT accuracy in a constrained part of the phase space.
The efficiency of such an optimized reference potential for the solid phase was shown in Ref. [38] within the framework of the original TU-TILD method. Here we concentrate on assessing the efficiency of the reference potential specifically fitted to the DFT liquid system, the "ref2" potential. To quantify the efficiency we investigate the standard deviation with respect to the DFT energies defined as where as needed in Eq. (7) for the thermodynamic integration from the reference potential to DFT. The smaller σ is the less MD steps involving computationally expensive DFT calculations are needed to statistically converge the thermodynamic integration. In fact 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 four, i.e., only 1/4 MD steps are required to get the same standard error.
The standard deviation for our optimized "ref2" potential is shown for the example of LDA by the black line in Fig. 5(a) as a function of the coupling parameter λ for the highest investigated temperature of 1600 K. The standard deviation increases strongly with temperature and investigating σ at such a high temperature constitutes therefore a critical test, regardless of the fact that this temperature was used for fitting. Despite this severe testing condition, the standard deviation for our optimized potential is only about 1.2 meV/atom and independent of λ. This is a very small value and thus the thermal average of the integrand in Eq. (7) converges quickly with the number of MD steps as exemplified in Fig. 5(b) by the black line. If we define a target value of for example 1 meV/atom for the standard error σ n , the required CPU time on a single core is less than one day (∼17 hours) as shown in Fig. 5(c) (i.e., 17 CPU hours).
In order to put the achieved standard deviation and CPU timings into perspective we have investigated two alternative reference potentials. One of these potentials is an EAM potential fitted to DFT MD trajectories of the solid phase, i.e., the "ref1" potential of the previous sections. It was fitted to the DFT solid at two volumes and five temperatures as indicated in Table III. Results for this potential are shown TABLE IV. Parametrization of the "ref1" (fit to the DFT solid) and "ref2" (fit to the DFT liquid) potentials for LDA and PBE. For the meaning of the parameters we refer to Ref. [38].   Example for the efficiency of the present TOR-TILD method. The figure represents results of a thermodynamic integration for the liquid phase at 1600 K from different reference potentials to DFT (low converged values and LDA functional). The "ref1" potential has been fitted to DFT solid energies, "ref2" to DFT liquid energies, and "lit" is a Cu potential available from literature [65]. in Fig. 5 by the blue lines. Although this potential provides reasonable efficiency, it is apparent that it is less efficient than the "ref2" potential, about a factor of three difference in the CPU timing. The other reference potential chosen for the thermodynamic integration to DFT is an EAM potential available from the literature [65] (labeled "lit"). Using this potential the efficiency further decreases (red lines in Fig. 5), about an order of magnitude in the CPU timing difference to our highly optimized reference potential.
We note that our method is fully independent of the availability of literature based empirical potentials for the system under study, the test being done here only for comparison purposes. This point is of practical importance because for other material systems empirical potentials might not be readily available. Table V compile the results of the TOR-TILD approach applied to Cu. Figures 6(a) and 6(b) show the Gibbs energies referenced with respect to the internal energy of fcc Cu at T = 0 K. The melting temperature is determined by the crossing point of the liquid and solid Gibbs energies. Our PBE calculations yield a value of 1251 ± 15 K. This result is in good agreement with that of Vočadlo et al. who predicted the melting temperature of Cu to be 1176 ± 100 K from the phase coexistence approach using GGA-PW91 [13]. Considering our LDA calculations, we find a substantially higher melting temperature of 1494 ± 5 K. This finding can be rationalized by recalling that LDA gives smaller lattice constants, stiffer bulk moduli, and phonon frequencies [42] due to its overbinding property. The stiffer LDA system is more resistant to melting than the softer GGA system. The stiffness of LDA also explains why the melting point of the LDA system is less sensitive to the statistical error in the MD simulations (smaller error bar for LDA than PBE in Table V).

Figure 6 and
Comparison of the DFT results to the experimental melting temperature of 1358 K [66] leads us to an interesting conclusion. Neither one of the functionals performs better: PBE underestimates by −107 K and LDA overestimates by 136 K. Thus, the experimental value resides in between the two functionals, by coincidence quite close to the middle of the interval. A similar behavior was observed before for properties of the solid phase of Cu and several other fcc elements [42]. The conclusion of Ref. [42] was that, although one cannot clearly assign a better performance to GGA nor LDA, using both functionals for the calculations provides an ab initio based confidence interval for the prediction of experimental data. Based on the present results, we expect that LDA and GGA provide also an approximate confidence interval for the melting temperature.
The virtue of the proposed approach is that it not only provides access to the melting temperature but to other      is observed from the solid to the liquid [Figs. 6(g) and 6(h)]. Measuring this volume increase is difficult but easily accessible within our approach. PBE gives a larger increase of 0.86Å 3 /atom and LDA a smaller one of 0.66Å 3 /atom as might have been anticipated due to the overall stiffer response of LDA. Considering the absolute volumes of solid and liquid copper at the melting point for which distinct experiments are available [67][68][69], we observe that all of the experiments fall in between our PBE and LDA results (see Table V). However, for the differences between the absolute volumes of solid and liquid this cannot be expected since it is highly sensitive to even small errors in the absolute volumes. Indeed, the volume increase from experiment is not inside the LDA/PBE bound but still very close (see Table V). The Gibbs energy of vacancy formation, G f , in fcc Cu was intensively studied in Ref. [53] including the influence of anharmonicity. We have used the data from Ref. [53] to derive the vacancy contribution to G solid DFT using [70]: Figure 7(a) shows the results for G vac (T ) for PBE and LDA. At the respective melting point, the vacancy contribution to the Gibbs energy is negligibly small (≈−0.03 meV/atom) for both functionals. Figure 7(b) shows the corresponding vacancy contribution to the entropy of the fcc bulk, and the conclusion is the same: There is only a negligible contribution (≈0.003 k B /atom) at the respective melting point of the two functionals. To provide an idea about the dependence of the vacancy contribution above the melting point, we have extrapolated the data from Ref. [53] using a third-order polynomial (dashed lines in Fig. 7). Even for PBE, for which the extrapolation is larger in the shown temperature window, the contribution of the vacancies to bulk properties remains very small. The vacancy concentration can be extracted from the Gibbs energy of vacancy formation [53]. At the experimental melting point (1358 K), the vacancy concentrations are 8.9 × 10 −4 for PBE (which agrees well with the experimental result of 7.6 ± 0.3 × 10 −4 [71,72]) and 6 × 10 −5 for LDA, i.e., significantly smaller for LDA. However, at the respective theoretical melting points of PBE (1251 K) and LDA (1494 K), the vacancy concentrations are similar with 3.4 × 10 −4 for PBE and 1.8 × 10 −4 for LDA. Thus, similarly as for the vacancy contribution to the Gibbs energy and entropy discussed above, we obtain a consistent result for the vacancy concentration for both exchange-correlation functionals at their respective theoretical melting points.

IV. CONCLUSIONS
We have developed a general, accurate, and highly efficient methodology, the two-optimized references thermodynamic integration using Langevin dynamics (TOR-TILD) method, to fully ab initio calculate melting properties of materials, such as the melting point, entropy of fusion, enthalpy of fusion, volume change, and vacancy concentration at the melting point. By extending the TU-TILD method from solid to liquid free energy calculations we have overcome the two main challenges of the standard FE approach, high computational cost and reliance on a good reference system. The key aspects of the new approach are tailored classical potentials fitted to the relevant part of the DFT phase space. Transferability is lost, but also not needed, because the potentials serve only as highly optimized references for the thermodynamic integration to DFT.
In the present work, TOR-TILD has been successfully applied to fcc Cu using the two standard exchange-correlation functionals PBE and LDA. Using both functionals provides a reliable ab initio confidence interval for the experimental melting properties.
Our methodology is not limited to calculate the melting properties of nonmagnetic phases. For magnetic materials, Eqs. (2) and (3) need to be extended by terms related to magnetic excitations. Corresponding calculations are nontrivial and under intensive research presently [73][74][75][76]. Assuming that close to the melting temperature the system is sufficiently above its Curie or Néel temperature, one could for example use the high temperature limit for the magnetic entropy [73,74] as an approximation. Additionally, the DFT calculations entering the other contributions in Eqs. (2) and (3) would require the spin-polarized extension of DFT and a proper representation of magnetic disorder, using for example special quasirandom structures in the spin space [76][77][78]. Given a method to treat magnetic excitations, our approach can be applied to magnetic materials.
To extend our methodology to binaries and higher order systems, additional terms related to configurational entropy might enter Eqs. (2) and (3) depending on the degree of disorder. Finally, our methodology is not limited to calculate the melting properties of stable phases, but can be likewise applied to metastable phases, which is desirable for updating, completing, and unifying the databases of CALPHAD based approaches.