Manual 03 for MD++
Equilibrium Lattice Constant and Bulk Modulus

Keonwook Kang and Wei Cai

Nov 26 , 2007

Reading a Log File

In the example scripts of Manual 02, the command setnolog was activated and the simulation results were printed out to the standard output, which is usually your terminal screen. If you want to keep the record of this information, you comment out the setnolog command by adding # before it. This will create a log ﬁle A.log in the directory speciﬁed by dirname. For example, if you run mo.script in Manual 02 with setnolog commented out, then you will ﬁnd a A.log ﬁle in the runs/mo-example directory. To read this ﬁle, type

$ more runs/mo-example/A.log

In the log ﬁle, the line starting with ASSIGN shows that a variable is assigned to the values as speciﬁed in the script ﬁle. The line with EXEC shows the execution of certain command. After the eval command is executed, properties such as number of atoms (NP), potential energy (EPOT), kinetic energy (KATOM), pressure (PRESSURE) and stress (Stress) are also printed. Sometimes, the log ﬁle may be zipped as A.log.gz, in which case you can use

$ gzip -cd runs/mo-example/A.log.gz | more

to see the content. To search for a speciﬁc property (e.g. EPOT) in the log ﬁle, you may use

$ grep EPOT runs/mo-example/A.log

or

$ gzip -cd runs/mo-example/A.log.gz | grep EPOT

depending on whether the log ﬁle is zipped or not. In case you want to “grep” several lines before and after your keyword (e.g. Stress), you can use

$ grep -3 Stress runs/mo-example/A.log

This will show three lines above and below any line containing the keyword “Stress”.

Obviously, this way of extracting the values of MD++ variables is very tedious. There is a better way to do this, which is to use the Tcl scripting language, as described in Manual 08 (will be available soon). The Tcl script for MD++ calculations described in this section is attaced at the end of this manual.

Equilibrium Lattice Constant and Cohesive Energy

Under ambient condition silicon takes the diamond-cubic (DC) structure. But in a computer simulation, we can create Si crystals with diﬀerent hypothet- ical structures, such as face-centered-cubic (FCC), body-centered-cubic (BCC), or simple-cubic (SC). A good potential model should be able to tell us that the DC structure is the one with the lowest energy, hence it is the most favorable structure for silicon. In this section, we will perform such calculations using MD++.

Let us deﬁne the lattice energy $Φ$ and the number density $ρ$ as

$\Phi = \frac{E_{pot}}{N}$ .

$\rho = \frac{N}{V}$ .

where $E p o t$ is the potential energy of the crystal, $N$ is total number of atoms in the simulation cell (corresponding to variable NP in MD++) and $V$ is the volume of the simulation cell. Run MD++ with following command line to use the Stillinger-Weber (SW) potential model<ref name="SW">F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985).</ref><ref name="BH">H. Balamane, T. Halicioglu, and W. A. Tiler, Phys. Rev. B 46, 2250 (1992)</ref>.

$ bin/sw_gpp scripts/ME346/si_polytype.script

Here is the content of the si_polytype.script.

# -*-shell-script-*-
#setnolog
setoverwrite
dirname = runs/si_polytype
#------------------------------------------------------------
#Create Perfect Lattice Configuration
#
crystalstructure = diamond-cubic
latticeconst = 5.4309529817532409 #(A) for Si
latticesize = [ 1 0 0 3
                0 1 0 3
                0 0 1 3]
makecrystal eval
latticeconst = 4.850 makecrystal eval
latticeconst = 4.950 makecrystal eval
latticeconst = 5.050 makecrystal eval
                :
(many lines omitted here for brevity)
                :
latticeconst = 5.900 makecrystal eval
latticeconst = 6.000 makecrystal eval
latticeconst = 6.100 makecrystal eval
crystalstructure = face-centered-cubic
latticeconst = 4.105 makecrystal eval
latticeconst = 4.110 makecrystal eval
                :
                :
latticeconst = 4.205 makecrystal eval
latticeconst = 4.215 makecrystal eval
crystalstructure = body-centered-cubic
latticeconst = 3.210 makecrystal eval
latticeconst = 3.220 makecrystal eval
                :
                :
latticeconst = 3.320 makecrystal eval
latticeconst = 3.340 makecrystal eval
crystalstructure = simple-cubic
latticeconst = 2.550 makecrystal eval
latticeconst = 2.600 makecrystal eval
                :
                :
latticeconst = 2.640 makecrystal eval
latticeconst = 2.650 makecrystal eval
quit

From the log ﬁle, you can ﬁnd number of atoms NP for the 3 X 3 X 3 DC cell to be 216. This number can also be obtained by calculating $8 \times 3^3 = 216$ , since there are eight atoms in the DC unit cell. For other crystal structures, the number of atoms in the 3X3X3 cell is given in the following table.

Table 1. Number of atoms in different crystal structures
Crystal Structure	No. of atoms in the unit cell	Total No. of atoms
SC	1	27
BCC	2	54
FCC	4	108
DC	8	216

The potential energy at each diﬀerent lattice constant can also be read from the log ﬁle by typing

$ grep EPOT runs/si_polytype/A.log

The volume of a simulation cell $V$ can be obtained from the determinant of the matrix $\mathbf{H}$ .<ref>The matrix $\mathbf{H}$ deﬁnes size and shape of the simulation cell. $\mathbf{H} = [\mathbf{c}_1 \; \mathbf{c}_2 \; \mathbf{c}_3]$ where $\mathbf{c}_i$ 's are three periodicity vectors. In MD++, the matrix $\mathbf{H}$ becomes a diagonal matrix when the supercell is a rectangular box (after reorienting the coordinate system). </ref> When $\mathbf{H}$ is a diagonal matrix (the same is true for an upper triangular matrix),

$V = \det(\mathbf{H}) = \mathbf{H}(1,1) \mathbf{H}(2,2) \mathbf{H}(3,3)$ .

the volume is the product of the entries in the main diagonal of $\mathbf{H}$ . The unit of length in MD++ is Å and the unit of volume is Å³ . From these, we can calculate the lattice energy $Φ$ (in eV) of silicon as a function of the number density $ρ$ (in 1/Å³), for diﬀerent crystal structures, as shown in Fig.1. The structure with the lowest energy is DC. The equilibrium lattice constant $a 0$ corresponds to the number density that gives the minimum of the $Φ(ρ)$ curve. The minimum of $Φ$ is also called the cohesive energy $E c o h$ . For DC silicon, $a 0 =$ 5.431 Å and $E c o h$ = −4.63 eV.

Fig.1 Lattice energy of Si in different structures. DC Si has the lowest energy well. The inset figure is the zoom-in image of the energy curve of DC Si.

Fig.2 Lattice energy vs lattice constant of DC Si. The curve is fitted to the quadratic polynomials. The equilibrium lattice constant is estimated to be 5.43095(Å).

Bulk Modulus

The curvature of the $Φ(ρ)$ curve near the minimum also tells us the bulk modulus of the crystal. The bulk modulus $B$ is deﬁned as<ref>This is because $B = -V \partial P / \partial V$ and $P = - \partial \Phi / \partial V$ .</ref>

$B = V \frac{\partial^2 \Phi}{\partial V^2}_{V=V_0}$ .

where $V 0$ is the equilibrium atomic volume (corresponding to the energy minimum). To compute the second derivative, plot $Φ$ as a function of $V = 1 / ρ$ and ﬁt the $Φ(V)$ curve by a quadratic function, i.e. $Φ = c 0 + c 1 V + c 2 V 2$ , as shown in Fig.3. This can be done by the Matlab command polyfit and the result is c₂ = 0.016879 (in unit of eV/Å⁶ ). Then, with $V 0 =$ 20.023 (Å³), the bulk modulus becomes

$B = 2\times 20.023 \mathrm{(\AA^3)}\times 0.016879 \mathrm{(eV/\AA^6)}\times 160.2 \mathrm{(GPa/(eV/\AA^3))}= 108.3 \mathrm{(GPa)}$ .

This result matches with an earlier report in the literature<ref name="BH"/>. In the elasticity theory<ref>J. P. Hirth and J. Lothe, Theory of Dislocations, (Wiley, New York, 1982)</ref>, the bulk modulus $B$ is related to the other elastic constants through

$B = \lambda + \frac{2}{3}\mu$	for isotropic material
$B = \frac{1}{3}(C_{11} + 2C_{12})$	for cubic material

where $λ$ and $μ$ are Lamé's constants and $C 11$ and $C 12$ are cubic elastic constants. From this relation, the bulk modulus would be 98.4 GPa if we calculate it from the experimental value of $C 11 =$ 165.7(GPa) and $C 12 =$ 63.9(GPa) for silicon. The discrepancy between the simulation and experimental results in the bulk modulus is partly because the simulation results corresponds to the ideal case of T = 0 K, while the experimental result is obtained at room temperature T = 300 K. The elastic constants generally decreases with increasing temperature. The discrepancy may also come from the fact that the potential model has empirical nature. Diﬀerent parameterization of the potential model generally leads to diﬀerent predictions of the elastic constants.

Fig.3 Lattice energy vs atomic volume of DC Si. The equilibrium atomic volume is estimated to be 20.023(Å³).

Tcl version of MD++ script

For anyone that wants to reproduce my data, I uploaded the Tcl version of MD++ script si_polytype.tcl and the matlab script for the post-process of the simulation data.

si_polytype.tcl

# -*-shell-script-*-
#*******************************************
# Definition of procedures
#*******************************************
proc initmd { } { MD++ {
#setnolog
setoverwrite
dirname = runs/si_polytype
zipfiles = 1   # zip output files
NIC = 200 NNM = 200
}}

proc makecrystal { } {
#--------------------------------------------
# Create Perfect Lattice Configuration
  MD++ element0 = Si
  MD++ latticesize = \[ 1 0 0 4 0 1 0 4 0 0 1 4 \]
  MD++ makecrystal
} 

proc readnwrite { a fileID } {
#--------------------------------------------
# read and write properties
for {set i -50} {$i <= 50} {incr i} {
  set latt [expr $a + $i/10000.0 ]
  MD++ latticeconst = $latt
  makecrystal

  MD++ eval    
  set EPOT     [MD++_Get EPOT]
  set N        [MD++_Get NP]
  set VOL      [MD++_Get OMEGA]
  set ELAT     [expr $EPOT / $N]
  set rho      [expr $N / double($VOL)]
  set AtomVol  [expr double($VOL) / $N]

  puts $fileID "[format %18.12e $latt]\t \
                [format %18.12e $rho] \t \
                [format %18.12e $AtomVol]\t \
                [format %18.12e $EPOT]\t \
                [format %18.12e $ELAT]"
}}

#*******************************************
# Main program starts here
#*******************************************
initmd

MD++ crystalstructure = diamond-cubic
set a 5.4309 
set fileID [open "EvsVol_DC.dat" w]
readnwrite $a $fileID
close $fileID

MD++ crystalstructure = face-centered-cubic
set a 4.1465 
set fileID [open "EvsVol_FCC.dat" w]
readnwrite $a $fileID
close $fileID

MD++ crystalstructure = body-centered-cubic
set a 3.245 
set fileID [open "EvsVol_BCC.dat" w]
readnwrite $a $fileID
close $fileID

MD++ crystalstructure = simple-cubic
set a 2.6120 
set fileID [open "EvsVol_SC.dat" w]
readnwrite $a $fileID
close $fileID

MD++ quit

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab script to plot the lattice energy of Si
% in different structures (DC, FCC, BCC and SC)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf

DC.data = load('EvsVol_DC.dat');
DC.latt = DC.data(:,1);
DC.rho  = DC.data(:,2);
DC.atomvol = DC.data(:,3);
DC.EPOT = DC.data(:,4);
DC.ELAT = DC.data(:,5);

FCC.data = load('EvsVol_FCC.dat');
FCC.latt = FCC.data(:,1);
FCC.rho  = FCC.data(:,2);
FCC.atomvol = FCC.data(:,3);
FCC.EPOT = FCC.data(:,4);
FCC.ELAT = FCC.data(:,5);

BCC.data = load('EvsVol_BCC.dat');
BCC.latt = BCC.data(:,1);
BCC.rho  = BCC.data(:,2);
BCC.atomvol = BCC.data(:,3);
BCC.EPOT = BCC.data(:,4);
BCC.ELAT = BCC.data(:,5);

SC.data = load('EvsVol_SC.dat');
SC.latt = SC.data(:,1);
SC.rho  = SC.data(:,2);
SC.atomvol = SC.data(:,3);
SC.EPOT = SC.data(:,4);
SC.ELAT = SC.data(:,5);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot the lattice energy vs number density
% with inset figure
h1 = figure(1);
set(gca,'fontsize',14)
plot(DC.rho, DC.ELAT, '.-',FCC.rho, FCC.ELAT, '.-', ...
     BCC.rho, BCC.ELAT, '.-',SC.rho, SC.ELAT, '.-')
xlabel('Number Density, \rho=N/V (1/Angstrom^3)','fontsize',16), 
ylabel('Binding Energy per atom (eV)','fontsize',16)
legend('DC','FCC','BCC','SC',2)
h2 = get(h1,'CurrentAxes')
h3 = axes('pos',[.5 .2 .35 .35]);

hold on; plot(DC.rho, DC.ELAT, '.-')
yticklong = get(h3,'YTick');
set(h3,'xlim',[0.0498 0.0501],'YTickLabel',yticklong,'box','on');
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot the lattice energy vs lattice constant
% and find the equilibrium lattice constant, a_0
[P, S] = polyfit(DC.latt,DC.ELAT,2);
latt_X = linspace(5.4305,5.4315);
ELAT_Y = polyval(P,latt_X);

figure(2)
set(gca,'fontsize',14)
plot(DC.latt, DC.ELAT, '.-',latt_X,ELAT_Y,'.-r')
xlabel('Lattice Constant, \it a \rm (Angstrom)','fontsize',16), 
ylabel('Binding Energy per atom (eV)','fontsize',16)
yticklong = get(gca,'YTick');
set(gca,'YTickLabel',yticklong,'Xlim',[5.426 5.436]);
legend('data','fitting curve')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot the lattice energy vs atomic volume
% and find the bulk modulus, B
[P2, S2] = polyfit(DC.atomvol,DC.ELAT,2);
vol_X = linspace(20.01,20.04);
ELAT_Y2 = polyval(P2,vol_X);

figure(3)
set(gca,'fontsize',14)
plot(DC.atomvol, DC.ELAT, '.-',vol_X,ELAT_Y2,'.-r')
xlabel('Atomic Volume, \it V \rm (Angstrom^3)','fontsize',16), 
ylabel('Binding Energy per atom (eV)','fontsize',16)
yticklong = get(gca,'YTick');
set(gca,'YTickLabel',yticklong);
legend('data','fitting curve')

LAMMPS script

For LAMMPS users, I present the LAMMPS script, in.equillatt_bulk_mod to carry out identical computation using the original SW potential model.

in.equillatt_bulk_mod

log log.Si_test
# log file name

shell rm out.data
variable latt0 equal 5.43
variable max_iter equal 101
variable shift_i  equal ceil(${max_iter}/2.0)
label loop 
variable i loop ${max_iter}
  if "$i > ${max_iter}" then "jump in.equillatt_bulk_mod"
  variable j equal $i-${shift_i}
  variable latt equal ${latt0}+$j/10000.0
  
  clear
  units metal
  dimension 3
  atom_style atomic
  boundary p p p
  newton on

  # number of cpus in x, y, and z dimensions
  processors 1 1 1

  lattice diamond ${latt} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
  region  box block 0 4 0 4 0 4 units lattice
  create_box 1 box
  create_atoms 1 region box

  mass 1 28.085
  
  # Stillinger and Weber, PRB 31, 5262 (1985)
  pair_style sw
  pair_coeff * * ../../potentials/Si.sw Si

  neighbor 0.31 bin
  neigh_modify delay 0 every 1 check yes

  thermo 10
  thermo_style custom step lx ly lz press pxx pyy pzz pe temp
  dump    1 all custom 10 Si_DC.dump id type xs ys zs
  dump_modify 1 format "%d %d %25.17E %25.17E %25.17E"

  variable EPOT equal pe
  variable NP equal atoms
  variable ELAT equal pe/${NP}
  variable VOL equal vol
  variable rho equal ${NP}/${VOL}
  variable AtomVol equal ${VOL}/${NP}

  min_style cg
  minimize 0.0 1e-25 1000 1000

  shell echo ${latt} ${rho} ${AtomVol} ${EPOT} ${ELAT} >> out.data 
  next i
  jump in.equillatt_bulk_mod loop
label break
variable i delete

Notes

Manual 03 for MD++ Equilibrium Lattice Constant and Bulk Modulus

Contents