MDS

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”); }