M03 Equilibrium Lattice Constant and Bulk Modulus

Manual 03 for MD++
Computing Equilibrium Lattice Constant

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 file A.log in the directory specified by dirname. For example, if you run mo.script in Manual 02 with setnolog commented out, then you will find a A.log file in the runs/mo-example directory. To read this file, type

$ more runs/mo-example/A.log

In the log file, the line starting with ASSIGN shows that a variable is assigned to the values as specified in the script file. 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 file 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 specific property (e.g. EPOT) in the log file, 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 file 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”.

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 different 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 define the lattice energy $\Phi$ and the number density 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 \rho} as

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 \Phi = \frac{E_{pot}}{N} } .

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 \rho = \frac{N}{V} } .

where 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_{pot}} is the potential energy of the crystal, 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 N} is total number of atoms in the simulation cell (corresponding to variable NP in MD++) 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 V} is the volume of the simulation cell. Run MD++ with following command line to use the Stillinger-Weber (SW) potential model.

$ 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 file, you can find number of atoms NP for the 3 X 3 X 3 DC cell to be 216. This number can also be obtained by calculating 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 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 different lattice constant can also be read from the log file by typing

$ grep EPOT runs/si_polytype/A.log

The volume of a simulation cell 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 V} can be obtained from the determinant of the matrix 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 \mathbf{H}} . When 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 \mathbf{H}} is a diagonal matrix (the same is true for an upper triangular matrix),

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 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 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 \Phi} (in eV) of silicon as a function of the number density 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 \rho} (in 1/ų), for different crystal structures, as shown in Fig.1. The structure with the lowest energy is DC. The equilibrium lattice constant 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} corresponds to the number density that gives the minimum of the 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 \Phi(\rho)} curve. The minimum of 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 \Phi} is also called the cohesive energy 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_{coh}} . For DC silicon, 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=} 5.431 Å 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 E_{coh}} = −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 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 \Phi(\rho)} curve near the minimum also tells us the bulk modulus of the crystal. The bulk modulus 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} is defined as

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 = V \frac{\partial^2 \Phi}{\partial V^2}_{V=V_0} } .

where 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 V_0} is the equilibrium atomic volume (corresponding to the energy minimum). To compute the second derivative, plot 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 \Phi} as a function of 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 V = 1/\rho} and fit the 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 \Phi(V)} curve by a quadratic function, i.e. 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 \Phi = 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 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 V_0=} 20.025 (Å⁶), the bulk modulus becomes

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 = 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. In the elasticity theory, the bulk modulus 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} is related to the other elastic constants through

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 = \lambda + \frac{2}{3}\mu }	for isotropic material
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 = \frac{1}{3}(C_{11} + 2C_{12})}	for cubic material

where 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 \lambda} 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 \mu} are Lamé's constants 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 C_{11}} 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 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 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}=} 165.7(GPa) 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 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. Different parameterization of the potential model generally leads to different 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')