Manual 04 for MD++
An Example Script File
Adrian Buganza, William Kuykendall and Wei Cai
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.in each line are an axis 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.in each line are an axis 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 setup
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<tt>, indicates that the simulation box (i.e. the <mathbf>H</mathbf>) 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 <tt>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 setup
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\rm \, sets the target temperature to 300 Kelvin; that will be used to initialize velocities. The second line sets total number of steps and the time step as a fraction of pico-seconds. Command \tt atommass = 28.0855 \rm specifies atomic mass in grams/mole. \\ Commands \tt atomTcpl = 200.0 boxTcpl = 20.0 \rm are variables to help control the thermostat and bariostat for certain ensembles. \\ In the nest line, If the variable \tt DOUBLE\_T \rm is set to 0 then MD++ will not double T\_OBJ and thus at equilibrium, T will be about T\_OBJ/2; if it is set to 1 MD++ will double T\_OBJ at the start of the simulation such that at equilibrium the temperature gets to the desired value T\_OBJ. \\ In the sixth line, \tt srand48bytime \rm seeds the random number generator by the computer's clock time and the next line initializes velocities of all the atoms. \\ Then, if \tt saveprop = 0 \rm, the property file will not be saved, if it is equal to 1, it will be saved; \tt savepropfreq \rm is how many time steps before you save the property and \tt openpropfile \rm opens the property file.\\ Lastly, \tt ensemble\_type = NVE \rm defines the ensemble. Typical ensembles are NVE and NVT. \tt integrator\_type = VVerlet \rm sets the integrator type. In this case, the integrator is chosen to be the \emph{Velocity Verlet} algorithm, but \emph{Gear6} can also be used.
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.
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 quitThe command sleep puts the main program to sleep for 600 seconds. You can change the default value by setting, e.g. sleepseconds = 30 .