M03 Equilibrium Lattice Constant and Bulk Modulus: Difference between revisions

From Micro and Nano Mechanics Group
Jump to navigation Jump to search
(Fig.2 added)
No edit summary
 
(11 intermediate revisions by 3 users not shown)
Line 1: Line 1:
<H1 ALIGN="CENTER"><FONT SIZE="-1">Manual 03 for MD++ </FONT>
<H1 ALIGN="CENTER"><FONT SIZE="-1">Manual 03 for MD++ </FONT>
<BR>
<BR>
Computing Equilibrium Lattice Constant</H1>
Equilibrium Lattice Constant and Bulk Modulus</H1>
<DIV>
<DIV>


Line 49: Line 49:
“Stress”.
“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 [[M08_MD++_Powered_by_Tcl_Language |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 ==
== Equilibrium Lattice Constant and Cohesive Energy ==
Line 70: Line 71:
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>.
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
$ bin/sw_gpp scripts/ME346/si_polytype.script


Here is the content of the '''si_polytype.script'''.
Here is the content of the '''si_polytype.script'''.
Line 153: Line 154:
|}
|}
the volume is the product of the entries in the main diagonal of \mathbf{H}. The unit
the volume is the product of the entries in the main diagonal of <math>\mathbf{H}</math>. The unit
of length in MD++ is &Aring; and the unit of volume is &Aring;<sup>3</sup> . From these, we can
of length in MD++ is &Aring; and the unit of volume is &Aring;<sup>3</sup> . From these, we can
calculate the lattice energy <math>\Phi</math> (in eV) of silicon as a function of the number
calculate the lattice energy <math>\Phi</math> (in eV) of silicon as a function of the number
Line 175: Line 176:
where <math>V_0</math> is the equilibrium atomic volume (corresponding to the energy minimum).
where <math>V_0</math> is the equilibrium atomic volume (corresponding to the energy minimum).
To compute the second derivative, plot <math>\Phi</math> as a function of <math>V = 1/\rho</math> and
To compute the second derivative, plot <math>\Phi</math> as a function of <math>V = 1/\rho</math> and
fit the <math>\Phi(V)</math> curve in the neighborhood of the minimum by a quadratic function,
fit the <math>\Phi(V)</math> curve by a quadratic function,
''i.e.'' <math>\Phi = c_0 + c_1 V + c_2 V^2</math>, as shown in Fig.1(b). This can be done by the Matlab
''i.e.'' <math>\Phi = c_0 + c_1 V + c_2 V^2</math>, as shown in Fig.3. This can be done by the Matlab
command '''polyfit''' and the result is c<sub>2</sub> = 0.0169 (in unit of eV/&Aring;<sup>6</sup> ).
command '''polyfit''' and the result is c<sub>2</sub> = 0.016879 (in unit of eV/&Aring;<sup>6</sup> ).
Then, with <math>V_0=</math>20.025 (&Aring;<sup>6</sup>), the bulk modulus becomes
Then, with <math>V_0=</math>20.023 (&Aring;<sup>3</sup>), the bulk modulus becomes


{|border="0" align="center"
{|border="0" align="center"
|<math> B = 2\times 20.025 \mathrm{(\AA^3)}\times 0.0169 \mathrm{(eV/\AA^6)}\times 160.2 \mathrm{(GPa/(eV/\AA^3))}= 108.5 \mathrm{(GPa)}</math>.
|<math> 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)}</math>.
|}
|}


Line 201: Line 202:
parameterization of the potential model generally leads to different predictions
parameterization of the potential model generally leads to different predictions
of the elastic constants.
of the elastic constants.
[[Image:E_vol_Si.jpg | frame | center | Fig.3 Lattice energy vs atomic volume of DC Si. The equilibrium atomic volume is estimated to be 20.023(&Aring;<sup>3</sup>).]]

== 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
<pre>
# -*-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
</pre>

<pre>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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')
</pre>

== 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'''
<pre>
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
</pre>


== Notes ==
== Notes ==

Latest revision as of 12:29, 3 July 2013

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 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”.

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 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 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} 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, is total number of atoms in the simulation cell (corresponding to variable NP in MD++) and is the volume of the simulation cell. Run MD++ with following command line to use the Stillinger-Weber (SW) potential model[1][2]. 

$ 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 , 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}} .[3] 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 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}} . The unit of length in MD++ is Å and the unit of volume is Å3 . 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/Å3), 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[4]

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 c2 = 0.016879 (in unit of eV/Å6 ). 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.023 (Å3), 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[2]. In the elasticity theory[5], 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(Å3).

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

  1. ↑ F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985).
  2. ↑ 2.0 2.1 H. Balamane, T. Halicioglu, and W. A. Tiler, Phys. Rev. B 46, 2250 (1992)
  3. ↑ 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}} defines size and shape of the 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 \mathbf{H} = [\mathbf{c}_1 \; \mathbf{c}_2 \; \mathbf{c}_3]} 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 \mathbf{c}_i} 's are three periodicity vectors. In MD++, 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}} becomes a diagonal matrix when the supercell is a rectangular box (after reorienting the coordinate system).
  4. ↑ This is because 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 \partial P / \partial V } 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 P = - \partial \Phi / \partial V} .
  5. ↑ J. P. Hirth and J. Lothe, Theory of Dislocations, (Wiley, New York, 1982)
Retrieved from "http://micro.stanford.edu/mediawiki/index.php?title=M03_Equilibrium_Lattice_Constant_and_Bulk_Modulus&oldid=5921"