\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 $$ for $0 \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 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} ={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}