\def\rr{\rho ({\bf r})}
\def\pmb#1{\setbox0=\hbox{$#1$}%
\kern-0.025em\copy0\kern-\wd0
\kern0.05em\copy0\kern-\wd0
\kern-0.025em\raise0.0433em\box0} 
\documentstyle[preprint,eqsecnum,aps]{revtex}

\begin{document}
%\draft
	
\title{The Dynamic Structure of Liquid Sodium from {\em Ab-Initio} Simulation}
\author{M. Foley, E.Smargiassi and P.A. Madden}
\address{
Physical Chemistry Laboratory, Oxford University, \\
South Parks Road, Oxford OX1 3QZ, U.K.}
\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}
%
%  Abstract goes in here.
%
 Data which pertain to experimental X-ray and inelastic neutron 
 scattering studies of liquid sodium are obtained from an {\it ab-initio}
 molecular dynamics 
 simulation.  A  
 long run was used, so that the statistical precision of
 the calculated quantities is as good as obtained experimentally. The
 long run is made possible by  the use of a recently 
 developed modification of the Car-Parrinello method which uses 
 an orbital-free density functional to prescribe the electron energetics.
 Agreement with experiment is excellent. 

\end{abstract}
 
\narrowtext

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}
 {\em Ab-initio} molecular dynamics simulations, as 
pioneered by Car and Parrinello \cite{Car_Par}, have greatly expanded the
 range of physical systems which can be accurately represented in 
 computer simulations \cite{Galli}. 
 In this method the nuclear and
electronic degrees of freedom describing the system are treated on the same 
 footing. In most practical calculations 
 the electronic degrees of freedom are the coefficients in the expansion of occupied orbitals in some basis set within the Kohn-Sham realisation of density
 functional theory \cite{Kohn}
  and evolve simultaneously in time with the nuclear motion \cite{Remler90}. 
However,  because of the amount of computer time required for the computations,
   calculations of  properties on systems of 100 atoms or more
  and which require long runs, either 
 for statistical averaging (such as free energies) or to study relatively
 slow dynamical processes, have not been made.

 The computational cost 
 arises from two sources. Firstly, the appropriate timestep in the C-P
 integration is dictated by the fast dynamics of the orbital coefficients
 \cite{Pastore}. Secondly, because of the requirement that the occupied orbitals
 be kept orthogonal, the cost of a single timestep scales like $N^2M$ for
 large sytems, where $M$ is the basis size and N the number of electrons.
 These problems are not circumvented by performing calculations with
 a timestep dictated by the nuclear motion and re-optimising the 
 wavefunction at each step \cite{Payne}. The cost of a single iteration has 
 the same poor scaling behaviour and the number of iterations required to
 achieve convergence depends on the nature and size of the system.

 Various ways of overcoming the scaling
problem  have been suggested (so-called ``order $N$" methods).
 Most \cite{Galli2,Baroni,Vanderbilt}
  depend on schemes which are appropriate when  the electronic density
 matrix is short-ranged, which is true for systems in which there is
 an energy gap between the occupied and virtual Kohn-Sham orbitals.
  We have followed a different approach \cite{Pearson93},
 which for the purpose of modelling the
electronic behaviour of the system, uses an orbital-free density
functional rather than the usual Kohn-Sham methodology. 
 Such an approach is explicitly of order $N$ for each electronic step.
 We have found suitable functionals for metals 
\cite{Smargiassi}, for which the density matrix decays algebraically,
  so that the method is, at present,
  complementary to the other
 order $N$ methods which have been introduced. A further
 advantage over the Kohn-Sham scheme for metals is that the 
 electronic degrees of freedom behave adiabatically \cite{Pearson93},
 a point which we further illustrate below. With the density as the
 basic variable we have shown that the electronic timestep is 
 significantly longer than with a Kohn-Sham scheme and can be made {\em  
 independent of system size}, so that the method is a {\em true}
 order $N$ method for dynamics. 

This opens the prospect of obtaining types of data which have
hitherto been unavailable from {\em ab-initio} molecular
dynamics simulations.
  In this
paper, we demonstrate that {\em ab-initio} methods 
can be used to generate  dynamical information at the level of 
 experimental accuracy and precision. Elsewhere, we illustrate the application to free energy
 calculations \cite{Smargiassi2}.
 
\section{Simulation Method}

The method used, which is exactly as described in reference
\cite{Pearson93}, will be outlined here.  
We represent the system being studied as a
set of ions in a cloud
of valence electrons, with the spatial density distribution 
of this cloud, $\bf \rho(r)$, being expanded in a Fourier 
series
\begin{eqnarray}
\rho({\bf r}) = \sum_{\bf g} \rho_{\bf g} e^{i{\bf g}.{\bf r}},
\end{eqnarray}
where the `$\bf g$'s are reciprocal-lattice vectors of the periodically
 replicated simulation sytem. 
Velocities and non-physical 
``masses'' $\bf \mu_g$ are allocated to
the Fourier coefficients.  Coupled equations of motion
are then derived for both the ion positions $\bf \{R_\alpha\}$
and the electronic coefficients $\bf \rho_g$:  
\begin{eqnarray}
\mu_{\bf g}\ddot {\rho}_{\bf g} =-{\delta E \over
\delta \rho_{\bf g}}, 
~~~~~ \mid {\bf g} \mid ^2 ~ \leq ~2E_c
\end{eqnarray}
and
\begin{eqnarray}
M_I \ddot
{\bf R}_{\alpha} = - \nabla _{\alpha} E ~ - \nabla _{\alpha}
V_{II},~~~~~{\alpha}=1,\ldots,N_I.
\end{eqnarray}
where $V_{II}$ is the potential energy of the direct ion-ion
interaction and the electronic energy 
 $E$ is a functional of $\bf \rho(r)$ for a
given $\bf \{R_\alpha\}$.  $N_I$ is the number of ions, and $E_c$ the energy cutoff for the planewave expansion.
 Molecular dynamics simulations are performed by
simultaneously integrating these equations of motion with the $\rho_{\bf g}$ 
 appropriate to the electronic ground state.

 As in the majority of Car-Parrinello
style simulations, $E[\rho]$ is provided by
density functional theory. However, for the reasons given above, the
 kinetic energy functional is constructed without the introduction of
 orbitals. How this is done, following work of Perrot \cite{Perrot}
 and similar developments by Teter and Wang \cite{Teter_Wang},
 is described in detail in reference \cite{Smargiassi}.
   The remainder of the energy functional
comprises the interaction of the electron density with
the ion potential (represented by the Topp-Hopfield pseudopotential \cite{Topp}), the coulombic electron-electron interaction and
the local exchange-correlation functional of Ceperley and Alder\cite{Ceperley}.

A simulation was performed on a system of 128 sodium atoms 
at $388 \pm 18$K with number 
density $0.0241{\rm\AA}^{-1}$ (as used in \cite{Kinell85}), in a
cubic cell with periodic boundary conditions, and starting from a 
 previously equilibrated liquid
configuration.  
The run was for 96ps, taking around 250 hours of CPU time on a
Convex C2 \cite{Comment}.  The molecular dynamics timestep
was $4.8\times 10^{-4}$ps, half that used in simulating the
solid \cite{Pearson93}.  The plane-wave cutoff  was 11.7
Rydberg,  larger than that used in the solid-state work.  These conservative 
 choices of parameters
were  found to be necessary to 
 avoid energy drift over the 
extremely long run length.  A
shortcoming of simulations using the Car-Parrinello algorithm in the Kohn-Sham 
scheme  is that energy can leak from the electronic to
the ionic degrees of freedom in metallic systems, for reasons 
 discussed in reference \cite{Pastore}. This can be cured by the
 introduction of ``thermostats" which control the ``temperatures" of
 the two subsytems \cite{Fois,Blochl}, at a cost of, at least, some elegance.
   Reference \cite{Pearson93} describes
how the present technique is able to minimize this ``non-adiabaticity" 
 problem   by adjusting the fake masses
assigned to the electronic degrees of freedom so as to widely
separate the frequencies of coefficient and ionic motions (``conditioning" the
 algorithm). The
target value for the coefficient frequencies was around
$580{\rm ps}^{-1}$.  The combination of conservative step length, more
accurate representation of the electrons and conditioning
 meant that no thermostatting was required:
  we did  requench the electronic degrees
of freedom three times during the course of the run, but
 this had no discernible effect on the ionic trajectories.

\section{Calculation of Observable Properties}

The properties of the simulated fluid were studied by means of
various correlation functions, some of which are also obtainable
from neutron and X-ray scattering experiments.  The calculated quantities
 were fourier coefficients of certain densities and
were obtained  for a range of
wavevectors (${\bf k}$) consistent with the periodic boundary conditions.
The correlation functions were averaged
over vectors of equal magnitude $k$ (as the
liquid is isotropic). 

\subsection{Static Structure Factor and X-ray Diffraction.}

The static structure factor 
\begin{equation} 
S(k) = \frac{1}{N_I}
\left< \left( \sum_{\alpha=1}^{N_I}
{ e^{-i {\bf k}.{\bf R}_\alpha(t_0)}}
\right) \left( \sum_{\alpha=1}^{N_I} { e^{i {\bf k}.{\bf R}_\alpha(t_0)}}
\right) \right>,
\end{equation}
is shown in figure \ref{static} in 
 comparison with data from an X-ray scattering experiment
conducted on liquid sodium at 373K \cite{Greenfield71}. The simulation values
 are calculated directly as indicated in the equation, rather than by transformation of the radial distribution function, so as to avoid truncation 
 errors.
  It can be seen that the overall 
agreement is excellent, except for the height of the first peak,
which is overestimated. 
 The figure shows results for $S(k)$ averaged over a relatively short 
 interval at the start of the run, over the first half of the run
 and over the full run. It can be seen 
that the peak has been growing and/or
 shifting to slightly lower $k$ as the run progresses.
 One $k$ point at the second peak becomes slightly
more pronounced as well but the remaining points are unchanged.  The increasing sharpness of the first peak implies that the simulated fluid is 
 developing longer range order at this wavevector than is present in the
 real fluid.
  We believe that we a seeing
an artifact of the small system size, exacerbated by the proximity of the
 simulated fluid to its solid phase.  Given that 128 atoms can form a
body-centred cubic lattice in a cubic box, that the density of sodium
changes by only around 2\% on freezing, and that we are in the
vicinity of the freezing point of sodium, we may be observing incipient
  solid-like order induced by the periodic boundary conditions.
   The $k$-vector which is affected by
this is the one at the position where the first Bragg peak
would appear for a lattice of solid sodium at our system's density.
However, the Bragg peak for a perfect crystal would have a height of
128, so the extent of the enhancement we see (of about 0.5)
suggests that any such effect is only small. We do not believe that
 this is artefact is a consequence 
  of the {\it ab-initio} nature of the simulation.
 Indeed simulation studies of liquid metals close to freezing which
 use effective interionic interactions have often shown an enhancement 
 of $S(k)$ at its first maximum over experiment \cite{Watabe94}. If our
 reasoning is correct, the simplest way to avoid 
 the artefact would be to use a particle number which could not 
 crystallise in the simulation cell and to use more particles.

\subsection{Dynamic Structure Factors and Neutron Scattering.}
 
Foremost amongst the correlation functions used to examine the 
dynamic properties of
the system is the intermediate scattering function
\begin{equation} 
F(k,t) = \frac{1}{N_I} \left< \left( \sum_{\alpha=1}^{N_I}
{ e^{-i {\bf k}.{\bf R}_\alpha(t+t_0)}}
\right) \left( \sum_{\alpha=1}^{N_I} 
{ e^{i {\bf k}.{\bf R}_\alpha(t_0)}} \right) \right> 
\end{equation}
which on Fourier transformation yields the dynamic structure
factor $S(k,\omega)$.

The time transform was taken over a 
  correlation  length of 1.50ps.  Multiplication
with a Blackman windowing function \cite{Allen&Tildesley} was performed
so as to remove effects due to the
finite length of the correlation function
and then Fourier transformation was carried out numerically by
the Filon method \cite{Allen&Tildesley}. Although the functions were
 calculated out to much longer times, in the hope of obtaining high
 resolution in frequency space, we observed that they became quite
 noisy at times longer than about 1.5ps, despite the long run length. 
 We attribute this to 
 the propagation of disturbances between periodic images. Using such a
 short correlation length means that the frequency resolution in
 the dynamic structure factor is only of order 1-2$ps^{-1}$. 
  The resulting
dynamic structure factor is displayed in Figure \ref{dynamic}: it is clear that
 any significant structure in this function occurs at frequencies higher than
 2-3$ps^{-1}$ and is thus unaffected by the windowing process.

There is a clear Brillouin peak at around $10{\rm ps}^{-1}$ for
the smallest ${\bf k}$-vector.  This persists almost up to the
principal peak of the structure factor, showing the presence of
collective motions in the fluid at these wavevectors, but it is
not observed beyond the principal peak. The persistence of a Brillouin peak
 up to much higher wavevectors  than observed in ionic or atomic fluids 
 is a well known property  of metallic liquids \cite{jacucci}.

In principle \cite{Hansen&McDonald}, $S(k,\omega)$ is
experimentally obtainable from coherent neutron 
scattering. In practice, the neutron
scattering experiments on sodium have not been able to separate
the    coherent  from the incoherent scattering as 
cross-sections for sodium, $\sigma_c$ and $\sigma_i$, are of the
same order ($\sigma_c$=1.66 barns, $\sigma_i$=1.74 barns) \cite{soderstrom}.
  What is immediately accessible is the {\em
effective} dynamic structure factor
\begin{equation} 
S_e(k,\omega) = \frac {\sigma_c S(
k,\omega) + \sigma_i S_s(k,\omega)} {\sigma_c +
\sigma_i}.
\label{eq:eff}
\end{equation}
where $S_s(k,\omega)$ is the the self dynamic structure factor arising from
the  incoherent scattering and 
  is the Fourier transform of the self intermediate scattering function:
\begin{equation}
F_s(k,t) = \frac{1}{N_I} \sum_{\alpha=1}^{N_I}
{ \left< e^{-i {\bf k}.{\bf R}_\alpha(t+t_0)} e^{i {\bf k}.{\bf
R}_\alpha(t_0)}\right>}.
\end{equation}
Such data have been presented in Reference \cite{Kinell85}.

We have constructed the effective dynamic structure factor 
from our simulation  data using Equation \ref{eq:eff}, where 
the $S_s(k,\omega)$ was  obtained by numerical Fourier
transform of the 12ps long self correlation function after
multiplication with a Blackman window.  (No problems with
interference from periodic images are encountered in the self
correlation.) The comparison with experiment over the range of $k$
 probed experimentally is shown in figure \ref{Figure:Seffective}. In practice,
 because of the fact that the amplitude of $S(k)$ is low except at the largest 
 $k$ value, and because $S(k,\omega)$ is much broader that $S_s(k,\omega)$, the
 shape of these functions is dominated by the incoherent scattering.

 Agreement with experiment is very good for most
of the range, except  at the largest value of $k$. Even at this $k$ value,
 if the amplitude of the
 simulated function is suitably scaled, the experimental and simulated curves are brought into excellent agreement; the disagreement apparent in figure
\ref{Figure:Seffective}e is therefore one of amplitude. The relevant $k$
 value lies
on a part of $S(k,\omega)$ which is very steep with respect to
$k$ because of the proximity to the principal peak of
$S(k)$, and so the coherent contribution to $S_e(k,\omega)$
  would be expected to be very sensitive to small
discrepancies in $S(k)$. However, to give the error in amplitude apparent
 in \ref{Figure:Seffective}e, the error in S(k) at this value of $k$
 would have to be of order 100\%, rather than the $\sim$5\% apparent from
figure \ref{static}, and changing the amplitude of the coherent relative to the incoherent would change the spectral shape for which the agreement is good. 
 It seems to us that the experimental amplitude is unreasonably low.
  At $k$=1.53, for which $S_e(k,\omega)$ is
 shown in figure \ref{Figure:Seffective}d, the effective scattering is
 dominated by $S_s(k,\omega)$. At $\omega$=0 the amplitude of $S_s(k,\omega)$
 should vary as $k^{-2}$. The experimental value of $S_e(k,0)$ at $k$=1.80
 is {\em lower} than would be expected from the  scaling $k$=1.53 data, i.e.
 lower than the expected incoherent scattering {\em alone}, without allowing 
 for the coherent contribution. 

Even though $S_s(k,\omega)$
has not been separately resolved from experiment, attempts have
been made to recover it indirectly at very low energy transfer through
the use of an estimate of the coherent contribution to
 the measured effective dynamic structure factor
\cite{Morkel86,Stangl93}.  Such data from a
very recent neutron scattering experiment \cite{Stangl93} are presented in
Figure \ref{selfdynamic} and compared with simulation data.
Since the neutron data are reported at very low frequency
  the simulated $S_s(k,\omega)$ 
was obtained by fitting the calculated $F_s(k,t)$ to  
  an exponential function
 for times in the range 3ps to 8ps and using this exponential to extend the 
 range of integration in calculating the transform.
  Over the range (1 -- 12 $ps$)  this fit was
 excellent and could not be systematically improved by choice of any other function.
 The transform was calculated by combining the analytic transform of the
 exponential with the numerical transform of the residue of the fit. 
 Doing the transform in this way allowed comparison 
to be made with the
very small values of $\omega$ where experimental data is
available. In figure \ref{selfdynamic}e, the value
 of $S_s(k,\omega)$ obtained from a direct transform of $F_s(k,t)$
 Blackman-windowed over 12ps range is shown, to illustrate the
 $\omega$ value at which the extrapolation procedure becomes significant.

 Good agreement is found between simulation and
experiment on the shape, width and amplitude of $S_s(k,\omega)$ is found,
 showing that the diffusive behaviour indicated by the exponential 
 representation of  the simulation data is realistic.
  The greatest discrepancy is found at $
k\approx 1$\AA (figure \ref{selfdynamic}e), although this can be partially attributed to 
the difference in $k$-vector used in simulation and
experiment.  Notice that our agreement with the  experimental data, for
$S_e(k,\omega)$, is good around this value of $k$.

\subsection{Current correlation functions}

Also calculated from the simulation were the longitudinal and
transverse current correlation functions
\begin{equation} C^l(q,t) = \left< \left( \sum_{\alpha=1}^{N_I}
{{\bf k}.{\bf v}_I(t+t_0) \over q} e^{-i {\bf k}.{\bf
R}_\alpha(t+t_0)}
\right) \left( \sum_{\alpha=1}^{N_I} {{\bf k}.{\bf v}_I(t_0) \over q}
e^{i {\bf
q}.{\bf R}_\alpha(t_0)} \right) \right> \end{equation}
and \begin{eqnarray} C^t(q,t) = \left< \left( \sum_{\alpha=1}^{N_I}
\left( {\bf v}_I(t+t_0) - {{\bf k}.{\bf v}_I(t+t_0) \over {\bf
q}.{\bf
q}} {\bf k} \right) e^{-i {\bf k}.{\bf R}_\alpha(t+t_0)} \right)
\right.
\nonumber\\ \left. ~~~~~~~~~~~~~~~~~~~~~
. \left(
\sum_{\alpha=1}^{N_I} \left( {\bf v}_I(t_0) - {{\bf k}.{\bf v}_I(t_0)
\over {\bf
q}.{\bf k}} {\bf k} \right) e^{i {\bf k}.{\bf R}_\alpha(t_0)} \right)
\right>. \end{eqnarray}
 The Fourier
transformation was performed as for $S(k,\omega)$.
The positions of the maxima of the peaks $\omega_l^{max}$ and $\omega_t^{max}$
in the Fourier transforms of these
functions provide curves analogous to phonon dispersion
curves for the solid \cite{Hansen&McDonald}. 
These are shown in Figures \ref{long} and \ref{trans}. These 
 curves have the same
form as seen in other simulations of alkali metals \cite{Kahl92}
  and, for $\omega_l^{max}$,
 with values which have been deduced from the peaks seen 
 in $\omega^2 \times S(k,\omega)$ from experimental scattering data
 \cite{Winter93}. These
 peaks have been attributed to the persistence of propagating sound modes
 even at these very high values of $k$, as anticipated from certain theories
 \cite{deSchepper}.

  

 One feature of interest, which might be seen in our data, is a ``Kohn anomaly"
 \cite{Ashcroft}  in the ``dispersion curves" 
  due to the singularity in the electronic reponse function caused 
 by  Fermi surface effects. This occurs at a wavevector equal to twice the 
 Fermi wavevector, or $2k_F = 1.79 \AA^{-1}$ at the density we have studied.
 Because of the {\it ab-initio} nature of our simulation, the origins of the anomaly are properly represented. Examination of Figure \ref{long} shows
 that any Kohn anomaly in the Brillouin peak position would be very difficult to
  observe, since $2k_F$ is very close to the position of the first peak in 
 the dynamic structure factor and hence the Brillouin peak is at very low
 frequency and very strongly damped. Our data is suggestive of a kink in 
$\omega_l^{max}$ at $2k_F$,
  but it should be remembered that $\omega_l^{max}$ has
 been extracted from the data as the position of a peak in a fairly
 broad spectral function. Figure \ref{trans} shows that $\omega_t^{max}$ does exhibit a sharp peak at
 $2k_F$: unfortunately,
  this function is not experimentally observable in a fluid. It would be 
 interetsting to compare this behaviour with that found in a sodium 
 simulation with an effective pair potential, which should not
 show a Kohn anomaly.
  
\section{Discussion}

 The primary purpose of this paper has been to show that {\it ab-initio}
 methods  can be used to simulate the dynamical properties of liquids
in the same 
manner as in conventional simulations based on effective internuclear 
potentials and 
 at the level of statistical precision achieved in experiments.
This is the first time {\it ab-initio} techniques
have been used to obtain comprehensive dynamical information about 
this---or indeed, any---liquid system.
 For the most part, agreement with experiment is excellent. 
 We have carefully discussed
 all apparent discrepancies between the simulated and experimentally observed behaviour: only in the case of the anomalously high peak in $S(k)$
 do we believe the discrepancy to be due to simulation methodology and
 here we anticipate that by working with a different number of 
 particles we  would eliminate even that.

  Ideally, the only 
 input to  a calculation of the type   we have described 
 is  the atomic number of sodium. However, we have fallen short of that ideal in
 two regards. We made use of the Topp-Hopfield pseudopotential
 because of its desireable computational 
 properties and because it does give a good account
 of the properties of the solid; however, this  is
 not a true {\it ab-initio} potential.
  Secondly, we performed the calculations 
 at the experimental liquid density at the melting point.
  In practice our simulated system is under a
 pressure of about 3 Kbar at this density. The source of this error
 could be associated with the pseudopotential, and further work is in progress
 with {\it ab-initio } potentials  to examine this possibility.
  However, the error corresponds to
 an overestimate of the lattice parameter of the solid of about 2\%, which
 is typical of calculations which use the local density approximation to
 to the exchange and correlation functionals and it may be that this 
 represents the real limitation on accuracy. 

 In order to accomplish the calculations 
 in  a reasonable amount of time, we have made use of a kinetic energy
 functional which has a restricted range of applicability (but does give
 an excellent representation of solid sodium).  The
 calculations otherwise employ exactly the methods laid down by Car and 
 Parrinello in their classic paper. Similar dynamical 
 calculations will become
 possible for more ``interesting" fluids than sodium both by
 further developments in the orbital-free functional scheme, which
 will expand its domain of applicability without sacrificing the
 advantages which we have sought to illustrate, and by improvements in
 Kohn-Sham methodology and computer technology.

The reason for regarding liquid sodium close to its freezing point
 as not ``interesting" is that it 
 has already been successfully treated using
conventional simulation methods \cite{Kinell85,Balucani93,Watabe94} with
 effective pair potentials derived from empirical pseudopotentials
 \cite{Hafner}.  The relationship of our methodology to that used in deriving effective potentials for metals was discussed in reference \cite{Smargiassi}.
 The quality of agreement with experiment reached in the best of
 these calculations is quite comparable to that we have demonstrated here.
   
\acknowledgments % % SERC,RSRE,IBM et al.  %
 We thank Prof. C. Morkel for sending us the
 incoherent neutron scattering data before publication.  
The work was supported by SERC and the MOD (UK) under grant
 GR/H38256.  It has also benefitted from further support from SERC
 under  grant no.
  GR/H10276. Michael Foley is  grateful to the Association of
 Commonwealth Universities for the award
 of a Commonwealth Scholarship.

\newpage
\begin{thebibliography}{00}

\bibitem{Car_Par} R. Car and M. Parrinello,
{\it Phys. Rev. Lett.} {\bf 55}, 2471 (1985).
\bibitem{Galli} G. Galli and A. Pasquarello, {\it Computer Simulation in
 Chemical Physics} Eds. M.P. Allen and D.J. Tildesley (Kluwer, Dordrecht, 1993).
 \bibitem{Kohn} W. Kohn and L.~J. Sham,
{\it Phys. Rev.} {\bf 140}, A1133 (1965).
\bibitem{Remler90} D. K. Remler and P. A. Madden, Mol. Phys. {\bf 70}, 921 (1990).
\bibitem{Pastore} G. Pastore, E. Smargiassi, and F. Buda,
{\it Phys. Rev. A}, {\bf 44}, 6334 (1991).
\bibitem{Payne} M.C. Payne, M.P. Teter, D.C. Allen, T.A. Arias and J.D. Joannopoulos, {\it Rev. Mod. Phys.} {\bf 64}, 1045 (1992).
\bibitem{Galli2} G. Galli, M. Parrinello,
{\it Phys. Rev. Lett.} {\bf }, (199?).
\bibitem{Baroni} S. Baroni and P. Gianozzi, {\it Europhys. Lett.  } {\bf
 17 }, 547 (1992). 
\bibitem{Vanderbilt} X.P. Li, W. Nu\~nes and  D. Vanderbilt, {\it
 Phys. Rev.  } {\bf B47 }, 10891 (1993).

\bibitem{Pearson93} M. Pearson, E. Smargiassi and P. A. Madden, J. Phys.: Condens. Matter
{\bf 5}, 3221 (1993).
\bibitem{Smargiassi} E. Smargiassi and P.A. Madden, {\it Phys. Rev.}{\bf B49},
 5220 (1994).
\bibitem{Smargiassi2} E. Smargiassi and P.A. Madden -- to be published.
\bibitem{Perrot} F. Perrot, to be published.
 \bibitem{Teter_Wang} L.-W. Wang and M.~P. Teter, 
{\it Phys. Rev. B} {\bf 45}, 13197 (1992).
 \bibitem{Topp} W.~C. Topp and J.~J. Hopfield,
{\it Phys. Rev. B} {\bf 7}, 1295 (1973).

\bibitem{Ceperley} D.~M. Ceperley and B.~J. Alder,
{\it Phys. Rev. Lett.} {\bf 45}, 566 (1980); J.~P. Perdew and A. Zunger,
{\it Phys. Rev. B} {\bf 23}, 5048 (1981).

\bibitem{Kinell85} T. Kinell, U. Dahlborg, O. S\"{o}derstr\"{o}m and
I. Ebbsj\"{o}, J. Phys. F: Met. Phys. {\bf 15}, 1033 (1985).

\bibitem{Comment} Further improvements to the code since this run was made
 would make the run time
 90 hours on an IBM RS6000 model 320H, using the ESSL library for the 
 Fourier transforms.

  \bibitem{Fois} Fois, E.S., Penman, J.N. and Madden, P.A., {\it
J. Chem. Phys.} {\bf 98}, 6361 (1993).
\bibitem{Blochl} P.E. Bl\"ochl and M. Parrinello, {\it Phys. Rev.} {\bf B45} 9413 (1992).

 \bibitem{Parr&Yang} R. G. Parr and W. Yang, {\em Density-Functional Theory of Atoms
and Molecules} (Oxford University Press, New York, 1989).
\bibitem{jacucci} G. Jacucci and I.R. McDonald, {\it Liquid and Amorphous Metals} Edited by  E. L\"uscher and H. Coufal (Sijthoff and Noordhoff, 1980). 

\bibitem{Hansen&McDonald} J. P. Hansen and I. R. McDonald, {\em Theory of Simple Liquids} (Academic, London, 1986).
\bibitem{soderstrom} O. S\"oderstr\"om and U. Dahlborg, {\it J. Phys. F: Metal Physics} {\bf 14}, 2297 (1984).
\bibitem{Greenfield71} A. J. Greenfield, J. Wellendorf and N. Wiser,
Phys. Rev. A {\bf 4}, 1607 (1971).
\bibitem{Watabe94} F. Shimojo, K. Hoshino and M. Watabe,
submitted for publication (1994).

\bibitem{Allen&Tildesley} M. P. Allen and D. J. Tildesley, {\em
Computer Simulation of Liquids} (Clarendon, Oxford, 1987).
 \bibitem{Morkel86} C. Morkel and W. Gl\"{a}ser, Phys. Rev. A {\bf 33}, 3383 (1986).
\bibitem{Stangl93} A. Stangl, Ph.D. thesis, Technishche
Universit\"{a}t M\"{u}nchen (1993).
\bibitem{Kahl92} S. Kambayashi and G. Kahl, Phys. Rev. A {\bf 46},
3255 (1992).

\bibitem{Winter93} R. Winter et al., J. Non-Cryst. Solid. {\bf 156--158},
9 (1993).

\bibitem{Ashcroft} N.~W. Ashcroft and N.~D. Mermin,
{\it Solid State Physics} (Holt-Saunders, Philadelphia, 1976); W.~A. Harrison,
{\it Solid State Theory}, (Mc Graw-Hill, 1970)

\bibitem{deSchepper} I.M. de Schepper and E.G.D. Cohen. {\it J. Stat. Phys.} {\bf 27}, 123 (1982).

\bibitem{Balucani93} U. Balucani, A. Torcini and R. Vallauri, Phys.
Rev. B {\bf 47}, 3011 (1993).
  

\bibitem{Hafner} J. Hafner,
{\it From Hamiltonian to Phase Diagrams} (Springer-Verlag, Berlin, 1987).


\end{thebibliography}

  
\newpage

%\figure{ caption
%\label{name}}

\begin{figure}
\caption{ The experimental static structure factor, from X-ray measurements
 (ref. 25) is compared with the simulation results. The symbols denote values
 calculated by averaging over different portions of the run, $\bigcirc $ --
 first 25 $ps$, + -- first 50 $ps$, $\times$ -- full run (96 $ps$).
\label{static}}
\end{figure}


\begin{figure}
\caption{ The simulated dynamic structure factor, $S(k,\omega)$, in a 3-d
 representation and with a superimposed contour plot. 
\label{dynamic}}
\end{figure}


\begin{figure}
\caption{ The effective dynamic structure factor, $S_e(k,\omega)$, obtained
 from neutron scattering (crosses, ref. 17) is compared with the simulation results
 (full lines) at several k values. From top to bottom the (experimental, simulation) k values ($ \AA^{-1}$) are:- (0.69,0.72), (0.94,1.02),
(1.25,1.25), (1.53,1.53) and (1.80,1.80).
\label{Figure:Seffective}}
\end{figure}


\begin{figure}
\caption{ Comparison of experimental self dynamic structure factors (refs.
 28, 29) shown as 
 crosses,
 with simulation values at several values of k.  From top to bottom the (experimental, simulation) k values ($ \AA^{-1}$) are:- (0.383,0.360), 
(0.518,0.509), (0.648,0.624), (0.740,0.720), and (0.966,1.02). The bottom figure also shows the simulation data transformed by truncation and Blackman 
 windowing at 12$ps$, dashed line.
\label{selfdynamic}}
\end{figure}


\begin{figure}
\caption{The position of the peak in the spectrum of the longitudinal current 
 correlation function. Note the close similarity to the behaviour of the contours in figure 2.
\label{long}}
\end{figure}


\begin{figure}
\caption{  The position of the peak in the spectrum of the transverse current 
 correlation function.
\label{trans}}
\end{figure}


   
\end{document}
