Manual 09 for MD++
Computing Ideal Strength
Keonwook Kang and Wei Cai
Ideal strength is the theoretical limit up to which a defect-free material can sustain external load before it breaks. We define ideal strength as the maximum derivative of the excess energy (per area) that comes from rigid translation of part of material.<ref>There is different definition of ideal strength. Roudy and Cohen calculated it from the excess energy of uniform deformation of material.(D. Roudy and M. L. Cohen, Phys. Rev. B 64 212103 (2001))</ref> In this manual, it is explained how the ideal strengths of diamond cubic structure of Si are calculated with MD++ codes. This task was part of work covered in our paper<ref>Keonwook Kang and Wei Cai, Philosophical Magazine, 87 2169–2189 (2007)</ref>.
Contents |
Ideal Tensile Strength
Let’s assume that we want to calculate the ideal tensile strength along the [1 1 0] direction of Si. First, a simulation cell is created with the three periodicity vectors , and , each of which is aligned along x, y and z direction, respectively. Make sure that the cell is subjected to periodic boundary condition (PBC) in all three directions. Next, the potential energy change is recorded as the function of the opening distance d when the (1 1 0) surface opens rigidly (See Fig. 1(a).) or
γ = γ(d). |
where γ is the excess energy to create two free surfaces. In order to rigidly separate a bulk crystal, we could collectively translate half of atoms but instead we used a simple trick. When we redefine the periodicity vector as
. |
while fixing the coordinates of all the atoms, it gives the same effect of separating the bulk crystal into two by d. We change the box length by redefining the periodicity vector from its initial definition with certain ratio and it will create a gap of width d in x direction between the periodic image (2) and the crystal (1) in the primary cell as shown in Fig. 1(a).
In MD++ Codes, there is a command changeH_keepR which changes one of periodicity vectors with atoms fixed. The input to this command has three numbers.
input = [ i j frac ] changeH_keepR
The first number is the index of the periodicity vector that you want to change. For example, if i=1, vector will be changed. (2 for and 3 for ) The second number is the index of direction. If j=1, the fraction frac of the vector is added to itself and can be longer or shorter depending on the sign of frac. If j=2, the vector will be tilted by the fraction frac of vector.
Do not forget to refresh the neighbor list whenever you change while atoms’ positions are fixed. MD++ automatically refreshes the neighbor list based on the change of the atoms’ positions. If atoms are fixed, the list will not be updated even if there could be a virtual displacement of atoms by changing . In the script (given at the end of this manual), you will see the command clearR0, right after changeH_keepR, which enforces the neighbor list to be refreshed. There are more MD++ commands to be explained: saveH stores the current 3×3 matrix H to H0 and restoreH restores H from H0. RHtoS converts the real coordinate of an atom i to the scaled coordinate by .
If we plot γ vs. d, the curve looks like Fig. 2(a). As d increases, the energy also increases up to a plateau 2γs, where γs is the (unrelaxed) surface energy of (1 1 0) plane of Si. Once the opening distance becomes longer than certain cutoff distance, the interaction between two parts no longer plays a role and the energy becomes constant and irrelative to d. You may be aware that the slope of the γ(d) curve has the unit of stress. We define the maximum value of the slope as the ideal tensile strength along the separation direction, e.g. [1 1 0]. The SW potential model gives =41.7 (GPa) for Si.
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 c1 = Nx[1 0 ¯1], c2 = Ny[1 1 1] and c3 = Nz[1 ¯2 1] to calculate the ideal shear strength � th (111)[10¯1]. The simulation cell is subject to PBC in all three directions. With the coordinates of all the atoms fixed, redefining c2 vector as
(GPa) |
Ideal Shear Strength in Non-tilting Cell
In particular,
Script
MD++ Script
# -*-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). 7 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 # enfore refresh neighborlist 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 8 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] # 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 9 # 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 }
Matlab Script
%#! /usr/local/bin/octave -qf % Matlab script, cal_idealstr.m % to calculate ideal tensile strength "sigma" and shear strength "tau" clear %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 10 % 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); 11 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