Manual 04 for MD++
An Example Script File

Apr 23 , 2012

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.

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

This can be done in two ways, by importing position data from a \emph{*.cn} file or by using the command \tt{makecrystal} \rm. For the former use the following:\\ \\ \fcolorbox{black}{gray!20}{\parbox{\linewidth}{ \tt{incnfile = relaxed.cn readcn} }}\\ \\ This is an easy way to get a setup you've already used. This reads in a file ("relaxed.cn") to set up your crystal. \\There is also \tt{saveintercn} \rm. This will save an intermediate configuration file. It also has a frequency variable \tt{plotfreq} \rm . \\ To create a perfect lattice configuration the following commands can be used\\ \\ \fcolorbox{black}{gray!20}{\parbox{\linewidth}{ \tt{crystalstructure = diamond-cubic \\ latticeconst = 5.4309529817532409 \begin{tabbing} latticesize = [\= 1 0 0 4\\ \> 0 1 0 4\\ \> 0 0 1 4 ] \end{tabbing} makecrystal\\ writecn } }}\\ \\ Starting from the first line, \tt{crystalstructure = diamond-cubic} \rm sets up the crystal structure to make a new crystal. \tt{latticeconst = 5.4309529817532409} \rm sets up the lattice constant for the crystal you want to make. Then \tt{latticesize =...} \rm Creates a lattice. The first three values in each line are an axis direction, while the fourth number is the number of repeats of this vector that exists. For example, [1 0 0 4] will produce the same thing as [2 0 0 2]. The atoms will be labeled differently, but the structure will be the same.\\ Again, to emphasize, the crystal structure and lattice constant determine the structure. The lattice size will change the size and orientation of the structure. \\ It is also worth commenting that Periodic Boundary Conditions (PBC) are implied in MD++. Everything is PBC. However, you can create gaps such that the next atoms are outside the cutoff radius. This will effectively create a free surface, for example.\\ Moreover, this code does not create half-atoms or fractional atoms. A simple cubic will have a single full atom at one corner. Then, using PBC, the other 7 corners will also have a single atom because of the repetition.\\ Finally, \tt{makecrystal} \rm is the command that actually creates the crystal structure and the command \tt{writecn} \rm writes the configuration to a file. Every line in the file is an atom. If other options are set, it may also save other things, such as velocities.\\

Plotting setup

This section if pretty standard, so we will be brief. The usual lines of code you can have in your code are:\\ \fcolorbox{black}{gray!20}{\parbox{\linewidth}{ \tt{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. The \tt bondlegth \rm parameter depends on the atom being modeled. The third line, \tt openwin \rm Opens the x window to display the configuration. The following line allocates the colors according to your specifications (atom colors, highlighting colors, bond colors, etc., then rotates the cell and saves the rotation as the default rotation. \\ The \tt eval \rm command is not a plotting command, but it is usually used before plotting the results. This command evaluates the total potential energy.\\ The final line, \tt plot \rm is the actual plotting command.

Relaxation

Relaxation consists of finding a minimum of the potential energy. This is done with the Conjugate Gradient iteration. To find the relaxed configuration use the following lines of code: \fcolorbox{black}{gray!20}{\parbox{\linewidth}{ \tt{ 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 \tt conj\_fixbox = 1\rm , indicates that the box is not changing. 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. \\ After the various variables for the CGR algorithm have been set in the first two lines, the third line relaxes the crystal, then it saves the relaxed configuration file.

MD simulation setup

In this part of the code you define parameters such as ensemble type or integrator to be used. Generic lines of code that might be useful are:\\ \fcolorbox{black}{gray!20}{\parbox{\linewidth}{ \tt{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 openpropfil \\ integrator\_type = 1\, ensemble\_type = NVE} }}\\ \\ \tt 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 a simulation there is only one command needed:\\ \fcolorbox{black}{gray!20}{\parbox{\linewidth}{ \tt{ run } }}\\ \\ To put MD++ to sleep and then quit the simulation use:\\ \fcolorbox{black}{gray!20}{\parbox{\linewidth}{ \tt{ sleep quit }