View source for M09 Computing Ideal Strength

<H1 ALIGN="CENTER"><FONT SIZE="-1">Manual 09 for MD++ </FONT> 
<BR>
Computing Ideal Strength</H1>
<DIV>

<P ALIGN="CENTER"><STRONG>Keonwook Kang and Wei Cai</STRONG></P>
</DIV>
<P ALIGN="CENTER">[[Dec. 8]], [[2007]]</P>
<P>
<BR><HR>

Ideal strength is the maximum stress that can exist in a perfect crystal.  When the applied stress exceeds the ideal strength, the crystal structure will collapse even at zero temperature.  The mode of failure depends on the type of applied stress.  When the tensile stress exceeds the ideal tensile strength, the crystal is unstable against cleavage fracture.  When the shear stress exceeds the ideal shear strength, the crysal is unstable against slip on crystallographic planes.  Because a real crystal is not perfect but contains defects (such as surfaces, grain boundaries, dislocations, vacancies), which lowers the mechanical strength of the crystal, the ideal strength provides an upper bound of the strength of a real crystal.  When the defect density is made sufficiently small (such as in whiskers), the strength of a real crystal can approach the ideal strength.  The ideal strength is also called the ''theoretical'' strength.

There are two alternative definitions of the ideal strength depending on the way that the deformation is imposed on the perfect crystal structure.  In the first definition, the deformation strain can be applied uniformly to the simulation cell (usually under periodic boundary conditions) resulting in an affine deformation of the entire crystal<ref>D. Roudy and M. L. Cohen, Phys. Rev. B '''64''' 212103 (2001))</ref>.  In the second definition, the deformation strain can be localized between two crystallographic planes<ref name=kangcai2007>Keonwook Kang and Wei Cai, Philosophical Magazine, '''87''' 2169–2189 (2007)</ref>.  The ideal strength values obtained from these two approaches are somewhat different but are also expected to co-relate with each other.  In this manual, we describe how to compute the ideal strength according to the second definition, using diamond-cubic structure of Si as an example.

<BR><HR>

== Ideal Tensile Strength ==
Suppose we want to compute the ideal tensile strength of Si along the [1 1 0] direction, i.e., <math>\sigma_{[110]}^{\mathrm{th}}</math>.  First, we create a simulation cell under periodic boundary condition (PBC) with the three periodicity vectors <math>\mathbf{c}_1 =
N_x [1 1 0]</math>, <math>\mathbf{c}_2 = N_y[\bar{1} 1 0]</math> and <math>\mathbf{c}_3 = N_z[0 0 1]</math>.  The coordinate system is defined such that <math>\mathbf{c}_1</math>, <math>\mathbf{c}_2</math> and <math>\mathbf{c}_3</math> are aligned along <math>x</math>, <math>y</math> and <math>z</math> directions, respectively. The crystal just created is surrounded by its own images in all three directions. 

Next, we change the length of the repeat vector <math>\mathbf{c}_1</math> while keeping the position of all atoms in the simulation cell fixed.  Let <math>\mathbf{c}_1^0</math> be the original value of the repeat vector along the <math>x</math> direction and <math>L_x^0 = |\mathbf{c}_1^0|</math>.  Suppose we assign the new repeat vector to be

{|border="0" align="center"
|<math>\mathbf{c}_1 = \left( 1 + \frac{d}{L_x^0}\right) \mathbf{c}_1^0 </math>.
|}

This will create a gap of width <math>d</math> between the crystal in the simulation cell (primary cell) with its images in the <math>x</math> direction, as shown in Fig. 1(a).  For each value of <math>d</math>, we compute the potential energy of the crystal (without relaxing the positions of the atoms).  The excess energy with respect to the case of <math>d=0</math>, dividied by the cross section area of the simulation cell perpendicular to the <math>x</math> direction, is defined as function of <math>d</math>

{|border="0" align="center"
|<math> \gamma = \gamma ( d ) </math>.
|}

and is plotted in Fig. 2(a).

In MD++, the transformation of the simulation cell as described above can be performed conveniently using the command '''changeH_keepR'''.  <math>\mathbf{H}</math> is a 3 by 3 matrix whose three columns correspond to the coordinates of the three repeat vectors, i.e. <math>\mathbf{H} = [\mathbf{c}_1 | \mathbf{c}_2 | \mathbf{c}_3]</math>.  The syntax for this command is

 input = [ ''i'' ''j'' ''frac'' ] changeH_keepR

This will lead to the following transformation of the repeat vectors

{|border="0" align="center"
|<math>\mathbf{c}_i := \mathbf{c}_i + frac \cdot \mathbf{c}_j </math>.
|}

Since we want to change the length of <math>\mathbf{c}_1</math> without changing its direction, both ''i'' and ''j'' should be 1 for computing the ideal tensile strength along <math>x</math> direction.  (''i'' and ''j'' will be different when computing the ideal shear strength as described in the next section.)  A related command '''changeH_keepS''' performs identical transformations to the repeat vectors, but the scaled coordinates, instead of the real coordinates of the atoms are kept fixed.  This will lead to an affine transformation of the simulation cell and is related to the first definition of the ideal strength as described above.

After the '''changeH_keepR''', the potential energy is computed by the command '''eval'''.  However, before computing the potential energy, we need to refresh the neighbor list.  MD++ uses the Verlet list algorithm, in which the neighbor list is automatically refreshed whenever any atom is displaced by more than half the skin of the Verlet list since the last time the neighbor list is constructed.  However, since here the atom positions do not change when we change the repeat vectors, the neighbor list will not be refreshed automatically.  We can force the neighbor list to refresh by calling '''clearR0''' right after the '''changeH_keepR''' command.

The MD++ script for this calculation is attached at the end of this manual.  A few more commands are used there to facilitate the calculation.  The '''saveH''' command saves the current value of the <math>\mathbf{H}</math> matrix to <math>\mathbf{H}_0</math>.  The '''restoreH''' command copies <math>\mathbf{H}_0</math> back to <math>\mathbf{H}</math>.  This commands are used to make the matrix <math>\mathbf{H}</math> is always restored to the original value whenever the '''changeH_keepR''' command is called.  The '''RHtoS''' command is always called before '''eval''' which converts the real coordinates <math>\mathbf{r}_i</math> of the atoms to the scaled coordinates <math>\mathbf{s}_i</math>, <math>\mathbf{s}_i = \mathbf{H}^{-1} \mathbf{r}_i</math>. This is needed because the calculation of potential energy in MD++ is based on the scaled coordinates.


Fig. 2(a) plots the excess energy per unit area of the simulation cell, <math>\gamma</math>, as a function of separation distance <math>d</math> for the Stillinger-Weber (SW) model of Si.  For large enough <math>d</math>, <math>\gamma(d)<\math> reaches a plateau, whose value is <math>2\gamma_s</math>, where <math>\gamma_s</math> is the (unrelaxed) surface energy of (1 1 0) plane of Si.  The plateau is reached when <math>d</math> becomes greater than the cut-off distance of the interatomic potential.  The slope of the <math>\gamma(d)</math> has the unit of stress.  The maximum slope of this curve can be interpreted as the maximum stress required to separate two rigid blocks of Si along the [1 1 0] direction, i.e. the ideal tensile strength.  For the SW model of Si, <math>\sigma_{[110]}^{\mathrm{th}}</math> = 41.7 GPa.

<gallery caption="Fig.1" widths="300px" heights="300px" perrow="2">
Image:Cut-surfL.jpg|(a) Separating a bulk crystal across a (1 1 0) plane by distance <math>d</math>.
Image:Cut-slip.jpg|(b) The generalized stacking fault energy for slip on a (1 1 1) plane can be computed by tilting the repeat vector <math>\mathbf{c}_2</math> of the simulation cell while keeping the position of all atoms fixed. These figures are excerpted from our paper.<ref name="kangcai2007"/> 
</gallery>

== Ideal Shear Strength ==
The procedure of calculating ideal shear strength is pretty much similar to that for
ideal tensile strength except that we now have to deal with 2-dimensional array
of energy. In Fig. 1(b), it is illustrated how we slip part of a bulk crystal to the
rest. For example, the periodicity vectors are chosen to be <math>\mathbf{c}_1 = N_x[1 0 \bar{1}]</math>, <math>\mathbf{c}_2 = N_y[1 1 1]</math> and <math>\mathbf{c}_3 = N_z[1 \bar{2} 1]</math> to calculate the ideal shear strength <math>\sigma^{\mathrm{th}}_{(111)[10\bar{1}]}</math>. The simulation cell is subject to PBC in all three directions. With the coordinates
of all the atoms fixed, redefining <math>\mathbf{c}_2</math> vector as

{|border="0" align="center"
|<math> \begin{align} \mathbf{c}_2 &= \left(1+ \frac{\delta}{L_z^0}\mathbf{c}_2^0 + \frac{d}{L_x^0}\mathbf{c}_1^0 \right) \\
                &= (L_z^0 + \delta)\mathbf{n} + d\mathbf{m}
        \end{align} </math>
|}

makes <math>\mathbf{c}_2</math> vector tilted and as a result the periodic image (2) will be shifted and moved from the crystal (1) in the primary cell by <math>d</math> and <math>\delta</math>, respectively. Notice that we do not simply slip but also allow the crystal move normal to slip plane and the energy <math>\gamma</math> is a function of <math>d</math> and <math>\delta</math> or

{|border="0" align="center"
|<math> \gamma = \gamma ( d, \delta ) </math>.
|}

Because the crystal will in reality adjust the gap <math>\delta</math> to reach energy-minimum state for a given slipped distance <math>d</math>, what we really want is <math>\min_{\delta}(d)</math> and denote it as <math>\Gamma</math> or

{|border="0" align="center"
|<math> \Gamma = \displaystyle\min_{\{\delta\} }\gamma ( d ) </math>.
|}

The curve <math>\Gamma(d)</math> is shown in Fig. 2(b). The shifted distance <math>d</math> is normalized by the magnitude of Burgers vector and we observe that the misfit energy is
periodic after one Burgers vector. The ideal shear strength is also defined as the
maximum slope of the misfit energy curve and it is 9.5 GPa from SW potential model.

You have to keep it in mind that the shuffle set plane of (1 1 1) is your slip
plane. Fig. 3 shows two configurations which would be same unless you tilt
the periodicity vector. After the vector <math>\mathbf{c}_2</math> is tilted, one gives 
<math>\Gamma</math> of shuffle-set plane and the other gives that of glide-set plane. And the ideal shear strength of glide-set plane (''e.g.'' 77.3GPa from SW) is much higher than that of shuffle-set plane. In MD++ script, this is handled by the command '''pbcshiftatom''', which shifts whole atoms periodically by '''dx''', '''dy''', and '''dz''' in all three directions.

 input = [dx dy dz] pbcshiftatom

The numbers '''dx''', '''dy''' and '''dz''' are given in scaled coordinate.

In MD++, there is a command '''calmisfit''' calculating the misfit energy <math>\gamma</math> as <math>d</math> and <math>\delta</math> change. 

<pre>
 input = [ surfacenormal
           x0 dx x1
           y0 dy y1
           z0 dz z1 ]
 calmisfit
</pre>

The input variable needs seven elements, the first of which indicates the surface normal index. If it is 2, the slip plane normal is in <math>y</math>-direction (1 for <math>x</math> and 3 for <math>z</math>-direction) and the periodicity vector <math>\mathbf{c}_2</math> is going to be tilted. The rest six numbers are basically grid points to designate to which
direction and how much the periodicity vector is tilted (or stretched). The numbers are normalized by the length of one lattice vector in each direction. If '''dx'''=0.001, the increment of tilting the vector in <math>x</math>-direction would be <math>0.001 \times a \times |\mathbf{c}_1|/N_x</math>, where <math>a</math> is the lattice constant. If '''dy'''=0.01, the vector <math>\mathbf{c}_2</math> increases by the amount of <math>0.01 \times a \times |\mathbf{c}_2|/N_y</math> in <math>y</math> direction. In our example, the last
three numbers are all zeros or '''z0''' = '''dz''' = '''z1''' = 0.

After running '''calmisfit''', a text file “Emisfit.out” will be dumped out and it
will be used to calculate the ideal shear strength by the Matlab script attached
at the end.

<gallery caption="Fig.2" widths="300px" heights="300px" perrow="2">
Image:Sisw_surfE.jpg|(a) The excess energy per area is plotted as a function of the opening
distance <math>d</math>. After certain cut-off distance, the energy curve becomes flat and the value corresponds to <math>2\gamma_s</math> where <math>\gamma_s</math> is the surface energy.
Image:Sisw_slipE.jpg|(b) Misfit energy <math>\Gamma</math> vs. the shifted distance <math>d</math>. SW potential model is used for atomistic simulation in (a) and (b). 
</gallery>

<gallery caption="Fig.3" widths="300px" heights="300px" perrow="2">
Image:Glideset_e.jpg |(a) 
Image:Shuffleset.jpg |(b) Two configurations in (a) and (b) are essentially identical before the periodicity vector <math>\mathbf{c}_2</math> is tilted. After <math>\mathbf{c}_2</math> is tilted along <math>\mathbf{c}_1</math> direction, the crystal is shifted on glide-set plane in (a) and shuffle-set plane in (b).
</gallery>

== Ideal Shear Strength in Non-tilting Cell ==
Sometimes we simply cannot tilt the periodicity vectors. For example, the
original Baskes’ MEAM codes are implemented only for a rectangular cell. In
this case, the procedures explained above can not be applied and we have to
move a chunk of atoms literally. Fig. 4(a) shows the initial configuration in
which a few layers at top and bottom are removed in order to have room for
atoms to be shifted as shown in Fig. 4(b). Other than that, the rest of process
to get the ideal shear strength is identical.

To remove top and bottom atoms, two MD++ commands '''fixatoms_by_position'''
and '''removefixedatoms''' are used. The command '''fixatoms_by_position''' fixes
atoms within a hexahedral volume defined by six numbers ('''x0''', '''x1''', '''y0''', '''y1''', '''z0''' and '''z1''' in scaled coordinate) in input

 input = [enable x0 x1 y0 y1 z0 z1 ] fixatoms_by_position

when '''enable'''=1. Then the command '''removefixedatoms''' removes those fixed
atoms.

MD++ codes have another command '''calmisfit2''' to calculate the misfit energy
in non-tilting cell. The input variable for '''calmisfit2''' has six more numbers
('''xmin''', '''xmax''', '''ymin''', '''ymax''', '''zmin''' and '''zmax''') than '''calmisfit''' to define a hexahedral volume which will be moved by the command. Those six numbers should be given in scaled coordinate. The first 10 numbers in the input variable are explained in the command '''calmisfit'''.

<pre>
 input = [ surfacenormal
           x0 dx x1
           y0 dy y1
           z0 dz z1
           xmin xmax ymin ymax zmin zmax]
 calmisfit2
</pre>

The command '''calmisfit2''' also produces a text file, “Emisfit.out” of the same
format and the same Matlab script will read it to calculate the ideal shear
strength.

<gallery caption="Fig.4" widths="300px" heights="300px" perrow="2">
Image:Beforeshift_e.jpg |(a) Configuration to calculate the ideal shear strength in non-tilting simulation cell: (a) the initial configuration and (b) the shifted configuration
Image:Aftershift.jpg |(b) 
</gallery>

== Script ==
The MD++ script and the Matlab script are given below for those who want
to reproduce the results or try to calculate the ideal strength of other crystal
materials. Let’s say MD++ script file name is “idealstr.tcl”. If you type

 $ bin/sw_gpp scripts/work/idealstr.tcl 0

it calculates the surface energy of (1 1 0) plane in terms of the opening distance
<math>d</math>. If you type

 $ bin/sw_gpp scripts/work/idealstr.tcl 1

it calculates the misfit energy of (1 1 1) plane along <math>[1\, 0\, \bar{1}]</math> direction by tilting the periodicity vector <math>\mathbf{c}_2</math>. If you type

 $ bin/sw_gpp scripts/work/idealstr.tcl 2

it will calculate the same misfit energy in different method explained in Sec. 3.


=== MD++ Script ===
<pre>
# -*-shell-script-*-
source "scripts/Examples/Tcl/startup.tcl"
#*******************************************
# Definition of procedures
#*******************************************
proc initmd { } {
  MD++ setnolog
  MD++ setoverwrite
  MD++ dirname = runs/si_idealstr
  MD++ NIC = 200 NNM = 200
}

proc openwindow { } { MD++ {
# Plot Configuration
  atomradius = 0.67 bondradius = 0.3 bondlength = 0
  atomcolor = orange bondcolor = red backgroundcolor = gray70
  win_width = 300 win_height = 300
  plotfreq = 10 rotateangles = [ 0 0 0 1.25 ]
  plot_atom_info = 1 # reduced coordinates of atoms
#  plot_atom_info = 2 # real coordinates of atoms
#  plot_atom_info = 3 # energy of atoms
openwin alloccolors rotate saverot eval plot
} }

#*******************************************
# 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 surface opening energy on (1 1 0).
  initmd
  MD++ {
    element0 = Si latticeconst = 5.43095
    crystalstructure = diamond-cubic
    latticesize = [ 1 1 0 5
                   -1 1 0 5
                    0 0 1 7 ]
    makecrystal }
  MD++ saveH eval
  set Lx0 [MD++_Get H_11]
  set Ly0 [MD++_Get H_22]
  set Lz0 [MD++_Get H_33]
  set EPOT0 [MD++_Get EPOT]
  set Area [expr $Ly0 * $Lz0]
  set fileID [open "Esurf.dat" w]
  for {set i -10} {$i <= 100} {incr i} {
    set strain [expr $i/1000.0]
    MD++ restoreH RHtoS
    MD++ input = \[ 1 1 $strain \] changeH_keepR
    MD++ clearR0 # enforce refresh neighbor list
    MD++ eval
    set Lx [MD++_Get H_11]
    set dist [expr $Lx - $Lx0]
    set EPOT [MD++_Get EPOT]
    # Esurf2 : excessive energy per area to create 2 free surfaces.
    set Esurf2 [expr ($EPOT - $EPOT0)/$Area]
    puts $fileID "[format %18.14e $dist]\t \
                  [format %18.14e $Esurf2]"
  }
  close $fileID
  MD++ quit

} elseif { $status == 1 } {
# calculate generalized stacking fault energy on (1 1 1) along [1 0 -1]
  initmd
# surface normal direction 1/x 2/y 3/z
  set surfacenormal 2
  set lattconst 5.43095
  set c1 "1 0 -1"; set c2 "1 1 1"; set c3 "1 -2 1"
  set Nx 5; set Ny 4; set Nz 3

  set gridx " 0.000 0.005 0.500"
  set gridy "-0.050 0.010 0.150"
  set gridz " 0.000 0.050 0.000"

# write the matlab script
  set fileID [open "parameter.m" w]
  puts $fileID "normaldir = $surfacenormal"
  puts $fileID "latticeconst = $lattconst"
  puts $fileID "latticesize = \[ $c1 $Nx "
  puts $fileID "                 $c2 $Ny "
  puts $fileID "                 $c3 $Nz \]"
  puts $fileID "dxyz = \[ $gridx "
  puts $fileID "          $gridy "
  puts $fileID "          $gridz \]"
  close $fileID

  MD++ {
    element0 = Si
    crystalstructure = diamond-cubic
  }
  MD++ latticeconst = $lattconst
  MD++ latticesize = \[$c1 $Nx $c2 $Ny $c3 $Nz\]
  MD++ makecrystal

# To make shuffle set of (111) plane exposed
  MD++ { input = [0 0.05 0 ] pbcshiftatom }
  MD++ input = \[ $surfacenormal $gridx $gridy $gridz \]
  MD++ calmisfit
  MD++ quit

} elseif { $status == 2 } {
# calculate generalized stacking fault energy on (1 1 1) along [1 0 -1]
  initmd
# surface normal direction 1/x 2/y 3/z
  set surfacenormal 2
  set lattconst 5.431
  set c1 "1 0 -1"; set c2 "1 1 1"; set c3 "1 -2 1"
  set Nx 5; set Ny 6; set Nz 3
  set gridx " 0.000 0.005 0.500"
  set gridy "-0.050 0.010 0.150"
  set gridz " 0.000 0.050 0.000"
  set xmin -1.00; set xmax 1.00
  set ymin  0.01; set ymax 1.00
  set zmin -1.00; set zmax 1.00

# write the matlab script
  set fileID [open "parameter.m" w]
  puts $fileID "normaldir = $surfacenormal"
  puts $fileID "latticeconst = $lattconst"
  puts $fileID "latticesize = \[ $c1 $Nx "
  puts $fileID "                 $c2 $Ny "
  puts $fileID "                 $c3 $Nz \]"
  puts $fileID "dxyz = \[ $gridx "
  puts $fileID "          $gridy "
  puts $fileID "          $gridz \]"
  close $fileID

  MD++ {
    element0 = Si
    crystalstructure = diamond-cubic
  }
  MD++ latticeconst = $lattconst
  MD++ latticesize = \[$c1 $Nx $c2 $Ny $c3 $Nz\]
  MD++ makecrystal

  MD++ input = \[1 -1 1 0.4 1 -1 1 \] fixatoms_by_position removefixedatoms
  MD++ input = \[1 -1 1 -1 -0.42 -1 1 \] fixatoms_by_position removefixedatoms
  MD++ input = \[ $surfacenormal $gridx $gridy $gridz $xmin $xmax $ymin $ymax $zmin $zmax \]
  MD++ calmisfit2

  openwindow
  MD++ sleep quit
  MD++ quit

} else {
puts "unknown status = $status"
MD++ quit
}
</pre>

=== Matlab Script ===
After running MD++, you just run the following matlab script in the directory
where the output files exist. Then the matlab script will plot the energy curve
and calculate ideal strengths.

<pre>
%#! /usr/local/bin/octave -qf
% Matlab script, cal_idealstr.m
% to calculate ideal tensile strength "sigma" and shear strength "tau"
clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate ideal tensile strength along [110] direction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data = load(’Esurf.dat’);
d = data(:,1);     % openning distance in A
gamma = data(:,2); % excessive energy per area creating 2 free surfaces (eV/A^2)

sigma = max(gamma(2:end)-gamma(1:end-1))/(d(2) - d(1)) * 160.2;
disp(sprintf(’Ideal tensile stregnth = %e(GPa)’,sigma));

figure(1)
set(gca,’fontsize’,14)
plot(d, gamma, ’.-’), grid
xlabel(’Opening Distance \it d \rm(Angstrom)’,’fontsize’,16),
ylabel(’\gamma (eV/Angstrom^2)’,’fontsize’,16)
line([-1 1.1],[gamma(end) gamma(end)],’linewidth’,2,’color’,’k’)
text(-1.3, gamma(end),’2\gamma_s’,’fontsize’,16)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate ideal shear strength on (111) along [110] direction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (~exist(’parameter.m’,’file’))
    error(’No parameter.m found’)
else
    parameter

    dx = dxyz(1,2); dy = dxyz(2,2); dz = dxyz(3,2);
    vx = latticesize(1,1:3); Lx = latticeconst*norm(vx);
    vy = latticesize(2,1:3); Ly = latticeconst*norm(vy);
    vz = latticesize(3,1:3); Lz = latticeconst*norm(vz);

    data = load(’Emisfit.out’);
    Mx=max(data(:,1)); My=max(data(:,2)); Mz=max(data(:,3));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert data to relaxed energy
    MM=-sort(-[Mx+1,My+1,Mz+1]);         % Then MM(1) > MM(2) > MM(3)
    E = reshape(data(:,4),MM(1),MM(2));

    Erlx = zeros(MM(1),1);
    Zrlx = zeros(MM(1),1);

    normalgrid = dxyz(normaldir,1):dxyz(normaldir,2):dxyz(normaldir,3);
    normalfinegrid = dxyz(normaldir,1):dxyz(normaldir,2)/1000:dxyz(normaldir,3);

    for i = 1:MM(1),
        a = E(i,:);
        a_interp = interp1(normalgrid,a,normalfinegrid,’spline’);
        [Erlx(i), ind] = min(a_interp);
        Zrlx(i) = normalfinegrid(ind);
    end

% retrieve unrelaxed energy
    Eunrlx=E(:,find(abs(normalgrid)<1e-8));

% calculate ideal shear strength
    tau = max(Erlx(2:end)-Erlx(1:end-1))/(dx*Lx) * 160.2;
    disp(sprintf(’Ideal shear stregnth = %e(GPa)’,tau));

    figure(2); set(gca,’fontsize’,14)
    plot([0:Mx]/Mx, Erlx-Erlx(1), ’b.-’,...
         [0:Mx]/Mx, Eunrlx-Eunrlx(1),’r-’);
    legend(’relaxed’,’unrelaxed’);
    xlabel(’Sheared Distance, d/|b| where b = [1 0 1]/2’,’fontsize’,16)
    ylabel(’\Gamma (eV/Angstrom^2)’,’fontsize’,16)

    figure(3); set(gca,’fontsize’,14)
    plot([0:Mx]/Mx,Zrlx, ’r.-’, ...
         [0 1],[1 1]*dxyz(normaldir,1),’b-’,...
         [0 1],[1 1]*dxyz(normaldir,3),’b-’);
    xlabel(’Sheared Distance, d/|b| where b = [1 0 1]/2’,’fontsize’,16)
    ylabel({’Lowest energy height,’; ’\delta/|c_2| where c_2 = [1 1 1]’},’fontsize’,16)
end
</pre>

== Notes ==

<references/>