Manual 05 for MD++
MD++ Tour with Sample Scripts
Keonwook Kang and Wei Cai
[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
![]() |