A Polygonal Dislocation Loop in MD++

Manual ?? for MD++
Creating a Polygonal Dislocation Loop

Keonwook Kang and Wei Cai

Jan 19 , 2009

This short summary explains how we create a polygonal dislocation loop in the atomistic simulation, using the displacement field of a triagonal dislocation loop.

---

Displacement Field of a Triagonal Dislocation Loop

Fig.1.

While studying a dislocation nucleation problem by atomistic simulations, the authors needed to create a polygonal dislocation loop, which requires the knowledge of the displacement field of a polygonal dislocation loop in the infinite isotropic medium.

In elasticity, the displacement field of a polygonal dislocation loop is merely the sum of fields of triangular loops that constitute the polygon loop. The coordinate-free form of analytic expression for the elastic displacement field u of a triangular loop ABC is given as

$\mathbf {u(x)} =-{\frac {\Omega }{4\pi }}\mathbf {b} -{\frac {1-2\nu }{8\pi (1-\nu )}}\left[\mathbf {f} _{AB}+\mathbf {f} _{BC}+\mathbf {f} _{CA}\right]+{\frac {1}{8\pi (1-\nu )}}\left[\mathbf {g} _{AB}+\mathbf {g} _{BC}+\mathbf {g} _{CA}\right]$

(1)

where

$\mathbf {f} _{AB}=\left(\mathbf {b} \times \mathbf {t} _{AB}\right)\ln \left({\frac {R_{B}}{R_{A}}}\cdot {\frac {1+{\boldsymbol {\lambda }}_{B}\cdot \mathbf {t} _{AB}}{1+{\boldsymbol {\lambda }}_{A}\cdot \mathbf {t} _{AB}}}\right)$

(2)

and

$\mathbf {g} _{AB}=({\boldsymbol {\lambda }}_{A}+{\boldsymbol {\lambda }}_{B}){\frac {\mathbf {b} \cdot ({\boldsymbol {\lambda }}_{A}\times {\boldsymbol {\lambda }}_{B})}{1+{\boldsymbol {\lambda }}_{A}\cdot {\boldsymbol {\lambda }}_{B}}}.$

(3)

f_BC and f_CA have the same expression as Eqn.(2) except cyclic permutation of A, B and C, and so do g_BC and g_CA with Eqn.(3). R_A, R_B and R_C are the magnitudes of R_A, R_B and R_C, respectively. (See Fig.1.) ${\boldsymbol {\lambda }}_{A}$, ${\boldsymbol {\lambda }}_{B}$ and ${\boldsymbol {\lambda }}_{C}$ are the unit vectors of R_A, R_B and R_C, respectively, e.g. ${\boldsymbol {\lambda }}_{A}=\mathbf {R} _{A}/R_{A}$. And t_AB, t_BC and t_CA are the unit tangent vectors of the segments AB, BC and CA, respectively. b is the loop Burgers vector.

The solid angle $\Omega$ in Eqn.(1) is expressed as

$\Omega =-2\tan ^{-1}{\frac {\mathbf {R} _{\textrm {A}}\cdot (\mathbf {R} _{\textrm {B}}\times \mathbf {R} _{\textrm {C}})}{R_{\textrm {A}}R_{\textrm {B}}R_{\textrm {C}}+(\mathbf {R} _{\textrm {A}}\cdot \mathbf {R} _{\textrm {B}})R_{\textrm {C}}+(\mathbf {R} _{\textrm {A}}\cdot \mathbf {R} _{\textrm {C}})R_{\textrm {B}}+(\mathbf {R} _{\textrm {B}}\cdot \mathbf {R} _{\textrm {C}})R_{\textrm {A}}}}.$

(4)

The sign of $\Omega$ is reversed from the original expression to be in accordance with the convention of Barnett (1985).

By displacing all the atoms according to Eqn.(1), we can create a triangular dislocation loop in the atomistic simulation as shown in Fig.2(a).

(a)	(b)
Fig.2 (a) A triangular and (b) a tetragonal dislocation loop in DC Si. $\mathbf {C} _{1}=8a[1\,1\,{\bar {2}}]$, $\mathbf {C} _{2}=20a[1\,1\,1]$ and $\mathbf {C} _{3}=12a[1\,{\bar {1}}\,0]$, where a is the lattice constant and $a=5.431\mathrm {\AA}$ in Si. The loop Bergurs vector $\mathbf {b} =a/2[1\,{\bar {1}}\,0]$ in both (a) and (b).

---

Construction of a Polygonal Dislocation Loop

For the displcement field of a polygonal dislocation loop, we first decompse the polygonal loop into a sum of triangular loops, calculate the fields of each loops according to Eqn.(1)-(4) and finally sum them up. For example, Fig.3 shows how a tetragonal dislocation loop can be decomposed of multiple triangular loops.

Note that the displacement given in Eqn.~(1) changes its sign depending on the circular direction of the loop, or $\mathbf {u} _{\triangle 123}=-\mathbf {u} _{\triangle 321}$. Hence, when we sum fields of two or more triangular loops, we need to carefully decide the circular direction of each loop such that the displacement fields due to the shared segment of neighboring triangular loops can be canceled each other. (See the circular arrows in Fig.3.)

---

Using the potential energy

Sometimes, the simulation code does not implement the Virial stress calculation. In ab initio simulations, the Virial stress is usually not very accurate. In these cases, we can compute the elastic constants by monitoring the potential energy of the system, because the strain energy density $W$ is defined as

$W={\frac {1}{2}}\sigma _{ij}\epsilon _{ij}={\frac {1}{2}}\epsilon _{ij}C_{ijkl}\epsilon _{kl}$.

If the only nonzero strain is $\epsilon _{11}$, then the streain energy per volume becomes $W=C_{11}\epsilon ^{2}/2$. We can determine $C_{11}$ by fitting 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 W(\epsilon_{11})} to a parabola as shown in Fig.2(a). The elastic constnat 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}} can be determined in a similar way as shown in Fig.2(b). The elastic constants obtained in this way 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 C_{11}=} 160.5(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_{44}=} 60.2(GPa), which are almost the same as the results obtained by using Virial stress.

The elastic 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 C_{12}} is easily determined by the elastic relation

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}) }

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 B} is the bulk modulus, which we have already calculated in the manual 03. Since the bulk modulus of silicon is 108.3 (GPa) with SW potential, the elastic 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 C_{12}} is

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} = \frac{1}{2}(3B - C_{11}) = 82.2 } (GPa)

• Fig.2
• 

  (a)

• 

  (b)

Finding shear elastic constant for a non-tilting simulation cell

In the above sections, we compute 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}} by shearing the simulation cell. The simulation cell was rectangular before the shear, but becomes tilted after the shear. However, sometimes a simulation code allows to use only a rectangular simulation box and we can not tilt the box at all. In this case, how can we find the shear elastic 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 C_{44}} ? Notice that 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}} can be always obtained with a rectangular cell, when the repeat vectors of the cell are aligned with the cubic axes of the crystal. 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 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}} C12 known, we can determine 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}} using the orthogonal transformation rule for the 4th order tensor. When the repeat vectors are not aligned with the cubic axes, we should be able to compute linear combination 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}} , 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}} 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_{44}} while keeping the simulation cell rectangular. Remember that 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}} , 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}} 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_{44}} are components of a 4th order tensor which transforms 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 C_{ijkl} = T_{mi} T_{nj} T_{sk} T_{tl} C_{mnst} }

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 T_{ij} = \mathbf{e}_i \cdot \mathbf{e}'_j} is relation between two coordinate systems {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{e}_i} } 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 \mathbf{e}'_j} }.

For example, if we have one coordinate system 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{e}_1} = [1 0 0], 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{e}_2} = [0 1 0] 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 \mathbf{e}_3} = [0 0 1], and another coordinate system 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{e}'_1} = [1 1 0]/ \sqrt{2}, 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{e}'_2} = [−1 1 0]/ \sqrt{2} 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 \mathbf{e}'_3} = [0 0 1], then the matrix [T] is

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 [T] = \begin{pmatrix} 1/\sqrt{2} & -1/\sqrt{2} & 0 \\ 1/\sqrt{2} & 1/\sqrt{2} & 0 \\ 0 & 0 & 1 \end{pmatrix}}

In particular,

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 \begin{align} C_{1111} &= T_{m1} T_{n1} T_{s1} T_{t1} C_{mnst} \\ &= \sum_{i=1}^3 T_{i1} T_{i1} T_{i1} T_{i1} C_{11} \\ & + 2(T_{11}T_{11}T_{21}T_{21} + T_{11}T_{11}T_{31}T_{31} + T_{21}T_{21}T_{31}T_{31})C_{12} \\ & + 4(T_{21}T_{31}T_{21}T_{31} + T_{11}T_{31}T_{11}T_{31} + T_{11}T_{21}T_{11}T_{21})C_{44}. \end{align}}

Suppose we create a perfect crystal such that x, y, and z axes of the simulation cell are aligned 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 \mathbf{e}'_1} , 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{e}'_2} 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 \mathbf{e}'_3} . Then by straining the crystal along x axis, we can obtain <maht>C'_{1111}</math> = 180.4 (GPa) using the SW potential. Given 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}} = 160.6, 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}} = 80.5, the matlab script below gives 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}} = 59.9(GPa) for silicon.

% With given C11, C12, and C1111’, find C44 in cubic materials.
% C11, C12, and C44 are the elastic constants in Cartesian
% coordinate system.
% C1111’ is the elastic constant in the new coordinate.

syms C11 C12 C44
C11n = 160.6; C12n = 80.5; C44n = 60.2;
C1111prime = 180.40;
xyz0 = [ 1 0 0; % original coordinate
         0 1 0;
         0 0 1 ];
xyz1 = [ 1 1 0; % new coordinate
        -1 1 0;
         0 0 1 ];
T = zeros(3,3);
for i=1:3
    for j=1:3
        T(i,j) = dot(xyz0(i,:),xyz1(j,:)/norm(xyz1(j,:)));
    end
end

p=1; q=1; r=1; s=1;
A = T(1,p)*T(1,q)*T(1,r)*T(1,s) + T(2,p)*T(2,q)*T(2,r)*T(2,s) + ...
                                  T(3,p)*T(3,q)*T(3,r)*T(3,s);
B = T(1,p)*T(1,q)*T(2,r)*T(2,s) + T(1,p)*T(1,q)*T(3,r)*T(3,s) + ...
                                  T(2,p)*T(2,q)*T(3,r)*T(3,s);
C = T(2,p)*T(3,q)*T(2,r)*T(3,s) + T(1,p)*T(3,q)*T(1,r)*T(3,s) + ...
                                  T(1,p)*T(2,q)*T(1,r)*T(2,s);
CC = A*C11 + 2*B*C12 + 4*C*C44;
cc = subs(CC, C11, C11n); cc = subs(cc, C12, C12n);

disp(sprintf(’C44 = %f’,double(solve(cc-C1111prime,C44))))

Tcl version of MD++ script

For those who want to reproduce my results, I upload the Tcl version of MD++ script and the matlab script.
si_elastconst.tcl

# -*-shell-script-*-
#*******************************************
# Definition of procedures
#*******************************************
proc initmd { } { 
#MD++ setnolog
 MD++ setoverwrite
 MD++ dirname = runs/si_elastconst
 MD++ NIC = 200 NNM = 200
}

proc makecrystal { } { MD++ {
#--------------------------------------------
# Create Perfect Lattice Configuration
  element0 = Si   latticeconst = 5.43095
  crystalstructure = diamond-cubic  
  makecrystal
}} 

proc relax_fixbox { } { MD++ {
# Conjugate-Gradient relaxation
conj_ftol = 1e-6 conj_itmax = 1000 conj_fevalmax = 1000
conj_fixbox = 1  relax
} }

#*******************************************
# Main program starts here
#*******************************************
if { $argc == 0 } {
 set status 0
} elseif { $argc > 0 } {
 set status [lindex $argv 0]
}
puts "status = $status"

if { $status == 0 } {
# calculate elastic constants C_{11}, C_{12} and C_{44}
  initmd 

  MD++ latticesize = \[ 1 0 0 5 0 1 0 5 0 0 1 5 \]
  makecrystal 

  MD++ saveH eval
  set vol0     [MD++_Get OMEGA]
  
  set fileID [open "C11_C12.dat" w]
  for {set i -10} {$i <= 10} {incr i} {    
    set strain [expr $i/10000.0]
    MD++ restoreH input = \[ 1 1 $strain \] changeH_keepS
    relax_fixbox
    MD++ eval
  
    set sigma11 [expr (-1)*[MD++_Get TSTRESSinMPa_xx]]
    set sigma22 [expr (-1)*[MD++_Get TSTRESSinMPa_yy]]
    set EPOT    [MD++_Get EPOT]

    puts $fileID "[format %18.14e $strain]\t \
                  [format %18.14e $sigma11]\t \
                  [format %18.14e $sigma22]\t \
                  [format %18.14e [expr $EPOT/$vol0]]"
  }
  close $fileID

  set fileID [open "C44.dat" w]
  for {set i -10} {$i <= 10} {incr i} {    
    set strain [expr $i/10000.0]
    MD++ restoreH input = \[ 1 2 $strain \] changeH_keepS
    relax_fixbox
    MD++ eval
  
    set sigma12 [expr (-1)*[MD++_Get TSTRESSinMPa_xy]]
    set EPOT    [MD++_Get EPOT]

    puts $fileID "[format %18.14e $strain]\t \
                  [format %18.14e $sigma12]\t \
                  [format %18.14e [expr $EPOT/$vol0]]"
  }
  close $fileID  
  MD++ quit
  
} elseif { $status == 1 } {
# calculate elastic constants C'_{1111} in different coordinate
  initmd 

  MD++ latticesize = \[ 1 1 0 5 -1 1 0 5 0 0 1 7 \]
  makecrystal
  MD++ saveH eval
  set vol0     [MD++_Get OMEGA]
  
  set fileID [open "C1111prime.dat" w]
  for {set i -10} {$i <= 10} {incr i} {    
    set strain [expr $i/10000.0]
    MD++ restoreH input = \[ 1 1 $strain \] changeH_keepS
    relax_fixbox
    MD++ eval
  
    set sigma1111 [expr (-1)*[MD++_Get TSTRESSinMPa_xx]]
    set EPOT    [MD++_Get EPOT]
    
    puts $fileID "[format %18.14e $strain]\t \
                  [format %18.14e $sigma1111]\t \
                  [format %18.14e [expr $EPOT/$vol0]]"
  }
  close $fileID
  MD++ quit
  
} else {

  puts "unknown status = $status"
  MD++ quit
}

% Matlab script to calculate the elastic constants C11, C12 and C44
% of DC Si
clear, clf

B = 108.3e3;         % Bulk modulus of Si (in MPa) obtained from M02
data = load('C11_C12.dat');
strain11 = data(:,1);
sigma11 = data(:,2); % stress component sigma_{11} in MPa
sigma22 = data(:,3); % stress component sigma_{22} in MPa
Wn = data(:,4);       % strain energy density in eV/A^3

data = load('C44.dat');
strain12 = data(:,1);
sigma12 = data(:,2);
Ws = data(:,3);

if(exist('C1111prime.dat','file'))
    data = load('C1111prime.dat');
    strain1111 = data(:,1);
    sigma1111 = data(:,2);
    W1111 = data(:,3);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot the lattice energy vs lattice constant
% and find the equilibrium lattice constant, a_0
[P, S] = polyfit(strain11,sigma11,1);
C11 = P(1);
[P, S] = polyfit(strain11,sigma22,1);
C12 = P(1);
[P, S] = polyfit(strain11,Wn,2);
C11_e = P(1)*160.2e3*2;   
[P, S] = polyfit(strain12,sigma12,1);
C44 = P(1);
[P, S] = polyfit(strain12,Ws,2);
C44_e = P(1)*160.2e3*2; 
if(exist('C1111prime.dat','file'))
    [P, S] = polyfit(strain1111,sigma1111,1);
    C1111 = P(1);
    [P, S] = polyfit(strain1111,W1111,2);
    C1111_e = P(1)*160.2e3*2;    
end

disp(sprintf('-----------------------------------------------------'))
disp(sprintf('From Virial stress calculation,'))
disp(sprintf(['the elastic constants C11 = %10.8e, C12 = %10.8e, ',...
             'C44 = %10.8e in MPa.'],C11,C12,C44))
disp(sprintf('Then bulk modulus B = (C11 + 2*C12)/3 = %10.8e in MPa.',(C11+2*C12)/3))         
disp(sprintf('-----------------------------------------------------'))
disp(sprintf('From energy calculation,'))
disp(sprintf('the elastic constants C11 = %10.8e, C44 = %10.8ein MPa.',...
             C11_e, C44_e))
disp(sprintf('Then C12 = (3*B - C11)/2 = %10.8e in MPa', (3*B - C11_e)/2))
disp(sprintf('-----------------------------------------------------'))
if(exist('C1111prime.dat','file'))
disp(sprintf('-----------------------------------------------------'))
disp(sprintf('In e1=[1 1 0], e2=[-1 1 0] and e3=[0 0 1] coordinate, '))
disp(sprintf('the elastic constants C1111 = %10.8e (MPa) from stress calculation',C1111))
disp(sprintf('                      C1111 = %10.8e (MPa) from energy calcualtion',C1111_e))
disp(sprintf('-----------------------------------------------------'))
end

figure(1)
set(gca,'fontsize',14)
plot(strain11, sigma11, '.-',strain11, sigma22,'.-r',...
     strain12, sigma12, '.-g'), grid
xlabel('Strain','fontsize',16), 
ylabel('Stress (MPa)','fontsize',16)
legend('\epsilon_{11} vs. \sigma_{11}','\epsilon_{11} vs. \sigma_{22}',...
       '\epsilon_{12} vs. \sigma_{12}',2)

figure(2)
set(gca,'fontsize',14)
plot(strain11, Wn, '.-', strain12, Ws, '.-r'), grid
xlabel('Strain','fontsize',16), 
ylabel('Strain Energy Density, \it W \rm (eV/Angstrom^3)','fontsize',16)
legend('\epsilon_{11}','\epsilon_{12}',2)