Comprehensive Nuclear Materials MD Case Studies

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.

---

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 $a_{0}$, cohesive energy $E_{\rm {coh}}$ 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 $a_{0}$, $E_{\rm {coh}}$ and $B$. When using the EPOT_a_coarse.dat file, the fitted values are: Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a_0 = 3.3058 \AA} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle E_{\rm coh} = -8.1000} eV and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle B = 197.19} GPa. When using the EPOT_a_fine.dat file, the fitted values are: Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a_0 = 3.3058 \AA} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle E_{\rm coh} = -8.1000} eV and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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 0

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: Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{11} = 266.00} GPa, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{12} = 161.20} GPa, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{44} = 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.