MD++ Tutorial
An Example Script File
Adrian Buganza, William Kuykendall and Wei Cai
This tutorial walks you through a complete .script file that performs an MD simulation.
To run a simulation using MD++ the basic approach consists on writing a script file. The script file consists of a series of variable assignments (e.g. latticeconst = 5.430) and commands (e.g. openwin) that MD++ uses to setup and run the simulation. The disadvantage of script files is that they do not allow for more complex programming such as for loops or similar programming tools. To incorporate this type of features the user can write a *.tcl file. More information on *.tcl files can be found in a different manual.
Contents |
Overall Structure of *.script files
In general, a script file will have the following parts:
- MD setup
- Initialize atoms positions
- Plotting setup
- Relaxation
- MD simulation settings
- MD Simulation
The last three items do not always appear because they actually depend on the specific problem to be solved. For example, it is a common practice to relax the structure (find a configuration with minimum potential energy) before running a simulation, but the initial configuration (e.g. read in from a file) might already satisfy such condition. Furthermore, there are some cases in which simulations at zero temperature are done in order to determine characteristic properties of a material like the ideal strength, in such a case, time integration is not needed and the last item on the list won't appear in the code. In summary, what operations are performed by MD++ depend on the commands listed in the script file. Nonetheless, in the following, we will briefly explain all the items listed above.
Suppose you have installed MD++ in your $HOME/Codes/MD++ folder.
export MDPP=$HOME/Codes/MD++ cd $MDPP
In this example, we will use the follwing script file:
$MDPP/scripts/Examples/example02a-si-md.script
All characters following the # (pound) sign are comments. They will be ignored by the MD++ program and we will not discuss them here either.
Output set up
Nearly every script file starts with the following three lines.
setnolog setoverwrite dirname = runs/si-example
The first two commands (setnolog and setoverwrite) are optional, but the line that specifies dirname must precede all other commands. Otherwise MD++ will exit with error.
The first line, setnolog indicates that everything will be printed to screen and the log file (A.log) created will be empty. If this line is commented out (by adding # in front), MD++ will redirect all screen outputs to the A.log file in the directory specified by dirname.
Command setoverwrite allows the new output files (including A.log file) to overwrite the old files in the directory (specified by dirname). If this line is commented out, MD++ will exit with error if the directory specified by dirname already exists.
Finally, dirname = runs/si-example specifies the output directory. If it does not exist, it will be created.
Initialize atoms positions
Before doing any simulations (e.g. commands eval, relax, run) or visualization (e.g. command plot), an atomic configuration must be loaded (by readcn) or created (e.g. by makecrystal) into MD++.
For the former use the following, e.g.
incnfile = relaxed.cn readcn
This is an easy way to read in a configuration (.cn) file you have created before. This reads in a file (relaxed.cn from the runs/si-example folder) that contains the positions of your atoms. Of course, if this is the first time you run this script file, your runs/si-example folder will not contain any .cn files. So this line will not work. You will have to comment it out and try to create an atomic configuration file from scratch.
A common way to create an atomic configuration is to create a perfect crystal structure, e.g. using the following commands.
crystalstructure = diamond-cubic
latticeconst = 5.4309529817532409
latticesize = [ 1 0 0 4
0 1 0 4
0 0 1 4 ]
makecrystal writecn
crystalstructure = diamond-cubic specifies the crystal structure. In this example, it specifies diamond-cubic, which is the crystal structure for Silicon. latticeconst = 5.4309529817532409 specifies the lattice constant (a0) in Angstrom. The lines latticesize = ... specifies the size of the lattice you wish to create. The first three values 1 0 0 specify the orientation of the first simulation cell vector (in Miller indices notation), followed by the number of repeats (4) in that direction. The next three values 0 1 0 specify the orientation of the second simulation cell vector, followed by the number of repeats (4) in that direction. The next three values 0 0 1 specify the orientation of the third simulation cell vector, followed by the number of repeats (4) in that direction.
For example, 1 0 0 4 will produce the crystal structure as 2 0 0 2. The atoms will be labeled differently, but the structure will be the same.
To emphasize, the crystalstructure and latticeconst determine the structure. The latticesize specifies the size and orientation of the crystal to be created.
It is also worth mentioning that Periodic Boundary Conditions (PBC) are always implied in MD++. PBC are always used. However, you can create gaps such that the atoms in the periodic images are outside the cutoff radius. This will effectively create a free surface. This is called the supercell approach.
When creating a crystal structuree, MD++ never creates half-atoms or fractional atoms (as might be seen in some solid state physics books). A simple cubic will have a single full atom at one corner (other corners have zero atoms). The other 7 corners will also have a single atom because of PBC.
Finally, makecrystal is the command that actually creates the crystal structure. It first computes the number of atoms to be allocated and then specifies their coordinates. The command writecn writes the configuration to a file in the .cn format. By default, the file name is final.cn, but you can change it by specifying variable, e.g. finalcnfile = perf.cn before calling writecn.
The .cn file starts with the number of atoms, then followed by that number of lines, each line giving the scaled coordinates of one atom. If other options are set (e.g. writevelocity = 1 or writeall = 1), it may also save other things, such as velocities, in each line.
Plotting set up
This section if pretty standard, so we will be brief. The usual lines of code you can have in your .script file are:
atomradius = 0.67 bondradius = 0.3 bondlength = 2.8285 atomcolor = orange highlightcolor = purple bondcolor = red backgroundcolor = gray70 plotfreq = 10 rotateangles = [ 0 0 0 1.25 ] openwin alloccolors rotate saverot eval plot
The first two lines are parameters for plotting. atomradius specifies how large a ball you want to draw for each atom. bondradius = 0.3 specifies how thick a line you want to draw for each bond. bondradius = 0.3 means any two atoms whose distance is less than 0.3 Angstrom will have a bond drawn between them.
plotfreq = 10 specifies that during your (future) MD simulation, the atomic structure in the X-window will be updated every 10 MD time steps. rotateangles specifies the HOME orientation (the first three values) and scaling (the fourth value) of your view. By default, your HOME viewing orientation has the x-axis pointing to the right, y-axis pointing upward, and z-axis pointing out of the plane. After rotating your viewing angle by mouse drag or arrow keys, you can return to this HOME viewing angle by pressing the HOME key.
openwin is the command that opens a new X-window (i.e. the canvas) where atoms will be plotted. You are better off if you always put the other commands alloccolors rotate saverot right after it, which help initialize the X-window.
The eval command is not a plotting command, but it is usually used before plotting the results. For example, sometimes we want to color atoms according to their local potential energy, which requires the potential function to be called before plotting. We can do so by calling eval. The last command plot puts the atoms onto the X-window. You can rotate your viewing angle by mouse drag. Press F1 in the X-window to print a list of commands you can use to change your view or save a snapshot.
Relaxation
Relaxation means moving the atomic coordinates to a local minimum of the potential energy function. This is done with the conjugate gradient (CG) algorithm (an iterative search algorithm). To find the relaxed configuration use the following lines.
conj_ftol = 1e-7 conj_itmax = 1000 conj_fevalmax = 10000 conj_fixbox = 1 conj_summary = 1 relax finalcnfile = relaxed.cn writecn
The first line sets up the function tolerance, maximum iterations and maximum number of times the code is allowed to call the potential function. In the second line conj_fixbox = 1, indicates that the simulation box (i.e. the 
) is not allowed to change. If this parameter is set to 0, the box size will go to the equilibrium size. Thus, for however many repeat cells you have, the box will go to the equilibrium size for that number of repeat cells. This is a convenient way to determine the equilibrium lattice constant of your crystal structure.
After the various variables for the CGR algorithm have been set in the first two lines, the relax calls the CG algorithm to relax the structure. The relaxed structure is saved into a configuration file (specified as relaxed.cn here).
MD simulation set up
In this part of the code you define parameters such as ensemble type or integrator to be used in your subsequent MD simulation. Generic lines of code that might be useful are:
T_OBJ = 300 equilsteps = 0 totalsteps = 100 timestep = 0.0001 atommass = 28.0855 atomTcpl = 200.0 boxTcpl = 20.0 DOUBLE_T = 0 srand48bytime initvelocity saveprop = 0 savepropfreq = 10 openpropfile ensemble_type = NVE
T_OBJ = 300 sets the target temperature to 300 Kelvin, which will be used to initialize velocities. (It will also be used by the thermostat if the NVT ensemble were used, though not in this example.) The second line sets total number of steps and the time step (in picosecond). atommass = 28.0855 specifies atomic mass in grams/mole. atomTcpl = 200.0 boxTcpl = 20.0 are variables pertaining to the thermostat and bariostat for certain ensembles (not needed for NVE).
In the next line, if the variable DOUBLE_T is set to 0 then MD++ will not double T_OBJ when initializing the velocity. When equilibrium is reached for an NVE simulation for solids, the temperature will usually be about half of your initial T_OBJ. If T_OBJ is set to 1, MD++ will double T_OBJ when initializing the velocity, so that at equilibrium the temperature gets close to the desired value T_OBJ.
srand48bytime seeds the random number generator by the computer's clock time, so that everytime you run the simulation, your atomic trajectories will be different (because initvelocity uses pseudo-random numbers when assigning velocities). For debugging purposes, you may want your simulations to be exactly repeatable. You can do so by always assigning the same random seed, e.g.
randseed = 12345 srand48
If saveprop = 0, the property file will not be saved during your MD simulation. If it were set to 1, a property file will be saved. savepropfreq = 10 specifies a line will be added to the property file every 10 MD time steps (if saveprop = 1). Command openpropfile opens the property file to prepare for subsequent MD simulation.
Lastly, ensemble_type = NVE defines the simulation ensemble. Typical ensembles are NVE and NVT. integrator_type = VVerlet sets the integrator type to be Velocity Verlet. You can also try Gear6 to see the difference.
Suppose you do want to output simulation properties into a file, then you can do so by setting, e.g.
saveprop = 1 savepropfreq = 10 openpropfile output_fmt = "curstep EPOT KATOM Tinst" outpropfile = thermo.out
These lines instruct MD++ to print out a line into the thermo.out file every 10 simulation steps. Each line will contain four entries: current simulation step, potential energy (in eV), kinetic energy (in eV), and instantaneous temperature (in K).
MD simulation
To run an MD simulation there is only one command needed:
run
The parameters of this MD simulation has been defined in the previous section. It is a good idea to put the command
closepropfile
before exiting the simulation. This makes sure that the output file on your hard disk gets updated before the simulation terminates.
The quit command instructs MD++ to exit without reading the rest of the script file. MD++ will also exit when it reaches the end of the script file. Sometimes, you want to view the atomic structure at the end of your simulation before MD++ exits. You can do so by putting the main program to sleep (the X-window is still awake) before exiting,
sleep quit
The command sleep puts the main program to sleep for 600 seconds. You can change the default value by setting, e.g. sleepseconds = 30 .