\documentclass[12pt]{article}
%\documentstyle[preprint,aps]{revtex}
\textwidth 16.5cm
\textheight 24cm
\topmargin -1.5cm
\evensidemargin 0cm
\oddsidemargin 0cm
\parindent 2.5em
\begin{document}

\raggedbottom

\begin{center}
\Large \bf
Computational Physics \\
\large 
\vspace{0.2cm}
21-may-2001
\end{center}

\noindent 
{\bf Ex.11  Molecular dynamics of Lennard-Jones atoms.} 


The basic program based on the 
velocity Verlet algorithm for a system of particles interacting via
a Lennard-Jones pair potential is  given in the fortran77 program $md.f$
which needs $md.inc$. Besides, a third file $md.old$
contains a starting configuration for the positions and velocities of the 64
particles considered. The program makes use of periodic boundary conditions for the evaluation of the forces.
The program has been written for didactical purposes by Dr. M. P. Allen, 
University of Bristol. It is well documented in itself and the lecture notes 
provide all the necessary background. A similar program in fortran90 by
Dr. F. Ercolessi (SISSA, Trieste) can also be used as guideline in rewriting 
your own program.

Notice that the calculation is performed 
in reduced units so that it applies to any system of particles interacting via 
a LJ potential. This potential may fairly well represent argon atoms if
$\epsilon/k_B=120 K$ and $\sigma$=3.4 \AA~. The mass of an Ar atom is 
39.95~X~1.66~X~10$^{-27}$ Kg.

The original purpose of this exercise 
was to calculate the diffusion coefficient
of a particle in a fluid as described at point 5. Therefore, in the program  
atom 1 has a mass of 0.5 units into the subroutines kinetic, movea and moveb
since this makes its motion into the liquid more noticeable.

Points, 1, 2 and 4 are compulsory, point 3 concerns visualization and is just 
given for information. You can choose whether to calculate the 
the diffusion coefficient as described at point 5.
Read all points anyway before deciding.
\begin{itemize}

\item[1)] read the program and understand how it works. Note that atom
1 has a mass of 0.5 units into the kinetic and move subroutines. 
We can think of it as 
a Ne atom in an Ar gas. If you do not intend to carry out the  
calculation at point 5 of the diffusion coefficient change the program so 
that also atom 1 has the same mass as the others.

\item[2)] Run the program beginning with the  
starting configuration $md.old$: use a trial timestep  
$\delta t$=0.005, a cut-off radius of 2.25 and run 200 blocks of 10 
time steps. Write the output in a file instead that 
on the screen.

Transform all average calculated values from adimensional to physical units,
using the values appropriate to argon 
(temperature in Kelvin, energy in eV,  ..).

\item[3)] 
Visualization of the motion can 
be obtained by use of the program molden (courtesy of Gijs Schaftenaar, 
University of Nijmegen) which needs a file $md.conf$ containing the 
configurations after each block. For each configuration, the file must contain,
in free format, in the first line the number of particles $n$, in the second 
line the step number $nblock \times nstep$, from the third to the $(n+2)th$ 
line a string with the atom type (ar or ne) and the three coordinates x, y, z 
in \AA.
An example of this file is $md.conf$.
Once you have run the same simulation as at point 2) and created the file 
md.conf, you type :

molden md.conf
 
and a screen will appear. Click on Solid to see the non-bonded atoms.
By clicking on the apple you can change the colors of atoms or change the Van 
der Waals radius (a bond is drawn between atoms which 
dist less than twice this value from each other). Change the colour of Ne to 
red so that you can follow its motion. Use a value of 1.8~\AA for the VDW 
radius.

If everything works, you see  the Ne atom moving around and you will also see 
the effect of the periodic boundary conditions. 

\item[4)] The choice of the time step does not have strict rules. It should be 
as large as possible to allow long simulations (by the way have you calculated 
your time steps in seconds?) and small enough as to ensure energy conservation
(the simulation is conducted in a microcanonical ensemble). A rule of the thumb is that fluctuations in the total energy should not exceed a few per cent of 
the fluctuations in the potential energy. Here there two options. One is just 
to try other (larger values) of the time step and see what happens, the second 
is to do the same while monitoring the value of the potential energy as already
done in the programs for the other quantities (energy, temperature etc).
On the basis of what you find can you confirm the above rule of the thumb?

\item[5)] Calculate the mean square displacement of atom 1 along the 
$x$-direction $<[x(t)-x(0)]^2>$ and its velocity autocorrelation function
$<v_x(t) v_x(0)>$ for $0<t<\Delta t$, with $\Delta t=nstep \times \delta t$. 
From these quantities the  
diffusion coefficient $D$ of atom 1 
can be calculated by use of the Einstein relation which relates it 
to the mean square displacement :
\begin{equation}
D= lim_{t \rightarrow \infty} {<[x(t)-x(0)]^2> \over 2t}
\end{equation}
and by use of the Green-Kubo relations which relates it to the integral of the 
velocity autocorrelation function:
\begin{equation}
D=lim_{t \rightarrow \infty} \int_0^t <v_x(\tau) v_x(0)> d \tau
\end{equation}

To calculate these quantities you can perform  simulations over a time 
interval of 200 steps and average over at least 100 blocks. 
The method is described in 
Sect. 7.1.2 of the book of Haile given in the notes. 
The idea is that you perform a simulation over a time interval 
$\Delta t=nstep \times \delta t$ and do an average of the behaviour over
$nblock$
simulations, where the time at the beginning of each block can be thought of as a new time origin. If we call, the time at the beginning of each block 
$t_k=k \Delta t$, with $k=1,nblock$, the average velocity autocorrelation 
function at time $t$,
can be calculated as
\begin{equation}
<v_x(t) v_x(0)>={1 \over nblock} \Sigma_{k=1}^{nblock} v_x(t_k) 
v_x(t_k+t)
\end{equation}
At the beginning of the first block, record $v_x(t=0)$ into a variable V0
and store
$v_x(0)^2$  as the first component of a vector VVX($nstep$). At each successive 
step, store in the successive elements of the vector the quantity
$v_x(o)v_x(t)$. In the successive blocks, set V0=$v(t_k)$ 
and add to each component of the 
vector VVX the value V0$\times v_x(t_k+t)$. At the end of the $nblock$ blocks
just divide each component of VVX by the number of blocks and plot it
as a function of time. Then calculate $D$ according to eq.(2).
 


\end{itemize}
\end{document}



