Temperature dependent effective potential method for accurate free energy calculations of solids

We have developed a thorough and accurate method of determining anharmonic free energies, the temperature dependent effective potential technique (TDEP). It is based on ab initio molecular dynamics followed by a mapping onto a model Hamiltonian that describes the lattice dynamics. The formalism and the numerical aspects of the technique are described in detail. A number of practical examples are given, and results are presented, which confirm the usefulness of TDEP within ab initio and classical molecular dynamics frameworks. In particular, we examine from first principles the behavior of force constants upon the dynamical stabilization of the body centered phase of Zr, and show that they become more localized. We also calculate the phase diagram for ${}^4$He modeled with the Aziz et al. potential and obtain results which are in favorable agreement both with respect to experiment and established techniques.

Figure 1

Illustration of the different terms included in free energy calculations using different approaches. The solid boxes denote the terms that are included, and the striped areas what is omitted. The main message is to point out that the internal energy has a nontrivial temperature dependence, something that is omitted in the quasiharmonic approximation. HOT indicate the higher order terms that are missing in the TDEP free energy, Eqs. (22) and (18).

Figure 2

Convergence of the free energy correction from thermodynamic integration [$\Delta F_{\text{AH}}$ in Eq. (23)] with respect to system size. At system sizes smaller than $\sim$700 atoms the uncertainty is about the same order as the correction. This particular case is for fcc Cu modeled with an embedded atom potential (Ref. 20).

Figure 3

Convergence of the free energy from TDEP with respect to simulation size. The system in question is fcc Cu with a classical embedded atom potential (Ref. 20). At sizes of about 100 atoms it is converged within 1 meV/atom.

Figure 4

Convergence of the free energy of bcc Zr at 1300 K with respect to the number of time steps. In (a) symmetry is treated numerically and in (b) it is treated in the different, analytical way. The same input data is used in both cases, and it converges to the same value, but we have improved the performance by several orders of magnitude.

Figure 5

Comparison of performance of TDEP and the harmonic Taylor expansion for the potential described by Eq. (24). Three examples are shown when the conventional harmonic approximation fails to describe the potential while TDEP succeeds. (a), (c), and (e) show the potentials and (b), (d), and (f) show the forces. $\alpha$ and $x_0$ are the parameters of Eq. (24). In (a) the harmonic approximation corresponds to a dynamically unstable system, whereas TDEP provides a dynamically stable solution. In (b) the harmonic approximation provides an inaccurate potential. (c) shows how TDEP finds the high temperature equilibrium position $x_0$.

Figure 6

Comparison of the effects of including higher order terms between the harmonic approximation and TDEP. When extending the Taylor expansion of an anharmonic potential (dashed blue line) in the Born–von Karman ansatz to higher order terms, we end up with the series of lines depicted in (b). (a) shows the same extension for TDEP. Even limiting oneself to the second order term ($n=2$), the fit will implicitly contain anharmonism to an arbitrary degree in the range that is thermally accessible. Extending TDEP to a higher order converges towards the true potential faster than when higher order terms are added to the harmonic approximation.

Figure 7

Comparison of the TDEP and quasiharmonic force constant matrices. We have plotted the Frobenius norm of the force constants vs coordination shell. In the inset we show the corresponding phonon dispersion relations. The open circles are 0 K harmonic values and the solid circles are TDEP values extracted at 1300 K. At high temperature the interactions at close distances are stronger, and fall off faster with increasing distance.

Figure 8

Equation of state for bcc Zr calculated with the TDEP method at temperatures 500 K (dashed line), 1000 K (dotted-dashed line), and 1500 K (dotted line). Equation of state obtained by means of conventional DFT calculations at $T=0$ K is also shown with a solid line, and the zero temperature equilibrium volume $V_0=22.83$ Å${}^3$ is chosen as a reference point at all the temperatures. Note that bcc Zr is dynamically unstable at $T=0$ K (see the bottom inset in Fig. 7).

Figure 9

The calculated phase diagram for ${}^4$He modeled with the Aziz et al. potential (Ref. 22). The red line indicates the experimental melting curve. The observation of the stabilization of the bcc phase before the melting demonstrates that TDEP treats this system accurately and in agreement with other approaches, even in such a strongly anharmonic system as ${}^4$He.