22.53 Statistical Processes and Atomistic Simulations

Greetings, class of 22.53

Lecture A:

  • Q/A about the Haile code (a better looking version)
  • Reduced units and dimensional analysis
  • Read section I, II and appendices of the manuscript for background information on
  • An experiment deals with one property, but an MD simulation deals with everything
  • Highly accurate methods exist (LDA, tight binding)/expensive
  • Semi-empirical potentials are less reliable/much cheaper

    Lecture B:

    The new MD code (version 1.1) for Argon features

  • con file input interface
  • LJ6-12 or Woon's ab initio pair potential (Chem. Phys. Lett. 204, 29, 1993)
  • NEV or NTP (Parrinello-Rahman) ensemble
  • g(r) and static structure factor S(k)
  • micro-canonical ensemble fluctuation formulas for
  • self-diffusion coefficient in liquids
  • Green-Kubo/spectral method to calculate the thermal conductivity
  • same for the shear viscosity

    For details on thermal conductivity calculations read this paper (text and figures) and this manuscript.

    To install the program, download 3D.tar to your home directory and then

    Ignore error messages like "rm: smile: No such file or directory". If make does not work on your system, try gmake. The executable after compilation is argon. It reads running conditions from standard input, so will feed an input file con to argon. The merit of doing that is that one can edit con at leisure before committing it to the run, and it is easier for later revision. con file has the format of questions and answers; you are supposed to only edit the answers. All questions are self-explanatory, and you are strongly encouraged to explore and revise the code.

    Details include:

  • The simulation cell can be triclinic, the three sides specified by the three rows of a certain H matrix

  • Particle positions are stored in dimensionless s, ranging from 0 and 1, related to the Cartesian coordinates by (x,y,z) = (s1, s2, s3) H

  • The MD configuration is always saved at the end of one run (config.woon.110), which can be used by later runs as starting points

  • If starting from scratch, the initial structure is cubic fcc cells (owning 4 atoms each) extended in three dimensions according to NC(1,2,3). The initial configuration is given twice as much the necessary kinetic energy to compensate for the potential energy

  • TAUT specifies how "gentle" the temperature rescaling should be; theoretically, the actual temperature converges to the desired temperature exponentially with that time constant

  • In NTP mode, H matrix is incorporated into the Hamiltonian and is treated on equal footing as particle coordinates. Thus it also has a "mass" and a "temperature". The larger the wall "mass", the slower it responds to stress inequilibrium

  • Short-term averages are printed on the screen

  • After kstart steps for the initial transients and equilibration to finish, the code starts to collect long-term averages that will appear at the end of temp.out as final results

  • If NEH mode is selected, then the total energy should be conserved to at least 5 significant digits after the temperature rescaling is turned off
  • If NTP mode is selected, temperature rescaling is never turned off except if thermal conductivity is demanded to be calculated: in that case, NTP is automatically switched to NEH mode after kstart, using the average H of previous steps

    After making sure the con file is correct, you can submit it to the run. Try a short one (con_liquid_short) for 110K liquid, make run is in "interactive mode" where one does not leave the terminal. It should finish in 5-10 minutes, and you may find the following files under your directory,

  • temp.out: the numeric part is the instantaneous conditions during the run - watch out for total energy conservation after kstart; the text part is a summary of long-term averaged properties

  • gr.out: pair correlation function g(r), which can be plotted out using grs.m in Matlab

  • struc.out: it records the time evolution of the H matrix, static structure factor etc.; at present it is of no great value but is useful in the future for studying solid-state instabilities

  • config: MD configuration at the end of the run

  • Data/: all outputs are automatically dated and backupped in that directory

    What the above short run cannot provide are the transport coefficients: thermal conductivity and shear viscosity, which can only be accurately calculated with longer runs (1-2 hrs). If you could access private machines, try the following (con_liquid_long), If everything turns out as planned (see Makefile), you should be able to log out of that machine and let the job run in background. And it will automatically time the job, save screen output to a file output, and notify you by e-mail as soon as it is done. To check if it is really running in the background, log in the same computer and use When it is finished, you should find additional files under your directory:

  • jx.out,jy.out, jz.out: heat current fluctuations

  • freq.out: heat current power spectrum

  • corr.out: heat current correlation function

  • sx.out,sy.out, sz.out: shear stress fluctuations

  • shear_freq.out: shear stress power spectrum

  • shear_corr.out: shear stress correlation function

    You do not have to fully understand these files. The final results will be clearly presented in temp.out. But in case you are interested, look up the thermal conductivity papers referenced above and play with the Matlab programs heat_spectra.m and heat_corr.m, which plot out the data.

    Our results compare well with the published liquid simulation results of Ermakova et al, J. Chem. Phys. 102, 4942 (1995), using the Woon's potential. The slight differences should be due to details in implementation, such as how the potential cutoff is treated. In argon, for instance, we do not carry out the analytical integration for the tail contribution since some properties such as thermal conductivity cannot be similarly treated.

    Now, how about some solid-state calculations? After all, the code's original purpose is to study these properties in great accuracy, especially the thermal conductivity. Try a 30K solid (con_solid_NEV) running under NEV mode, It turns out in the end that the average pressure for this condition is approximately zero. That is of course not a coincidence. One can determine the equilibrium (zero-pressure) volume of a structure for a given temperature by trying out many different volumes (lattice constants) under NEH mode and see which leads to P=0. But a much simpler way is to run the simulation under NTP mode with P=0, starting from some arbitrary volume (con_solid_NTP); the volume relaxation is then done automatically, and property averaging whose results appear in temp.out only begins after P is approximately zero.

    A general Parrinello-Rahman scheme not only relaxes the volume, but also the shape of the simulation box, if a non-hydrostatic constant external stress is applied. Our code is able to do that after a tiny revision, but is not implemented here for aesthetic reasons, since under non-hydrostatic constant external stress one cannot define a thermodynamic potential for the system as in the hydrostatic case, and there is no Hamiltonian dynamics with conservation laws.

    Food for thought:

  • The micro-canonical fluctuation formulas provide us a means to calculate the heat capacities and thermal expansion coefficient in one shot. But how can one check independently that these results are actually correct, that we/others have not made a mistake in those ghastly derivations?

  • In the same light, how can one calculate the individual elastic constant, like C11, C12 and C44, through their definitions?

  • The Verlet list algorithm in argon 1.1 take O(N2) effort to maintain, where N is the total number of particles. If N is large, say 100,000, which one often encounters in dislocation and fracture simulations, the cost is going to be prohibitive. Can anyone come up with a good O(N) design?

  • How can one improve the user-interface of this code, which has also to be very robust and portable?

    Fun Links and References:

  • Modeling from a physicist's view, a chemist's view, a material scientist's view
  • Moldy
  • MDS
  • PMD
  • Rasmol
  • Simulated Annealing
  • MIT 3.320 (Spring)
  • NIST Physical Constants
  • Web Elements
  • Numerical methods
  • People, research and links in the Multiscale Materials Modeling Group
  • Order-N MD project
  • Astronomy Picture of the Day
    Email: [email protected]. This page has been accessed at least several times since Oct.6, 1998.