\documentstyle[preprint,aps]{revtex}
\newtheorem{theorem}{Theorem}

\begin{document}

\title{Density Functional Methods at Finite Temperature}

\author{Jean-Pierre Hansen} 

\address{Laboratoire de Physique, Unit\'e de Recherche Associ\'ee 1325
du CNRS, Ecole Normale Sup\'erieure de Lyon, 69364 Lyon CEDEX 07
(France)}

\author{Enrico Smargiassi}

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

%\date{\today}

\maketitle

\begin{abstract}

After a brief summary of Density Functional Theory (DFT) for both
classical and quantum inhomogeneous systems, the theory is applied to
the simulation of systems having widely different time and/or length
scales. These include ``ab-initio'' Molecular Dynamics and Monte Carlo
simulations of colloidal suspensions, polyelectrolytes and metallic
systems.

\end{abstract}

\section{Inhomogeneous systems and effective interactions}

Density Functional Theory (DFT) is the art of dealing with {\it
inhomogeneous} many-particle systems of Condensed Matter Science, in
which translational and/or rotational invariance is broken either
spontaneously (as in crystals) or by an external field. The main
objective of these notes is to show how DFT can be put to use in the
simulations of highly asymmetric systems, made up of at least two
species of particles with widely different masses and sizes. The
resulting difference in time and length scales does not allow the
``large'' and ``small'' particles to be described on an equal footing
within traditional Monte Carlo (MC) or Molecular Dynamics (MD)
algorithms. It will often prove more efficient to adopt a discrete
(``particle'') representation of the large particles, referred to
hereafter as species 1, which will always be assumed to obey classical
Statistical Mechanics, while a continuum (``field'') representation is
used for the small particles (species 2), which may be either
classical or partially degenerate. The basic field for the latter is
the one-particle density (or density profile) $\rho ({\bf r})$.

\par

Two examples of such asymmetric systems are polyelectrolytes, made up
of mesoscopic polyions and microscopic counterions (both classical),
and ion-electron systems, including semiconductors, metals and
plasmas, where the ions can usually be treated as classical, while the
electrons are almost invariably degenerate.

\par

A formal reduction of the initial two-component system to an effective
one-component system, involving only the degrees of freedom of the
large particles, may be carried out as follows; it is worth stressing
that this derivation does not depend in any crucial way on the
asymmetry between the two two species. For the sake of simplicity, all
particles will be assumed to have translational degrees of freedom
only. Let $\{ {\bf R}_i, {\bf P}_i \},\; 1\le i \le N_1$ and $\{ {\bf
r}_j, {\bf p}_j \}, \; 1\le j \le N_2$ denote the sets of coordinates and
momenta of the large and small particles respectively. The total
Hamiltonian of the system may be split into three terms:

\begin{equation}
H = H_1(\{ {\bf R}_i, {\bf P}_i \}) + H_2 (\{ {\bf r}_j,{\bf p}_j \})
  + V_{12}(\{ {\bf R}_i, {\bf r}_j \}).
\label{totham}
\end{equation}

The parts $H_\alpha$ ($\alpha=1,2$), made up of kinetic and potential
energy contributions $(H_\alpha = K_\alpha+V_{\alpha\alpha})$, involve
the coordinates and momenta of species $\alpha$ only, while $V_{12}$
is the potential energy of interaction between particles of different
species. The canonical partition function of the total system is:

\begin{equation}
Q_{N_1,N_2} = \mbox{Tr} \, e^{ -\beta (H_1 + H_2 + V_{12})}
\label{total_partition}
\end{equation}

where $\beta = 1/k_B T $ is the inverse temperature ($k_B$ is the
Boltzmann's constant), and the trace is a sum over a complete set of
quantum states, which reduces to a phase-space integration for degrees
of freedom which may be treated classically. Since this is the case
for the large particles (species 1), the trace operator factorizes
according to:

\begin{equation}
Q_{N_1,N_2} = \mbox{Tr}_1 \left[ e^{-\beta H_1}
 \mbox{Tr}_2 e^{-\beta (H_2+V_{12})} \right].
\label{factored_partition}
\end{equation}

The second trace, over the degrees of freedom (or quantum states) of
species 2, is just the partition function $ Q_{N_2} (\{ {\bf R}_i \})$
of the {\it inhomogeneous} system of small particles in the
``external'' field of the large particles frozen at positions $\{ {\bf
R}_i\}$. The corresponding Helmholtz free energy $F_2 (\{ {\bf R}_i\})
= - k_B T \log (Q_{N_2} (\{ {\bf R}_i\}) )$ depends parametrically on
the $\{ {\bf R}_i\}$, and the total partition function now reduces to:

\begin{equation}
Q_{N_1,N_2} = \mbox{Tr}_1 e^{-\beta [H_1 + F_2(\{ {\bf R}_i\})]} =
\frac{1}{N_1!\Lambda^{3N_1}} \int d^{N_1} {\bf R} e^{-\beta
V^{eff}_{11}(\{ {\bf R}_i\})}
\label{effective_partition}
\end{equation}

where the factor $\Lambda^{-3N_1}$ ($\Lambda$ denotes the de Broglie
wavelength of a large particle) results from the trivial integration
over the (classical) momenta of species 1, while the {\it effective}
potential energy of interaction between particles of that species is
the sum of the direct term $V_{11}$ and an induced contribution $F_2$:

\begin{equation}
V_{11}^{eff} (\{ {\bf R}_i\}) = V_{11} (\{ {\bf R}_i \}) + F_2
(\{ {\bf R}_i \}).
\label{pot_eff}
\end{equation}

It is important to stress that, while $V_{11}$ may be essentially
pair-wise additive (as is the case, e.g., for direct coulombic
interactions), the induced contribution $F_2$ is state-dependent and
generally many-body in nature.

\par

The above reduction is formally exact, but an effective calculation
scheme hinges on the formidable problem of obtaining accurate
estimates of the free energy of species 2, for any configuration of
species 1. In practice this is achieved within the density functional
formalism, which relies in an essential way on a uniqueness theorem
and a variational principle, to be explicitly formulated in the next
section. In connection with the present reduction scheme, these
theorems imply that the free energy $F_2$ is a {\it functional} of the
density profile $\rho_2({\bf r})$ of species 2. This functional
reaches its minimum when its profile reaches its equilibrium value for
any configuration $\{{\bf R}_i \}$ of large particles [which produces
the external potential via $V_{12}$; cf. Eq.(\ref{totham})].

\par

The reduction process may hence be implemented, at least in principle,
by standard minimization procedures, within MC or MD simulations. 

\par

In MC simulations, configurations of a system of large particles may
be sampled from the Boltzmann distribution $P_{N_1}(\{ {\bf R}_i \})
\sim e^{-\beta V_{11}^{eff}(\{ {\bf R}_i \})}$ using the usual
Metropolis algorithm, provided $F_2[\rho({\bf r})]$ can be minimized
efficiently, for each configuration $\{ {\bf R}_i \}$, by a steepest
descent or conjugate gradient method. An illustration of this
procedure will be reported in section V.

\par

An efficient MD implementation of the reduction process relies
crucially on a clear separation of time scales\cite{pastore91}, which
requires $m_2/m_1 <\!\!< 1$. Such a separation validates adiabatic
dynamics, whereby the physical dynamics of the large (heavy) particles
and the dynamic minimization of the free energy of the small (light)
particles are carried out in parallel, using some standard finite
time-step solver of coupled (classical) equations of motion
(Car-Parrinello method\cite{carpar85}). Illustrations of this
``ab-initio'' MD methodology will be given in sections IV and VI. The
effective force acting on particle $i$ $(1\leq i \leq N_1)$ follows,
at each step, from an implicit application of the Hellman-Feynman
theorem, provided particles of species 2 constantly remain
``sufficiently'' close to their ``Born-Oppenheimer'' (BO) surface
(i.e. the surface of lowest free-energy $F_2$). According to this
theorem:

\begin{equation}
{\bf F}_i = - \nabla_i V_{11}^{eff} (\{ {\bf R}_j \}) = - \nabla_i
V_{11} (\{ {\bf R}_j \}) - \langle \nabla_i V_{12}(\{ {\bf R}_j
\}, \{ {\bf r}_l \}) \rangle_{\{ {\bf R}_i \}},
\label{hell_feynm}
\end{equation}
where the angular bracket denotes an equilibrium average over the
degrees of freedom of species 2, for a fixed configuration $\{ {\bf
R}_i \}$ of species 1.

\section{Density functional theory: basic principles and standard 
approximations}

The reduction process of asymmetric systems, sketched in the previous
section, may be effectively carried out, provided an efficient
theoretical scheme is available for dealing with inhomogeneous systems
characterized by spatially modulated density profiles $\rho({\bf r})$.
Familiar examples of systems with broken translational symmetry are
the liquid-vapour interface, studied as far back as 1893 by van der
Waals\cite{vanderwaals93}, an electrolyte solution in the vicinity of
a planar electrode, a problem tackled in 1910 by Gouy\cite{gouy10}, or
a degenerate Fermi gas of valence electrons in the periodic field of
an ionic lattice. The common theoretical framework to describe such
inhomogeneous systems is DFT, pioneered by Percus\cite{percus62} for
classical fluids and by Hohenberg, Kohn and
Sham\cite{hohenberg64,kohn65} for the degenerate electron
gas. Excellent reviews are available, both for classical\cite{evans79} and
quantum\cite{parr89} DFT, and this section only provides a brief summary of
of the basic principles involved, and of some approximation schemes
which will be put to use in the perspective of computer
simulations. The general formulation will be given for a system at
equilibrium in an external field at finite temperature, and holds both
for classical and quantum statistics.

\par

The Hamiltonian of a system of $N$ identical particles subject to an
external potential $\phi({\bf r})$ reads:

\begin{equation}
H = K_N(\{ {\bf p}_i \}) + V_N (\{ {\bf r}_i\}) +
  \Phi_N(\{ {\bf r}_i \}).
\label{ham_ext}
\end{equation}

The contribution of the external field to the potential energy may be
cast in the form:

\begin{equation}
\Phi_N = \sum_{i=1}^N \phi({\bf r}_i) = \int_V d{\bf r}
 \tilde{\rho}({\bf r}) \phi({\bf r}),
\label{pot_en}
\end{equation}
where $\tilde{\rho}({\bf r}) = \sum_{i=1}^N \delta({\bf r}-{\bf r}_i)$
is the microscopic density operator.

Note that within the reduction scheme of a two-component system to an
effective one-component system, the external potential acting on the
``small'' particles is provided by the ``large'' particles, i.e.:

\begin{equation}
\phi({\bf r}) = \sum_{i=1}^{N_1} v_{12} ({\bf r}-{\bf R}_i)
\label{eff_phi}
\end{equation}
where $v_{12}$ is the interaction between one particle of species 1
and the particles of species 2. The considerations of the present
section are independent of such an identification. First consider the
inhomogeneous system at fixed canonical potential $\mu$, volume $V$
and temperature $T$. The grand partition function involves a summation
over all possible values of $N$:

\begin{equation}
\Xi (\mu,V,T) = \sum_{N=0}^\infty \mbox{Tr}_N\, e^{-\beta (H_N-N\mu)}
\label{grandpart}
\end{equation}
where $\mbox{Tr}_N$ refers to a classical or quantum-mechanical trace
over all states with a fixed value of $N$ (canonical trace).
The corresponding grand potential:

\begin{equation}
\Omega (\mu,V,T) = -k_B T \log [\Xi (\mu,V,T) ]
\label{grandpot}
\end{equation}
is clearly a {\it functional} of the external potential $\phi({\bf
r})$ or, more conveniently, of the ``intrinsic'' local chemical
potential $\Psi({\bf r}) = \mu - \phi ({\bf r})$. Two successive
functional differentiations of $\Omega$ with respect to $\Psi$ yield
the one-particle density and the two-point linear-response function
(or static susceptibility):

\begin{eqnarray}
-\frac{\delta \Omega [\Psi]}{\delta \Psi({\bf r})} & = & \frac{1}{\Xi}
\sum_N \mbox{Tr}_N \{ e^{-\beta (H_N -N\mu)} \tilde{\rho} ({\bf r})
\} = \langle \tilde{\rho} ({\bf r}) \rangle \equiv \rho ({\bf r})
\label{static_susc_1} \\
-\frac{\delta^2 \Omega[\Psi]}{\delta \Psi ({\bf r}) \delta \Psi ({\bf
r}^\prime)} & = & \frac{\delta \rho ({\bf r})}{\delta \Psi({\bf
r}^\prime)} = \chi({\bf r},{\bf r}^\prime).
\label{static_susc_2}
\end{eqnarray}

Eq. (\ref{static_susc_1}) defines the one-particle density (or density
profile) as the equilibrium statistical average of the microscopic
density, while Eq. (\ref{static_susc_2}) defines the susceptibility
$\chi$. Its physical significance is easily understood by considering
the effect of a small change in the external potential ($\phi({\bf r})
\to \phi ({\bf r}) + \delta \phi ({\bf r})$) on the density profile
($\rho({\bf r}) \to \rho ({\bf r}) + \delta \rho ({\bf r})$), at
constant chemical potential $\mu$. A functional Taylor expansion,
limited to linear order, then leads, according to
Eq. (\ref{static_susc_2}), to the static linear response result:

\begin{equation}
\Delta \rho({\bf r}) = \int d {\bf r}^\prime \chi({\bf r},{\bf
r}^\prime) \delta \Psi ({\bf r}^\prime) = - \int d {\bf r}^\prime
\chi({\bf r},{\bf r}^\prime) \delta \phi ({\bf r}^\prime)
\label{taylor_linear}
\end{equation}

For a homogeneous system ($\phi({\bf r}) = 0$), transitional
invariance implies that the density is uniform, i.e. $\rho({\bf
r})=n=N/V$, and that the susceptibility depends only on the relative
distance, i.e. $\chi ({\bf r},{\bf r}^\prime) = \chi ({\bf r}-{\bf
r}^\prime)$.

A standard perturbation calculation leads to the following expression
for $\chi$\cite{forster75}:

\begin{equation}
\chi ({\bf r},{\bf r}^\prime) = \frac{1}{k_B T} \int^{\beta\hbar}_0 \frac{d
\lambda}{\beta\hbar} \left\langle \Delta\tilde{\rho}({\bf r},-i\lambda)
\Delta\tilde{\rho}({\bf r}^\prime)\right\rangle
\label{chi_perturb}
\end{equation}
where $\tilde{\rho}({\bf r},i\lambda)$ denotes the Heisenberg
time-displaced density-fluctuation operator ($\delta\tilde{\rho} =
\tilde{\rho} - \rho$) taken at the imaginary time $t=i\lambda$. In the
classical ($\hbar \to 0$) limit, Eq. (\ref{chi_perturb}) reduces to:

\begin{equation}
k_B T \chi ({\bf r},{\bf r}^\prime) = \left\langle
\Delta\tilde{\rho}({\bf r}) \Delta\tilde{\rho}({\bf
r}^\prime)\right\rangle \equiv \rho({\bf r}) \left[ \delta({\bf r}-{\bf
r}^\prime) + \rho({\bf r}^\prime) h^{(2)}({\bf r},{\bf r}^\prime)
\right]
\label{chi_class}
\end{equation}
which defines the two-point correlation function
$h^{(2)}({\bf r},{\bf r}^\prime)$ \cite{hansen86}.

\par

It turns out to be much more fruitful to consider the density profile
$\rho({\bf r})$, rather than $\Psi({\bf r})$, to be the fundamental
field variable. This is exemplified by the fact that a system may be
spontaneously inhomogeneous at the microscopic level and thus be
characterized by a spatially modulated $\rho({\bf r})$, even in the
absence of an external field $\phi ({\bf r}) \to 0$, as is the case
for crystalline solids, where the local atomic density $\rho({\bf r})$
is periodic rather than uniform. Switching from $\Psi$ to $\rho$ as
independent variable ($V$ and $T$ being held constant) entails a
corresponding change from the thermodynamic potential $\Omega [\Psi]$
to to the potential ${\cal F} [\rho]$. This is achieved by a
generalized Legendre transformation:

\begin{equation}
{\cal F} [\rho] = \Omega [\Psi] + \int d {\bf r} \Psi({\bf r}) \rho ({\bf r})
\label{legendre_trasf}
\end{equation}
where the new thermodynamic potential ${\cal F}[\rho]$ is the
intrinsic free-energy which is a functional of $\rho ({\bf r})$. The
total Helmholtz free energy $F$ is recovered by adding the
contribution from the external field:

\begin{equation}
F = {\cal F} [\rho] + \int_V d {\bf r} \rho ({\bf r}) \phi ({\bf r})
\label{helmholtz_tot}
\end{equation}

Functional differentiation of (\ref{legendre_trasf}) leads to:

\begin{equation}
\frac{\delta {\cal F}[\rho]}{\delta \rho ({\bf r})} = \Psi ({\bf r})
\label{helmholtz_diff}
\end{equation}
while a further differentiation yields, according to
(\ref{static_susc_2}):

\begin{equation}
\frac{\delta^2 {\cal F}}{\delta \rho({\bf r})\delta \rho({\bf
r}^\prime)} = \frac{\delta \Psi({\bf r}^\prime)}{\delta\rho({\bf r})}
= \chi^{-1} ({\bf r},{\bf r}^\prime)
\label{inverse_susc}
\end{equation}
where $\chi^{-1}$ denotes the functional inverse of the susceptibility
$\chi({\bf r},{\bf r}^\prime)$, i.e.:

\begin{equation}
\int d {\bf r}^{\prime\prime} \chi^{-1} ({\bf r},{\bf
r}^{\prime\prime}) \chi ({\bf r}^{\prime\prime},{\bf r}^{\prime}) =
\delta ({\bf r}-{\bf r}^\prime).
\label{inverse_susc_def}
\end{equation}

It often proves convenient to split the functional into an ``ideal''
part, corresponding to a system of non-interacting particles (with the
same density profile and temperature), and an ``excess'' part stemming
from the interaction between particles:

\begin{equation}
{\cal F} [\rho] = {\cal F}_{id} [\rho] +{\cal F}_{ex} [\rho].
\label{split_helm}
\end{equation}

The inverse response function defined by Eq. (\ref{inverse_susc})
splits accordingly into two terms:

\begin{equation}
\chi^{-1} ({\bf r},{\bf r}^\prime) = \chi^{-1}_{id} ({\bf r},{\bf
r}^\prime) - k_B T c^{(2)} ({\bf r},{\bf r}^\prime)
\label{split_chi}
\end{equation}
where the second term, resulting from the second functional derivative
of ${\cal F}_{ex}$, is referred to as the two-particle direct
correlation function \cite{hansen86,chihara78}.

\par

The functional ${\cal F}[\rho]$, which reduces to the ground-state energy
of the inhomogeneous system at zero temperature, obeys two fundamental
theorems, which were first stated and proved by Hohenberg and Kohn
\cite{hohenberg64} for the ground state, and generalized to finite
temperature by Mermin \cite{mermin65}. Both theorems follow easily
from a lemma stating that there is a one-to-one correspondence between
the fields $\rho({\bf r})$ and $\phi({\bf r})$.

\begin{theorem}[Hohenberg-Kohn]
For any given system, characterized by the unperturbed Hamiltonian
$H_N = K_N + V_N$, the intrinsic free-energy ${\cal F}$ is a {\it unique}
functional of the equilibrium density $\rho({\bf r})$. In other words, the
functional form of ${\cal F}[\rho]$ is independent of the external
potential $\phi({\bf r})$.
\label{hk_theorem}
\end{theorem}

\begin{theorem}[Variational Principle]
Let $\rho^*({\bf r})$ be a properly normalized trial density profile. The
functional
\begin{equation}
\Omega_\phi [\rho^*] = {\cal F} [\rho^*] - \int d {\bf r} \rho^*({\bf
r}) \Psi({\bf r})
\end{equation}
reaches its minimum when $\rho^*({\bf r})$ coincides with the
equilibrium density $\rho({\bf r})$, i.e.
\begin{equation}
\left. \frac{\delta\Omega_\phi [\rho^*]}{\delta
\rho^*({\bf r})}\right|_{\rho^*=\rho} = 0
\end{equation}
and, moreover, $\Omega_\phi[\rho]$ is equal to the equilibrium grand
potential $\Omega$.
\label{var_theorem}
\end{theorem}

In terms of the functional ${\cal F}[\rho]$, the variational principle reads:

\begin{equation}
\left. \frac{\delta {\cal F} [\rho^*]}{\delta
\rho^*({\bf r})}\right|_{\rho^*=\rho} = \Psi({\bf r}) = \mu - \phi({\bf r}).
\label{variational}
\end{equation}

These two theorems are instrumental in all practical applications of
DFT.

\par

Henceforth, up to section VI, restriction will be made to classical
DFT ($\hbar=0$ limit). For classical systems, the ideal part of the
free-energy functional takes on the particularly simple form
appropriate for a classical ideal gas, namely:

\begin{equation}
{\cal F}_{id} [\rho] = k_B T \int_V d {\bf r}\rho({\bf r}) \{
\log[\Lambda^3\rho({\bf r})]-1\}
\label{idgas_classical}
\end{equation}
where $\Lambda$ is some irrelevant length scale, generally identified
with the De Broglie thermal length.

\par

No general expression exists for the non-trivial ``excess'' part
${\cal F}_{ex}$, for which {\it approximate} functionals have to be
constructed, which are tailored for given physical conditions. In the
following, four of the most commonly used approximations are briefly
described. Each is suited for fairly well-defined physical situations,
and has an analog in quantum-mechanical DFT.

\subsection{Linear response}

Linear response theory applies when the external perturbation is of
small amplitude compared to the interactions between particles. If the
unperturbed system ($\phi=0$) is homogeneous, the density profile
induced by a weak perturbation $\phi({\bf r})$ is given by
Eq. (\ref{taylor_linear}),

\begin{equation}
\Delta\rho({\bf r}) = \rho({\bf r})-n_0=-\int d {\bf r}^\prime
\chi_0({\bf r}-{\bf r}^\prime) \phi({\bf r}^\prime)
\label{linear_resp}
\end{equation}
where the index 0 refers to the homogeneous ``reference'' system. For
a spatially periodic perturbation $\phi({\bf r}) = \phi_{\bf k}
e^{i{\bf kr}}$ the Fourier components of the density profile are
easily expressed with the help of the convolution theorem. Using the
result (\ref{chi_class}), valid for classical systems, as adapted to a
homogeneous reference system, the final linear response expression
reads\cite{yvon58}:

\begin{equation}
\Delta \rho_{\bf k} = - \hat{\chi}({\bf k})\phi_{\bf k} = -\beta n_0 [n_0
\hat{h}^{(2)}({\bf k}) +1]\phi_{\bf k} = -\beta n_0 S_0({\bf k})\phi_{\bf k}
\label{class_deltarho}
\end{equation}
where $\hat{f}$ denotes the Fourier transform of $f$, and $S({\bf k})$
is the static structure factor of the homogeneous reference system
\cite{hansen86}.

\subsection{Gradient expansion}

If the inhomogeneity is slowly varying in space, on a characteristic
scale $\xi$ much larger than the mean spacing between particles, the
free-energy density $f[\rho]$ may be formally expanded in powers of
$1/\xi$. For an isotropic fluid, rotational invariance then leads to
the following gradient expansion\cite{evans79}:

\begin{equation}
{\cal F} [\rho] = \int_V d {\bf r} f[\rho({\bf r})] = \int_V d {\bf r} \{
f_0(\rho({\bf r}))+f_2(\rho({\bf r}))|\nabla\rho({\bf r})|^2\} + O(\xi^{-4}).
\label{grad_exp}
\end{equation}

Keeping only the zeroth order term amounts to making the {\it local
density approximation} (LDA); $f_0(\rho({\bf r}))$ is the free-energy
per unit volume of a homogeneous fluid of density $n=\rho({\bf r})$, a
function (as opposed to a functional) of $\rho({\bf r})$. Note that
the ideal contribution (\ref{idgas_classical}) to ${\cal F}$ is
precisely of the LDA form, so that an approximation is only made for
${\cal F}_{ex}$.

\par

The quadratic gradient correction to the LDA involves the function
$f_2(\rho({\bf r}))$ which is expressible in terms of the direct
correlation function [cf. Eq. (\ref{split_chi})] of the homogeneous
reference system\cite{evans79}:

\begin{equation}
f_2(\rho({\bf r})) = \frac{k_BT}{12} \int d {\bf r}^\prime {\bf r}^{\prime 2}
c_0^{(2)}(\rho({\bf r});|{\bf r}-{\bf r}^\prime|).
\label{grad_corr_classical}
\end{equation}

The gradient expansion (\ref{grad_exp}), limited to second order, is
the basis of van der Waals' theory of the liquid-vapour interface
\cite{vanderwaals93,rowlinson82} and of Cahn's theory of
wetting\cite{cahn77}. Higher-order corrections involve higher moments
of $c_0^{(2)}$ (as well as more-than-two-body direct correlation
functions). Since $c_0^{(2)}$ is known to behave asymptotically (as
$r\to \infty$) like the pair potential\cite{hansen86}, the integrals defining
higher-order coefficients $f_{2n}$ diverge for long-range
interactions, so that the expansion is ill-defined, and should wisely
be truncated at second order. In fact, for coulombic interactions,
$c_0^{(2)} (r) \sim 1/r$ at large $r$, so that even $f_2$ does not
exist. An entirely different approximation scheme is required to
handle ionic systems.

\subsection{Mean-field approximation}

For systems with long-range interactions, like ionic systems, {\it
non-local} effects may be expected to become important, which cannot
be properly coped with by a gradient expansion. In the spirit of
thermodynamic perturbation theory\cite{hansen86}, it proves convenient
to split the pair interaction potential between particles into a
short-range part $v_0$, usually repulsive, and a long-range
contribution $w$, the Coulomb potential in the case of ionic systems:

\begin{equation}
v({\bf r},{\bf r}^\prime) = v_0({\bf r},{\bf r}^\prime) + w({\bf
r},{\bf r}^\prime).
\label{split_pot}
\end{equation}

If the long-range part $w$ is treated within the van der Waals
mean-field approximation, the free-energy functional may be cast in
the form\cite{evans79,hansen86}:

\begin{equation}
{\cal F} [\rho] = {\cal F}_0[\rho] + \frac{1}{2} \int d {\bf r} \int
d {\bf r}^\prime \rho({\bf r}) w({\bf r},{\bf r}^\prime)\rho({\bf
r}^\prime)
\label{meanfield_classical}
\end{equation}
where ${\cal F}_0 [\rho]$ is the intrinsic free-energy functional of
the reference system, i.e. a system of particles interacting via $v_0
({\bf r},{\bf r}^\prime)$ only. Since the latter is short-ranged, the
gradient expansion may be appropriate to estimate ${\cal
F}_0[\rho]$. If the LDA is adopted for the latter, the variational
principle (\ref{variational}) leads to the following Euler-Lagrange
equation for $\rho({\bf r})$:

\begin{equation}
f^\prime_0(\rho({\bf r})) = \mu - \phi({\bf r}) - \int d {\bf r} w
({\bf r},{\bf r}^\prime) \rho({\bf r}^\prime)
\label{euler_meanfield}
\end{equation}
where $f_0$ is the free-energy density of the homogeneous reference
system, and the prime denotes a derivative with respect to the
argument. Eq. (\ref{euler_meanfield}) shows that the long-range
interactions act like a self-consistent external (or Hartree)
potential; it provides the basis for the study of electric
double-layers (cf. Section IV).

\subsection{Functional Taylor expansion}

If, contrary to the situation in section 2B, the inhomogeneity varies
rapidly in space, while remaining of small amplitude, it is natural to
seek an expansion of ${\cal F}_{ex}[\rho]$ around its value for a
homogeneous reference system, of density $n$, in powers of the local
``order parameter''

\begin{equation}
\Delta \rho({\bf r}) = \rho({\bf r}) - n.
\label{order_param}
\end{equation}

Using the results (\ref{inverse_susc}) and (\ref{split_chi}), one
easily arrives at:

\begin{eqnarray}
{\cal F}_{ex} [\rho] & = & {\cal F}_{ex} (n) + \mu_{ex} (n) \int_V d
{\bf r} \Delta \rho ({\bf r}) \nonumber \\ & - & \frac{k_B T }{2} \int
d {\bf r} \int d {\bf r}^\prime \Delta \rho({\bf r}) c_0^{(2)}
(n_0,{\bf r}-{\bf r}^\prime) \Delta \rho({\bf r}^\prime) + O((\Delta
\rho)^3).
\label{func_taylor_classical}
\end{eqnarray}

The linear term vanishes whenever the mean density of the
inhomogeneous system equals that of the reference system. The
expansion (\ref{func_taylor_classical}) holds even for a vanishing
external potential $\phi({\bf r})$, i.e. in particular for systems
with spontaneously broken translational invariance. It provides the
basis of DFT of first-order phase transitions like
freezing\cite{baus90}.

\section{A gentle application of DFT: sedimentation equilibrium
and osmotic equation-of-state}

As a simple illustration of how DFT can be efficiently combined with
computer simulation, consider a suspension of mesoscopic colloidal
particles under the action of gravity, which supplies the external
potential $\phi(z) = mgz$ where $m$ is the buoyant mass of the
particles, $g$ is the gravitational acceleration and $z$ the vertical
coordinate. Gravity breaks translational invariance in that direction
and induces a density profile $\rho(z)$. The overall colloid
concentration is characterized by the average number of particles per
unit area at the bottom:

\begin{equation}
n_s = \int_0^h d z \rho(z)
\label{av_bott}
\end{equation}
where $h$ is the height of the sample. The details of the z-dependence
of the density profile depend on the resolution of the observational
apparatus.

\par

At {\it high resolution} (finer than the particle diameter $\sigma$),
$\rho(z)$ is found to exhibit strong oscillatory behaviour near the
bottom wall (say $z=0$), characteristic of {\it layering}, as shown by
MC simulations\cite{biben93}. At sufficiently high concentrations
$n_s$, the lowest layers exhibit 2D crystalline order (triangular
lattice in the x-y plane\cite{biben94}). Sophisticated versions of DFT
(the so-called ``weighted density approximation''\cite{curtin85}) are
required to reproduce such refinements\cite{biben94}.

\par

At {\it much lower resolution} (i.e. on a scale significantly larger
than $\sigma$), as is achieved in depolarized light scattering
experiments\cite{piazza93}, one is looking at a coarse-grained density
profile:

\begin{equation}
\overline{\rho} (z) = \int_0^\infty d z^\prime \rho(z^\prime)
w(z-z^\prime)
\label{profile_highres}
\end{equation}
where $w(z)$ is a suitably normalized resolution function,
characteristic of the experimental set-up, or of the bin-width of
computer-generated histograms. Sufficient coarse-graining (resolution
of typically 2 -- 3 $\sigma$) smoothes out the low-z oscillations in
$\rho(z)$ and leads to a monotonically decreasing
profile\cite{biben93,piazza93}. Such a smoothly varying
$\overline{\rho}(z)$ may be calculated directly from DFT at the LDA
level. The Euler-Lagrange equation associated with the variational
principle (\ref{variational}) is particularly simple within the LDA,
namely (dropping the index 0):

\begin{equation}
f^\prime(\overline{\rho}(z)) \equiv
\frac{\Pi(\overline{\rho}(z))}{\overline{\rho}(z)} +
\frac{f(\overline{\rho}(z))}{\overline{\rho}(z)} = \mu -mgz
\label{euler_sedim}
\end{equation}
where $f(\overline{\rho})$ is the free-energy density of a homogeneous
suspension of concentration $\overline{\rho}$, and
$\Pi(\overline{\rho})$ is the corresponding osmotic pressure. Taking
the gradient of both sides of Eq. (\ref{euler_sedim}) leads to the
macroscopic osmotic equilibrium condition:

\begin{equation}
\frac{d \Pi(z)}{dz} = -mg\overline{\rho}(z).
\label{osmot_eq}
\end{equation}

Along an isotherm, $\Pi$ varies with altitude $z$ via its
$\overline{\rho}$-dependence (this is the local thermodynamic
equilibrium assumption implicit in LDA). Hence Eq. (\ref{osmot_eq})
can be re-cast into the following non-linear differential equation for
$\overline{\rho}(z)$:

\begin{equation}
\frac{d \overline{\rho}(z)}{dz} = - \chi_T (\overline{\rho}(z))
\frac{\overline{\rho}(z)}{\xi}
\label{nonlin_eq}
\end{equation}
where $\xi= k_B T /mg$ is the gravitational length and
$\chi_T=(\partial \beta \Pi/\partial \rho)^{-1}$ denotes the
dimensionless osmotic compressibility. Since $\chi_T>0$, the LDA can
only lead to a monotonously decreasing $\overline{\rho}(z)$. For an
ideal suspension of non-interacting particles, $\chi_T=1$, and the
density profile decays exponentially. At finite concentrations,
excluded-volume effects between the colloidal particles lead to a much
more slowly decaying $\overline{\rho}(z)$ at low
altitudes\cite{biben93,piazza93,vrij80}.

\par

An interesting application of sedimentation equilibrium is the {\it
inverse problem}, which amounts to extracting the osmotic
equation-of-state $\Pi(\rho)$ from experimental
\cite{piazza93,hachisu82} or computer-generated \cite{biben93} density
profiles of concentrated suspensions.

\par

Integration of the osmotic equilibrium condition (\ref{osmot_eq}) from
$z$ to an altitude $h$ sufficiently high for the osmotic pressure to
vanish yields:

\begin{equation}
\beta\Pi(z) = \frac{1}{\xi} \int_z^h dz^\prime \overline{\rho}(z^\prime).
\label{inverse_sedim}
\end{equation}

At each altitude, $\overline{\rho}(z)$ and $\Pi(z)$ are thus known
simultaneously, and hence, upon elimination of $z$, the whole osmotic
equation-of-state $\Pi(\rho)$ is determined, for $\rho \leq \rho_0 =
\overline{\rho}(z=0)$, by a single simulation or
measurement\cite{piazza93}. The resulting $\Pi(\rho)$ should be
independent of the intensity of the gravitational field, as
characterized by the dimensionless ratio $\alpha=\sigma/\xi$; this
independence provides a test for the validity of the procedure based
on the LDA. The method has proved its efficiency by reproducing quite
accurately the well-documented equation-of-state of a hard sphere
liquid, by inversion of density profiles determined by MC
simulations\cite{biben93} and by depolarized light scattering
measurements on a sterically stabilized colloidal
suspension\cite{piazza93}.

\par

In a recent application, a system of infinitely thin disks, carrying a
quadrupole moment $\stackrel{\leftrightarrow}{Q}$, was used as a model
for a clay colloid suspension\cite{dijkstra95}. The quadrupole
accounts for the electrostatic interactions between the double-layers
formed by the charged clay platelets and their associated counterion
charge distribution. If $\sigma$ denotes the disc diameter, the
relevant dimensionless state variables are $n^*=N\sigma^3/V$ and
$\beta^* = Q^2/(\varepsilon k_B T \sigma^5)$, where $Q$ is is the
magnitude of the quadrupole moment and $\varepsilon$ the dielectric
constant of the suspending fluid (water). The density profiles
$\overline{\rho}(z)$ computed by MC simulations for several values of
$\beta^*$ show a marked change in behavior around $\beta^*\simeq0.1$
(cf. Fig. \ref{fig_clay_a}); inversion of these density profiles leads
to the osmotic equations-of-state shown in Fig. \ref{fig_clay_b}; the
dramatic change at $\beta^* = 0.1$ points to a first-order phase
transition between a dilute sol-phase and a concentrated
gel-phase\cite{dijkstra95}. The nearly-incompressible gel results from
a competition between entropy, which favors stacking at high
concentrations\cite{eppenga84} (the loss of orientational entropy
being outweighted by the gain in translational entropy), and the
quadrupolar interactions favoring ``edge-to-face'' pair
configurations, which leads to a characteristic ``house-of-cards''
structure.

\section{A density functional approach to polyelectrolytes}

Polyelectrolytes are made up of multiply charged, mesoscopic {\it
polyions} and microscopic {\it counterions} of opposite sign. If salt
is added, microscopic coions (carrying charges of the same sign as the
polyions) are also present, but for the sake of simplicity this
physically more realistic situation will not be considered here.
These charged particles are in solution in a solvent (generally
water), which will be assumed to act as an inert continuum of
dielectric constant $\varepsilon$. Taking the counterions to be
monovalent, and the polyions to carry $Z$ elementary charges,
electrical neutrality requires that $N_1Z = N_2.$ Polyions may be
linear objects (flexible charged polymer chains), rod-like (a good
example is the tobacco mosaic virus), lamellar (charged clay platelets
or membranes), or globular (spherical colloidal particles). The case
of spherical polyions will be considered first.

\par

Apart from overlap repulsion at short distances (which will be assumed
to be of the hard sphere type in the following), the {\it bare}
interactions between ions are purely coulombic; thus, for instance,
returning to the notations of section I:

\begin{eqnarray}
V_{11} (\{ {\bf R}_I\}) &=& \sum_I \sum_{J>I} v_{11} (|{\bf R}_I-{\bf
R}_J|) \label{interaction_colloids} \\
v_{11} ({\bf R}) &=& \infty ~~~~~~~~~~R < \sigma \nonumber \\
 \phantom {v_{11} ({\bf R})}
&=& \frac{Z^2e^2}{\varepsilon R} ~~~~~~~R > \sigma. \nonumber
\end{eqnarray}

Counterions may be regarded as point charges, attracted by the
oppositely charged polyions, which they cannot penetrate, so that an
{\it electric double layer} builds up at the surface of of the
polyions; this double layer is characterized by a counterion density
profile $\rho({\bf r})$ which takes on large values whenever $|{\bf r} - {\bf
R}_I|$ is close to $\sigma$.

\par

From the analysis of section I it follows that, after integration over
the counterion coordinates, the total {\it effective} potential energy
of interaction of the macroions is the sum of a direct term, $V_{11}$,
which is pairwise additive, and an indirect (or induced) term, a
complicated function of the polyion configuration, which is
intrinsically many-body in nature:

\begin{equation}
V^{eff}_{11} (\{{\bf R}_I\}) = V_{11} (\{{\bf R}_I\}) + F_2
[\rho({\bf r}),\{{\bf R}_I\}]
\label{colloid_interaction}
\end{equation}
where $F_2$ is the free-energy of the inhomogeneous fluid of
counterions, which is a functional of $\rho({\bf r})$ and depends
parametrically on the polyion coordinates. It is convenient to split
the non-trivial excess part into a mean-field Coulomb (also called
Hartree) term (cf. section IIc) and a correction due to
correlations. Dropping the index 2:

\begin{equation}
F [\rho({\bf r}),\{{\bf R}_I\}] = {\cal F}_{id} + {\cal F}_{coul} +
{\cal F}_{corr} + {\cal F}_{ext}
\label{split_excess}
\end{equation}
where ${\cal F}_{id}$ is given by Eq. (\ref{idgas_classical}),

\begin{eqnarray}
{\cal F}_{coul} &=& \frac{e^2}{2\varepsilon} \int d {\bf r} \int d {\bf
r}^\prime \frac{\rho({\bf r})\rho({\bf r}^\prime)}{|{\bf r}-{\bf
r}^\prime|} \\
{\cal F}_{ext} &=& \int d {\bf r} \rho({\bf r})
\phi({\bf r},\{{\bf R}_I\}) = \sum_{I=1}^{N_1} \int d {\bf r}
\rho({\bf r}) v_{12} (|{\bf r}-{\bf R}_I|).
\label{fcoul_fext}
\end{eqnarray}

The (normally small) correlation part can be approximated by the LDA
ansatz:

\begin{equation}
{\cal F}_{corr} = \int_V d {\bf r} f^{ex}_{OCP} (\rho({\bf r})) 
\label{corr_lda}
\end{equation}
where $f^{ex}_{OCP}$ is the excess free-energy density of a
homogeneous system of point charges in a uniform background of
opposite charge. Such a system is referred to as a ``one-component
plasma'' (OCP), where $f^{ex}_{OCP}$ is well documented\cite{baus80}.

\par

For any macroion configuration $\{{\bf R}_I\}$, $\rho({\bf r})$ is
determined by the variational principle (\ref{variational}). The
resulting minimal free-energy is then substituted into
Eq. (\ref{colloid_interaction}) to determine the effective many-body
potential between polyions.

\par

An analytic solution of the variational problem can be obtained
\cite{lowen93}, provided two approximations are made: a) ${\cal
F}_{corr}$ is neglected [which is akin to Poisson-Boltzmann (PB)
theory in electrochemistry]; b) ${\cal F}_{id}$ is expanded to second
order in $\rho({\bf r})-n$, where $n=N_2/V$ is the mean counterion
density. With these approximations the free-energy reduces to a
quadratic functional of $\rho({\bf r})$ and the resulting
$V_{11}^{eff}$ turns out to be pairwise additive; the effective
polyion pair potential is of the ``DLVO'' form\cite{vervey48}:

\begin{equation}
v_{11}^{eff} ({\bf r}) = \frac{Z^2_{eff}}{\varepsilon r} e^{-k_{D}r};
~~~~~~~ r > \sigma
\label{dlvo}
\end{equation}
where $Z_{eff} = Z e^{k_D\sigma/2}/(1+k_D\sigma/2)$, and $k_D$ is the
inverse of the Debye screening length associated with the counterions.

\par

If the full functional (\ref{split_excess} -- \ref{corr_lda}) is used
instead, minimization must be carried out numerically. In principle
this could be achieved, for each polyion configuration, by steepest
descent or conjugate gradient methods. The resulting effective
interaction between polyions could then be used in a regular
Metropolis MC algorithm to obtain macroion configuration space
averages. This procedure will be used in section V for linear
polyelectrolytes. Alternatively, the effective forces between polyions
could be calculated from the classical version of the Hellman-Feynman
theorem (\ref{hell_feynm}), which results in:

\begin{equation}
{\bf F}_I = - \nabla_I V_{11} (\{ {\bf R}_J \}) -
\int_V d {\bf r} \rho({\bf r}) \nabla_I v_{12}({\bf r}-{\bf R}_I).
\label{hell_feynm_classical}
\end{equation}

These forces could then be used to evolve polyion positions in an MD
or Brownian dynamics simulation. More efficiently, along the lines
pioneered by Car and Parrinello for ion-electron systems
\cite{carpar85}, successive minimization and polyion updating steps
are replaced by a single dynamical evolution, which involves physical
motions of the polyions and fictitious dynamics of the counterion
density, arranged such that the latter readjusts adiabatically to the
much slower motion of the polyions. In practice, the time-dependent
density profile $\rho({\bf r},t)$ is parameterized by a plane-wave
expansion

\begin{equation}
\rho({\bf r},t) = \sum_{{\bf k}} \rho_{\bf k}(t) e^{i{\bf k}{\bf r}}
\label{pwexp}
\end{equation}
where the Fourier components $\rho_{{\bf k}}(t)$ are considered to be
dynamical variables. A fictitious mass $m_{{\bf k}}$ is associated
with the corresponding ``velocities'' $\dot\rho_{{\bf k}}(t)$ and the
following generalized Lagrangian is introduced:

\begin{equation}
{\cal L} = \frac{1}{2} M \sum_{I=1}^{N_1} \dot {\bf R}_I^2 +
\frac{1}{2} \sum_{{\bf k}} m_{{\bf k}} \dot\rho_{{\bf k}}^2 -
V_{11}(\{{\bf R}_I\})-F[\rho,\{{\bf R}_I\}]
\label{carparlag}
\end{equation}

The resulting equations of motion couple the physical time evolution
of the $\{{\bf R}_I\}$, and the fictitious dynamics of the
$\rho_{{\bf k}}$:

\begin{eqnarray}
 M\ddot {\bf R}_I &=& - \nabla_I V_{11} (\{{\bf R}_I\}) - \nabla_I
F[\rho,\{{\bf R}_I\}]\\
 m_{{\bf k}}\ddot {\rho}_{{\bf k}} &=&-\frac{\delta F [\rho,\{{\bf
R}_I\}]}{\delta \rho_{{\bf k}}} + \alpha \int_V d{\bf r}
[\rho({\bf r})-n]
\label{eom_carpar}
\end{eqnarray}
where the last term implements the constraint that the number of
counterions is conserved: $\rho_{\bf k=0}=n$. To achieve adiabaticity,
the ``masses'' $m_{{\bf k}}$ must be chosen sufficiently small
(compared to $M$, in some suitable units) to ensure a clear-cut
separation of time scales and thus to prevent equipartition between the
degrees of freedom associated with the polyions and the counterion
density. With these precautions, the counterion density remains close
to the optimum density (which minimizes the free-energy functional $F$
for each polyion configuration), all along the dynamical evolution of
the system.

\par

A technical difficulty arises from the discontinuous variation of the
density profile $\rho({\bf r})$ near the surfaces of the the polyions,
which are impenetrable to the counterions. This would make it
necessary to include a large number, in principle infinite, of Fourier
components in the representation (\ref{pwexp}). Fortunately, this
difficulty may be overcome via the introduction of a continuous
polyion-counterion pseudopotential, reminiscent of the
pseudopotentials used in quantum ion-electron systems \cite{lowen93}.

\par

Results obtained by such ``ab-initio'' simulations for various
properties of concentrated suspensions of moderately charged colloidal
particles ($Z \geq 100$) differ significantly from the data obtained
from standard MD or Brownian Dynamics simulations based on the DLVO
effective pair potential, Eq. (\ref{dlvo}). In view of the high
computational cost of the ``ab-initio'' simulations, it has been
suggested to fit the many-body ``ab-initio'' forces by forces derived
from the pairwise additive DLVO potential for a representative set of
polyion configurations\cite{kramposthuber93}. In this scheme the
effective valence $Z_{eff}$ of the polyions and the inverse screening
length $k_D$ are considered as free parameters to be determined by a
least-squares fitting procedure.

\par

The ``ab-initio'' method for polyelectrolytes with spherical polyions
has been successfully extended to the case of rod-like
polyions\cite{lowen94}, and is at present being generalized to
treat lamellar clay-colloid suspensions\cite{dijkstra95}.

\section{``Ab-initio'' Monte Carlo sampling}

Linear polyelectrolytes, consisting of charged polymer chains and
counterions, pose methodological problems quite different from those
encountered in the simulation of globular polyelectrolytes. In order
to explore the phase space of a polymer chain efficiently, large
conformational changes must be mapped out in a simulation. Such
changes can only be generated by allowing ``non-physical''
displacements of the individual monomers, and hence MC sampling of
configuration space is generally more efficient for polymers than
``deterministic'' MD.

\par

Again, for the sake of simplicity, only the case of a single
polyelectrolyte chain, together with an adequate number of
counterions, will be considered hereafter. The following model will be
adopted: the chain consists of $\nu$ identical segments of length $l$,
carrying a charge $Ze$. A segment is generally made up of several
monomers. The center of a segment and its orientation will be denoted
by ${\bf R}_i$ and the unit vector ${\bf u}_i,\; (1\le i\le \nu)$,
respectively. The segments are assumed to be freely-jointed, but
intrinsic rigidity is introduced via a harmonic coupling between the
orientations of two successive segments. This rigidity is enhanced by
the Coulomb repulsion between segments, as epitomized by the delicate
concept of electrostatic persistence length\cite{odijk77}. The total
potential energy of a chain thus reads:

\begin{equation}
V_{11} = \frac{k}{2} \sum_{i=1}^{\nu-1} \theta_i^2 +
\sum_i\sum_{j>i}\frac{Z^2e^2}{\varepsilon|{\bf R}_i-{\bf R}_j|}
\label{pot_chain}
\end{equation}
where $\cos(\theta_i) = {\bf u}_i\cdot{\bf u}_{i+1}$. The dielectric
constant $\epsilon$ will be set to 1.

\par

The counterions are once more described in terms of the density
profile $\rho({\bf r})$, and a free-energy functional in the form of
Eqs. (\ref{split_excess} -- \ref{corr_lda}). However, in contrast to
the case of the colloidal particles considered in section IV, the size
$l$ of the individual segments is comparable to that of the
microscopic counterions, so that excluded volume effects are not
expected to play an important role; in particular $l$ is generally
significantly smaller than the Debye screening length $\lambda_D =
k_D^{-1}$, which determines the scale of the counterion charge
distribution around a segment. Under these conditions it is reasonable
to replace the bare Coulomb attraction between a segment and a
counterion by the following smooth ``pseudopotential'':

\begin{equation}
v_{12} ({\bf r}) = - \frac{Ze^2}{\varepsilon r} \mbox{erf} (r/r_c)
\label{pseudopot}
\end{equation}
with $r_c \simeq l$; this functional form will subsequently prove very
favorable for a cut-off of the plane-wave representation (\ref{pwexp})
at a low wavenumber. Segments are prevented from overlapping by the
bare Coulomb repulsion in Eq. (\ref{pot_chain}). An MC simulation of
the polyelectrolyte then proceeds via Metropolis sampling of the
Boltzmann distribution:

\begin{equation}
P \sim \exp \{ - \beta V_{11} (\{{\bf R}_i\}) - \beta F_2[\rho,\{{\bf
R}_i\}] \} = P_{harm} (\{{\bf u}_i\}) P_{coul} (\{{\bf R}_i\})
\label{metropolis}
\end{equation}
where the harmonic factor corresponds to the intrinsic rigidity term
in the potential energy (\ref{pot_chain}).

\par

Each MC trial move comprises a) a conformational change of the
polymer, and b) a subsequent minimization of the counterion
free-energy $F_2$, which amounts to readjusting $\rho({\bf r})$ to the
new polymer configuration.

\par

a) The ``reptation'' algorithm \cite{wall75} provides a simple and
reasonably efficient way of displacing the polymer segments. A trial
move consists in detaching one segment from one of the ends (chosen at
random and designated as ``tail'') of the chain and reattaching it to
the opposite end (the ``head''), with an orientation sampled from
$P_{harm}$. The polymer chain thus ``slithers'' like a snake. In
approaches based on a discrete description of counterions, care must
be taken when defining an attempted MC move, as any move leading to an
overlap between monomers and counterions will be rejected, leading to
a slow sampling of phase-space. On the contrary, in a DFT approach
like the present one this problem does not exist, since the
short-range repulsion is implicit in $v_{12}$, and therefore even
complex trial moves of the polymer, which may lead to very efficient
algorithms, can be easily used.

\par

b) The simplest algorithm to minimize $F_2[\rho]$ is steepest
descent. The density profile is parameterized as in Eq. (\ref{pwexp}),
a fictitious time variable $s$ is introduced, and the coupled
equations for ``overdamped dynamics'' of the form

\begin{equation}
m_{{\bf k}} \frac{d}{ds} \rho_{{\bf k}} = - \frac{\delta
F_2[\rho]}{\delta \rho_{{\bf k}}}
\label{steepest}
\end{equation}
are solved by a finite ``time-step'' algorithm. The $\rho_{{\bf k}}$
become stationary when the free-energy minimum is reached. The choice
of the coefficients $m_{{\bf k}}$ is crucial since they determine the
convergence rate of the procedure without affecting the location of
the free-energy minimum. To see why, consider what happens close to
that minimum. There the free-energy is a quadratic function of the
deviations $\delta \rho_{{\bf k}}$ of the Fourier coefficients from
their equilibrium values:

\begin{eqnarray}
F_2[\rho] & = & F_2^{(0)} + \left. \frac{1}{2} \sum_{{\bf k}}
\sum_{{\bf k}^\prime} \frac{\delta^2F_2}{\delta\rho_{{\bf
k}}\delta\rho_{{\bf k}^\prime}} \right|_{\rho_{eq}} \delta\rho_{{\bf
k}} \delta \rho_{{\bf k}^\prime} + O(\delta \rho^3_{{\bf k}})
\nonumber \\ & = & F_2^{(0)} + \frac{1}{2} \sum_{{\bf k}} \sum_{{\bf
k}^\prime} \hat{\chi}^{-1}({\bf k},{\bf k}^\prime) \; \delta\rho_{{\bf
k}} \delta\rho_{{\bf k}^\prime} + O(\delta \rho^3_{{\bf k}})
\label{quadratic}
\end{eqnarray}
according to the Fourier-space version of Eq. (\ref{inverse_susc}).
Assume that the counterion density in the equilibrium state is nearly
uniform; as a consequence of translational invariance,
$\hat{\chi}^{-1}$ is then diagonal, $\hat{\chi}^{-1} ({\bf k},{\bf
k}^\prime) \simeq \hat{\chi}^{-1} ({\bf k}) \delta_{{\bf k},{\bf
k}^\prime}$. The quadratic free-energy expansion is then substituted
into the steepest descent equations (\ref{steepest}):

\begin{equation}
m_{{\bf k}} \frac{d}{ds} \rho_{{\bf k}} = - \frac{\delta}{\delta
\rho_{{\bf k}}} \{ F_2^{(0)} + \frac{1}{2} \sum_{{\bf k}} \sum_{{\bf
k}^\prime} \hat{\chi}^{-1}({\bf k},{\bf k}^\prime) \delta\rho_{{\bf
k}} \delta\rho_{{\bf k}^\prime} \} = -\sum_{{\bf k}^\prime}
\hat{\chi}^{-1} ({\bf k},{\bf k}^\prime)\delta\rho_{{\bf k}^\prime}
\simeq - \hat{\chi}^{-1} ({\bf k}) \delta\rho_{{\bf k}}.
\label{pre_precondition}
\end{equation}

The resulting linear equation is trivially solved in the form 

\begin{equation}
\delta\rho_{{\bf k}} \sim e^{-\gamma({\bf k})s}
\label{solve_sd}
\end{equation}
where the ``friction'' coefficient is

\begin{equation}
\gamma({\bf k}) = \frac{\hat{\chi}^{-1}({\bf k})}{m_{{\bf k}}}.
\label{friction}
\end{equation}

The rate of convergence is clearly dictated by the smallest
$\gamma({\bf k})$, while the stability of the numerical integration is
governed by the largest $\gamma({\bf k})$. The optimal situation would be
to have all $\gamma({\bf k})$ equal. This may be achieved by choosing

\begin{equation}
m_{{\bf k}} \propto \hat{\chi}^{-1} ({\bf k}).
\label{precondition}
\end{equation}

In order to obtain a rough estimate of $\hat{\chi}^{-1}$, one can
start from the general expression for a diagonal $\chi$

\begin{equation}
\hat{\chi}^{-1} ({\bf k}) \delta_{{\bf k},{\bf k}^\prime} =
\left[\hat{\chi}^{-1}_{id} ({\bf k}) + \frac{4\pi e^2}{\varepsilon{\bf
k}^2} - \hat{c}^{(2)}({\bf k}) \right] \delta_{{\bf k},{\bf k}^\prime}
\label{chidiag}
\end{equation}
where the susceptibility has been split into its ideal, mean-field (or
Hartree) and correlation contributions. In classical Statistical
Mechanics $\hat{\chi}^{-1}_{id}$ reduces to the constant $k_B T /
n$. For low counterion concentration the correlation term
$\hat{c}^{(2)}({\bf k})$ may be neglected (PB theory), so that

\begin{equation}
\hat{\chi}^{-1} ({\bf k},{\bf k}^\prime) \simeq \left[ k_B T / n +
\frac{4\pi e^2}{\varepsilon{\bf k}^2} \right] \delta_{{\bf k},{\bf
k}^\prime}.
\label{chi_pb}
\end{equation}
and a sensible choice for the $m_{\bf k}$ becomes

\begin{equation}
m_{\bf k} \propto \frac{k_B T}{n} + \frac{4\pi e^2}{\varepsilon{\bf
k}^2}.
\label{precondition_1}
\end{equation}

This amounts to the ``preconditioning'' proposed by Pearson {\sl et
al.} \cite{pearson93} to achieve optimal convergence of the steepest
descent minimization, and optimal adiabaticity in the dynamical
runs. They noted that, since steepest descent is most efficient if the
parabolic potential well (\ref{quadratic}) is as nearly as possible
axially symmetric, the choice (\ref{precondition}) for the $m_{\bf k}$
optimizes the efficiency of the algorithm by transforming an initially
``long and narrow valley'' into a ``circular'' one.

\par

Preliminary MC runs have been made on a thirty-segment chain in the PB
approximation, using both the full free energy functional
(\ref{split_excess}), neglecting ${\cal F}_{corr}$ and the linearized
version alluded to in Section IV. The effective charge of the monomers
is $Z=1$, and the counterion density is such that the Debye length
$\lambda_D=1/k_D=6l$. The simulations were performed in a simple cubic
cell. This of course introduces a spurious periodicity and thus
unphysical interactions between the charges inside the box and their
periodic images. To minimize this effect, a cell edge $L=46l$ was
used; with this value of $L$, even in the worst possible case (a fully
stretched chain parallel to one of the box edges) the distance between
a monomer and any one of the periodic images is about $3 \lambda_D$,
which gives a spurious interaction energy of less than 0.01 $k_B
T$. The results obtained with the use a pseudopotential of the form
(\ref{pseudopot}), $r_c=l$, were indistinguishable from those obtained
using the bare coulombic potential.

\par

The results show that excellent convergence to a relative accuracy
$\Delta F_2/ k_B T\sim 10^{-3}$ is achieved in about 10 steps for each
chain configuration. It is worth noting that in simulations using a
linearized energy functional [which is equivalent to using an
interaction of the DLVO form, Eq. (\ref{dlvo})], the steepest-descent
procedure fully converges after a single step, showing the
preconditioning at its best. At each attempted MC move, the trial
density used to start the energy minimization is the converged density
at the previous polymer configuration. However this is probably not
the best startup, since a large charge imbalance is created by the
reptation algorithm when a ``tail-to-head'' move is attempted
\cite{smargi96}. The new head in fact is by necessity placed in a
region of nearly uniform density, while leaving behind a peak in the
counterion density. More suitable choices for the starting density are
possible and should make the simulation appreciably faster.

\par

The simulation runs were 10000 MC steps (i.e., attempted head-to-tail
moves) long. All the calculations were carried out on a desktop
workstation. The acceptance ratio was 65\% for the linear-screening
run and 76\% for the full PB simulation, rather higher than the
``optimal'' 50\% ratio. This suggests that the reptation trial moves
produce changes in the polymer configuration that are unnecessarily
small. More complex MC algorithms, involving the simultaneous
displacement of many segments, should lead to more efficient sampling
of phase space.

\par

The physical problems which can be addressed in these simulations are
Manning ``condensation'' of counterions\cite{manning69}, which takes
place when the ionic strength is high enough to cause a massive
accumulation of counterions on the polyelectrolyte, and the
electrostatic persistence length\cite{odijk77}. The latter has
attracted much attention recently since the theory predicts the
existence of two different regimes, one in which the persistence
length is proportional to $1/k_D$ and another where it is proportional
to $1/k_D^2$. The parameters chosen for the present simulation
correspond roughly to those appropriate to the transition regime, a
case where non-linear effects on screening begin to be important. This
can be seen by the computed persistence lengths in the two cases of
linear screening and of PB. Fig. \ref{fig_persis} shows the average
correlation function between the orientations of the monomers for the
zero-charge case (an exponential with persistence length of 1) and for
the two runs. The difference between charged and ideal cases is
considerable, showing the expected charge-induced stiffening of the
chain. Longer chain lengths are required to obtain significant results
and are presently being investigated\cite{smargi96}.

\par

Refinement of the physical model are straightforward and are currently
being implemented. In particular, the important case of non-dialized
solutions (i.e. systems with added salt), which have never been
addressed numerically, is simply dealt with by adding a second density
field which represents the coionic density. It is also easily possible
to go beyond PB theory by adding a correlation term to the energy
functional, for example within the LDA.

\section{The quantum regime: orbital-free density functionals}

The basic theorems and most of the approximation schemes introduced in
the previous sections apply to classical and quantum systems
alike. When the ``light'' particles are valence electrons while the
``heavy'' particles are ions, the mass ratio $m_2/m_1 < 10^{-3}$, so
that the BO approximation is, under most circumstances, quite
adequate.

\par

However, under room-temperature conditions, the electrons are almost
invariably degenerate. Indeed, if $T_F$ denotes the Fermi temperature
of the valence electron gas of average density $n$, it is customary to
define the degeneracy parameter $\theta=T/T_F=2(9/4\pi)^{2/3}r_s^2
T^*$ where $r_s=(3/4\pi n)^{1/3}$ and $T^*=k_B T$ expressed in atomic
(Hartree) units: $e=\hbar=m_{el}=1$, which will be adopted throughout
in this section. In metals, $\theta<10^{-2}$, so that the valence
electrons may be considered as fully degenerate, and the free energy
$F$ reduces for all practical purposes to the ground state energy $E$,
which, according to the Hohenberg-Kohn theorem, is a unique functional
of the local electron density $\rho({\bf r})$, which reaches its minimum
when $\rho$ coincides with the equilibrium density. One is thus left
with the problem of determining the ground state energy functional of
an inhomogeneous electron gas in its ground state, subjected to the
instantaneous ``external field'' of the (classical) ions.

\par

This functional can be decomposed, in analogy with
(\ref{split_excess}), into a kinetic, Hartree, exchange-correlation
and external part:

\begin{eqnarray}
E[\rho,\{{\bf R}_I\}] & = & K + E_{hartree} + E_{xc} + E_{ext} \\
E_{hartree} [\rho] & = & \int d {\bf r} \int d {\bf r}^\prime
\frac{\rho({\bf r})\rho({\bf r}^\prime)}{|{\bf r}-{\bf r}^\prime|} \\
E_{ext} [\rho] & = & \sum_{I=1}^{N_1} \int d {\bf r} \rho({\bf
r}) v_{12}({\bf r}-{\bf R}_I)
\label{split_quantum}
\end{eqnarray}
The ion-electron interaction potential $v_{12}$ is often modeled using
a pseudopotential, which eliminates core electrons from the
calculation.

\par

The main problem, which does not arise in classical DFT, is to write
down an accurate kinetic energy functional $K$. Neglecting spin
variables for simplicity, if $\psi({\bf r}_1,\cdots,{\bf r}_N)$
denotes the many-body ground state wavefunction, $K$ is given exactly
by:

\begin{eqnarray}
K &=& -\frac{1}{2} \sum_{I=1}^{N} \int d {\bf r}^N \psi^*({\bf
r}_1,\cdots,{\bf r}_N) \frac{\partial^2 \psi({\bf r}_1,\cdots,{\bf
r}_N)}{\partial {\bf r}^2_I} \nonumber \\
 &=& -\frac{N}{2} \int d^N {\bf r}
\psi^*({\bf r}_1,\cdots,{\bf r}_N) \frac{\partial^2 \psi({\bf
r}_1,\cdots,{\bf r}_N)}{\partial {\bf r}^2_1} \nonumber \\
 &=& -\frac{N}{2} \int d
{\bf r} \int d {\bf r}^\prime \delta ({\bf r}-{\bf r}^\prime)
\psi^*({\bf r},{\bf r}_2,\cdots,{\bf r}_N) \frac{\partial^2 \psi({\bf
r},{\bf r}_2,\cdots,{\bf r}_N)}{\partial {\bf r}^2} \nonumber \\
 &=&
-\frac{1}{2} \int d {\bf r} \left.  \frac{\partial^2 \rho^{(1)}({\bf
r},{\bf r}^\prime)}{\partial {\bf r}^2} \right|_{{\bf r}^\prime={\bf
r}}
\label{ke_exact}
\end{eqnarray}
where $\rho^{(1)}({\bf r},{\bf r}^\prime)=N \int d {\bf r}_2 \cdots d
{\bf r}_N \psi^*({\bf r},{\bf r}_2,\cdots,{\bf r}_N)\psi({\bf r},{\bf
r}_2,\cdots,{\bf r}_N)$ is the density-matrix in real space. The usual
density is the diagonal part of $\rho^{(1)}$:

\begin{equation}
\rho({\bf r}) = \int d^N {\bf r} \tilde{\rho} ({\bf r}) |\psi|^2 =
\rho^{(1)}({\bf r},{\bf r}).
\label{density_diag}
\end{equation}

The usual DFT representations of $K$ are of one of the following
types:

\begin{itemize}
\item[a)] Write $K$ as the kinetic energy of a system of
non-interacting electrons moving in a self-consistent, effective
potential. That is, $\rho({\bf r})=\sum_i|\phi_i({\bf r})|^2$ where
the $\phi$'s are single particle wavefunctions of electrons moving in
the mean-field generated by the others (Kohn-Sham
orbitals\cite{kohn65}). In this case $K$ assumes the familiar form
$K=\frac{1}{2}\sum_i \int d {\bf r} \phi^*_i \nabla^2\phi_i$. This
approach is extensively discussed in the lectures of R. Car and
W. Kohn included in this volume.
\item[b)] Use approximations of $K$ which do not use orbitals but
directly $\rho$ as the basic variable. This is closer in spirit to the
Hohenberg-Kohn approach but, on the other hand, it is not clear a
priori whether good approximations for $K$ can be found.
\end{itemize}

The simplest approach of type (b) is to use the LDA approximation for
$K[\rho]$: this leads to the Thomas-Fermi (TF) and Thomas-Fermi-Dirac
(TFD) theories. Gradient-corrected functionals\cite{parr89} can also
be used to improve upon them. Another straightforward road for
improvement is is to use the linear response theory, as in
Eq. (\ref{linear_resp}) (Hohenberg-Kohn
functional\cite{hohenberg64}). These approaches, together with the
more recent ``hybrid'' schemes, are briefly discussed in the
following.

\subsection{Thomas-Fermi Molecular Dynamics}

The kinetic energy per particle of a uniform electron gas of density
$n$ is

\begin{equation}
E_K^0 = c_k n^{2/3} \sim 1/r^2_s, ~~~~~~ c_k = \frac{3}{10}
(3\pi^2)^{2/3}.
\label{ke_uniform}
\end{equation}

If one assumes that the kinetic energy of an inhomogeneous electron
gas is a function of the local density only, the result is the
well-known TF kinetic energy:

\begin{equation}
K_0 = c_k \int d {\bf r} [\rho({\bf r})]^{5/3}.
\label{ke_tf}
\end{equation}

This can be inserted in the energy functional (\ref{split_quantum})
and $\rho$ can be used as variable in a Car-Parrinello scheme, in
complete analogy with the dynamical treatment of polyelectrolytes in
Section IV.

\par

The obvious way of proceeding is to write $\rho$ on a discrete grid in
real space and to use the values of $\rho$ at the grid points as the
additional degrees of freedom. This leads to the following Lagrangian:

\begin{equation}
{\cal L} = \frac{1}{2} M \sum_{I=1}^{N_1} \dot {\bf R}_I^2 +
\frac{1}{2} m \int d {\bf r} \dot\rho^2({\bf r}) - E[\rho,\{{\bf R}_I\}]
- \sum_J \sum_{I<J} \frac{Z^2}{|{\bf R}_I-{\bf R}_J|}
+ \mu\left\{ \int d {\bf r} \rho({\bf r}) - N_2\right\}.
\label{carparlag_elec}
\end{equation}
The last term ensures the conservation of the number of the electrons
via the Lagrange multiplier $\mu$, the chemical potential of the
electrons. As discussed in Sec. IV, the ``fake'' mass $m$ must be
small enough to ensure adiabatic separation between ionic and (fake)
electronic degrees of freedom. The coupled equations of motion derived
from Eq. (\ref{carparlag_elec}) are

\begin{eqnarray}
M\ddot {\bf R}_I &=& Z^2 \sum_{I=1}^{N_1} {\bf R}_{IJ}/R^3_{IJ} -
\frac{\partial E[\rho,\{{\bf R}_I\}]}{{\bf R}_I}\\
 m \ddot {\rho}({\bf
r}) &=&-\frac{\delta E [\rho,\{{\bf R}_I\}]}{\delta \rho} + \mu =
-\frac{5}{3} c_k \rho^{2/3}({\bf r})-\phi({\bf r}) - \int d {\bf
r}^\prime \frac{\rho({\bf r}^\prime)}{|{\bf r}-{\bf r}^\prime|} + \mu.
\label{eom_carpar_elec}
\end{eqnarray}

This real-space approach has been successfully applied to the study of
dense hydrogen plasmas \cite{clerouin92}. However it has the drawback
that no effective preconditioning procedure is known and as a
consequence the efficiency of the algorithm is sensibly reduced. Thus
a plane-wave expansion of $\rho$ [Eq. (\ref{pwexp})], using fast
Fourier transform techniques to switch back and forth from real to
reciprocal space, seems preferable.

\par

Other, more basic problems lie with the functional itself. TF is known
to have several unphysical features: for instance, it leads to
divergence of the density at the atomic nucleus and, perhaps even
worse, cannot give chemical binding. An effortless improvement is to
include the LDA exchange functional \cite{kohn65}:

\begin{equation}
E_{xc} = c_x \int d {\bf r} [\rho({\bf r})]^{4/3}
\label{lda_xc}
\end{equation}
which leads to the TFD theory. Unfortunately this only make matters
worse: the non-binding theorem still holds and the quantitative
(dis-)agreement with both experiments and other, more accurate
theoretical results worsens.

\subsection{Gradient correction to TF}

The leading gradient correction to $K_0$ turns out to have the same
form for both the cases of rapidly and slowly varying density:

\begin{equation}
K_2[\rho] = K_0 + \frac{c_2}{8} \int d {\bf r}
\frac{(\nabla\rho({\bf r}))^2}{\rho({\bf r})}
\label{gradex_el}
\end{equation}
where $c_2=c_g=1/9$ if $\rho$ is slowly varying (second-order gradient
expansion) and $c_2=c_w=1$ if $\rho$ is rapidly varying
(Thomas-Fermi-von Weizs\"acker functional). In both cases the curses
of the divergence at the atomic nucleus and of the absence of chemical
binding are nullified.

\par

Still, this functional has the problem of not being able to yield the
correct linear-response behaviour. It is instructive to examine this
point in detail. Let us consider a system of non-interacting electrons
subject to a weak external potential $\phi$. Using the above
gradient-corrected functional,

\begin{equation}
E [\rho] = c_k \int d {\bf r} \rho^{5/3}({\bf r}) + \frac{c_2}{8} \int
d {\bf r} \frac{|\nabla \rho({\bf r})|^2}{\rho({\bf r})} + \int d {\bf
r} \phi({\bf r}) \rho({\bf r})
\label{vwfunc}
\end{equation}

The minimum condition reads

\begin{equation}
\frac{\delta E}{\delta\rho({\bf r})} = \frac{5}{3} c_k \rho^{2/3}({\bf r}) -
\frac{c_2}{8} \frac{|\nabla\rho({\bf r})|^2}{\rho^2({\bf r})} - \frac{c_2}{4}
\frac{\nabla^2\rho({\bf r})}{\rho({\bf r})} + \phi({\bf r}) = \mu
\label{eul_vw}
\end{equation}
where the presence of $\mu$, the chemical potential, ensures the
conservation of the number of electrons. Writing now $\rho({\bf
r})=n+\delta\rho({\bf r})$, where $\delta\rho({\bf r})=O(\phi)$, and
expanding $\frac{\delta E}{\delta\rho}$ to first order in $\delta\rho$
and $\phi$, one obtains upon Fourier transformation :

\begin{equation}
\delta\rho_{\bf k} = \chi_2({\bf k})\phi_{\bf k}; ~~~~~~ \chi_2({\bf
k}) = - \frac{k_F}{\pi^2}\frac{1}{1+3c_2q^2}
\label{lrt_vw}
\end{equation}
where $q=k/2k_F$ and $k_F=(9\pi/4)^{1/3} / r_s$ is the Fermi
wavevector (in a.u.). On the other hand the exact linear response
function for non-interacting electrons is the Lindhard
function\cite{lindhard54} and it is easy to check that the gradient
corrected kinetic energy gives the correct answer in the limit of
large $q$ (if $c_2=c_w=1$, von Weizs\"acker), or in the limit of small
$q$ (if $c_2=c_g=1/9$, gradient expansion).

\subsection{The Hohenberg-Kohn functional}

The correct linear response behaviour can, on the other hand, be
directly built in the functional as suggested by Hohenberg and
Kohn\cite{hohenberg64}. In this case one may write

\begin{equation}
K_L[\rho] = 
K_0(n) + \left. \frac{1}{2} \sum_{{\bf k}} \sum_{{\bf
k}^\prime} \frac{\delta^2 K}{\delta\rho_{{\bf k}}\delta\rho_{{\bf
k}^\prime}} \right|_n \delta\rho_{{\bf k}} \delta \rho_{{\bf
k}^\prime} 
= K_0(n) + \frac{1}{2} \sum_{{\bf k}} \chi_{id}^{-1}({\bf
k}) \delta\rho_{\bf k} \delta\rho_{-\bf k}
\label{quadratic_hk}
\end{equation}
where again $\chi_{id}$ is the Lindhard function. For weak density
modulations this functional is superior to TF, since the latter
implies a constant $\chi_{id} = k^2_F/\pi^2$. The Hohenberg-Kohn
functional has been used to derive effective ion-ion
potentials\cite{hafner} for simple metals. In this approach, a weak
ion-electron pseudopotential is constructed and the corresponding
response of a uniform electron gas at the appropriate $n$ is
calculated in the linear response approximation (which is equivalent
to using Eq. (\ref{quadratic_hk}) plus possibly a linearized
exchange-correlation potential). The resulting energy $E_2$ is then
used to calculate the effective potential $V_{11}^{eff}$ via
Eq. (\ref{pot_eff}). If the assumption of weak pseudopotentials holds,
it is easy to show that $V_{11}^{eff}$ is pairwise additive, as in the
classical case of polyelectrolytes [cf. Eq. (\ref{dlvo})]. A
disadvantage of the Hohenberg-Kohn functional is that it can lead to
unphysical negative values of $\rho({\bf r})$, even in conjunction
with weak pseudopotentials like those for Na\cite{pearson93}.

\subsection{Hybrid functionals}

Given that the approaches described in the last two paragraphs run
into complementary difficulties, it is quite natural to ask whether it
is possible to combine the above functionals into a new one which is
free from these problems\cite{pearson93,wang92,perrot94}. In general,
one may write a family of functionals of the form:

\begin{equation}
K_\alpha = K_w + \frac{1}{2} \sum_{\bf k} \rho^{\alpha *}_{\bf k}
[\chi^{-1}_{id}-\chi^{-1}_2] \rho^\alpha_{\bf k}.
\label{alpha_func}
\end{equation}

By $K_w$ it is meant the functional (\ref{gradex_el}) with $c_2=c_w$,
and $\chi_2$ is the linear response function derived from it
[Eq. (\ref{lrt_vw})], while $\chi_{id}$ is the Lindhard function.
$\rho^\alpha_{\bf k}$ is the Fourier transform of $\rho^\alpha({\bf
r})$, and the exponent $\alpha$ is a free parameter which has been set
to 1 \cite{perrot94}, 1/2 \cite{pearson93}, and 5/6 \cite{wang92}.

\par

By construction, $K_\alpha$ is exact in the linear regime and is
superior to Hohenberg and Kohn's in that it includes non-linearity via
the TF-von Weizs\"acker term. This term gives, for instance, the
correct asymptotic behavior in the case of rapidly varying density
even for large density variations.

\par

These functionals give good results for alkali
metals\cite{pearson93,smargi94} and plasmas \cite{penman}. Extensions
and improvements are under way\cite{foley}.

\subsection{Kohn-Sham vs. orbital-free functionals}

To conclude this section, it seems appropriate to summarize the
comparative advantages and disadvantages of the Kohn-Sham and
orbital-free approaches.

\par

The Kohn-Sham functional certainly gives a more accurate
representation of the kinetic energy, in particular for semiconductors
and insulators. (One should be careful not to consider it as
``exact'', though: it is simply the best possible $K$ in a
single-particle, mean-field approach). Another advantage of an
orbital-based scheme is that the density is, by construction, positive
($\rho = \sum_i |\phi_i|^2$).

\par

On the other hand orbital-free schemes offer great practical advantages:

\begin{itemize}
\item The basic variable being $\rho$, the number of degrees of
freedom is $O(N_2)$ times smaller, with savings both in computer time
and memory.
\item There is no need for Brillouin zone sampling. The use of an
insufficient number of sampling points is a major source of error in
orbital-based calculations, occasionally making them {\sl less}
accurate than orbital-free ones \cite{smargi94}.
\item No orthonormality has to be imposed (the charge conservation
constraint is trivial). Orthonormalization causes the computer time
required by orbital-based schemes to scale as $N_2^3$; orbital-free
methods are instead intrinsically $O(N_2)$.
\item Very good preconditioning (along lines similar to those
illustrated in section V) is possible, giving excellent adiabaticity
in the MD runs even for metallic systems, where Kohn-Sham-based
Car-Parrinello methods often run into difficulties\cite{pastore91}.
\item Inclusion of finite-$T$ effects is more straightforward. All one
has to do is to modify the functional as follows\cite{perrot79}.
First, change the Thomas-Fermi-von Weizs\"acker part:

\begin{eqnarray}
K_w &=& \sqrt{2}/\pi^2 (k_B T)^{5/2} [-\frac{2}{3} I_{3/2}(\eta)+\eta
I_{1/2}(\eta)] \nonumber \\
& + & \frac{c_w}{8} \int d {\bf r}
\left[-\frac{3\sqrt{2}\pi^2}{8(k_B T)^{3/2}} \frac{d}{d\eta}
(1/I_{-1/2}(\eta))\right] (\nabla\rho({\bf r}))^2
\label{finite_dft_quantum}
\end{eqnarray}
where $I_\nu (\eta) = \int_0^\infty dx \frac{x^\nu}{e^{x-\eta}+1}$ and
$\eta$, the reduced chemical potential, is defined via $\rho=
\sqrt{2}/\pi^2 (k_B T)^{3/2} I_{1/2}(\eta)$. Secondly, replace the
linear response function in Eq. (\ref{alpha_func}) with its
finite-temperature counterpart:

\begin{equation}
\chi_{id} (k;T) = -\frac{k_F^2}{\pi^2 k} \int_0^\infty dy
\frac{y}{e^{(1/2 k_F^2 y^2-\mu)/k_B T}+1} \ln \left| \frac{k + 2 k_F
y}{k - 2 k_F y} \right| \nonumber
\end{equation}

\end{itemize}

\vskip 0.5truecm

On balance, orbital-free functionals appear to be better adapted to
the simulation of metals and plasmas, especially if dynamical
properties, requiring long trajectories in phase space, are to be
investigated.

\acknowledgements

The work on polyelectrolytes has been done in collaboration with
J.-L. Barrat, and that on the quantum orbital-free energy functionals
with P.~A. Madden. We thank M. Dijkstra for providing us with the
figures on clays. Additional funding has been provided by the EEC
network ``Monte Carlo and Molecular Dynamics Simulations of Classical
and Quantum Systems''.

\begin{thebibliography}{99}

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

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

\bibitem{vanderwaals93} J.~D. van der Waals,
Verhandel. Konink. Akad. Weten., Amsterdam (Sect.1), {\bf 1}, 56
(1893).

\bibitem{gouy10} G. Gouy, J. Physique {\bf 9}, 457 (1910).

\bibitem{percus62} J.~K. Percus, Phys. Rev. Lett. {\bf 8}, 462 (1962).

\bibitem{hohenberg64} P. Hohenberg and W. Kohn, Phys. Rev. {\bf 136},
B864 (1964).

\bibitem{kohn65} W. Kohn and L.~J. Sham, Phys. Rev. {\bf 140}, A1133
(1965).

\bibitem{evans79} R. Evans, Adv. Phys. {\bf 28}, 143 (1979); {\sl
ibid}, in {\it Fundamentals of of Inhomogeneous Fluids}, edited by
D. Henderson (Marcel Dekker, New York, 1992).

\bibitem{parr89} R.~G. Parr and W. Yang, {\sl Density Functional
Theory of Atoms and Molecules} (Oxford University Press, Oxford,
1989); R.~M. Dreizler and E.~K.~U. Gross, {\sl Density Functional
Theory} (Springer, Berlin, 1990); R.~O. Jones and O. Gunnarson,
Rev. Mod. Phys. {\bf 61}, 689 (1989).

\bibitem{forster75} See, for instance, D. Forster, {\sl Hydrodynamic
Fluctuations, Broken Symmetry and Correlation Functions} (Benjamin,
Reading (MA), 1975).

\bibitem{hansen86} See, for instance, J.-P. Hansen and
I.~J. MacDonald, {\sl Theory of simple liquids}, 2nd edition (Academic
Press, London, 1986).

\bibitem{chihara78} J. Chihara, Prog. Theor. Phys. {\bf 60}, 1640
(1978); J. Phys. -- Cond. Matt. {\bf 3}, 8715 (1991).

\bibitem{mermin65} N.~D. Mermin, Phys. Rev. {\bf 137}, A1441 (1965).

\bibitem{yvon58} J. Yvon, Suppl. Nuovo Cim. {\bf 9}, 144 (1958).

\bibitem{rowlinson82} J.~S. Rowlinson and B. Widom, {\sl Molecular
theory of capillarity} (Oxford University Press, Oxford, 1982).

\bibitem{cahn77} J.~W. Cahn, J. Chem. Phys. {\bf 66}, 3667 (1977).

\bibitem{baus90} See, for instance, M. Baus, J. Phys. Cond. Matt. {\bf 2},
2111 (1990).

\bibitem{biben93} T. Biben, J.-P. Hansen, and J.-L. Barrat,
J. Chem. Phys {\bf 98}, 7330 (1993).

\bibitem{biben94} T. Biben, R. Ohnesorge, and H. L\"owen,
Europhys. Lett. {\bf 28}, 665 (1994).

\bibitem{curtin85} W.~G. Curtin and N.~W. Ashcroft, Phys. Rev. A {\bf
32}, 2909 (1985).

\bibitem{piazza93} R. Piazza, T. Bellini, and V. Degiorgio,
Phys. Rev. Lett. {\bf 71}, 4267 (1993).

\bibitem{vrij80} A. Vrij, J. Chem. Phys. {\bf 72}, 3735 (1980).

\bibitem{hachisu82} S. Hachisu and K. Tokano, Adv. Colloid. Interface
Sci. {\bf 16}, 233 (1982).

\bibitem{dijkstra95} M. Dijkstra, J.-P. Hansen, and P.~A. Madden,
Phys. Rev. Lett. {\bf 75}, in press (1995).

\bibitem{eppenga84} R. Eppenga and D. Frenkel, Mol. Phys. {\bf 52},
1303 (1984).

\bibitem{baus80} For reviews, see M. Baus and J.~P. Hansen,
Phys. Reports {\bf 59}, 1 (1980), or S. Ichimaru, H. Iyetomi, and
S. Tanaka, {\sl ibid.} {\bf 149}, 91 (1987).

\bibitem{lowen93} H. L\"owen, J.-P. Hansen, and P.~A. Madden,
J. Chem. Phys. {\bf 98}, 3275 (1993).

\bibitem{vervey48} E.~J.~W. Vervey and J.~T.~G Overbeek, {\sl Theory
of the stability of lypophobic colloids} (Elsevier, Amsterdam, 1948).

\bibitem{kramposthuber93} H. L\"owen and G. Kramposthuber,
Europhys. Lett. {\bf 23}, 673 (1993).

\bibitem{lowen94} H. L\"owen, J. Chem. Phys. {\bf 100}, 6738 (1994).

\bibitem{odijk77} T. Odijk, J. Polym. Sci. {\bf 15}, 477 (1977);
J.-L. Barrat and J.-F. Joanny, Europhys. Lett. {\bf 24}, 333 91993).

\bibitem{wall75} F.~T. Wall and F. Mandel, J. Chem. Phys. {\bf 63},
4592 (1975).

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

\bibitem{smargi96} E. Smargiassi, J.-L. Barrat, and J.-P. Hansen, to
be published.

\bibitem{manning69} G.~S. Manning, J. Chem. Phys. {\bf 51}, 924
(1969).

\bibitem{clerouin92} J. Cl\'erouin, G. Z\'erah, D. B\'enisti, and
J.-P. Hansen, Europhys. Lett. {\bf 13}, 685 (1990); G. Z\'erah,
J. Cl\'erouin, and E.~L. Pollock, Phys. Rev. Lett. {\bf 69}, 446
(1992); J. Cl\'erouin, E.~L. Pollock, and G. Z\'erah, Phys. Rev. A
{\bf 46}, 5130 (1992).

\bibitem{lindhard54} J. Lindhard, Kgl. Danske Videnskab. Selskab.,
Mat.-Fys.  Medd., {\bf 28}, 8 (1954).

\bibitem{hafner} J. Hafner, {\sl From Hamiltonians to Phase Diagrams}
(Springer, Berlin, 1987).

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

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

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

\bibitem{penman} J.~I. Penman and G. Zerah, to be published.

\bibitem{foley} M. Foley and P.~A. Madden, to be published.

\bibitem{perrot79} F. Perrot, Phys. Rev. A {\bf 20}, 586 (1979).

\end{thebibliography}

\begin{figure}
\caption{Density profiles for the model clay system described in the
text, for several values of the dimensionless parameter $\beta^* =
Q^2/(\varepsilon k_B T \sigma^5)$.}
\label{fig_clay_a}
\end{figure}

\begin{figure}
\caption{Osmotic equation of state for clay, extracted from the
profiles of Fig. \ref{fig_clay_a}. The stars are results obtained from
direct NVT MC simulations.}
\label{fig_clay_b}
\end{figure}

\begin{figure}
\caption{Averaged monomer-monomer correlation function. Dashed line:
PB simulation. Dotted line: simulation using the linearized
functional. Full line: theoretical result in the absence of charges.}
\label{fig_persis}
\end{figure}

\end{document}



