\tolerance=10000
\documentstyle[preprint,aps]{revtex}
\begin{document}
\draft

\title{Free-energy calculations in solids from first-principles
molecular dynamics: vacancy formation in sodium}

\author{Enrico Smargiassi\cite{present} and Paul A. Madden}

\address{Physical Chemistry Laboratory,
University of Oxford, Oxford OX1 3QZ (GB)}

\maketitle

\begin{abstract}
 The application of thermodynamic integration methods to
  free-energy calculations via {\sl ab-initio} computer simulations
 of condensed matter
systems is discussed, with particular reference to the 
 calculation  of the formation free-energy of
defects and of the absolute free-energy of crystals.
  We find that
 specific techniques have
to be used to ensure canonical sampling, particularly 
for the case of the perfect crystal. We illustrate the method with a
calculation of the free-energy of formation of the sodium vacancy and
of the free-energy of the sodium crystal.
\end{abstract}

\pacs{PACS numbers: 02.70.Ns, 61.70.Ji, 65.50.+m}

\section{Introduction}
\label{sec_intro}
 The calculation of free energies in classical computer simulations, in which
 effective internuclear potentials are used to describe the interactions, 
 has become a standard practice \cite{frenkel86}. Since the entropy cannot be 
calculated as a statistical average at a single thermodynamic 
 state-point, thermodynamic
 integration methods must be used to relate the free energy of the system of 
interest to that of some reference system for which it is known or can
 be calculated analytically \cite{bennett75}. 

Following Frenkel and Ladd \cite{frenkel84},
 the standard method of calculating free energies in the solid state has 
 been to construct a reversible path between the system of interest and
 an Einstein crystal with the same structure. This is achieved by 
 decoupling the interparticle interactions whilst simultaneously switching
 on a harmonic potential which binds each nucleus to its lattice site. The
 decoupling is achieved by performing a series of simulations with potentials 
 intermediate between that of the system of interest and the harmonic limit. 
 The switching is controlled by a parameter $\lambda$ and the derivative 
 of the free energy with respect to $\lambda$ ($
 \frac{\partial F(\lambda)}{\partial \lambda}$) is obtained
 from a  canonical average in each simulation, by using either a Monte-Carlo 
(MC) or Molecular Dynamics (MD) method: finally, the free energy is calculated 
by integrating this derivative over $\lambda$ and adding the free 
 energy of the reference state \cite{squire69,broughton83}.

Concurrent with the development of such free energy techniques has been the
 rapid evolution of ``{\it ab-initio}" or ``first-principle" molecular
 dynamics methods (FPMD), following their introduction by Car and Parrinello
 \cite{carpar85}. In such simulations the interactions are derived from an 
 explicit representation of the electronic ground state and the introduction of
 effective potentials is avoided. The FPMD methodology thereby greatly expands 
 the range of physical systems which may be simulated reliably \cite{galli93}.
 In particular, problems in which a non-trivial rearrangement of the 
 electrons may accompany the different nuclear configurations visited in
 the course of the simulation (``itinerant electron" problems) become 
 accessible. 

 Application of the thermodynamic integration methodology in the context
 of FPMD presents a number of problems which we address in this paper. In
 order to implement the Frenkel-Ladd procedure we must not only replace 
 a direct internuclear potential by the harmonic one in a series of 
 $\lambda$ steps, but also simultaneously decouple the nuclei from
 the valence electrons in such a way as to allow 
$ \frac{\partial F(\lambda)}{\partial \lambda}$ to be calculated. We have 
examined several ways to achieve this and present one of them below. 
 The second problem is that, in order to calculate
$ \frac{\partial F(\lambda)}{\partial \lambda}$  accurately, relatively 
 long simulation runs are necessary in order to achieve acceptable statistical
 precision in the canonical average. For the most part, FPMD simulations
 are very  time consuming and this has doubtlessly inhibited their use in 
 free energy calculations hitherto. In the present work we partially 
 circumvent the problem by the use of a recently introduced,
  extremely efficient {\it ab-initio}
 scheme which avoids the introduction of orbitals to represent the kinetic
 energy functional \cite{our}. This technique is particularly appropriate for 
metallic 
 systems. Its use is not an essential ingredient of the free-energy development;
 the same scheme could be implemented (at substantially greater computational
 cost) within a general Kohn-Sham (KS) 
 density functional framework and, indeed, this has has already been
 done \cite{tesi}. 

 Besides these obvious problems for {\it ab-initio} free energy simulations,
 we have spent a considerable amount of time investigating a difficulty
 of more subtle origin. The tractability of FPMD methods is derived from
 the fact that they are {\em dynamics} schemes, the adiabatic electronic state
 is predicted from one step to the next and not re-computed {\it a priori}
 for each new nuclear configuration. However, molecular dynamics methods
 become very poor at sampling phase space when the underlying interactions are 
 harmonic, which is precisely the situation which pertains near the 
 Einstein crystal limit of the thermodynamic integration procedure. This 
 problem does not normally arise in the classical free energy simulations
 where MC methods may be used for the partially decoupled calculations.
 Hence we have examined various extended dynamics schemes for the nuclear 
 motion to examine the extent to which they render the motion ergodic, without
 wrecking the adiabatic dynamics which makes FPMD so efficient.

 In this paper we illustrate the application of the methodology by
 calculating the free energy of vacancy formation $F_f$ in sodium at various
 temperatures. 
This is an important quantity because it controls the vacancy
concentration, $ c = e^{-F_f/k_B T} $, which, in turn, is relevant
for a description of many dynamical processes, such as diffusion
processes in the solid phase. 
 The calculation of  defect properties of metals is a 
 non-trivial application of FPMD techniques. The atomic and electronic 
environment around an ion at a defect is very different from that in the bulk
 and the use of the same effective potential to describe interactions in
 both locations is fraught with conceptual and practical limitations. A
 more extensive discussion of the free energies of various point defects in 
sodium is presented elsewhere.


\par
We stress that the applicability of the method is not restricted to
the field of defects. For example, if the free energies of suitable
reference systems are available, e.g. by calculations using classical
empirical potentials, it is possible to calculate the phase diagram of
materials from first principles. Work is in progress along these
lines.

\section{Free-energy calculations}
\label{sec_free}

\subsection{Free-energy calculations via MD}

We begin by recapitulating the main points of the  thermodynamic integration 
methods. We write the energy of the system of $N$  ions in volume $\Omega$
as $E (\{ {\bf  R}_I \},\{ {\bf \dot R}_I \};\lambda) = 
 K_I(\{ {\bf \dot R}_I \}) + U (\{ {\bf  R}_I \}; \lambda)$, where $K_I$ is the
ionic kinetic energy, $U$ the potential energy and $\lambda$ is a
parameter such that $\lambda \in [0,1]$. 
 In the FPMD scheme the potential energy 
is obtained from an explicit representation of the electron density
 but this does not alter the fact that the potential energy is simply a 
 function of the nuclear positions $\{ {\bf  R}_I \}_{I=1,N}$ since the {\em adiabatic}
 electron density is itself determined by the instantaneous $\{ {\bf  R}_I \}_{I=1,N}$.
 Formally then, the calculation of the free energy in FPMD is just like in
 a classical simulation, in which the dependence of $U$ on $\{ {\bf  R}_I \}_{I=1,N}$ is given by an 
 explicit function.

 The thermodynamic integration method is based on two facts:
\begin{itemize}
  \item The derivative of the free energy with respect to $\lambda$
        has the form of a canonical average, since
        \begin{eqnarray}
             \frac{\partial F(\lambda)}{\partial \lambda}
            & = & - k_B T \frac{\partial}{\partial \lambda}
            \ln \int_{\Omega^N} {\rm d}^{3N}R ~ e^{-U (\lambda)/k_B T}
            \label{dfdl} \\
            & = & \frac{1}{{\cal Z}_q (\lambda)} \int_{\Omega^N} {\rm d}^{3N}R
             ~ \Phi (\lambda) \, e^{-U (\lambda)/k_B T} =
             \langle \Phi (\lambda) \rangle_\Omega
            \nonumber
        \end{eqnarray}
        where ${\cal Z}_q$ is the configurational integral entering
        the partition function,
        $\Phi = \frac{\partial U}{\partial \lambda}|_\Omega$, and 
	the angular brackets indicate canonical ensemble average. 
  \item The free energy is a function of state so that, irrespective of the
        functional form of $U (\lambda)$ at intermediate values of
        $\lambda$ (provided it is a smooth function), we have
        \begin{equation}
          F(\lambda=1) - F(\lambda=0) = \int^1_0 {\rm d}
          \lambda \left\langle \Phi (\lambda) \right\rangle_\Omega ,
          \label{f1f0}
        \end{equation}
\end{itemize}
where $F(\lambda)$ is the free-energy of a system with energy
$E(\lambda)$.  In other words, one may use the parameter $\lambda$ to
drive the system via a reversible thermodynamic transformation from
the state corresponding to the energy functional $E(\lambda=0)$ to
another corresponding to $E(\lambda=1)$ (or vice versa) and eq.
(\ref{f1f0}) will give the desired free energy difference between the
two states. It is not necessary that this transformation be physically
realizable, but it is desirable that the states involved be as similar
as possible. In this case $\Phi$ will be a slowly varying function of
$\lambda$ and the integration in eq.  (\ref{f1f0}), by necessity
numerical, will require as few $\lambda$ points as possible.
\par
In general, one has to compare free energies of systems with different
total volumes, hence the need for a slightly more general form than
eq.  (\ref{dfdl}). We have then to allow $\Omega$ to be a function of
$\lambda$. We stress that even in this case $F(\lambda)$ is a
Helmholtz free-energy (albeit of an unphysical system), since the
volume of the system is kept constant as the configurational integral
is performed and it is only after that $F(\lambda)$ has been
calculated that the variation of $\Omega(\lambda)$ takes place.
\par
To derive this more general equation we first introduce new scaled
variables in the integrals in eq. (\ref{dfdl}): $s_i = R_i / \Omega^{1/3}
\Rightarrow {\rm d}^{3N} s = {\rm d}^{3N} R / \Omega^N $. Hence the
integral is now performed over the unit hypercube ${\cal Q}^N =
[0,\ldots,0]\times [1,\ldots,1]$ and 
\begin{eqnarray}
  \frac{\partial F(\lambda)}{\partial \lambda}
  & = & - k_B T \frac{\partial}{\partial \lambda}
  \ln \left[ \Omega^N(\lambda) \int_{{\cal Q}^N} {\rm d}^{3N}s~
  e^{-U^\prime (\lambda)/k_B T} \right]
  \label{dfscal} \\
  & = & - \frac{N k_B T}{\Omega} \frac{\partial \Omega}{\partial \lambda} +
  \int_{{\cal Q}^N} {\rm d}^{3N} s~  e^{-U^\prime (\lambda)/k_B T}
  \frac{\partial U^\prime}{\partial\lambda} /
  \int_{{\cal Q}^N} {\rm d}^{3N} s~ e^{-U^\prime (\lambda)/k_B T} 
  \nonumber
\end{eqnarray}
where $U^\prime(\{ s\};\Omega,\lambda) =  U(\{\Omega^{1/3} s\};
\Omega,\lambda)$. Scaling back to the old variables: 
\begin{eqnarray}
  \frac{\partial F(\lambda)}{\partial \lambda}
  & = & - \frac{N k_B T}{\Omega} \frac{\partial \Omega}{\partial \lambda} +
  \frac{1}{\cal Z}_q\int_{\Omega^N} {\rm d}^{3N} R~
  e^{-U(\lambda)/k_B T}
  \frac{\partial U}{\partial\lambda} = - P_{ideal~gas}
  \frac{\partial \Omega}{\partial\lambda} +
  \left\langle \frac{\partial U}{\partial\lambda} \right\rangle
  \nonumber
\end{eqnarray}
In taking the derivative of $U$ its full $\lambda$ dependence has to
be taken into account, including the implicit dependence via $\Omega$.
Therefore,
\begin{equation}
  \partial F/ \partial\lambda =
  - P_{ideal~gas} \frac{\partial\Omega(\lambda)}{\partial \lambda} +
  \left.\left\langle \partial U/ \partial \Omega \right|_\lambda\right\rangle
  \partial \Omega /\partial \lambda +
  \left.\left\langle\partial U/\partial \lambda \right|_\Omega\right\rangle
  = - P(\lambda) \frac{\partial\Omega(\lambda)}{\partial \lambda}
     + \left\langle \Phi(\lambda) \right\rangle_\Omega
  \label{dfdlp}
\end{equation}
where $P(\lambda)$ is the thermodynamic pressure exerted by the
system at volume $\Omega(\lambda)$ (see Appendix A).
\par
Alternative, equivalent formulae can be obtained by dividing the
thermodynamic integration into two parts: first, one performs the
integration over $\lambda$ at fixed $\Omega=\Omega(\lambda=0)$, and
then a conventional integration of the pressure over the volume at
fixed $\lambda=1$ (or equivalently one may first perform the pressure
integration at $\lambda=0$ and then the $\lambda$-integration).
However, compared to the above method this requires an additional
thermodynamic integration and it is thus rather more expensive.

\subsection{The case of defects}
\label{subsec_defects}

We specialize now   to the case in which we have a crystal
represented by an $N$ atoms cell and we wish to decouple $L \leq N$
atoms from it. This is the case of the calculation of the free energy
of formation of an intrinsic defect: vacancy, divacancy, interstitial,
etc.  (generalization to more complicated cases such as alloys or
impurities is straightforward).
\par
The quantity of experimental interest is the concentration $c = n/N$
of defects, where $n$ is the number of defects in the sample and $N$
the number of crystal sites, $N>\!\!> n >\!\!> 1$. This can be
obtained by minimizing the Gibbs free energy of the system $G(n) =
G_0(n) - T S^{conf}$ with respect to $n$, where $G_0(n)$ is the free
energy of a crystal with a specified configuration of $n$ defects
(which includes electronic relaxation and vibrational contributions
for that configuration) and $T S^{conf} = k_B T \ln N^{conf}(n,N)$ is
the contribution stemming from the configurational disorder. (We have assumed 
that the external pressure is negligible).
\par
In the case of a truly pointlike defect like a vacancy with symmetric
relaxation,
$N^{conf}(n,N)$ is the number of ways in which $n$ sites can be chosen
out of $N$, which is $N!/[(N-n)!n!]$. In general however a defect
sitting on a specific lattice site may exist in several different but
equivalent states. For example, in a cubic crystal an ``interstitial
[100] dumbbell'' consists of a lattice site occupied by a pair of
atoms forming a ``rod'' pointing in the [100] direction; it is clear
that  the [010] and [001]
orientations are completely equivalent. Hence, per each lattice site
occupied by a dumbbell we may have one of these $m=3$ orientations at
random, and the number of configurations becomes $N^{conf}(n,N) = m^n
N!/[(N-n)!n!]$. Since usually $c=n/N <\!\!< 1$, the defects are indeed
not interacting, $G_0(n)$ is independent of the arrangement of the
defects, and $G_f=\partial G_0/ \partial n|_{n=0}$ is a well defined
quantity: the formation free energy of a single, isolated defect with
a specified orientation.
\par
One finally obtains, in the thermodynamic limit, $G_f - k_b T \ln
(N/n) - k_B T \ln m = \partial G/ \partial n = 0$, hence $c = m
e^{-G_f/k_B T}$. Equivalently, one may omit the preexponential factor
and include into $G_f$ an ``orientational entropy'' $S^{or} = k_B \ln
m$. This is not a negligible quantity, as often $m$ is of order 10,
giving $S^{or}\sim 2-3 \, k_B$, but as it is a purely geometrical
factor it can be easily calculated once and for all per each kind of
defect.
\par
In what follows we will deal with $G_f$ only, or, to be precise, with
its Helmholtz counterpart $F_f$. We will also neglect the electronic
contribution to the free energy of formation of defects, due to the
change in the electronic excitation energies induced by the defect.
This is unlikely to be of importance in simple metals, as it is
proportional to the change in the density of states at the Fermi
energy. In semiconductors and insulators, on the other hand, a defect
can introduce levels in the band gap, thus dramatically altering the
threshold of the electronic excitations. In this case, the electronic
contribution may be significant \cite{blochl93}.
\par
The formation free-energy $F^N_f$ as calculated by using an $N$ atoms
supercell is
\begin{equation}
  F^N_f (\Omega_C,T) =
   F^{N-L}_D (\Omega_D,T) - \frac{N-L}{N} F^N_C (\Omega_C,T)
  \label{free_nolambda}
\end{equation}
where $F^N_C(\Omega_C,T)$ is the free energy of a crystal with $N$
atoms, volume $\Omega_C$ and temperature $T$, while $F^{N-L}_D
(\Omega_D,T)$ is the free energy of a crystal with the same number of
 lattice sites but $L$ of them being vacant. We choose as the volume of the
supercell containing the defect $\Omega_D = \frac{N-L}{N}\Omega_C =
(N-L)\Omega_a$, $\Omega_a$ being the atomic volume in the perfect
crystal.  Note that according to this definition the atomic volume
stays constant: other choices are possible, but this should be the one
which gives the weakest $N$-dependence of $F_f^N$ \cite{devita91}.
  We start from the
perfect crystal and decouple $L$ atoms from it (say, the atoms with
indexes $I = 1, \ldots , L$); i.~e., we pick up $L$ atoms and put each
of them in reference states described by potential functions $u_I
({\bf r}_I)$, reversibly. We choose these reference states as states
in which the atoms have no interaction with any other. We choose
$\lambda$ as the coupling parameter of these atoms to the rest of the
crystal, including both the direct ion-ion and the ion-electron
interaction, so
\begin{itemize}
  \item $E^N(\lambda = 1) = E^N_C = $ energy functional of the perfect crystal,
  \item $E^N(\lambda = 0) = E^{N-L}_D + \sum_{I=1}^L \varepsilon_I = $ energy
    functional of the crystal with $N-L$ atoms plus the energy of the
    decoupled atoms.
\end{itemize}
Therefore the free energies of the bulk system and of the system with
decoupled atoms are respectively
\begin{eqnarray}
  F^N_C     & =  & F^N (\lambda = 1) \nonumber \\
  F^{N-L}_D & =  & F^N (\lambda = 0) - \sum_{I=1}^L \varphi_I
\end{eqnarray}
where $\varphi_I$ is the free-energy of an atom in a potential $u_I$,
at the appropriate temperature. This is a well defined quantity as
long as each potential $u_I$ is a single-body potential.  We impose a
linear $\lambda$ dependence of $\Omega$: $\Omega(\lambda) =
(1-\lambda) \Omega_D + \lambda \Omega_C$, which gives
$\frac{\partial\Omega(\lambda)}{\partial \lambda} = L \Omega_C / N = L
\Omega_a$. Then, using eq. (\ref{dfdlp}),
for the free energy $F^N_f$ we get
\begin{eqnarray}
  \label{effelambda}
   F^N_f & = & F^{N-L}_D  - \frac{N-L}{N} F^N_C =
   F^N (\lambda=0) - F^N (\lambda=1) +
   \frac{L}{N} F^N_C - \sum_{I=1}^L \varphi_I \nonumber \\
   & = & - \int^1_0  {\rm d}\lambda \, \left\langle \Phi (\lambda)
   \right\rangle_\Omega
   + L \Omega_a \int^1_0  {\rm d}\lambda \, P (\lambda) + L \Delta F.
\end{eqnarray}
\par
The two integrals in eq. (\ref{effelambda}) refer to the same
transformation and may be performed using the same set of $\lambda$
points, provided it gives the desired accuracy for both.  It is
apparent that the states described by the $u_I$ are arbitrary: their
free energy is implicitly added to $F^N_f$ when computing $F(0)$ and
subtracted in the term $\Delta F = \frac{1}{N} F^N_C - \frac{1}{L}
\sum_{I=1}^L \varphi_I $. A suitable choice for bulk systems is to
choose $u_I$ as identical harmonic potential wells, i.~e.
\begin{equation}
  \label{well}
  u_I = \varepsilon_0 + \frac{1}{2} M_I \omega^2 ({\bf r}_I - {\bf r}^0_I)^2
~~~~~~I=1,\ldots,L
\end{equation}
$M_I$ being the ionic masses, and ${\bf r}^0_I$, $\varepsilon_0$, and
$\omega$ arbitrary quantities. This is effective if ${\bf r}^0_I$ are
taken to be the position vector of the atomic sites to be vacated and
$\omega$ a typical crystal frequency: in this case the decoupled atoms
will behave roughly as bulk atoms and the integral over $\lambda$ will
run over fairly similar systems.
\par
The function $\Delta F$ is the free-energy per atom of the crystal
relative to an Einstein crystal having the same thermodynamic
parameters. There are several ways of calculating it: as $\frac{1}{N}
F^N_C$ is a bulk quantity, one might, for example, extract its value
from experimental data or from accurate quasi-harmonic \cite{jacucci84}
first-principles calculations of phonons.
 The free-energy of the
Einstein crystal, of course, is known analytically \cite{frenkel84}.
 The most
consistent way of proceeding, however, involves an additional
$\lambda$ integration: one uses the same method discussed above, but
with $\Omega=$ constant (the integration of the pressure then
disappears) and $L=N$. In other words, one uses the parameter
$\lambda$ to make all the atoms in a cell disappear and at the same
time reappear as independent harmonic oscillators. In this case the
knowledge of the free energy of the Einstein crystal is not even
necessary. Calculations of this kind were first been performed by
Frenkel and Ladd \cite{frenkel84}, using an MC method, and have become
 standard since then.


\subsection{Thermodynamic integration in FPMD}

To compute the free energy differences with the above method in a
FPMD  scheme not only $L$ ions have to be adiabatically
decoupled from the others and from the electrons, but also a
corresponding number of electrons must be removed.
  A simple way in which this can be
accomplished is as follows:
\begin{itemize}
        \item[i)]  Multiply by $\lambda$ the ionic charge $Z$ and the bare
		    pseudopotential $v_I^{ps}$ of each of the ions to be
                    decoupled.
        \item[ii)] Adjust the electronic degrees of freedom so as that
                    the total electronic charge is $Z(N-L(1 -\lambda))$.
                    The specific way in which this is accomplished depends
                    on the specific functional used to describe the
                    electrons. In the KS formalism, a good choice
                    is to multiply the occupation numbers of $ZL$
                    orbitals by $\lambda$ \cite{tesi}.
\end{itemize}
In other words, the interaction Hamiltonian between the ions being
decoupled and the rest of the crystal is made disappear with $\lambda$
and external harmonic springs are at the same time switched on, so as
to compensate, roughly, for the missing interaction. When the density is
used as basic variable, we have the following equations:
\begin{eqnarray}
 E &=& K_I [\{ {\bf \dot R}_I \}] + U [\rho,\{ {\bf R}_I \}; \lambda]
  + \sum_{I=1}^L \left( \frac{1}{2} M_I \dot {\bf R}^2_I + w({\bf r}_I)
  \right) 
  \nonumber \\
 \rho_{\bf G = 0} &=& Z (N - L(1 - \lambda)) / \Omega \nonumber \\
  Z_I &=& \lambda Z,
    \phantom{\lambda v^{ps} ({\bf r}),} I \leq L \\
  v^{ps}_I ({\bf r}) &=& \lambda v^{ps}_I ({\bf r}),
    \phantom{\lambda Z,} I \leq L \nonumber \\
  w ({\bf r}_I) &=& (1-\lambda) \, u_I ({\bf r}_I). \nonumber
\end{eqnarray}
 where $\rho_{\bf G}$ is a Fourier component of the electron density, $\rho$ at 
 wavevector ${\bf G}$.
It is easy to verify that this functional satisfies the conditions
posed in section \ref{subsec_defects}.
\par
Another possible choice is to keep all the electrons in the system,
and just switch off the pseudopotential interaction. This would
greatly simplify the equations for $\partial E/\partial \lambda$,
which would be a subset of ours. Since in this case the electrons no longer 
interact with the ions but are otherwise still present, in
order to preserve charge neutrality the ionic charges must be kept in
full. This unfortunately implies that the ``decoupled'' atom is, in
fact, still interacting with the others at $\lambda = 0$, via the
direct ion-ion interaction. This makes this choice unsuitable for
defect calculations. On the other hand it may be useful in the
calculation of the free energy of the crystal or of the liquid, the
reference state (the one at $\lambda=0$) being a one-component
classical plasma in this case.

\section{FPMD Methodology}

\subsection{The energy functional}

First-principles MD calculations are, in general,
 very expensive. In particular in
metals the need for Brillouin zone sampling is pressing and this involves a
huge computational effort. Worse still, due to the need to
orthogonalize the orbitals the computation time increases as
$\Omega^3$. These requirements limit practical calculations to small
cells ($\sim$ 50 atoms) and short times, leaving the doubt that size
effect may still be present, especially in the dynamics.

 To overcome
these problems we use a scheme which has been developed recently
 \cite{our}. We
write the total energy as follows:
\begin{eqnarray}
 E [&\rho&, \{ {\bf R}_I \}, \{ {\bf \dot R}_I \}]
  = K_I + U = \label{our_ke}\\
 &K_I& + T_{TF} + T_{vW} + \int\!\int d{\bf r} d{\bf r}^\prime \,
 \Delta\rho^\alpha ({\bf r}) K_\alpha ({\bf r}-{\bf r}^\prime) 
 \Delta\rho^\alpha ({\bf r}^\prime)  + E_H + E_{XC} + E_{ps} + E_{II} 
 \nonumber
\end{eqnarray}
where $E_H$ is the coulombic interaction term between electrons,
$E_{XC}$ the exchange-correlation energy (for which we use the Local
Density Approximation), $E_{ps}$ the ion-electron
pseudopotential energy (we describe the electron-ion interaction by
the Topp-Hopfield pseudopotential
\cite{topp73}), and $E_{II}$ is the potential energy of the  direct ion-ion
interaction. Here $\Delta \rho = \rho -\rho_0$, with $\rho_0 =
N_{el}/\Omega $  is the average electron density. These terms are well known and
are routinely used in the framework of KS calculations; they can be found
written explicitly in Ref.\cite{remler90}. The kinetic term, however, is
replaced by one suggested by Perrot~\cite{perrot94} (a closely related 
 functional was proposed by Wang and Teter~\cite{wang92}) and includes the
well-known Thomas-Fermi and von Weizs\"acker terms \cite{dreizler90}, plus a
quadratic term in  $\Delta \rho $ which is chosen such
as to give the correct linear response. This functional does not require the
explicit introduction of wavefunctions.
\par
This scheme has numerous computational advantages in FPMD with respect to the
conventional KS scheme. First, the electronic part of the calculation
scales linearly with the system volume. Second, it can be easily and
effectively ``preconditioned''\cite{our}, which means that it is
possible to use comparatively large timesteps while at the same time
maintaining a very good adiabatic separation between the ionic and
pseudo-electronic subsystems; this feature is independent of the size
of the system, which means that the scheme has linear scaling in
practice as well as in theory. By avoiding the use of
wavefunctions any problems connected with poor sampling of the Fermi
surface (which usually forces one to use k-point sampling) are
avoided. Finally,  possible instabilities connected with level
crossings \cite{pastore91}, particularly acute in metals \cite{blochl92},
do not arise.  All in all, the gain in speed is about two orders of magnitude
for cells of about 100 atoms, and arguably greater for larger system
sizes.
\par
Of course, there is a price to pay for this, and it is some loss of
generality. However, the Perrot functional has been shown to give rather
accurate results for a range of properties of Na \cite{our,our1}, and
we believe it is perfectly adequate for the system which is studied
here.
\par
We use the above scheme in a Car-Parrinello-type implementation (for a
full description see Ref.~\cite{our}) in which the density is used as
basic variable. The resulting expression for $ \frac{\partial
E}{\partial \lambda} $ in this scheme is reported in Appendix B.

 We use  Nos\`e-Hoover thermostats \cite{noseh} on the
electronic and ionic subsystems in order to control their temperatures.
For the ions we have 
\begin{eqnarray}
 M_I \ddot {\bf R}_I &=& - \partial E/\partial {\bf R}_I - M_I \dot {\bf R}_I
 P_\xi/Q_i \nonumber \\
 \dot P_\xi &=& \sum_I M_I \dot{\bf R}_I^2 - 3 (N-1) k_B T
 \label{ions}
\end{eqnarray}
The Nos\`e - Hoover degree of freedom is here $P_\xi$, its ``mass'' $Q_i$
and $T$ the target temperature (i.e. the physical temperature of the
 system of interest).
\par
After expanding the density in plane waves, the Car-Parrinello
equations of motion read
\begin{eqnarray}
 \mu_{\bf G}\ddot {\rho}_{\bf G} &=&-{\delta E \over \delta \rho_{\bf G}}
  - \mu_{\bf G} \dot\rho_{\bf G} p_\xi/Q_e
  ~~~~~~~~~|{\bf G}|\le G_c, \nonumber \\
  \dot p_\xi &=& \sum_{\bf G} \mu_{\bf G} \dot\rho_{\bf G}^2 - N_g k_B T_e
\label{fake}
\end{eqnarray}
with the constraint that the total number of electrons be conserved:
$\rho_{{\bf G} = {\bf 0}} = \rho_0$. The parameter $\mu_{\bf G}$ is a
``fake mass'' or inertia parameter and $G_c$ the plane-wave cutoff.
The function $\mu_{\bf G}$ was chosen, following Ref.~\cite{our}, as
$-\mu_0/F(G)$, with $\mu_0$ a constant and $F(G)$ the response
function of the free electron gas at the same average density of the
system being simulated, within the random phase approximation
\cite{dreizler90}. The Nos\`e - Hoover degree of freedom is here
$p_\xi$, its ``mass'' $Q_e$, and $N_g$ is the number of independent
electronic degrees of freedom.

 The target temperature for the
electrons, $T_e$, is set to be a small fraction of the target ionic temperature,
 as we have found that this provides the optimum way of keeping the electron
 density close to its adiabatic value \cite{fois93}. In practice 
 we  set  $T_e = T \frac{\mu_0}{M_I} f$, where the factor $f$ can be 
 shown to be proportional to 
 $\mbox{thermal average of } \int_\Omega {\rm d}^{3N} {\bf r} \left(
\nabla \rho \right)^2$.
The proportionality constant depends only on the function $\mu_{\bf
G}$ and, more weakly, on the thermodynamic parameters. Therefore $f$
is dependent on the kind of  simulated system (pseudopotential and
chemical state) but is virtually independent of the number of
electronic or ionic degrees of freedom.  This means that one can
estimate it by taking the ratio $\frac{T_e/\mu_0}{T/M_I}$ just after
the ``backstep'' procedure \cite{remler90} of a test run.
  

\subsection{Choice of the parameters}
\label{sc_par}

We used an energy cutoff ($E_c = \frac{1}{2} G_c^2$)
in the planewave expansion of the density of about 9 Ryd, which gives a total
energy for the perfect crystal at $T=0$ of -0.225155 Hartree/atom, to
be compared with a converged value of -0.225159 Hartree/atom.
Dynamical test runs with $E_c = 20$ Ryd showed no meaningful
differences in the ionic trajectories. The time per step on a
54 atom cell with 9 Ryd is about 1.5 sec on a Convex 200, but in our
code there is still scope for further optimization.
\par
We found that choosing $\mu_0 = 5\ 10^4$ au we were able to run
dynamics with a time step as long as 100 atu (atomic time units; 1 atu
= $2.4188\ 10^{-5}$ ps), with energy fluctuations below $10^{-4}$
Hartree/atom. No signs of drift in the total energy were found.
Moreover, we also found that the amount of energy transferred to the
fake electronic degrees of freedom was always less than 1\% of the
kinetic energy of the ions. We found that for solid sodium at 293 K a 
 suitable value for  $f$ is 0.018 a.u.. This is more or less independent of 
$\lambda$ in the
calculations of the vacancy, as the energy functional changes only on
the order $O(1/N)$ in that case. When calculating $\Delta F$, on the
other hand, the system changes substantially and it turns out that $f$
must be multiplied by roughly $\lambda$ for the thermostat on the
electron to work properly. The parameter $Q_e$ was set to 1000 a.u.,
following the usual recipe \cite{noseh}. In all the runs the system was
started with ionic velocities picked at random from a Maxwell-Boltzmann
distribution at the appropriate $T$, while the electronic velocities
were set using the backstep procedure of Ref. \cite{remler90}.
\par
All the runs were performed at constant atomic volume, which at every
temperature was chosen to be the experimental value. The volume at
$T=100~K$ was assumed to be the same that at $T=0$. We set
$\varepsilon_0$ equal to the energy per particle for the corresponding
bulk system at $T=0$, and $\omega$ equal to the empirical Debye
frequency ($4.75\:10^{-4}$ a.u.).
\par
We set ourselves an accuracy limit of 0.01 eV. Any difference
smaller than that was not regarded as significant.

\subsubsection{Convergence of the integrals over  $\lambda$}

The integrals over $\lambda$ obviously must be discretized. Every
$\lambda$-point used requires a statistically converged,
finite-temperature simulation, hence it is vital to choose an
integration scheme which minimizes their number. If the systems at
$\lambda<1$ are suitably chosen, one may expect that $\Phi(\lambda)$
and $P(\lambda)$ can be closely approximated by low-order polynomials.
If this is the case, the Gauss-Legendre formula is optimal, as it
allows the exact evaluation of any polynomial $\gamma(\lambda)$ of
order $2j-1$ using only $j$ points: $\int_0^1 \gamma(\lambda) {\rm d}
\lambda = I_j = \sum_{i=1}^j W_i^j \lambda_i^j \gamma(\lambda_i^j)$,
where $\{ \lambda_i^j \}$ are zeros of the Legendre polynomials of
order $j$ and $W_i^j$ appropriate weights\cite{numrecip}. If $\gamma$ is
an arbitrary smooth function then, $ \int_0^1 \gamma {\rm d}\lambda =
\lim_{j\to\infty} I_j$.
\par
In order to choose the smallest value of $j$ we performed a series of
tests at $T=0$ (but allowing full ionic relaxation). In this case, the
thermal averaging disappears, the free-energies become total energies
and thus $F_f (\Omega, 0)$ can be obtained directly and exactly from
eq.~(\ref{free_nolambda}). We can then calculate $I_j$ for increasing
values of $j$ and check which value is sufficient to achieve the
preset convergence.  We found that for the vacancy formation energy
$j=2$ suffices (the integral of $P$ being actually already converged for
$j=1$). The points to be used are then $\lambda_1^2 = 0.2113248654$
and $\lambda_2^2 = 1 - \lambda_1^2 = 0.7886751346$, with $W_1^2 =
W_2^2 = 1/2$. On the other hand, the same convergence level on $\Delta
F$ was only obtained by using $j=4$. This is similar to the findings
of other authors on different systems
\cite{frenkel84,broughton83}. The set of points in this case is
$\lambda_1^4 = 0.0694318442$, $\lambda_2^4 = 0.3300094782$, $\lambda_3^4 =
1 - \lambda_2^4 = 0.6699905218$, $\lambda_4^4 = 1 - \lambda_1^4 =
0.9305681558$, with weights $W_1^4 = W_4^4 = 0.1739274226$ and $W_2^4 =
W_3^4 = 0.3260725774$.

\section{ The problem of thermalization}

The calculation of the free energy $\Delta F$ of the real crystal
relative to the Einstein crystal may run into a fundamental problem,
which to our knowledge has yet to be explicitly addressed in this
context. As $\lambda \to 0$, the coupling of the atoms vanishes,
leaving one with a system of uncoupled harmonic oscillators. It is
well known that a single harmonic oscillator has a poor ergodic
behaviour. It is true that at the Gauss-Legendre quadrature points 
 ($\lambda \ne 0$)these oscillators are coupled, but the
coupling is very small and since the Fermi-Pasta-Ulam calculation
\cite{fermi55} it is known that in these circumstances an ergodic
behaviour is by no means guaranteed. It is easy to realize that in our
scheme the coupling is essentially proportional to $\lambda^2$. 
 Both the number of electrons and the pseudopotentials are
proportional to $\lambda$, hence the electron-mediated interaction is
proportional to $\lambda^2$. The same is true of the direct ion-ion
interaction. The problem is compounded by the fact that the value of
$\lambda$ required may be rather small: for example, we wish to use
$\lambda=\lambda_1^4$, and in this case the coupling is about $0.07^2
\approx 0.5\%$ that of the real crystal, which means truly weak coupling.
\par
This problem does not arise in an MC calculation, where the canonical
sampling is imposed, or in MD calculations with certain ``random''
thermostats like Andersen's \cite{andersen80}. In the framework of
FPMD however the last thing one wishes is some abrupt external
disturbance on the system, since in that case the adiabatic behaviour
of the electronic degrees of freedom may be easily disrupted.
Unfortunately, ``continuous'' thermostats like the Nos\`e-Hoover,
which has many useful characteristics, are not sufficient to guarantee
canonical sampling. Indeed, if the system is not ergodic the correct
sampling is not obtained\cite{noseh}.
\par
The problem has been tackled by Martyna {\sl et al.}\cite{martyna92} who
diagnosed the root of the problem to be the lack of control over the
fluctuations of the thermostat degree of freedom. To cure it, they
simply imposed an additional thermostat on the first; the
procedure can be repeated by imposing another thermostat on the
second, and so on as many times as needed. In fact, they found that a
chain of two thermostats is sufficient to give a good representation
of the canonical distribution even in the case of the harmonic
oscillator. The equations of motion in this case read as follows:
\begin{eqnarray}
 M_I \ddot {\bf R}_I &=& - \partial E/\partial {\bf R}_I - M_I \dot {\bf R}_I
 P_\xi/Q_i \nonumber \\
 \dot P_\xi &=& \sum_I M_I \dot{\bf R}_I^2 - 3 (N-1) k_B T -
  P_\xi P_\eta/Q_t
  \nonumber \\
 \dot P_\eta &=& P_\eta^2 / Q_i - k_B T
 \label{NHC}
\end{eqnarray}
where $P_\eta$ is the second thermostat of the chain and $Q_t$ its
``mass''.  In Figure \ref{fighar} we show the coordinate distribution
function of an harmonic oscillator without thermostatting, with a
Nos\`e-Hoover thermostat, and with a chain of two thermostats. The
target temperature was $200~K$, the oscillator mass was chosen to be
that of sodium and the frequency to be the Debye frequency of sodium.
The parameters $Q_i$ and $Q_t$ were chosen according to Martyna's
prescription \cite{martyna92}, which in this case gives $Q_i = Q_t =$
3000 a.u.. It can be seen that the coordinates distribution when a
chain of thermostats is used differs from the exact canonical
distribution by about the amount of the thermal noise, showing that if
the statistical error is under control the residual deviation from the
canonical distribution is under control as well. The other two
distributions are  qualitatively wrong.
\par
A less obvious problem arises  when one tries to simulate an
{\sl ensemble} of harmonic oscillators using this method. In Figure
\ref{fignh} we show the individual atomic energies of a system of 16
identical harmonic oscillators (masses, frequencies and temperature as
above), thermostatted by using a chain of two thermostats. It is clear
that the system is not going towards an equipartition of the energy.
Rather, some oscillators stay hotter than the mean. Even worse, there
is no sign of any evolution after 1 ps or so. A plot of the
distribution of coordinates reveals that again the  canonical
distribution has not been reached.
\par
It seems that the atom-atom coupling induced by the thermostat is wholly
insufficient to make the system canonical. A close inspection of Fig.
\ref{fignh} indicates  why this is so. Let us suppose that, as in the case of
the figure, the instantaneous temperature at the start is too high.
The thermostat will try to cool down the system by damping the atomic
velocities, but the value of the damping is determined by the {\em mean}
 departure of the oscillator kinetic energies from their target values.
 This means that  atoms which are hotter than the target $T$
are underdamped, while the cooler ones are overdamped.
   Later on, when the average kinetic energy of the
atoms is correct, the damping parameter relaxes to a small value. 
 If some atoms are outside the physical range of kinetic energy fluctuations,
 they will stay out of it forever, since the only thing that the thermostat
is actually doing after the first 1 or 2 ps is to synchronize the
oscillations of the individual atomic energies.
\par
It is clear that adding other thermostats to the chain is not going to
help, the coupling between the oscillators via the thermostat is 
 not effective in transferring energy between them.
 We  use another method, which  forces the canonical distribution; 
 this method too has been devised by Martyna {\sl et al.}
\cite{martyna92} who called it ``massive Nos\`e-Hoover chain'' (MNHC).
It simply consists of attaching an independent Nos\`e-Hoover chain to
each degree of freedom of the system. If the chain method works for a
single, isolated degree of freedom it must work also if the said
degree of freedom belongs to a larger system.  In Figure
\ref{figmnhc} we show the effect of the introduction of such a system
of chains over the same system of Fig. \ref{fignh}: after about 3 or 4
ps the system looks completely chaotic, as desired. The resulting
coordinate distribution function after a longer run is shown in Fig.
\ref{figdist}, alongside with the exact distribution. It is clear that the
system has been satisfactorily thermalized.

\subsection{Dynamics of the center of mass}

If the forces on the simulated system are internal only, the system
may translate without changing its energy. Therefore the number of
degree of freedom is decreased by 3, as the center of mass is
effectively ``frozen''.  When using external forces on the system,
like the harmonic springs in our case or a thermostat, this is not
true anymore and forces on the center of mass may appear.
 These may cause difficulties. It is clear that if the centre of mass
 starts to wander there will be uncontrolled contributions
 to the harmonic energy in the Frenkel-Ladd procedure. Furthermore, the
number of degrees of freedom may change and cause problems 
  in the definition of instantaneous temperature; in 
$ \frac{1}{2} N_d k_B T = K_I $ it is not clear if the number of degrees of
freedom $N_d$ should be $3N$ or $3(N-1)$.
\par
Sometimes the problem can be avoided quite simply. In a system described by
the equations of motion (\ref{ions}), the forces $ - \partial
E/\partial {\bf R}_I$ are the sum of internal forces, depending only
on the relative coordinates of the atoms, and of identical external
springs $u_I$ [eq.  (\ref{well})]. By summing over $I$ and remembering
that the sum of the internal forces is zero one obtains an equation of
motion for the center of mass:
\begin{equation}
 \ddot {\bf R}_{CM} = - \omega^2 \left( {\bf R}_{CM} - {\bf R}_{CM}^0 \right)
 - \dot {\bf R}_{CM} P_\xi/Q_i
\label{cofm}
\end{equation}
where ${\bf R}_{CM}^0 = \sum_I M_I {\bf r}^0_I / \sum_I M_I $. This
shows that, in general, the center of mass is acted upon by a spring
with the same natural frequency of the springs applied to the single
atoms, and will be thermalized by the same thermostat acting on the
ions. However, if one starts with ${\bf R}_{CM} = {\bf R}_{CM}^0$ and
$\dot {\bf R}_{CM} = 0$ the acceleration will be 0 and the center of
mass will stand still forever (apart from numerical and integration
errors). Notice that this is still true if three thermostats are used,
one per cartesian direction.
\par
This kind of quick fix is not possible when different thermostats are
used on different atoms or subsets of atoms, for example when a MNHC
is used. Each thermostat will exchange momentum independently with the
degree(s) of freedom to which is attached and there is nothing that
guarantees that these exchanges will add up to zero. 
 The thermalisation problems which led us to introduce the MNHC chain
 are only pronounced close to $\lambda$ equal to 1 for
 the perfect crystal calculation, i.e. close to the Einstein 
 crystal limit. We therefore carry out the calculations for other $\lambda$
 points with a single Nos\'e-Hoover chain and impose a stationary center
 of mass by the ``quick fix" described above. The
 temperature in these calculations is obtained from 
 $ T = K_I / \frac{3}{2} (N-1) k_B $. At $\lambda$ close to 1 we resort 
 to the MHNC with a moving center of mass and use  $ T = K_I / \frac{3}{2} 
 N k_B $ and subtract the potential energy of the center of mass 
 (equation (14)) from the total
 energy. Note that in this limit the center of mass is bound by a harmonic 
 spring and consequently does not make large excursions at thermal energies.
 We have checked that, on systems of 50 atoms or more, both calculations
 give the same results at intermediate values of $\lambda$ where thermalisation
 takes place as a result of the anharmonic interactions between the ions.
 
\section{Results}

\subsection{The perfect crystal}

A 54 atom cell was chosen for the perfect crystal calculations.  A
convergence study with cells up to 128 atoms showed that 54 atoms
gives the required convergence in the free energies.  All the values
$Q_i$ and $Q_t$ were set to 3000 a.u..
\par
The run were found to be stable and adiabatic, as described in sec.
\ref{sc_par}. The total length of the runs was
$\tau = $ 6 ps. The instantaneous values of $\Phi ( \lambda ) $ were
recorded from the start of the run; afterwards, we calculated their
running averages by starting from different instants $t_0$. This
amounts to allowing a time $t_0$ for thermalization and a time
$\tau-t_0$ for the statistical averages $\left< \phi(\lambda)
\right>_{t_0}$. We have consistently found that that when $t_0$ ranged
from 2 to 4 ps the dispersion in the values of $\left< \phi(\lambda)
\right>_{t_0}$ was steadily small [$O(10^{-2}$ eV)], but whenever
$t_0$ was 1 ps then $\left< \phi(\lambda)\right>_{t_0}$ lay
significantly outside this dispersion band. The same happened whenever
$t_0 = $ 5 ps. We conclude that a thermalization period of 2 ps is
necessary and probably sufficient for our purposes, and that a similar
time span is required in order to have converged statistical averages.
The running averages of $\left< \phi(\lambda) \right>_{t_0}$ in a
typical case for $t_0 = $ 2, 4, 6 ps are shown in Fig.
\ref{figdeltaf}.
\par
The results for $\Delta F$ are shown in Table \ref{tabfree}.

\subsection{Vacancy}
\label{res_vac}

For the vacancy calculation we chose a 54 atom cell again. Not only
 does this cell gives a good convergence with size for the vacancy formation
energy, but it also includes the $\frac{2}{3} \frac{2}{3}
\frac{2}{3}$ phonon which is believed to be crucial for many vacancy
properties in bcc crystals\cite{kohler88}.
\par
The thermalization problems which beset the calculation of $\Delta F$ 
 were less pronounced in the vacancy calculations. The critical
issue is  whether the decoupled atom thermalizes at all
the $\lambda$ values required. If thermalization happens for that
atom, one can safely assume that the other atoms thermalize as
well, since they experience the full anharmonic interatomic interaction. 
  The coupling of the
decoupled atom to the remaining (anharmonically interacting) atoms 
 is now only proportional to $\lambda$, whereas the
 couplings disappear as $\lambda^2 $ in the perfect crystal calculations.

We performed a test run over our 54 atom cell, without any thermostat,
at $\lambda = \lambda_1^2$ (which is arguably the most difficult case)
and $T=200~K$. The run was 12 ps long and the corresponding coordinate
distribution function for the decoupled atom in 
  the first and the second halves of the run
are shown in Fig. \ref{figlam}. The continuous line represents the
canonical distribution for a harmonic oscillator with a frequency
$\omega$, i.e. the Einstein oscillator to which the decoupled atom
reduces as $\lambda \to 0$. We expect this to be a good approximation
to the true, unknown canonical distribution function of the
coordinates of the partially decoupled atom
  in sodium, the main difference being that the
true distribution may have a tail in the large-$\Delta R$ region, due
to the presence of anharmonic effects and low-frequency phonons.
\par
The coordinate distribution functions calculated from the run  do indeed look
 canonical, the difference between them and the theory being
mostly within the statistical noise or  explicable with the
arguments above. We conclude that a 54 atom cell is, at least within
our model, close enough to the thermodynamic limit to give canonical
distributions, even if one atom is coupled only to 20\%.
\par
Running averages of $\Phi (\lambda)$ analogous to those of Fig.
\ref{figdeltaf} are shown in  Fig. \ref{figphi} for a typical case.
Again, we found that $t_0 =$ 2 ps and $\tau - t_0 = $ 2 ps were enough
for a statistically converged result. The calculated values of
$F_f(T)$ are shown in Fig. \ref{figfree}. The uncertainty on $F_f$ is
not small enough to allow us to calculate the formation entropy, which
is its derivative respect to $T$, however the downward trend of $F_f$
with increasing $T$ is clear: this is indeed unsurprising, as as the
anharmonic effects (which usually increases the entropy) should be
larger in the defective crystal as a consequence of the
vacancy-induced phonon softening. Also discernible is the downward
curvature of $F_f (T)$, indicating that, as expected, the anharmonic
effects increase with the temperature. This is particularly apparent
in the point at $T=370~K =$ experimental melting point of Na, where
also the statistical fluctuations on $\Phi$ increase sharply. We
interpret this as a precursory sign of the melting.

\section{Conclusions}

In this work we have shown that it is possible to apply the
thermodynamic integration technique to FPMD calculations. This is
straightforward if the integration parameter is a physical parameter;
in this cases the only difference with respect to the classical MD
calculations is the way the forces on the ions are calculated. We have
concentrated then on the non-trivial case in which the coupling among
atoms, among atoms and electrons, and among electrons is changed
during the integration, and the coupling parameter becomes the
integration parameter. In this case the the fact that the electrons
are explicitly present in the model introduces novel complications and
problems in the method. The details of course depend on the specific
model used to describe the electrons: we have chosen a method we
recently developed that allows one to perform easily the large-scale
computations required for a realistic calculation. (A complete
calculation of $F_f (T)$ for a single temperature required only 6 -- 7
hours of CPU time on a Convex 200 machine, including the calculation
of $\Delta F$ and all the necessary minimizations). We nevertheless
feel that many of our findings have a general validity.
\par
We have shown that, as long as the number of atoms to be decoupled is
a small fraction of the total number of atoms present in the cell, no
problem of thermalization arises. The system becomes thermally
equilibrated after about 2 ps and the averages yield statistically
converged results if performed over 2 ps. No instabilities,
non-adiabatic behaviour, or other numerical problems were encountered.
On the other hand, when most or all of the atoms become decoupled from
each other, one must be careful in order to get canonical averages.
The system may easily stay out of thermal equilibrium for a very long
time unless corrective action is taken. The method of Martyna {\sl et al.}
 \cite{martyna92} has been found to be effective in this respect. 

\par
Finally, we wish to stress that the method is not limited to the case
illustrated here (defects); other problems can be tackled in a similar
way. For example, one may compute the mixing free-energy of a solid or liquid
alloy by transmutating some of the atoms of one kind in atoms of the
other kind, in this case $\lambda$ is the parameter driving this
``alchemistic'' transformation. At present, we are using the
$\lambda$-integration method to calculate the melting point of caesium
from first-principles.

\acknowledgements

Thanks to C. Ursenbach for valuable technical help with the program.
The work was supported by SERC and the MOD (UK) under grant GR/H38256.
It has also benefited from further support from SERC under grant no.
GR/H10276.

\appendix

\section{Pressure in the FPMD scheme}

The pressure can be calculated by just taking the trace of the stress
tensor given in \cite{nielsen85}. Here we derive the same expression
by using a simple direct method. By definition,
\begin{eqnarray}
  P &=& \left. -\partial F/\partial \Omega \right|_T =
   \left. k_B T\partial \ln {\cal Z} /\partial \Omega\right|_T =
   N k_B T/ \Omega - \frac{1}{{\cal Z}_q} \int_{\Omega^N} {\rm d}^{3N} q \,
   e^{-U/k_B T} \frac{\partial U}{\partial \Omega} \nonumber \\
   &=& P_{ideal~gas} - \left \langle \partial U/\partial \Omega
   \right \rangle_T
  \label{pressure}
\end{eqnarray}
In the Born-Oppenheimer approximation, $U$ includes the direct
interaction energy between the ions $E_{II}$ plus the energy of the
electrons in the ground state: $E_k$, $E_H$, $E_{XC}$, and $E_{ps}$.
Using the Hellmann-Feynman theorem \cite{nielsen85}, we can explicitly
calculate the thermal average in (\ref{pressure}) as the thermal
average of the derivative respect to $\Omega$ of the {\sl functional}
$U$, if during the run the electronic system stays always very close
to the instantaneous ground state. We first note that the derivative
of $U$ involves a virtual expansion of the cell, which means that the
ionic coordinates and the electronic density have to be simply
rescaled with $\Omega^{1/3}$.  Therefore, any energy term $E_x$ which
scale with the volume as $V^\gamma$ gives a contribution $-\gamma E_x
/ \Omega$ to the pressure:
\begin{eqnarray}
-\partial E_k/\partial \Omega |_T &=& \frac{2}{3}E_k/\Omega \\
-\partial E_H/\partial \Omega |_T &=& \frac{1}{3}E_H/\Omega \\
 -\partial E_{II}/\partial \Omega |_T &=& \frac{1}{3}E_{II}/\Omega
\end{eqnarray}
We have assumed here that the kinetic energy scales correctly with
$\Omega$. Although this is not rigorously true for the quadratic term
in the Perrot functional, the error introduced in the case of Na is in
practice negligible (less than 1\%). If this turns out not to be the
case in other systems, or if greater accuracy is needed, it is
possible to calculate it exactly, as it is done in Appendix B for
$\partial E_k/\partial \lambda$.
\par
Noting that $\partial \left[ \rho({\bf r}){\rm d}^3 {\bf r}\right]/\partial
\Omega = \partial \left[ \rho({\bf r}) \Omega \right]/\partial \Omega = 0 $
it is then easy to show that
\begin{equation}
 -\partial E_{XC}/\partial \Omega |_T = \frac{1}{\Omega} \int_\Omega {\rm d}^3
{\bf r} \,  \rho({\bf r}) \left( \mu_{XC} - \varepsilon_{XC} \right)
\end{equation}
where $\mu_{XC} = \partial \left[ \rho({\bf r}) \varepsilon_{XC} (\rho)
\right] /\partial \rho$. This form is specific to the Local Density 
Approximation, which introduces
a characteristic length of its own (the radius of the
exchange-correlation hole), therefore concealing the scaling of the
underlying Coulombic interaction. If one uses, say, the Hartree-Fock
approximation, then $E_{XC}$ scales exactly as $\Omega^{-1/3}$ and the
contribution to the pressure is simply $\frac{1}{3}E_{XC}/\Omega$.
\par
Finally, the pseudopotential term is best calculated in reciprocal
space.  We assume that only one species is present. Using the
conventions of \cite{remler90}, the dependence of the pseudopotential
$v_I^{ps}({\bf G})$ on $\Omega$ is twofold: first, there is an
explicit dependence of $v_I^{ps}({\bf G})$ on $\Omega$, since it is
normalized to the atomic volume, and second, there is an implicit
dependence via the $G$ vectors. Therefore ${\rm d} v_I^{ps} ({\bf G})
/ {\rm d} \Omega = - v_I^{ps}({\bf G})/\Omega - G/3\Omega ~~
\partial v_I^{ps}({\bf G}) / \partial {\bf G}$, and
\begin{eqnarray}
 -\partial E_{ps}/\partial \Omega |_T
 &=& -\Omega \sum_{\bf G} \rho^*({\bf G}) S ({\bf G})
 {\rm d} v_I^{ps} ({\bf G})/ {\rm d} \Omega \nonumber \\
 &=& E_{ps}/\Omega +
 \frac{1}{3} \sum_{\bf G} G \rho^*({\bf G}) S ({\bf G})
 \partial v_I^{ps}({\bf G})/\partial {\bf G}
 \label{pseupres}
\end{eqnarray}
where we used the invariance of $\rho({\bf r}) \Omega $ and of the
structure factor $S ({\bf G}) = \sum_I e^{-i {\bf G}{\bf R}_I}$
respect to virtual changes in the volume. Notice that the contribution
stemming from the $G=0$, noncoulombic part of the pseudopotential,
which in Ref. \cite{nielsen85} is treated separately, is here included
in the first term in the right-hand side of (\ref{pseupres}).



\section{Formulae for $\partial E/\partial\lambda$}

As described in Appendix A for the pressure, the Hellman-Feynman
theorem can be used to calculate the derivatives of the energy with
respect to $\lambda$. For the case described in the text we have
\begin{eqnarray}
\partial \rho ({\bf r})/\partial \lambda &=& L Z / \Omega \Rightarrow 
\partial \Delta \rho ({\bf r})/\partial \lambda = 0 \\
\partial v_I^{ps}/\partial \lambda &=& v_I^{ps} /\lambda \mbox { for }
I\leq L, \mbox { 0 otherwise}\\ 
\partial Z_I/\partial \lambda &=& Z_I /\lambda \mbox { for } I\leq L,
\mbox { 0 otherwise} \\
\partial k_F/\partial \lambda &=& L Z k_F / 3 N_{el} \\
\partial \theta /\partial \lambda &=& - L Z \theta/3 N_{el}
\end{eqnarray}
with $\theta = G /2 k_F$. It is clear that $v_I^{ps} / \lambda$ and
$Z_I / \lambda$ are simply the values of the bare pseudopotential and
of the ionic charges prior to the scaling. One obtains:
\begin{eqnarray}
  \partial E_k /\partial \lambda &=& L Z \frac{5}{3}\frac{C_k}{\Omega}
     \int_\Omega
    {\rm d}^3{\bf r} \, \rho^\frac{2}{3} - \frac{LZ}{8\Omega}
     \int_\Omega {\rm d}^3
    {\bf r} \, \left( \nabla\rho/\rho \right)^2 \nonumber \\
   &+& \frac{\Omega^2}{N_{el}}
    \sum_{\bf G} |\rho^\alpha({\bf G})|^2 K_\alpha(G)
    \kappa_\alpha(\theta) \\
  \partial E_H /\partial \lambda &=& 0 \\
  \partial E_{XC} /\partial \lambda &=& \frac{LZ}{\Omega}
    \int_\Omega {\rm d}^3{\bf r} \, \mu_{XC} \\
  \partial E_{ps} /\partial \lambda &=& L Z \sum_{I=1}^N v_I^{ps} (G=0) +
    \Omega \sum_{\bf G} \rho^* \sum_{I=1}^L e^{i{\bf GR}_I}
    v_I^{ps}/\lambda\\
  \partial E_{II} /\partial \lambda &=& - \pi L Z_1 \sum_{I=1}^N
    Z_I/\lambda \Omega \eta^2 - 2 L \eta Z_1^2 /\sqrt\pi \lambda
    \nonumber \\
  &+& \frac{1}{\lambda} \sum_{I=1}^L
    \left[ 2\sum_{J=I+1}^L Z_J c_{IJ} + \sum_{J=L+1}^N Z_J c_{IJ} \right] \\
  \partial w({\bf r}_I)/\partial \lambda &=& w({\bf r}_I)/(\lambda-1)
\end{eqnarray}
having defined
\begin{eqnarray}
  \rho^\alpha ({\bf G}) &=& 
    \mbox{Fourier transform of } \rho^\alpha ({\bf r})\\
  \kappa_\alpha &=&
  \left( 1-2\alpha \right) + \left[ \frac{1}{2} - \frac{1+\theta^2}{4\theta}
  \ln\left| \frac{1+\theta}{1-\theta} \right| +
  6 \theta^2 \chi^2 \right] / \left[
  \chi \left( 1-\chi-3\chi\theta^2 \right) \right] \nonumber \\
  \chi &=& \frac{1}{2} + \frac{1-\theta^2}{4\theta}\ln\left|
  \frac{1+\theta}{1-\theta} \right| \nonumber \\
  c_{IJ} &=& \frac{Z_I}{2} \left\{ \sum_{\bf A}
   \left. \mbox{erfc} (\eta x)/x \right|_{{\bf x}={\bf A}+{\bf R}_I+{\bf R}_J}
   + \sum_{\bf G} \frac{4\pi}{\Omega G^2} \cos\left[ {\bf G}
    \left( {\bf R}_I - {\bf R}_J \right) \right] e^{-G^2/4\eta^2} \right\}.
\end{eqnarray}

\begin{references}

\bibitem[*]{present} Present address: Laboratoire de Physique,
Ecole Normale Sup\'erieure de Lyon, 46 All\'ee d' Italie, 69364 Lyon
C\'edex 07 (F).

\bibitem{frenkel86} D. Frenkel, in {\it
Proceedings of the International School of Physics ``Enrico Fermi'' ---
Course XCVII}, edited by G. Ciccotti and W.~G. Hoover  (North-Holland,
Amsterdam,1986).

\bibitem{bennett75} C.~H. Bennett, in {\it Diffusion in Solids},
edited by A. S.Nowick and J.~J. Burton  (Academic Press, New York, 1975).

\bibitem{frenkel84} D. Frenkel and A.~J.~C. Ladd,
J. Chem. Phys. {\bf 81}, 3188 (1984).

\bibitem{squire69} D.~R. Squire and W.~G. Hoover,
J. Chem. Phys. {\bf 50}, 701  (1969).

\bibitem{broughton83} J.~Q. Broughton and G.~H. Gilmer,
J. Chem. Phys. {\bf 79}, 5095  (1983).

\bibitem{carpar85} R. Car and M. Parrinello,
Phys. Rev. Lett. {\bf 55}, 2471 (1985).

\bibitem{galli93} G. Galli and A. Pasquarello, in {\it Computer Simulations in
Chemical Physics}, edited by M.~P. Allen and D.~J. Tildesley (Kluwer,
Dordrecht, 1993).

\bibitem{our} M. Pearson, E. Smargiassi, and P.~A. Madden,
J. Phys. -- Cond. Matter {\bf 5}, 3321 (1993).

\bibitem{tesi} E. Smargiassi, PhD thesis, Scuola Internazionale di Studi
Avanzati, Trieste (I), 1990; R. Car, P.~E. Bl\"ochl, and E.
Smargiassi, Mater. Sci. Forum {\bf 83-87}, 443 (1992).

\bibitem{blochl93} P.~E. Bl\"ochl, E. Smargiassi, R. Car, D.~B. Laks,
W. Andreoni, and S.~T. Pantelides, Phys. Rev. Lett. {\bf 70}, 2435
(1993).

\bibitem{devita91} A. Devita and M.~J. Gillan, J. Phys. -- Cond.
Matter {\bf 3}, 3929 (1991).

\bibitem{jacucci84} G. Jacucci, in {\it Diffusion in Crystalline Solids},
edited by G.~E. Murch and A.~S. Nowick (Academic Press, New York, 1984).

\bibitem{topp73} W.~C. Topp and J.~J. Hopfield,
Phys. Rev. B {\bf 7}, 1295 (1973).

\bibitem{remler90} D.~K. Remler and P.~A. Madden,
Mol. Phys. {\bf 70}, 921 (1990).

\bibitem{perrot94} F. Perrot, J. Phys. -- Cond. Matt. {\bf 6}, 431 (1994).

\bibitem{wang92} L.-W. Wang and M.~P. Teter, Phys. Rev. B {\bf 45}, 13197
(1992).

\bibitem{dreizler90} R.~M. Dreizler and E.~K.~U. Gross,
{\it Density Functional Theory} (Springer-Verlag, Berlin, 1990).

\bibitem{pastore91} G. Pastore, E. Smargiassi, and F. Buda, 
Phys. Rev. A {\bf 44}, 6334 (1991).

\bibitem{blochl92} P.~E. Bl\"ochl and M. Parrinello, Phys. Rev. B {\bf
45}, 9413 (1992).

\bibitem{our1} E. Smargiassi and P.~A. Madden, Phys. Rev. B {\bf 49}, 5220
(1994).

\bibitem{noseh} S. Nos\`e, J. Chem. Phys. {\bf 81}, 561 (1985); W.~G.
Hoover, Phys. Rev. A {\bf 31}, 1695 (1985).

\bibitem{fois93} E.~S. Fois, J.~I. Penman, and P.~A. Madden, J. Chem. 
Phys. {\bf 98}, 6361 (1993).

\bibitem{martyna92} G.~J. Martyna, M.~L. Klein, and M. Tuckerman,
J. Chem. Phys. {\bf 97}, 2635 (1992).

\bibitem{numrecip} W.~H. Press, B.~P. Flannery, S.~A. Teukolsky, and
W.~T. Vetterling, {\it Numerical Recipes} (Cambridge University Press,
Cambridge, 1989).

\bibitem{fermi55} E. Fermi, J.~G. Pasta, and S.~M. Ulam,
LASL Rep. LA-1940 (1955); G. Benettin, in {\it Proceedings of the
International School of Physics ``Enrico Fermi'' --- Course XCVII},
edited by G. Ciccotti and W.~G. Hoover (North-Holland, Amsterdam, 1986).

\bibitem{andersen80} H.~C. Andersen, J. Chem. Phys. {\bf 72}, 2384
(1980).

\bibitem{kohler88} U. Kohler and C. Herzig, Philos. Mag. A {\bf 58}, 769 
(1988).

\bibitem{nielsen85} O.~H. Nielsen and R.~M. Martin, Phys. Rev. B
{\bf 32}, 3780 (1985).

\end{references}

\begin{figure}
\caption[]{Dotted line: probability distribution of the atomic displacement
of a harmonic oscillator in the microcanonical ensemble, at $T=200~K$.
Dashed line: with Nos\`e--Hoover thermostat. Long dashed line: with
a Nos\`e--Hoover chain of two thermostats. Solid line: exact canonical
distribution.}
\label{fighar}
\end{figure}

\begin{figure}
\caption[]{Individual atomic energies as a function of time in a system
of 16 uncoupled harmonic oscillators. A Nos\`e-Hoover thermostat is
used, $T = 200~K$.}
\label{fignh}
\end{figure}

\begin{figure}
\caption[]{Individual atomic energies as a function of time in a system
of 16 uncoupled harmonic oscillators. A massive Nos\`e-Hoover chain of
2 thermostats per degree of freedom is used, $T= 200~K$.}
\label{figmnhc}
\end{figure}

\begin{figure}
\caption[]{Dotted line: probability distribution of the displacements of
the atoms for the system of Fig. \ref{figmnhc} (averaged over 24.2
ps).  Solid line: exact canonical distribution. $T=200~K$.}
\label{figdist}
\end{figure}

\begin{figure}
\caption[]{Cumulative average of $\Phi$ in a typical run for $\Delta F$
($T=200~K$, $\lambda=\lambda_2^4$). The three curves represent
averages starting at different times along the same run (0, 2, 4 ps).}
\label{figdeltaf}
\end{figure}

\begin{figure}
\caption[]{Dotted line: probability distribution of the atomic displacement
of the decoupled atom in a test run for the vacancy formation energy,
at $T=200~K,~\lambda=\lambda_1^2$ (first 6 ps). Dashed line: from 6 to
12 ps.  Solid line: theoretical distribution the corresponding
Einstein oscillator at the same $T$ (see text).}
\label{figlam}
\end{figure}

\begin{figure}
\caption[]{Cumulative average of $\Phi$ in a typical run for the vacancy:
$T=200~K$, $\lambda=\lambda_1^2$ (see Fig. \ref{figdeltaf}).}
\label{figphi}
\end{figure}

\begin{figure}
\caption[]{Free energy of formation for the sodium vacancy as a function of
$T$. The error bars reflect the estimated overall uncertainty of the
results ($\pm 0.02$ eV).}
\label{figfree}
\end{figure}

\narrowtext
\begin{table}
\caption[]{Excess free energy of the crystal with respect to the Einstein
crystal as a function of the temperature 
  using a massive Nos\`e--Hoover chain of two
thermostats per degree of freedom.}
\begin{tabular}{ccc}
 $T (K)$ &  {$\Delta F$ (eV)} \\
 \tableline
   100    & -0.019 \\
   200    & -0.031 \\
   300    & -0.042 \\
   370    & -0.052 \\
\end{tabular}
\label{tabfree}
\end{table}

\end{document}
 

