Metals; Free Energy Functional

In the free energy approach [53,54], the excited states are populated according to the Fermi-Dirac (finite-temperature equilibrium) distribution which is based on the assumption that the electrons ``equilibrate'' more rapidly than the timescale of the nuclear motion. This means that the set of electronic states evolves at a given temperature ``isothermally'' (rather than adiabatically) under the inclusion of incoherent electronic transitions at the nuclei move. Thus, instead of computing the force acting on the nuclei from the electronic ground-state energy it is obtained from the electronic free energy as defined in the canonical ensemble. By allowing such electronic transitions to occur the free energy approach transcends the usual Born-Oppenheimer approximation. However, the approximation of an instantaneous equilibration of the electronic subsystem implies that the electronic structure at a given nuclear configuration $\{ {\bf R}_I \}$ is completely independent from previous configurations along a molecular dynamics trajectory. Due to this assumption the notion "free energy Born-Oppenheimer approximation" was coined in Ref. [55] in a similar context. Certain non-equilibrium situations can also be modelled within the free energy approach by starting off with an initial orbital occupation pattern that does not correspond to any temperature in its thermodynamic meaning, see e.g. Refs. [56,57] for such applications.

The free energy functional as defined in Refs. [53,54] is introduced most elegantly by starting the discussion for the special case of non-interacting Fermions

$$H_{\rm s} =-{1\over 2} \nabla^2 -\sum_{I} {Z_I \over {\vert{\bf R}_I - {\bf r} \vert} }$$ (109)

in a fixed external potential due to a collection of nuclei at positions $\{ {\bf R}_I \}$ . The associated grand partition function and its thermodynamic potential (``grand free energy'') are given by

$$\Xi_{\rm s} (\mu V T ) = {\det}^2 \left( 1+ \exp \left[-\beta \left(H_{\rm s}-\mu \right) \right] \right)$$ (110)

$$\Omega_{\rm s} (\mu V T ) = - k_{\rm B} T \ln \Xi_{\rm s} (\mu V T ) \enspace ,$$ (111)

where $\mu$ is the chemical potential acting on the electrons and the square of the determinant stems from considering the spin-unpolarized special case only. This reduces to the well-known grand potential expression

$$\Omega_{\rm s} (\mu V T ) = - 2 k_{\rm B} T \ln \det \left( 1+ \exp \left[-\beta \left(H_{\rm s}-\mu \right) \right] \right)$$

$$= - 2 k_{\rm B} T \sum_i \ln \left( 1+ \exp \left[-\beta \left(\epsilon_{\rm s}^{(i)}-\mu \right) \right] \right)$$ (112)

for non-interacting spin-1/2 Fermions where $\{ \epsilon_{\rm s}^{(i)}\}$ are the eigenvalues of a one-particle Hamiltonian such as Eq. (109); here the standard identity $\ln \det {\bf M} =$   Tr$\ln {\bf M}$ was invoked for positive definite ${\bf M}$ .

According to thermodynamics the Helmholtz free energy ${\cal F} (NVT)$ associated to Eq. (111) can be obtained from a Legendre transformation of the grand free energy $\Omega (\mu VT)$

$${\cal F}_{\rm s} (NVT) = \Omega_{\rm s} (\mu VT) + \mu N + \sum_{I < J} { {Z_I Z_J} \over {\vert{\bf R}_I - {\bf R}_J \vert} }$$ (113)

by fixing the average number of electrons $N$ and determining $\mu$ from the conventional thermodynamic condition

$$N= - \left( {{\partial \Omega} \over {\partial \mu }} \right)_{VT} \enspace .$$ (114)

In addition, the internuclear Coulomb interactions between the classical nuclei were included at this stage in Eq. (113). Thus, derivatives of the free energy Eq. (113) with respect to ionic positions $-\nabla_I {\cal F}_{\rm s}$ define forces on the nuclei that could be used in a (hypothetical) molecular dynamics scheme using non-interacting electrons.

The interactions between the electrons can be ``switched on'' by resorting to Kohn-Sham density functional theory and the concept of a non-interacting reference system. Thus, instead of using the simple one-particle Hamiltonian Eq. (109) the effective Kohn-Sham Hamiltonian Eq. (25) has to be utilized. As a result, the grand free energy Eq. (110) can be written as

$$\Omega^{\rm KS} (\mu V T ) = - 2k_{\rm B}T \ln \left[ \det \left( 1+ \exp \left[-\beta \left(H^{\rm KS}-\mu \right) \right] \right) \right]$$ (115)

$$H^{\rm KS} = - {1\over 2} \nabla^2 -\sum_{I} {Z_I \over {\vert{\bf R}_I - {\bf... ..._{\rm H} ({\bf r}) + {{\delta \Omega_{\rm xc} [n]} \over {\delta n ({\bf r}) }}$$ (116)

$$H^{\rm KS} \phi_i = \epsilon_i \phi_i$$ (117)

where $\Omega_{\rm xc}$ is the exchange-correlation functional at finite temperature. By virtue of Eq. (112) one can immediately see that $\Omega^{\rm KS}$ is nothing else than the ``Fermi-Dirac weighted sum'' of the bare Kohn-Sham eigenvalues $\{ \epsilon_i\}$ . Whence, this term is the extension to finite temperatures of the "band-structure energy" contribution to the total electronic energy.

In order to obtain the correct total electronic free energy of the interacting electrons the corresponding extra terms (properly generalized to finite temperatures) have to be included in $\Omega^{\rm KS}$ . This finally allows one to write down the generalization of the Helmholtz free energy of the interacting many-electron case

$${\cal F}^{\rm KS} (NVT) = \Omega^{\rm KS} (\mu VT) + \mu \int d{\bf r} \enspace n({\bf r}) + \sum_{I < J} { {Z_I Z_J} \over {\vert{\bf R}_I - {\bf R}_J \vert} }$$

$$- {1\over 2} \int d{\bf r} \enspace V_{\rm H} ({\bf r}) \; n ({\b... ... \enspace {{\delta \Omega_{\rm xc}[n]}\over {\delta n({\bf r})}} \; n ({\bf r})$$ (118)

in the framework of a Kohn-Sham-like formulation. The corresponding one-particle density at the $\Gamma$ -point is given by

$$n ({\bf r}) = \sum_i f_i (\beta ) \> \left\vert \phi_i ({\bf r} )\right\vert^2$$ (119)

$$f_i (\beta ) = \left( 1+ \exp \left[\beta \left(\epsilon_i -\mu \right) \right] \right)^{-1} \enspace ,$$ (120)

where the fractional occupation numbers $\{ f_i \}$ are obtained from the Fermi-Dirac distribution at temperature $T$ in terms of the Kohn-Sham eigenvalues $\{ \epsilon_i\}$ . Finally, ab initio forces can be obtained as usual from the nuclear gradient of ${\cal F}^{\rm KS}$ , which makes molecular dynamics possible.

By construction, the total free energy Eq. (118) reduces to that of the non-interacting toy model Eq. (113) once the electron-electron interaction is switched off. Another useful limit is the ground-state limit $\beta \to \infty$ where the free energy ${\cal F}^{\rm KS} (NVT)$ yields the standard Kohn-Sham total energy expression $E^{\rm KS}$ after invoking the appropriate limit $\Omega_{\rm xc} \to E_{\rm xc}$ as $T\to 0$ . Most importantly, stability analysis [53,54] of Eq. (118) shows that this functional shares the same stationary point as the exact finite-temperature functional due to Mermin [58], see e.g. the textbooks [34,35] for introductions to density functional formalisms at finite temperatures. This implies that the self-consistent density, which defines the stationary point of ${\cal F}^{\rm KS}$ , is identical to the exact one. This analysis reveals furthermore that, unfortunately, this stationary point is not an extremum but a saddle point so that no variational principle and, numerically speaking, no direct minimization algorithms can be applied. For the same reason a Car-Parrinello fictitious dynamics approach to molecular dynamics is not a straightforward option, whereas Born-Oppenheimer dynamics based on diagonalization can be used directly.

The band-structure energy term can be evaluated by diagonalizing the Kohn-Sham Hamiltonian after a suitable ``preconditioning'' [53,54]. Specifically, a second-order Trotter approximation is used

$$\exp \left[ -\beta H^{\rm KS}\right] = \sum_i \exp \left[ - \beta \epsilon_i \right] = \sum_i \rho_{ii}(\beta )$$ (121)

$$= \Biggl( \Biggl\{ \exp \left[-{{\Delta \tau} \over 2} \left( - {1\... ...a^2 \right) \right] \exp \left[-\Delta \tau V^{\rm KS} \left[ n \right] \right]$$

$$\exp \left[-{{\Delta \tau} \over 2} \left( - {1\over 2} \nabla^2 \right) \right] \Biggr\} + {\cal O} \left( \Delta \tau^3 \right) \Biggr)^P$$ (122)

$$\approx \sum_i \left\{ \rho_{ii}(\Delta \tau ) \right\}^P = \sum_i \left\{ \exp \left[ - \Delta \tau \epsilon_i \right] \right\}^P$$ (123)

in order to compute first the diagonal elements $\rho_{ii}(\Delta \tau)$ of the "high-temperature" Boltzmann operator $\rho (\Delta \tau)$ ; here $\Delta \tau = \beta /P$ and $P$ is the Trotter ``time slice''. To this end, the kinetic and potential energies can be conveniently evaluated in reciprocal and real space, respectively, by using the split-operator / FFT technique [59]. The Kohn-Sham eigenvalues $\epsilon_i$ are finally obtained from the density matrix via $\epsilon_i = - (1/\Delta \tau) \ln \rho_{ii} (\Delta \tau)$ . They are used in order to compute the occupation numbers $\{ f_i \}$ , the density $n({\bf r})$ , the band-structure energy $\Omega^{\rm KS}$ , and thus the free energy Eq. (118).

In practice a diagonalization / density-mixing scheme is employed in order to compute the self-consistent density $n({\bf r})$ . A suitably constructed trial input density $n_{\rm in}$ is used in order to compute the potential $V^{\rm KS}[n_{\rm in}]$ . Then the lowest-order approximant to the Boltzmann operator Eq. (123) is diagonalized using an iterative Lanczos-type method. This yields an output density $n_{\rm out}$ and the corresponding free energy ${\cal F}^{\rm KS}[n_{\rm out}]$ . Finally, the densities are mixed and the former steps are iterated until a stationary solution $n_{\rm scf}$ of ${\cal F}^{\rm KS}[n_{\rm scf}]$ is achieved. Of course the most time-consuming part of the calculation is in the iterative diagonalization. In principle this is not required, and it should be possible to compute the output density directly from the Fermi-Dirac density matrix even in a linear scaling scheme [60], thus circumventing the explicit calculation of the Kohn-Sham eigenstates.

As a method, molecular dynamics with the free energy functional is most appropriate to use when the excitation gap is either small, or in cases where the gap might close during a chemical transformation. In the latter case no instabilities are encountered with this approach, which is not true for ground-state ab initio molecular dynamics methods. The price to pay is the quite demanding iterative computation of well-converged forces. Besides allowing such applications with physically relevant excitations this method can also be straightforwardly combined with ${\bf k}$ -point sampling and applied to metals at ``zero'' temperature. In this case, the electronic ``temperature'' is only used as a smearing parameter of the Fermi edge by introducing fractional occupation numbers, which is known to improve greatly the convergence of these ground-state electronic structure calculations [60,61,62,,64,65,66].

Finite-temperature expressions for the exchange-correlation functional $\Omega_{\rm xc}$ are available in the literature. However, for most temperatures of interest the corrections to the ground-state expression are small and it seems justified to use one of the various well-established parameterizations of the exchange-correlation energy $E_{\rm xc}$ at zero temperature.