Temperature dependent effective potential method for accurate free energy calculations of solids
Abstract
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 details. 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 firstprinciples the behavior of force constants upon the dynamical stabilization of body centered phase of Zr, and show that they become more localized. We also calculate phase diagram for He modeled with the Aziz et al. potential and obtain results which are in favorable agreement both with respect to experiment and established techniques.
I Introduction
One of the common goals for first principles calculations is the comparison of energies, such as configurational energies, surface energies, mixing enthalpies or lattice stabilities. Usually the comparison is limited to total energies. This is appropriate when the effects of temperature can be neglected. However, for many problems within physics, materials and Earth sciences they are not negligible and Gibbs free energy represents the proper thermodynamic potential when the temperature and pressure are external parameters. The quasiharmonic approximation can bridge the temperature gap, but there are cases where it falls short. Strongly anharmonic systems are not described wellGrimvall et al. (2012), especially dynamically unstable systems. Traditionally the problem of dynamical instability was addressed either by including more terms in the Taylor expansion of the potential energy or via a selfconsistent approach.Klein and Horton (1972); Born (1951); Born and Huang (1964)
HootonHooton (1955) realised that even though the second derivatives at the equilibrium positions are negative, the atoms in a solid rarely occupy these positions. They move in the effective potential of their nonstationary neighbours. By sampling the potential energy surface not at the equilibrium positions but at the most probable positions for a given temperature one can obtain a harmonic approximation that describes the system at elevated temperatures. The selfconsistent formalism employs an iterative procedureGillis et al. (1968) by creating a harmonic potential, which is used to describe the thermal excitations that again give a new harmonic potential. This is then repeated until selfconsistency.
The doubletime Green’s functions, developed by Choquard,Choquard (1967) use a cumulant expansion in the higher order terms. Although formally exact, this formalism requires knowledge of the higher order force constants. Obtaining these accurately for something other than the structurally simple systems is computationally very demanding from first principles.Chaput et al. (2011)
A recent implementation of the selfconsistent formalism by Souvatzis et al.Souvatzis et al. (2008, 2011) uses ab initio supercell calculations. A problem with this approach is that the excitations could only be done in the harmonic sense which means probing phase space with a limited basis set. Where the harmonic approximation works well this is not a problem. When the harmonic approximation fails it is due to a strong anharmonic contribution. Strongly anharmonic systems are by definition badly described with the harmonic approximation.
Several techniques use BornOppenheimer molecular dynamics to obtain anharmonic corrections to quasiharmonic free energies.Grabowski et al. (2009); Wu and Wentzcovitch (2009); Wu (2010) They focus on anharmonic corrections to materials that are well described in the quasiharmonic approximation, and the applicability to strongly anharmonic systems is questionable when the phonon renormalization due to increased temperature can not be described by a linear equation.
We have developed a new methodHellman et al. (2011) that is similar to HootonsHooton (1955) original idea, but with a foundation in (ab initio) molecular dynamics. In this paper we present a substantial refinement and generalization of the temperature dependent effective potential method (TDEP), showing how it deals with: a model onedimensional anharmonic oscillator, a strongly anharmonic system, bcc Zr, treated from firstprinciples, and He modelled with the Aziz et al. potential.
Ii TDEP formalism
The starting point of our method is to introduce a model Hamiltonian for a Bravais lattice in the harmonic form:
(1) 
which describes the lattice dynamics. Here and are the momentum and displacement of atom , bold symbols indicate vectors and doubly overlined symbols matrices respectively. The reference point for the displacements is the 0K relaxed lattice (initially, in Sec. IV we will revisit this). The interatomic force constants and the ground state energy are yet to be determined. Given atoms in this model system the forces acting on the atoms are given by
(2) 
As detailed in our previous paperHellman et al. (2011) we seek to determine the force constant matrices through minimization of the difference in forces from the model system a real system, simulated, for instance, by means of ab initio molecular dynamics (AIMD). An AIMD simulation will provide a set of displacements , forces and potential energies . We seek to minimize the difference in forces from AIMD and our harmonic form () at time step , summed over all time steps :
(3) 
This is realized with a a MoorePenrose pseudoinverse
(4) 
to obtain the linear least squares solution for that minimize . This is then mapped to the form
(5) 
where are indices to atoms in a unit cell with atoms and Cartesian indices. The pair vectors in the supercell are mapped to stars of lattice vectors connecting atoms of type and . From this form the phonon dispersion relations, free energy and all other quantities can be extracted. This direct implementation works well,Hellman et al. (2011) but the numerical efficiency can be improved, as is demonstrated below.
Iii Symmetry constrained force constant extraction
The form of the force constant matrices depends only on the supercell used in the AIMD and the crystal lattice. We begin by populating the force constant matrices with unknown variables
(6) 
and so on, including vectors up to a cutoff determined by the size of the supercell. Some of these are equivalent. To find out which, we look at the symmetry relations the force constant matrices have to obey. We haveMaradudin and Vosko (1968)
(7)  
(8)  
(9) 
that stem, in order, from the facts that there is no net translation of the crystal, all Bravais lattices have inversion symmetry, and that the order of the second derivatives does not matter. Each relation will give us a few equations for the unknowns , reducing their number. In addition to these fundamental properties of the force constant matrices we can benefit from the symmetry of the lattice. If symmetry relation , belonging to the point group of the lattice, relates two vectors we have the following relation:
(10) 
By applying equations (7)–(10) the number of unknown variables is massively reduced. For example, a bcc lattice modelled as a supercell (128 atoms) would have 147456 unknown variables in , if one does not consider symmetry arguments. Application of Eqs. (7)–(10) reduces the problem to 11 unknown variables. Having found the reduced problem with unknown variables, it can be substituted back into (2). The expression for the forces at timestep will schematically look like this:
(11) 
The actual distribution of the will depend on the problem at hand. Carrying out the matrix product gives us a new set of equations for the forces
(12) 
where second sum describes the coefficients for each contained in the expression for force component . These coefficients are linear combinations of the displacement components . The explicit form is determined by the lattice. In matrix form this is written as
(13) 
Eq. 13 is equivalent to Eq. 2, it is just rewritten in terms of the symmetry inequivalent interactions. This implementation symbolically reduces the number of unknowns, generates the function that gives the matrix from a set of displacements and the mapping from the set of back to the force constant matrix .^{1}^{1}1In practice, this is implemented in Mathematica
To find we minimize the difference in forces between AIMD simulations and the model hamiltonian
(14) 
This is again realised with the least squares solution for
(15) 
Having determined we can substitute the components into the force constant matrix and proceed to calculate thermodynamic properties of the original (real) system.
The suggested scheme is a superior way of using symmetry to improve the numerical accuracy. Most schemes involving symmetry revolve around determining the interaction in one direction and then using the symmetry relations to translate and rotate that interaction. Any numerical errors–which are always present–will propagate to all interactions, whereas in the present approach the errors will be averaged and should cancel each other to a large degree.
The advantage of this procedure will be illustrated below.
Iv Internal degrees of freedom
If the system has internal degrees of freedom for the structural relaxations, such as a crystal with point defects, an interface or a random alloy the atoms ideal positions and equilibrium positions do not coincide. While one could find the relaxed positions from 0K calculations, the equilibrium positions are by no means constant with respect to temperature. TDEP handles this in an elegant way. Note that we find the second order terms in (1) with a least squares fit of the slope of force versus displacement. Originally, the displacements could be calculated with respect to the wrong equilibrium positions that do not correspond to the temperature of the simulation. Still our experience shows that slope will be the well approximated. That allows for the following procedure:

Guess equilibrium positions, usually the ideal lattice positions.

Use these to calculate the displacements from the AIMD simulations.

Determine and from that the residual force
(16) where are the AIMD forces and are given by equation 13. is the number of timesteps, and subscript denote the corresponding forces at time step .

These forces are then used to move the atoms in a steepest descent scheme towards equilibrium positions. The whole process is repeated until convergence.
Our test shows that this procedure is remarkably stable. The second order force constants are weakly affected by the choice of equilibrium positions. The vibrational entropy and phonon dispersion relations are largely unaffected as well. Eliminating the first order term, however, is formally important and crucial when extracting higher order terms. Note that selfconsistent iterations are numerically efficient, because the most timeconsuming step for applications of TDEP is MD simulations, while their mapping on model Hamiltonian (1) represents postprocessing of the MD results with minimal computaional cost.
V Determining the free energy
We will begin by reiterating the traditional way how free energy is determined in the quasiharmonic approximation. Divided into parts it will be
(17) 
where a division of the right hand side into parts makes a clear distinction between the electronic contribution and the vibrational contribution . The electronic contribution is divided into the total energy of the lattice, , and the electronic entropy . The vibrational contribution is divided into average kinetic and potential energy of the ions, vibrational entropy and zero point energy . The lattice contribution is obtained from DFT calculations with the Mermin functional and the vibrational part from the harmonic approximation via
(18) 
Where is the phonon density of states. In this approach, all the vibrational contributions are calculated within the harmonic approximation.
Turning to AIMD, the free energy (in the canonical ensemble) is divided as
(19) 
were the potential energy is temperature dependent. Since the ions move as classical particles the zero point energy is missing. There is unfortunately no information about the entropy, but through the force constant matrices obtained using TDEP, the vibrational entropy and zero point energy can be estimated. For TDEP to have an accurate free energy the potential energy should, on average, be equal to that of AIMD. This would ensure that the full anharmonic is included. The problem is that is rapidly oscillating over time, requiring a long simulation time to converge. If we look at the TDEP potential energy
(20) 
and recognise that it should model the thermal fluctuations of the original system we can overcome the numerical issues. Setting the average potential energies equal, , gives us
(21) 
By removing the thermal excitations of the harmonic form the fluctuations can be decreased by roughly an order of magnitude and the accuracy of is thus increased be the same amount. Including higher order terms in the energy expansion would further serve to minimize these fluctuations.
The potential energy that was removed will be added again when the Helmholtz free energy is calculated:
(22) 
where is the phonon contribution given Eq. 18, with the phonon density of states determined in the TDEP formalism. It includes the kinetic and potential energy of the ions. In Fig. 2 we further illustrate the difference in methods of obtaining the free energy.
The formally exact method of thermodynamic integrationFrenkel and Smit (2002) can be used to determine the free energy. This method will determine the anharmonic correction to the free energy. If the TDEP model Hamiltonian is used as the reference point the full free energy will be
(23) 
The integral is over the Kirkwood coupling parameter , and the potential energy difference is between the model Hamiltonian and the molecular dynamics potential energy.
The model TDEP Hamiltonian is constructed to describe the system as accurately as possible while retaining the harmonic form. It is then easy to argue that the anharmonic correction term should be small. While it is difficult to make general statements regarding this, our experience so far is that this correction is very small in every system we have tested.
In addition, the thermodynamic integration technique can be numerically inefficient when high accuracy is needed. While one can accurately control the numerical accuracyGrabowski et al. (2009) the finite size effects are more difficult to control, especially in ab intio simulations. In Fig. 3 we show that the error due to the limited size is on the same order as the correction to the TDEP free energy in reasonable simulation sizes^{2}^{2}2The thermodynamic integrations were carried out from a TDEP potential extracted for fcc Cu modelled with an embedded atom potentialFinnis and Sinclair (1984). A Langevin thermostat was used to control the temperature and break the modelocking. The numerical integration over coupling parameter was carried out over 15 discrete steps, and the we ran the MD simulations for 300000 timesteps to ensure convergence within 0.1meV/atom.. It makes little use to add a correction to the TDEP free energy where the uncertainty is of the same order of magnitude as the correction itself.
The TDEP free energy, on the other hand, behaves well with respect to limited simulation cell. In figure 4 we see that at the reasonable system size of 100 atoms the free energy is converged within 1meV/atom.^{3}^{3}3The convergenve tests for TDEP used an an embedded atom potentialFinnis and Sinclair (1984) for fcc Co. We used A Langevin thermostat and ran the MD simulations for 150000 timesteps to ensure convergence within 0.1meV/atom. It is also easily converged in terms of simulation length: in Fig. 5 we illustrate the advantage of analytically treating symmetry. In our previous work,Hellman et al. (2011) we studied Zr in the bcc phase. There, we found convergence within 1meV/atom for the free energies after 25000 time steps. With the symmetry constraints, we converged to the same value using 50 time steps, an improvement by several orders of magnitude.
Vi Application of TDEP to a onedimensional anharmonic oscillator
To illustrate TDEP we first apply it to a onedimensional anharmonic oscillator. Consider the following onedimensional potential:
(24) 
Here is the position and , and are known parameters. The equilibrium position can depend on temperature and is assumed to be 0 at . The aim is to find the second degree polynomial fit to Eq. 24 that best describes the system. If this polynomial only consist of a quadratic and a constant term it will describe a harmonic oscillator with well defined free energy. If one applies the harmonic approximation to this potential it will not work well. The second derivative
(25) 
will determine the force constant . The temperature dependence of is omitted and we will end up with
(26) 
This, as seen in Fig. 6, will not be a particularly good model for the true potential. These issues arise from the fact that the potential energy surface is only probed at , the equilibrium positions. To work around this problem, let us apply TDEP and put a particle in the potential given by equation (24) and perform a MD simulation. When controlled by an appropriate thermostat, the particle will yield a set of forces, positions and energies, , one for each timestep. This data can now be used to fit a potential of the form
(27) 
a harmonic potential centered at . Let us begin by determining the second order term. As discussed in Sec. IV we guess a value for , use the forces from molecular dynamics and minimize
(28) 
This is easiest realized as a least squares fit of a straight line in forces, as demonstrated in the right panels of Fig. 6. Equation (28) determines the second order term. The residual force at , , can be used to find the equilibrium position. It is done in the following manner: a guess for gives us a and . This residual force is used to move to a new position, and the process is repeated until selfconsistency is reached. When we have found the equilibrium position we can safely assume that any first order term in our polynomial can be set to 0. As described in Sec. V The constant energy term, , can be determined from the energies obtained from molecular dynamics simulations:
(29) 
This is the best possible potential of the harmonic form at a given temperature that approximates the original potential in Eq. 24. In Fig. 6 the true potential and the fit are illustrated for different , and . The anharmonism of the potential is implicitly described by the polynomial fit. In Fig. 7 the expansion in Eq. 27 has been extended to higher orders for an anharmonic potential. TDEP, probing the effective potential at finite temperature, converges to the true potential rapidly whereas including more terms in the Taylor expansion in Eq. (1) does by no means guarantee numerical stability at finite temperature.
Vii Practical applications of TDEP
Let us summarize this and present the scheme used to calculate accurate Gibbs free energy surfaces in the TDEP formalism from first principles.

First calculate DFT total energy as a function of volume. This provides an equation of state and allows us to choose the volume interval, which covers pressure of interest.

If feasible, obtain approximate harmonic potentials for the systems at hand in this volume interval. These potentials are used to speed up the calculations as described in Steneteg et al. (submitted to PrB, BZ11707).

On the grid of volumes and temperatures, perform AIMD simulations in the canonical ensemble.

To increase further the accuracy of the calculation, it is recommended to select a subset of uncorrelated samples from the AIMD simulations upsample these to high accuracy, as described in Ref. Grabowski et al., 2009. A new free energy is calculated.

The equation of state is interpolated over the grid of temperatures and volumes providing the Gibbs free energy surface. This is then repeated for each structure, compound or composition of interest.
TDEP is a thorough and timeconsuming method, but the results are excellent. The phonon dispersion relations of a material that is dynamically unstable at zero temperature is a good example. When reevaluating the results for Zr obtained in Ref. Hellman et al., 2011, we observe a striking difference in harmonic and TDEP force constants. This is illustrated in Fig. 8. The effective TDEP force constants decrease faster with distance compared to the harmonic ones, a behavior that is expected. It is a vivid illustration of the temperature dependence of the potential energy landscape, and at the same time a confirmation that the TDEP technique describes this renormalisation well. From the free energy surface we can extract the finite temperature equation of state for bcc Zr, as illustrated in Fig. 9.
To test the performance of TDEP close to melting, we turn to solid He modelled with the Aziz et al. potentialAziz et al. (1979)^{4}^{4}4He was modelled using the Aziz et al. potential.Aziz et al. (1979) The bcc and fcc phases were both described by supercells (125 atoms), a 0.5fs timestep for a total simulation length of 20ps. Free energies were converged within 1meV/atom after about 5ps. Free energy surfaces were interpolated from a grid of six volumes and five temperatures in the appropriate range.. The melting curve and bccfcc transition just before melting has been extensively studied, see Ref. Belonoshko et al., 2012 and references therein. As demonstrated in Fig. 10 strongly anharmonic He pose no problem for the presented method. The stabilization of the bcc phase before melting is consistent with results from phasecoexistence simulations. This are at the moment considered the most accurate methods for determining phase stabilities at high temperature. They do, however, require simulation cells much larger than what is accessible to AIMD, and can only be used with classical potentials. We show here that with the presented method we can with simulation sizes of 125 atoms accurately reproduce the same transition temperatures. This verifies the accuracy of the method and opens up the applicability to high pressure high temperature studies of phase stabilities close to melting.
Viii Conclusions
We have presented detailed description of the Temperature Dependent Effective Potential method for the treatment of lattice dynamics of strongly anharmonic solids, including an extension and refinement to this accurate technique. Moreover, we have detailed how the temperature dependence of all components of the free energy should be taken into account, and presented several successful examples, including a model anharmonic potential, firstprinciples calculations of equation of state for bcc Zr, and classical moleculardynamics simulations of bcctofcc transition in He.
Ix Acknowledgements
Support from the Knut & Alice Wallenberg Foundation (KAW) project âIsotopic Control for Ultimate Material Propertiesâ, the Swedish Research Council (VR) projects 62120114426 and LiLiNFM, and the Swedish Foundation for Strategic Research (SSF) program SRL100026 is gratefully acknowledged. SSI acknowledges support from the Swedish Government Strategic Research Area Grant in Materials Science to AFM research environment at LiU. Supercomputer resources were provided by the Swedish National Infrastructure for Computing (SNIC). We would also like to thank Prof. Anatoly Belonoshko for suggesting the investigation of He phase transitions.
References
 Grimvall et al. (2012) G. Grimvall, B. MagyariKöpe, V. OzoliÅš, and K. Persson, Reviews of Modern Physics 84, 945 (2012).
 Klein and Horton (1972) M. L. Klein and G. K. Horton, Journal of Low Temperature Physics 9, 151 (1972).
 Born (1951) M. Born, Festschr. Akad. Wiss. Göttingen 1951, Math.Phys. Kl., 116 , 1 (1951).
 Born and Huang (1964) M. Born and K. Huang, Dynamical theory of crystal lattices (Oxford University Press, Oxford, 1964).
 Hooton (1955) D. J. Hooton, Zeitschrift für Physik 142, 42 (1955).
 Gillis et al. (1968) N. Gillis, N. Werthamer, and T. Koehler, Physical Review 165, 951 (1968).
 Choquard (1967) P. Choquard, The anharmonic crystal (W. A. Benjamin, INC., New York, 1967).
 Chaput et al. (2011) L. Chaput, A. Togo, I. Tanaka, and G. Hug, Physical Review B 84, 094302 (2011).
 Souvatzis et al. (2008) P. Souvatzis, O. Eriksson, M. I. Katsnelson, and S. Rudin, Physical Review Letters 100, 95901 (2008).
 Souvatzis et al. (2011) P. Souvatzis, S. Arapan, O. Eriksson, and M. I. Katsnelson, EPL (Europhysics Letters) 96, 66006 (2011).
 Wu and Wentzcovitch (2009) Z. Wu and R. Wentzcovitch, Physical Review B 79, 104304 (2009).
 Wu (2010) Z. Wu, Physical Review B 81, 172301 (2010).
 Hellman et al. (2011) O. Hellman, I. A. Abrikosov, and S. I. Simak, Physical Review B 84, 180301 (2011).
 Maradudin and Vosko (1968) A. A. Maradudin and S. Vosko, Reviews of Modern Physics 40, 1 (1968).
 (15) In practice, this is implemented in Mathematica.
 (16) Modelled as a 128 atom supercell using VASPKresse (1999); Kresse and Furthmüller (1996); Kresse and Hafner (1993); Kresse (1996). We used a 1 fs timestep and ran the simulations using the point for BZ integration for 30000 timesteps. A random subset of these configurations where then upsampled using a kpoint grid and increased cutoff.
 Frenkel and Smit (2002) D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications, Computational science (Academic Press, 2002).
 Grabowski et al. (2009) B. Grabowski, L. Ismer, T. Hickel, and J. Neugebauer, Physical Review B 79, 134106 (2009).
 (19) The thermodynamic integrations were carried out from a TDEP potential extracted for fcc Cu modelled with an embedded atom potentialFinnis and Sinclair (1984). A Langevin thermostat was used to control the temperature and break the modelocking. The numerical integration over coupling parameter was carried out over 15 discrete steps, and the we ran the MD simulations for 300000 timesteps to ensure convergence within 0.1meV/atom.
 (20) The convergenve tests for TDEP used an an embedded atom potentialFinnis and Sinclair (1984) for fcc Co. We used A Langevin thermostat and ran the MD simulations for 150000 timesteps to ensure convergence within 0.1meV/atom.
 Finnis and Sinclair (1984) M. W. Finnis and J. E. Sinclair, Philosophical Magazine A 50, 45 (1984).
 Aziz et al. (1979) R. a. Aziz, V. P. S. Nain, J. S. Carley, W. L. Taylor, and G. T. McConville, The Journal of Chemical Physics 70, 4330 (1979).
 (23) He was modelled using the Aziz potential.Aziz et al. (1979) The bcc and fcc phases were both described by supercells (125 atoms), a 0.5fs timestep for a total simulation length of 20ps. Free energies were converged within 1meV/atom after about 5ps. Free energy surfaces were interpolated from a grid of six volumes and five temperatures in the appropriate range.
 Belonoshko et al. (2012) A. B. Belonoshko, L. Koči, and A. Rosengren, Physical Review B 85, 012503 (2012).
 Kresse (1999) G. Kresse, Physical Review B 59, 1758 (1999).
 Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Physical Review B 54, 11169 (1996).
 Kresse and Hafner (1993) G. Kresse and J. Hafner, Physical Review B 48, 13115 (1993).
 Kresse (1996) G. Kresse, Computational Materials Science 6, 15 (1996).