Comprehensive Nuclear Materials - Modelling Section (R. Stoller, Editor)
Case Studies for Chapter 128 Molecular Dynamics
Wei Cai, Ju Li and Sidney Yip
This tutorial describes how to reproduce the case studies on crystals and dislocations described in Comprehensive Nuclear Materials, Chapter 128 Molecular Dynamics. We provide the MD++ script files and describe how to use them in detail.
Contents |
Download files
First, you need to install MD++ on your computer by following the instructions on MD++ Manuals.
Second, copy files ta_perfect_crystal.tcl and ta_edge_dislocation.tcl to your input file directory of MD++. You can do so by the following commands (assuming you have installed MD++ in ~/Codes/MD++).
export MDPP=~/Codes/MD++
mkdir -p ${MDPP}/scripts/work/cnm_chapter
cd ${MDPP}/scripts/work/cnm_chapter
wget http://micro.stanford.edu/mediawiki/images/0/08/Ta_perfect_crystal.tcl.txt -O ta_perfect_crystal.tcl
wget http://micro.stanford.edu/mediawiki/images/4/4a/Ta_edge_dislocation.tcl.txt -O ta_edge_dislocation.tcl
Both files contain a comment section in the beginning that describes in detail the various options to reproduce various steps of the case studies. In this tutorial, we describe how to run these scripts and what kind of results you should expect.
Compile executable file
Different executables in MD++ contains different interatomic potentials, for example,
sw - Stillinger-Weber Silicon (modified by Balamane et al.) sworig - Original version of SW Silicon swge - Stillinger-Weber Germanium tersoff - Tersoff potential for Silicon meam-lammps - MEAM (taken from lammps) meam-baskes - MEAM (taken from Baskes's code dynamo) meam - MEAM (directly implemented in MD++) eam - Embedded Atom Method fs - Finnis-Sinclair
In this tutorial, we compile the Finnis-Sinclair potential model fs by typing
cd ${MDPP}
make fs build=R SYS=gpp
If compilation is successful, this will create binary file fs_gpp file in your ${MDPP}/bin directory. If you would like to use the intel compiler, you should use SYS=intel in the make command. To compile MD++ on a Mac, use SYS=mac.
The parameters for the fs potential for several BCC metals are specified by the potential files in the ${MDPP}/potentials folder, e.g. ta_pot for Ta.
Perfect Crystal
Zero Temperature Properties
To create a perfect crystal of BCC Ta and compute its equilibrium lattice constant a0, cohesive energy Ecoh and bulk modulus B, run the following command.
cd ${MDPP}
bin/fs_gpp scripts/work/cnm_chap/ta_perfect_crystal.tcl 0
If you would like to visualize the crystal during this (short) simulation, uncomment the following lines (i.e. remove the # sign) in ta_perfect_crystal.tcl file.
#openwindow #MD++ sleep
The data files of this simulation are saved in the ${MDPP}/runs/ta-perf directory. To analyze the result, use the following commands.
cd ${MDPP}/runs/ta-perf
wget http://micro.stanford.edu/mediawiki/images/2/26/Fit_a0EB.m.txt -O fit_a0EB.m
octave
fit_a0EB('EPOT_a_coarse.dat',1,2)
fit_a0EB('EPOT_a_fine.dat', 1,2)
Octave should plot the potential energy of the crystal as a function of the lattice constant, and print out the fitted value for a0, Ecoh and B. When using the EPOT_a_coarse.dat file, the fitted values are: 
, Ecoh = − 8.1000eV and B = 197.19 GPa. When using the EPOT_a_fine.dat file, the fitted values are: , Ecoh = − 8.1000eV and B = 196.14 GPa. (This reproduces Fig.5 of the chapter.)
You can also run fit_a0EB.m in Matlab if it is installed on your computer.
To compute the other (cubic) elastic constants of BCC Ta, run the following command.
cd ${MDPP}
bin/fs_gpp scripts/work/cnm_chap/ta_perfect_crystal.tcl 1
The data are saved in the files C11_C12.dat, C44.dat in the ${MDPP}/runs/ta-perf folder. To analyze the result, use the following commands.
cd ${MDPP}/runs/ta-perf
wget http://micro.stanford.edu/mediawiki/images/1/11/Fit_C11C12C44.m.txt -O fit_C11C12C44.m
octave
fit_C11C12C44('C11_C12.dat','C44.dat')
Octave should plot the stress as a function of strain and print out the fitted value of the elastic constants: C11 = 266.00 GPa, C12 = 161.20 GPa, C44 = 82.41 GPa. (This reproduces Fig.6 of the chapter.)
Finite Temperature Properties
To perform MD simulations in the NVE ensemble (Fig.7a of the chapter), run the following command.
cd ${MDPP}
bin/fs_gpp scripts/work/cnm_chap/ta_perfect_crystal.tcl 2
The data are stored in the prop.out file in the ${MDPP}/runs/ta-perf folder. The format of this property file is specified in the output_fmt line in the ta_perfect_crystal.tcl file. For example, the first four columns are simulation step number, potential energy (in eV), kinetic energy (in eV) and instantaneous temperature (in K), respectively. To plot the data, do the following.
cd ${MDPP}/runs/ta-perf
octave
data = load('prop.out');
figure(1); plot(data(:,1)*1e-3, data(:,4));
figure(2); plot(data(:,1)*1e-3, mean(data(:,7:9)*1e-3,2));
Figure 1 plots instantaneous temperature (in K) as a function of time (in ps) and figure 2 plots the instantaneous pressure (in GPa) as a function of time (reproducing Fig.7a of the chapter).
To perform a series of MD simulations in the NVT ensemble and adjust the box sizes iteratively to relax the stress (Fig.7b of the chapter), run the following command.
cd ${MDPP}
bin/fs_gpp scripts/work/cnm_chap/ta_perfect_crystal.tcl 4
The instantaneous temperature and stress data are still stored in the prop.out file. The averaged pressure versus box size data for each MD simulation are stored in the P_L.dat file in the ${MDPP}/runs/ta-perf folder. To plot the data, do the following.
cd ${MDPP}/runs/ta-perf
octave
data = load('prop.out');
figure(1); plot(data(:,1)*1e-3, data(:,4));
figure(2); plot(data(:,1)*1e-3, mean(data(:,7:9)*1e-3,2));
PL = load('P_L.dat');
figure(3); plot(PL(:,1),PL(:,2));
Figure 1 plots instantaneous temperature (in K) as a function of time (in ps) and figure 2 plots the instantaneous pressure (in GPa) as a function of time (reproducing Fig.7b of the chapter). Figure 3 plots the averaged pressure during each MD simulation in the series.
Dislocation
Peierls Stress at Zero Temperature
To create an edge dislocation in a BCC Ta crystal and relax the structure to a local energy minimum, run the following command.
cd ${MDPP}
bin/fs_gpp scripts/work/cnm_chap/ta_edge_dislocation.tcl 0
An X-window will automatically appear, displaying the dislocation core atoms (in blue) and the surface atoms (in grey). If you would like to turn off the X-window (e.g. if this job is submitted to a queue), comment out (i.e. add the # sign in front of) the following line in ta_edge_dislocation.tcl file.
openwindow
The data files are stored in the ${MDPP}/runs/ta-edge-0 folder. File A.log stores the print-outs during the simulation. File edge_relaxed.cn contains the atomic positions of the relaxed structure. The format of .cn files is explained in [Output File Formats in MD++]. File edge_relaxed.cfg stores the dislocation core and surface atoms in Atomeye's cfg format and can be visualized by the following command (reproducing Fig.8b of the chapter).
cd ${MDPP}/runs/ta-edge-0
${MDPP}/Tools/AtomEye/A.i686 edge_relaxed.cfg
In the Atomeye window, press Alt-0 to color the atoms according to their central symmetry parameters.
To compute the Peierls stress, run the following command.
cd ${MDPP}
bin/fs_gpp scripts/work/cnm_chap/ta_edge_dislocation.tcl 1
The dislocation core atoms will be displayed in an X-window as the crystal is relaxed under increasing amount of shear stresses. The potential energies of the relaxed structures are saved in the EPOT.dat file. When the Peierls stress is exceeded, the dislocation moves over large distances during the relaxation and the potential energy decreases abruptly. The relaxation also takes a much longer time to terminate compared with the case when the applied stress is below the Peierls stress. All these can be used as indications for having reached the Peierls stress. Dislocation motion can also be directly observed in the X-window. The atomic configurations of relaxed structures at different applied stresses can be examined using Atomeye by the following commands, e.g.,
cd ${MDPP}/runs/ta-edge-0
${MDPP}/Tools/AtomEye/A.i686 ta-stress20.cfg
This plots the relaxed configuration under σxy = 20 MPa.
The Peierls stress predicted by this simulation is 28 MPa.
Mobility at Finite Temperature
To perform MD simulations, we first need to equilibrate the simulation cell at T = 300 K using the following command.
cd ${MDPP}
bin/fs_gpp scripts/work/cnm_chap/ta_edge_dislocation.tcl 2
This takes about 2 hours to finish. After this, we are ready to simulate dislocation motion at finite stress, using the following command, e.g.,
cd ${MDPP}
bin/fs_gpp scripts/work/cnm_chap/ta_edge_dislocation.tcl 4 100
This set the shear stress to be σxy = 100 MPa. The simulation data are stored in the ${MDPP}/runs/ta-edge-100 folder. The auto####.cfg files are the simulation snapshots that can be visualized by Atomeye. The inter####.cn files are the restart files containing all atomic positions and velocities.
Here is a PBS script for submitting MD simulations jobs under different stresses to a parallel cluster (mc2.stanford.edu).
The dislocation core position as a function of time and the average dislocation velocity can be obtained by running the find_disl_vel.m file in Octave or Matlab. Use the following commands.
cd ${MDPP}/runs/ta-edge-0
wget http://micro.stanford.edu/mediawiki/images/d/d6/Find_disl_vel.m.txt -O find_disl_vel.m
octave
filepath = '../../runs/ta-edge-100';
find_disl_vel
Modify the filepath line to analyze data from MD simulation under different stresses. This reproduces Fig.9 of the chapter.