\documentstyle[preprint,aps]{revtex}

\begin{document}

\draft

\title{First-principles free-energy calculations on condensed matter
systems: lattice vacancy in silicon}

\author{Enrico Smargiassi} 

\address{Centre Europ\'een de Calcul Atomique et Mol\'eculaire, ENS Lyon,
69364 Lyon CEDEX 07, France}

\author{Roberto Car}

\address{Departement de Physique, University of Geneva, Switzerland}

\date{\today}

\maketitle

\begin{abstract}
We illustrate a method for performing first-principles free-energy
calculations within the Kohn-Sham scheme. We show that the method can
be used in cases in which electrons have to be decoupled from the
system and illustrate it with results for the formation free-energy of
the Si vacancy. The results agree well with the available data from
experiments and {\sl ab-initio} calculations.
\end{abstract}

\pacs{61.72.Ji 02.70.Ns 61.72.Bb} % MD & particle methods; Pt defects
                                  % & def. clusters; Theory & models
                                  % of crystal def.

\section{Introduction}

In recent years there has been growing interest in the calculation of
thermal properties of condensed matter via computer simulation. The
most general and direct way is to calculate an appropriate
free-energy, from which all the equilibrium thermodynamic properties
can be obtained. For instance, the knowledge of the free-energy of
formation $F_f$ of a defect allows one to compute its equilibrium
concentration $c=e^{F_f/k_B T}$ which, in addition to being important
in itself, also enters the defect-induced diffusion coefficient and
other important quantities.

\par

Unfortunately absolute free-energies are not directly measurable,
neither experimentally or in a computer simulation, as they contain an
entropic part which cannot be computed in any conceivable amount of
time. Therefore, indirect methods must be used \cite{frenkel86}. A
popular choice is ``thermodynamic integration''. In this method one
takes the potential energy $U$ of the system studied to be a smooth
function of a parameter $\lambda$. It is then easy to show that
$\frac{\partial F}{\partial \lambda} = \langle \Phi (\lambda) \rangle
$ (where $\Phi (\lambda) = \frac{\partial U}{\partial \lambda}$ and
the angle parentheses indicate thermal averaging); the free-energy
difference between the states corresponding to the values $\lambda=1$
and $\lambda=0$ can then be expressed as $\int_0^1 \langle \Phi
\rangle \, d\lambda$. As $F$ is a function of state, the states at
intermediate values of $\lambda$ do not need to have any physical
meaning: only the smoothness of $U(\lambda)$ and its values at the
endpoints are relevant.

\par

This method has been used widely in the past on simple model systems
like Lennard-Jones \cite{jacucci84}, showing that a proper account of
anharmonic effects is essential for a full understanding of diffusion
processes. On the other hand, similar calculations on important
materials like semiconductors have been hindered to date by the lack
of sufficiently reliable interatomic potentials. In fact, even
well-established potentials may easily lead to unphysical results for
defects if not constructed and handled with great care
\cite{ghaisas93}. And even in this cases, the results are not
sufficiently accurate to allow one to draw reliable quantitative
conclusions, as it can be seen by comparing the results of
Ref. \cite{maroudas93} with those obtained {\sl
ab-initio} \cite{bloechl93}. The alternative route of using {\sl
ab-initio} methods has for long been confined to $T=0$ or, at best,
harmonic-like approximations \cite{clark94} due to the huge
requirements in terms of computer resources. To overcome these
problems, first-principles molecular dynamics (FPMD)
schemes\cite{carpar85} have been developed, which make calculations of
thermal properties of real materials from first principles
affordable. Therefore the extension to FPMD of the above methods for
calculating free-energies is now an important issue.

\par

The characteristics of the thermodynamic integration method depend in
an important way on the transformation which is chosen. The parameter
$\lambda$ may, in principle, be anything. It can be chosen as a
physical parameter of the system, in which case the system is driven
through a physically realizable transformation. In this case the
extension of the method to first-principles calculations is
straightforward: all one has to do is to compute energies and forces
using one of the existing {\sl ab-initio} schemes instead of using
classical interatomic potentials. The calculation of the thermal
average of $\Phi$ in a first-principles scheme is in this case no more
complicated than the calculation of the forces acting on the ions
(actually $\Phi$ can be considered a ``generalized force''): First,
one computes, analytically, the derivative of the {\sl functional}
$U(\{ {\bf r}_i\})$ with respect to $\lambda$. Second, one performs an
MD simulation at the chosen $\lambda$-points, taking care of keeping
the system as close as possible to the instantaneous Born-Oppenheimer
(BO) ground state. By invoking the Hellman-Feynman theorem and
assuming ergodicity in the subspace of the ionic coordinates one can
then show that the calculated $\langle \partial U / \partial \lambda
\rangle $ coincides with the desired ensemble average. In a recent
work \cite{milman93} for instance, $\lambda$ was taken to be the
reaction coordinate of an atom jumping across a potential barrier. In
this case $\Phi$ was the component of the force acting on the jumping
atom along the reaction path, thus obtaining the free-energy of
migration across the barrier.

\par

Many interesting application however require that the system be led
through an unphysical transformation. In several important cases one
or more atoms have to be decoupled from the others, thus $\lambda$ is
the coupling of one or more atoms to the rest of the crystal. This is
the case of the calculation of formation free-energies of defects, or
of free-energies of extended systems. When using {\sl ab-initio}
schemes, where electrons are treated explicitly, some or all of them
must then be decoupled, together with the corresponding ions. This
procedure has no counterpart in the approaches that make use of
effective interatomic potentials only (a discussion of this point in
general terms is given in Ref. \cite{smargi95}) and makes the whole
exercise quite more complicated. As it will be shown in this paper in
fact, the applicability of the thermodynamic integration method in the
framework of the Kohn-Sham formalism to these unphysical
transformations is not obvious. We had already shown that it is
nevertheless feasible \cite{tesi} and used the method to calculate the
free-energy of formation of the vacancy in Si \cite{sugino95}. The
results have been published previously \cite{car92}, therefore in this
paper we mainly focus the discussion on the technical aspects of our
work. We treat in detail the case of the vacancy, but modifications of
the method for use in other cases are straightforward. To simplify
matters further, we shall discuss only the case of simulations at
constant lattice parameter.

\section{Formalism}

The formation free-energy $F^N_f$ of a vacancy as calculated by using
an $N$ atoms supercell is
\begin{equation}
  F^N_f (\Omega_C,T) =
   F^{N-1}_V (\Omega_V,T) - \frac{N-1}{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-1}_V
(\Omega_V,T)$ is the free energy of a crystal with the same number of
lattice sites but one of them being vacant. After choosing a suitable
single-body reference state defined by a potential energy $u$ we
define now
\begin{itemize}
  \item $U^N(\lambda = 1) = U^N_C = $ potential energy of the perfect crystal,
  \item $U^N(\lambda = 0) = U^{N-1}_V + u = $ potential energy of the vacancy
                             crystal plus the energy of the decoupled atom.
\end{itemize}

In other words, one ion is adiabatically removed from the crystal and
put in the state described by $u$, while at the same time its valence
electrons are made disappear. The free energies of the bulk system and
of the vacancy system are thus respectively
\begin{eqnarray}
  F^N_C     & =  & F^N (\lambda = 1) \nonumber \\
  F^{N-1}_V & =  & F^N (\lambda = 0) - \varphi
\end{eqnarray}

where $\varphi$ is the free-energy of an atom in a potential $u$, at
the appropriate temperature. This is a well defined quantity as long
as the potential $u$ is a single-body potential. Clearly the form of
$u$ is arbitrary: its free-energy $\varphi$ is implicitly added in
$F^N (\lambda=0)$ and then explicitly subtracted. For calculations on
solids, $u$ is well chosen as a harmonic potential, which has the
additional advantage that in this case $\varphi$ is analytically
computable. We use the following form:
\begin{equation}
 u = \varepsilon_0 + \frac{1}{2} M_I \omega^2 \left( {\bf R - R_0} \right)^2
\end{equation}
where $\varepsilon_0$ is the energy of an atom in the perfect crystal
at $T=0$, ${\bf R}_0$ is the lattice site of the atom to be made
disappear (i.e. that in which the vacancy will appear), $M_I$ the
ionic mass and $\omega$ is a typical lattice frequency, which we chose
to be the experimental Debye frequency. In this case the decoupled
atom behaves approximately as a bulk atom.

\par

We now exploit the fact that $\Phi$ has the form of a canonical average
and can therefore be computed by MD:
\begin{eqnarray}
 \frac{\partial F(\lambda)}{\partial \lambda} & = & - k_B T
 \frac{\partial}{\partial \lambda} \ln \int {\rm d}^{3N}R ~ e^{-U
 (\lambda)/k_B T} \label{dfdl} \\
 & = & \frac{1}{{\cal Z}_q (\lambda)}
 \int {\rm d}^{3N}R ~ \Phi (\lambda) \, e^{-U (\lambda)/k_B T} =
 \langle \Phi (\lambda) \rangle \nonumber
\end{eqnarray}
where ${\cal Z}_q$ is the configurational integral entering the
partition function. Then,
\begin{equation}
 F(\lambda=1) - F(\lambda=0) = \int^1_0 {\rm d}
  \lambda \left\langle \Phi (\lambda) \right\rangle.
  \label{f1f0}
\end{equation}
The formation free-energy is finally obtained as
\begin{equation}
  \label{effelambda} F_f = F (\lambda=0) - F (\lambda=1) + \frac{1}{N}
   F_C - \varphi = - \int^1_0 {\rm d}\lambda \, \left\langle \Phi
   (\lambda) \right\rangle + \Delta F.
\end{equation}
The correction factor $\Delta F$ is the difference between $\varphi$,
which is known independently (here analytically), and the free energy
per atom in the perfect crystal. The latter is easily obtained from
accurate {\sl ab-initio} calculations of phonons or directly from
experimental data.

\par

The next step is now to find a potential energy $U(\lambda)$
satisfying the appropriate conditions at $\lambda=0$ and
$\lambda=1$. We chose it in the following way:
\begin{itemize}
 \item We multiply by $\lambda$ the occupation numbers $f_i$ of the
  states occupied by $Z$ electrons ($Z=4$ being the valence of the Si
  ion).
 \item We multiply by $\lambda$ the ionic charge $Ze$ of the ion to be
  decoupled and its bare pseudopotential.
\end{itemize}

\par

To compute the thermal averages needed, we use FPMD with the standard
Car-Parrinello (CP) method \cite{carpar85} with the following
Lagrangian:
\begin{equation}
 {\cal L} = \frac{1}{2} \sum_I^N M_I {\bf \ddot R}_I +
  \frac{1}{2} \sum_i \mu_i \left\langle \dot \psi_i | \dot \psi_i \right\rangle
  - E_{KS} [\{ {\bf R}_I \}, \{ \psi_i \}] - E_{II} + \mbox{constraints}
 \label{lagcp}
\end{equation}
where $E_{KS}$ is the Kohn-Sham energy functional, $E_{II}$ is the
direct ion-ion interaction and the ``constraint'' term imposes the
orthonormality condition on the wavefunctions.

The electron pseudo-masses $\mu_i$ must be chosen in such a way that
the pseudo-electronic frequencies are much higher than all ionic
frequencies\cite{pastore91}. When $Z$ electrons are decoupled however
the occupation numbers of the corresponding states become smaller by a
factor $\lambda$ and it is easy to see\cite{pastore91} that for
$\lambda \to 0$ this requirement will be violated. We therefore choose
to set the pseudo-electronic masses as follows: for the electrons
which stay fully coupled, $\mu_i=\mu=$ constant, for the electron
states corresponding to disappearing electrons $\mu_i=\lambda\mu$ (it
should be noticed that this affects the form of the orthonormalization
procedure: the appropriate form is given in ref. \cite{tesi}). With
this choice the less coupled electrons have correspondingly low masses
and thus the frequency of their motion remains high.

\section{The Si vacancy}

We applied the above scheme to the case of the Si vacancy. We used a
64 atom cell and the standard ``local density approximation plus
planewaves plus pseudopotential'' scheme. We used the $\Gamma$ point
to sample the Brillouin zone, and a plane-wave cutoff of 8
Rydberg. This values are sufficient to give phonon frequencies within
10\% of the converged value. We checked that using a sufficient number
of $\lambda$ points the results obtained from Eq. \ref{dfdl} at $T=0$
coincide with the ones obtained directly from
eq. \ref{free_nolambda}. This is indeed the case, within 0.05 eV, if
one uses a 3-point Gauss-Legendre formula. It is possible to show that
the states corresponding to the electrons being decoupled are
necessarily the ones with the highest orbital energies. In our case
these states become, at decreasing $\lambda$, the two highest states
of the three introduced by the vacancy in the band gap. We verified
that this is indeed remains true even during the dynamical runs.

\par

We verified that the CP procedure gives correct results for $\Phi$ by
performing a series of calculations as follows. First, we evolved the
system according to the CP lagrangian using the same parameters which
were subsequently used in the calculations: $\mu=875$ a.u. and $\Delta
t=10$ atu (atomic time units; 1 atu $= 2.4188~10^{-5}$ ps). Second, we
evolved the same system from the same initial conditions by minimizing
explicitly the electronic energy at every ionic step by using $n$
steepest descent steps and computing the forces on ions via the
Hellman-Feynman theorem; this must give the correct ionic evolution in
the BO approximation in the limit $n\to\infty$. The results are shown
in Fig. \ref{nstepismo}. The slow oscillation, which is present in all
curves, comes from the ionic dynamics. The fast oscillation present
only in the CP curve is due to the fast electronic motion around the
instantaneous ground state. Apart from this oscillation, which should
average to zero, it is seen that as $n$ increases the broken curves
converge to the full curve and for $n \ge 20$ they already lie within
the oscillation band of the fast motion.

\par

On longer runs we however observed that the usual criteria used to
check for drifts from the BO surface are not sufficiently tight to
guarantee accurate results for $\Phi$; drifts which could be tolerated
in other simulations are not acceptable here. Indeed, $\Phi$ is a
differential quantity and as such it can be expected to be prone to
large numerical errors. In this particular case the problem arises
because $\Phi$ mainly depends on the behavior of a small number of
electrons, i.e. those which are being decoupled, which being the
highest in energy are also the ones which are more likely to create
problems with the adiabatic behavior \cite{pastore91}. A good
strategy would be to force the electrons to stay close to the ground
state by imposing a thermostat on them\cite{nosele}. We used a simpler
approach and just quenched periodically the electronic degrees of
freedom on the instantaneous BO surface (one minimization every 100
steps). Of course this implies that the system of electrons plus ions
is not conservative anymore and, in fact, energy is gradually
transferred from the hot ions to the cold electrons. We kept the ionic
temperature constant by using a Nos\'e thermostat \cite{nose85},
coupled to the ionic degrees of freedom only. This also ensured that
under the assumption that the ionic subsystem is ergodic -- for a
discussion of this point see Ref. \cite{smargi95} -- the sampling of
the ionic phase space is indeed canonical. We found that our method is
sufficient to give accuracy on $\Phi$ better than 0.1 eV. We are then
satisfied that, if the system is evolved via the CP dynamics and care
is taken that the electrons do not drift away from the BO surface, the
values obtained for $\Phi$ in this way are correct.

\par

We run the program for two different temperatures, $T=500$ and
$T=1000$, and the 3 Gauss-Legendre points of the interval $[0,1]$. We
used a fixed cell whose lattice parameter $a$ was chosen to be the
experimental one at $T=0$. By performing a $T=0$ calculation of $F_f$
at the value of $a$ corresponding to a temperature of 1400 $K$ we
checked that thermal expansion would give a correction of the order of
0.1 eV or less. The value of $F_f$ at $T=0$, which gives the largest
contribution, was taken to be the one obtained by accurate calculation
of the formation energy \cite{bloechl93} and the dynamical runs were
only used to compute the temperature-dependent part. We checked that a
simulation time of 6 ps is sufficient to give reasonable
thermalization and statistics, i.e. formation entropies with a
statistical uncertainty of $ \sim 1~k_B$, this agrees with the
findings of Ref. \cite{smargi95}.

\par

The results have been previously published \cite{car92} and here we
only discuss them. The value at $T=500~K$ agrees with the value
obtained \cite{bloechl93} by using the same model of Si but using the
Local Harmonic Approximation \cite{najafabadi89}, which shows that the
anharmonic effects are still small at that temperature. Also, $F_f$
decreases with increasing $T$. This is to be expected, as it is
reasonable that the entropy of the defective crystal is larger than
the perfect crystal and thus its free energy decreases more
rapidly. It also looks like that $F_f$ is bent downward; this is an
expected effect too, since the anharmonic effects, which become more
and more important as $T$ increases, are arguably larger in the
defective crystal and so its entropy should increase faster than the
entropy of the perfect crystal. In contrast, values calculated by
using $T=0$ calculations in connection with standard
statistical-mechanics methods \cite{clark94} do not show any
temperature dependence of the entropy.

\par

It should be noted that the calculated values only include the ionic
free-energy. As the vacancy introduces closely spaced electronic
levels in the band gap, only one of whom is occupied at $T=0$, there
has to be a sizeable contribution to $F_f$ coming from electronic
excitations. To calculate this electronic contribution one must have a
good description of excited states, which is not the case within local
density approximation. This contribution however has been estimated to
be approximately 1$k_B$ \cite{bloechl93}. Even neglecting this, the
calculated values for $F_f$ agree well with the experimental estimates
\cite{dannefaer85}. Our data therefore show that the observed high
entropy of vacancies is a direct consequence of the vacancy-induced
phonon softening. This may possibly be true for other point defects as
well.

\section{Conclusions}

We have shown how the thermodynamic integration method can be used
within the Kohn-Sham scheme, including the cases in which the initial
and final states have a different number of electrons and thus some
electrons have to be decoupled from the rest of the crystal. The
implementation of the method requires some care to avoid spoiling the
adiabatic separation between the ionic and pseudo-electronic motion;
we showed how these problems can be successfully overcome and that
good results can be achieved.

\par

We benefited from stimulating discussions with P.~E. Bl\"ochl,
S.~T. Pantelides, G. De Lorenzi, A. Voter and M. Parrinello. This work
has been supported in part by the International School for Advanced
Studies (SISSA) - Centro Interuniversitario del Nord-Est per il
Calcolo Automatico (CINECA) colaborative project, under the
sponsorship of the Italian Ministry for Public Education, by the
European ARO grant N. DAJA45-89-C-0025 and by the Italian National
Research Council (CNR) grant N. 89.0006.69.

\begin{references}

\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{jacucci84} G. Jacucci, in {\it Diffusion in Crystalline
Solids}, edited by G.~E. Murch and A.~S. Nowick (Academic Press,
Orlando, 1984).

\bibitem{ghaisas93} S.~U. Ghaisas, Phys. Rev. B {\bf 43}, 1808 (1993).

\bibitem{maroudas93} D. Maroudas and R.~A. Brown, Phys. Rev. B {\bf
47}, 15562 (1993).

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

\bibitem{clark94} S.~J. Clark and G.~J. Ackland, Phys. Rev. B {\bf
48}, 10899 (1994).

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

\bibitem{milman93} V. Milman {\sl et al.}, Phys. Rev. Lett. {\bf 70},
2928 (1993).

\bibitem{smargi95} E. Smargiassi and P.~A. Madden, Phys. Rev. B {\bf 51},
117 (1995).

\bibitem{tesi} E. Smargiassi, PhD thesis, Scuola Internazionale di Studi
Avanzati, Trieste (I), 1990.

\bibitem{sugino95} A similar approach has also been used to calculate
the melting point of Si {\sl ab-initio} [O. Sugino and R. Car,
Phys. Rev. Lett. {\bf 74}, 1823 (1995)].

\bibitem{car92} R. Car, P.~E. Bl\"ochl, and E. Smargiassi, Mater. Sci.
Forum {\bf 83-87}, 443 (1992).

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

\bibitem{nosele} P.~E. Bl\"ochl and M. Parrinello, Phys. Rev. B {\bf 45},
9413 (1992); E.~S. Fois, J.~I. Penman, and P.~A. Madden, J. Chem. Phys.
{\bf 98}, 6361 (1993).

\bibitem{nose85} S. Nos\`e, J. Chem. Phys. {\bf 81}, 561 (1985).

\bibitem{najafabadi89} R. Najafabadi, R. LeSar, and D.~J. Srolovitz,
Phys. rev. Lett. {\bf 63}, 624 (1989).

\bibitem{dannefaer85} S. Dannefaer, P. Mascher, and D. Kerr,
Phys. Rev. Lett. {\bf 56}, 2195 (1985).

\end{references}

\begin{figure}

\caption{$\Phi$ as calculated using the CP dynamical procedure (full
line). The dashed, dotted, and dashed-dotted lines represent
calculations done by minimizing the total electronic energy at every
ionic step, using $n$ steepest descent steps per each ionic
step. $\Delta t =$ 10 a.u., $\mu = $ 875 a.u..}
\label{nstepismo}
\end{figure}

\end{document}
