Manual 05 for MD++
MD++ Tour with Sample Scripts

[First written, [ Jan 19 ]], 2007

[Last Modified, [ Feb 25 ]], 2008

Finite Temperature Simulation

In Manual 02, we learned how to create a perfect BCC, FCC or DC crystal by the MD++ command makecrystal, by which we assigned the initial position to each of the atoms. For nonzero finite temperature simulation or molecular dynamics (MD) simulation, we need to additionally give the initial velocities and set parameters such as numerical integrator, time step and so on. In this manual, we aim to be familiar with MD++ commands and variables related with MD simulations. MD simulations can be categorized as microcanonical (NVE), canonical (NVT), isoenthalpic-isobaric (NPH) or constant pressure/temperature (NPT) ensemble type according to which physical quantities are expected to be constant or controlled. For example, a system of microcanonical ensemble has the total energy (E), the system volume (V) and the number of atoms (N) to be conserved during MD simulation. Basically, you need to choose a proper ensemble set for your own purpose of simulation.

Microcanonical(NVE) Ensemble

It would be a good start explaining MD simulation of NVE ensemble to understand what you do with MD simulations. Imagine a system of N-particle gas inside of a rigid container which is well insulated that the system is isolated from the environment. There is no heat and mass transfer through the container wall. Then you expect that the number of particles, the container volume and the total energy will not change. In MD simulations, you do not see a visible wall but you may think of a closed system with fictitious walls.

Now I give the example script mo_NVE.tcl below, written in Tcl version, which shows how to run typical NVE MD simulation at 300 K from the perfect BCC molybdenum (Mo) structure. You can run the script by typing

$ bin/fs_gpp scripts/mo_NVE.tcl

At this point, you may raise a question how we can say that the temperature is 300 K in NVE ensemble in which there is no temperature we control. It is true. The temperature in NVE ensemble is used to generate atoms' initial velocities and that does not mean the temperature will be maintained at 300 K during the simulation. This will be covered more detail later.

# -*-shell-script -*-
# run NVE MD simulation of perfect crystal of Mo.
#
#*******************************************
# Definition of procedures
#*******************************************
proc initmd { } { MD++ {
#setnolog
setoverwrite
dirname = runs/mo-example
zipfiles = 1   # zip output files
#
#--------------------------------------------
# Create Perfect Crystal
element0 = Mo
crystalstructure = body-centered-cubic
latticeconst = 3.1472           # in Angstrom for Mo
latticesize = [ 1 0 0 5
                0 1 0 5
                0 0 1 5]
}}
#------------------------------------------------------------
proc readMoPot { } { MD++ {
# Read the potential file
potfile = ~/Codes/MD++/potentials/mo_pot readpot
}}
#-------------------------------------------------------------
proc openwindow { } { MD++ {
# Plot Configuration
atomradius = 1.0 bondradius = 0.3 bondlength = 2.725
atomcolor = orange highlightcolor = purple
bondcolor = red backgroundcolor = gray70
plotfreq = 10 rotateangles = [ 0 0 0 1.25 ]
openwin alloccolors rotate saverot eval plot
}}
#--------------------------------------------
proc setup_md { } { MD++ {
#MD settings
T_OBJ = 300        # (in Kelvin) Desired Temperature
atommass = 95.94   # (in g/mol)
timestep = 0.001   # (in ps)
totalsteps = 2000
saveprop = 1 savepropfreq = 10
savecn = 0 savecnfreq = 100
writeall = 1
DOUBLE_T = 1 randseed = 12345 srand48
ensemble_type = "NVE" integrator_type = "VVerlet"

#*******************************************
# Main program starts here
#*******************************************
initmd
readMoPot
MD++ makecrystal finalcnfile = perf.cn writecn
openwindow
#---------------------------------------------
# run MD
setup_md
MD++ initvelocity finalcnfile = init.cn writecn
MD++ {output_fmt = "curstep EPOT KATOM Tinst"}
MD++ outpropfile = thermo.out openpropfile
MD++ run closepropfile
MD++ finalcnfile = 300K-5X5X5.cn writecn eval
MD++ sleep quit
}}

I will skip some of commands in the above script that are already explained in previous manuals. zipfiles = 1 will compress output files. For example, the configuration file init.cn and the thermodynamic data file thermo.out will be zipped as init.cn.gz and thermo.out.gz, respectively. If we do not specify outputpropfile, the default file name is prop.out. The format of output file is determined by output_fmt and every line of the file thermo.out.gz has current step (curstep), potential energy (EPOT), kinetic energy (KATOM), and temperature (Tinst) as defined in output_fmt. If we do not specify output_fmt, the default array of thermodynamic data output is current step (curstep), kinetic energy (KATOM), potential energy (EPOT), pressure (PRESSURE), −σxy (TSTRESS_xy), −σyz (TSTRESS_yz), −σzx (TSTRESS_zx), HELM, HELMP, thermodynamic frction coefficient (zeta), dEdlambda and system volume (OMEGA). Units are given as eV for energy and eV/A˚3 for stress.

Atom species can be specifed like element0 = Mo. We define a Tcl function readMoPot to read the interatomic potential of Mo. Commands and variables related with visualization in openwindow will be explained in a separate manual in detail. A configuration .cn file always has information of number of atoms (NP) in the first line and 3 by 3 matrix H in the last three lines. Inbetween, writeall=1 stores NP lines of the following information: scaled coordinates, sacaled veloc- ities, individual potential energy, whether-or-not fixed, topology, atom species and atom groups. Without declearation of writeall=1, only scaled coordinates will be saved at each line. You can check this differences by comparing two con- figuration files perf.cn and init.cn. Because perf.cn is written before setting writeall = 1, it has only number of atoms, scaled coordinates of each atom and the box matrix H. Alternatively, writevelocity = 1 can be used, which stores scaled coordinates and velocities.

• Fig.1 Visualization of FCC structure
• 

  (a) Conventional image of FCC unit cell

• 

  (b) FCC unit cell in MD++. One corner atom and three face-centered atoms per cell

• 

  (c) FCC of 3X3X3 cells

.

Notes