Computer simulations in statistical physics
Homework 5: Molecular dynamics simulations
Usage and functionality of program MD LJ
First analysis of program output (default parameters)
Code of program MD LJ
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 1
Usage and functionality of program MD LJ
What can the program do, what are the parameters?
WS2006_Simulation/HW5_MD> ./MD_LJ -h
#************************************************************#
# MD_LJ: molecular dynamics program for Lennard-Jones atoms #
# Version: 0.6, 22.12.2006 by Nils Bluemer, compiled for D=3 #
# For help on usage and available options: run MD_LJ -h #
#************************************************************#
# options: -h this help #
# -v increase verbosity #
# (-n# number of particles - set automatically) #
# -l* lattice type for init: (h)ypercubic or (f)cc #
# -s# system size (# particles/row) | default: 8 #
# -r$ density | default: 0.600 #
# -e$ energy per particle | default: 0.500 #
# -w$ warmup time | default: 0.200 #
# -t$ simulation time | default: 2.000 #
# -d$ timestep | default: 0.002 #
# -c$ cutoff (of potential/force) | default: 2.500 #
# -b# number of bins in paircorr | default: 200 #
# symbol key for parameters: # integer, $ real, * character #
#************************************************************#
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 2
Explanation of input parameters
-v increase verbosity: print state vectors after initialization and at end of run
-l lattice type for initialization: hypercubic or (hyper)fcc (default: fcc)
hypercubic: atoms on cartesic lattice in D dimensions (linear, square, simple cubic, . . . )
(hyper)fcc: atoms on every second site of hypercubic lattice (even taxi-cab distance from origin)
-s system size: integer extent of base lattice in each direction
number of atoms is sD for hypercubic, sD/2 for fcc
determines finite-size effects
-r density (in natural units)
1st physical parameter; density and number of atoms unit distance
-e energy per particle
2nd physical parameter; enforced at initialization, but nearly constant during simulation
-c cutoff of potential and forces (in natural units)
determines accuracy of Hamiltonian; must be smaller than half linear box size
-d time step (in natural time units)
determines accuracy of MD integration
-t total simulation time (in natural time units, not MD steps!)
determines “statistical” error of averages; must be much larger than autocorrelation times
-w warm-up time (in natural time units, not MD steps!)
determines convergency error; must be larger than “melting” or reordering times
-b number of bins in pair correlation function
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 3
Explanation of output (default parameters)
When started, the program first prints a header identifying itself (version, date, author, and dimension
from compiler directives):
WS2006_Simulation/HW5_MD> ./MD_LJ | more
#************************************************************#
# MD_LJ: molecular dynamics program for Lennard-Jones atoms #
# Version: 0.6, 22.12.2006 by Nils Bluemer, compiled for D=3 #
# #
# For help on usage and available options: run MD_LJ -h #
#************************************************************#
It then states all parameters used in the simulation:
# numparts = 256, densisty = 0.600000, length = 7.528288
# twarmup = 0.200000, simtime = 2.000000, timestep = 0.002000
# E_N = 0.500000, cutoff = 2.50, corrbins = 200, lattice = f
Then follows (for usual run, i.e. without the option -v) the primary output: all observables are printed
at each time step in the simulation (negative times indicate warmup)
# time Epot/N Ekin/N Etot/N temperature pressure
-0.2000 -4.404346 4.904346 0.500000 3.269564 -2.266062
-0.1980 -4.404603 4.904603 0.500000 3.269735 -2.265327
-0.1960 -4.405370 4.905369 0.499999 3.270246 -2.263113
...
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 4
.Finally, the pair correlation function is printed (based on all positions observed after warm-up in the
simulation):
## # distance paircorr(distance)
## 0.009410 0.000000
## 0.028231 0.000000
...
## 0.818701 0.000000
## 0.837522 0.000862
## 0.856343 0.006825
## 0.875164 0.031667
## 0.893984 0.126551
...
In verbose mode (option -v), the full states [number of particle, position (D coordinates), velocity (D
coordinates), force=acceleration (D coordinates)] are given directly after initialization and at the end
of the simulation:
# time Epot/N Ekin/N Etot/N temperature pressure
-0.2000 -4.404346 4.904346 0.500000 3.269564 -2.266062
1 0.471 0.471 0.471 -1.094 -2.233 2.277 -0.000 -0.000 -0.000
2 2.353 0.471 0.471 2.312 0.220 -0.951 -0.000 -0.000 -0.000
3 4.235 0.471 0.471 -1.595 -0.542 1.868 0.000 -0.000 -0.000
...
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 5
First analysis of program output (default parameters)
Step 1: run simulation and save output
WS2006_Simulation/HW5_MD> ./MD_LJ > results_MD_LJ_default.dat
Step 2: analyze time-dependent observables
All non-comment lines show observables as a function of time. Let’s first check the total energy:
Etot stays constant
within about 10−3
Etot = 0.49958 ±
0.00001
Significant autocorrelation time (after warm-up): 7 time steps, i.e. τ ≈ 0.014 in natural units
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 6
How were the plot and the numbers obtained?
The plot was generated using gnuplot (http://www.gnuplot.info/):
WS2006_Simulation/HW5_MD> gnuplot results_MD_LJ_energy.gnu
with the input file
set term post eps enh color solid 12
set out "results_MD_LJ_energy.eps"
set size 0.7,0.5
set xlabel "t"
set ylabel "E_{tot}"
set nokey
plot [-0.2:2] "results_MD_LJ_default.dat" u 1:4 w lp lw 2 pt 7 ps 0.4
Here, the important part is plot [-0.2:2] ”results MD LJ default.dat” u 1:4 (u as abbreviation for
“using” 1st and 4th column are plotted)
The data was analyzed using the tool awk (http://www.gnu.org/software/gawk/gawk.html) and
“our” usual statistics tool:
HW5_MD> grep -v "#" results_MD_LJ_default.dat | awk ’{if ($1>=0) print $4}’ | stats -a
Average: 0.49959043, variance: 1.4889776e-08, error: 3.679152e-06
Korrelation time: 8.3545387, corrected error: 1.0634302e-05, transient: -8.4977902e-08
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 7
Now, let’s look at all energy observables:
Explain the lead-in behavior, in particular the initially negative pressure!
Is the warm-up time well chosen?
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 8
Code of program MD LJ
Program header
#define PROGNAME "MD LJ"
#define VERSION "0.6"
#define DATE "22.12.2006"
#define AUTHOR "Nils Bluemer"
/* project: molecular dynamics program for Lennard-Jones system */
/* new: 0.3: cutoff, full interface
0.4: pair correlation function
0.5: pressure, warm-up with negative times, tail corrections
for potential and pressure, code clean-up */
0.6: fcc lattice (alternative to hypercube), numlin->systemsize
/* todo: velocity distribution,
binning for pair correlation function (-> error analysis),
extra files for correlation function, velocity distribution etc.,
make dimension regular variable, implement D-general volume,
make computation of Epot more efficient (together with force)
next major version: neighbor lists */
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 9
#include
#include
#include
#include
#include
#include "pointer utils.c" /* vector-definitions */
#define BUFSIZE 100 /* used for reading command line */
#define EPSILON 0.000001 /* used for end time in time loop */
#define DIMENSION 3 /* todo: introduce dimension variable */
#define DEF verbosity 0
#define DEF systemsize 8
#define DEF simtime 2.0
#define DEF twarmup 0.2
#define DEF timestep 0.002
#define DEF density 0.6
#define DEF E N 0.5
#define DEF cutoff 2.5
#define DEF corrbins 200
#define DEF lattice ’f’
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 10
main program
int main (int argc, char *argv[]) {
char c, lattice;
double **position, **velocity, **force;
unsigned long *paircorr, ncollect;
double time,simtime,timestep,twarmup;
double density, length, length2, cutoff, cutoffsq;
double Epot, energy, E N, corrfac, virial;
unsigned long numparts;
unsigned int systemsize, corrbins;
unsigned int verbosity;
/* default values for parameters */
verbosity = DEF verbosity;
systemsize = DEF systemsize;
simtime = DEF simtime;
twarmup = DEF twarmup;
timestep = DEF timestep;
density = DEF density;
E N = DEF E N;
cutoff = DEF cutoff;
corrbins = DEF corrbins;
lattice = DEF lattice;
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 11
/* parse command line (read in parameters) */
while (--argc > 0 && (*++argv)[0] == ’-’)
while (c= *++argv[0])
switch (c) {
case ’h’: printhelp(); exit(0); break;
case ’n’: sscanf(++argv[0],"%ld\n",&numparts); error("change of numparts
not yet implemented"); break;
case ’r’: sscanf(++argv[0],"%lf\n",&density); break;
case ’e’: sscanf(++argv[0],"%lf\n",&E N); break;
case ’t’: sscanf(++argv[0],"%lf\n",&simtime); break;
case ’d’: sscanf(++argv[0],"%lf\n",×tep); break;
case ’w’: sscanf(++argv[0],"%lf\n",&twarmup); break;
case ’c’: sscanf(++argv[0],"%lf\n",&cutoff); break;
case ’s’: sscanf(++argv[0],"%d\n",&systemsize); break;
case ’l’: sscanf(++argv[0],"%c\n",&lattice); break;
case ’v’: verbosity++; break;
}
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 12
/* compute derived parameters */
if (lattice==’h’)
numparts=powD(systemsize);
else { /* lattice == ’f */
if (systemsize%2==1)
error (”fcc lattice only for even systemsize”);
numparts = powD(systemsize)/2;
}
length=pow(numparts/density,1.0/DIMENSION);
length2=0.5*length;
energy=E N*numparts;
cutoffsq=cutoff*cutoff;
if (cutoff>length2)
error(”box too small for cutoff; decrease density or increase particle number”);
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 13
/* first output */
print header ();
printf (”# numparts = %ld, densisty = %f, length = %f\n”, numparts, density, length);
printf (”# twarmup = %f, simtime = %f, timestep = %f\n”, twarmup, simtime, timestep);
printf (”# E N = %f, cutoff = %4.2f, corrbins = %d, lattice = %c\n”, E N, cutoff, corrbins,
lattice);
printf (”#\n”);
print obs step (numparts,velocity,virial,Epot,length,time,cutoff,density,0);
/* set up pointers, init binned observables */
position = dmatrix(1,numparts,1,DIMENSION);
velocity = dmatrix(1,numparts,1,DIMENSION);
force = dmatrix(1,numparts,1,DIMENSION);
paircorr = lvector(0,corrbins); /* component 0 holds number of bins */
init paircorr (paircorr,corrbins);
/* initialize positions and variables */
init positions(numparts,position,length,systemsize,lattice);
Epot=Epot fun(numparts,position,length2,cutoff,density);
if (energy negative Ekin”);
init velocities(numparts,velocity,energy-Epot);
ncollect=0;
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 14
/* initialize simulation */
time=-twarmup;
compute forces(numparts,position,force,length2,cutoffsq,&virial,paircorr,time);
/* compute and print observables for initial state */
print obs step (numparts,velocity,virial,Epot,length,time,cutoff,density,1);
if (verbosity>0)
print pos vel force(numparts,position,velocity,force);
/* start of simulation loop */
for (time=time+timestep;time<=simtime+EPSILON;time+=timestep){
if(time>=0)
ncollect++; /* number of measurements (e.g. for paircorr) */
/* advance system by one time step */
velocity verlet(numparts,position,velocity,force,timestep,length,cutoffsq,
&virial,paircorr,time);
Epot=Epot fun(numparts,position,length2,cutoff,density);
print obs step (numparts,velocity,virial,Epot,length,time,cutoff,density,1);
} /* end of simulation loop */
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 15
/* print collected observables (e.g. pair correlation) */
if (verbosity>0)
print pos vel force(numparts,position,velocity,force);
corrfac=2.0/3.0*M PI*ncollect*numparts*density*powD(length2);
if (DIMENSION!=3)
warning(”incorrect scale in pair correlation function!”);
print paircorr(paircorr,length2,corrfac);
/* finish up */
freevectors(position,velocity,force,paircorr,numparts);
return(0);
}
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 16
double powD(double x) /* computes xˆDIMENSION */
{
unsigned int d;
double powD;
powD=x;
for (d=2;d<=DIMENSION;d++)
powD*=x;
return(powD);
}
/* potential enters here!!!!!!! */
double pairpot (double r2)
{
double rinv6;
rinv6=1.0/(r2*r2*r2);
return (4.0*rinv6*(rinv6-1.0));
}
double forcefac (double r2)
{
double rinv6;
rinv6=1.0/(r2*r2*r2);
return (48.0*rinv6*(0.5-rinv6)/r2);
}
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 17
double Epot fun (unsigned long numparts, double **position, double length2,
double cutoff, double density)
{
unsigned long n,m;
double distsq,sum,*vec,cutoffsq,cutoffc,offset;
cutoffsq=cutoff*cutoff;
cutoffc=cutoffsq*cutoff;
sum=0.0;
vec = dvector(1,DIMENSION);
offset=pairpot(cutoffsq);
for (n=1;n<=numparts;n++)
for (m=n+1;m<=numparts;m++){
comp rel(length2,position[n],position[m],vec,&distsq);
if (distsqlength2)
dx-=2.0*length2;
relvec[d]=dx;
*distsq+=dx*dx;
}
}
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 19
double Ekin fun (int numparts, double **velocity)
{
unsigned int d,n;
double x;
x=0.0;
for (n=1;n<=numparts;n++)
for (d=1;d<=DIMENSION;d++)
x+=velocity[n][d]*velocity[n][d];
return(0.5*x);
}
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 20
/* place atoms on regular grid */
void init positions (int numparts, double **position, double length, int systemsize, char lattice) {
int n,m,k,d, *ivec, fac, sum;
if (lattice ==’h’)
fac=1;
else /* lattice == ’f’ */
fac=2;
n=0;
ivec=ivector(1,DIMENSION);
for (m=0;m0))
update paircorr(paircorr,distsq/length2sq);
}
free dvector(vec,1,DIMENSION);
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 26
void velocity verlet (unsigned long numparts, double **position, double **velocity,
double **force, double timestep, double length,
double cutoffsq, double *virial, unsigned long *paircorr,
double time)
{
int n,d;
double ts2;
ts2=timestep*timestep;
for (n=1;n<=numparts;n++)
for (d=1;d<=DIMENSION;d++){
position[n][d]+=velocity[n][d]*timestep+0.5*force[n][d]*ts2;
if (position[n][d]<0)
position[n][d]+=length;
else if (position[n][d]>length)
position[n][d]-=length;
velocity[n][d]+=0.5*timestep*force[n][d];
}
compute forces (numparts,position,force,0.5*length,cutoffsq,virial,paircorr,time);
for (n=1;n<=numparts;n++)
for (d=1;d<=DIMENSION;d++)
velocity[n][d]+=0.5*timestep*force[n][d];
}
Computer simulations in statistical physics - HW 5 · WS 2006/07 · Nils Blu¨mer (Univ. Mainz) C ←↩ 4 B 27
void printvector (double *vec)
{
unsigned int d;
for (d=1;d<=DIMENSION;d++)
printf(”%6.3lf ”,vec[d]);
}
void print pos vel force (unsigned long numparts, double **position,
double **velocity, double **force)
{
unsigned long n;
for (n=1;n<=numparts;n++){
printf(”%ld ”,n);
printvector(position[n]);
printf(” ”);
printvector(velocity[n]);
printf(” ”);
printvector(force[n]);
printf(”\n”);
}