|
| |
Calculation of strain in nano heterostructures
The physical foundations are described in
Stefan Hackenbuchner's PhD thesis
in the chapter about strain (see also there in the appendix).
Here the details of the program are explained.
Subroutine
set_strain determines in which way to treat strain
(homogeneous-strain (1D) or elasticity theory (2D/3D)).
For option 'strain-minimization',
subroutine
calculate_strain is called:
First the POINTER
uV which containes the displacements is allocated and set to zero.
This subroutine then calls
subroutine
define_strain_problem (setup of strain matrix) and
then solves the linear equation
with the BICGSTAB-method (bi-conjugate gradient stabilized) or with Kent's
eigenvalue solver (preconditioned conjugate gradient: PCG_Kent). The solution is
written to the displacement vector
uV.
The POINTER uM
that stores the displacements for each point (MODULE
strain_data) is allocated and set to the according values of
uV. uM
refers to the
calculation system.
Now the strain is provided by
subroutine
extract_strain_for_guenther(...,crystalL=.FALSE.,...), which provides the strain for
material
grid points in the
calculation system.
NOTE: The strain is
basically calculated for each octant (from view of Stefan's grid points).
So one takes the eight relevant octants and calculates the average to get the
strain on material
grid point.
What we do now:
- Strain tensor is
calculated on the physical grid and then interpolated to the material grid.
What we should
do instead:
- Interpolate
displacements to the material grid and then calculate strain tensor.
The subroutine
extract_strain_on_octant provides the strain tensor for a special
octant [(i,j,k) on the physical grid and boxnum=1,..,8
(1,...,4 in 2D].
The subroutine
extract_strain_for_stefan(...,crystalL=.FALSE.,...)
provides the strain tensor for a point (i,j,k)
on the physical grid in the
calculation system (i.e. takes the average of the relevant octants).
There are also routines that provide the strain tensor in the
crystal system:
-
extract_strain_for_guenther(...,crystalL=.TRUE.,...) (->
flag
crystalL=.TRUE.)
-
extract_strain_for_stefan (...,crystalL=.TRUE.,...) (->
flag
crystalL=.TRUE.)
The strain matrix is set up with MODULE
strain_matrix:
In order to calculate the strain several data from the input are necessary:
subroutine
get_reference_grid_constants |
As pointed out in
Stefan Hackenbuchner's PhD thesis
the dislocations are relative to a reference frame.
This subroutine gets the grid constants for the reference frame (a_ref and c_ref)
Note: For zincblende it must be c_ref=a_ref. |
subroutine
strain_dirichlet_to_zero |
At least one point must be
set to a Dirichlet value for uM.
This subroutine returns for a point (i,j,k)
if the dislocations there are set to zero: uM(i,j,k,:) = 0d0
Note: This only corresponds to zero
strain if the lattice constants there are identical to those
of the reference material. |
subroutine
get_regions_strained |
It is possible to exempt
regions from the strain calculation: This subroutine tells one for any point
if it is within the simulation region or not.
On boundaries strain zero is assumed.
Numbering of 'octants' in 2D:
3 | 4
-----
1 | 2
Simulation area: i,j
1,n_2
n_1,n_2
--------------------------------|
|2
1|2
1|
|
|
|
|
|4
3|
|-
-|
|2
1|
|
|
|
|
|4
3|4
3|
--------------------------------|
1,1
n_1,1
|
subroutine
get_forces |
On each point forces through
internal interfaces can be specified by this subroutine. |
As described in
Stefan Hackenbuchner's PhD thesis
for each point (i,j,k)
a three-dimensional vector
uM(i,j,k,1:3)
is defined
which contains the dislocation in x-,y- and
z-direction. They are grouped in one whole vector
uV
with dimension 3*n_1*n_2*n_3.
So the numbering is according to the following rule:
point (i,j,k), component
l (l=1,2,3)
<==> ( (k-1)*n_1*n_2 +
(j-1)*n_1 + i - 1
)*3 + l in
uV
The differential equation -
div sigma = f (sigma stress tensor,
f force vector)
is discretized --> StrainMatrixM * uV = bV
The StrainMatrixM can be accessed via
subroutine strain_matvec and subroutine psolve_strain.
In a first step the matrix elements different
from zero are counted, later on the matrix is allocated and finally defined.
Setup of Strain Matrix
The matrix assembler loops through all points by
do k = 1,n_3
do j = 1,n_2
do i = 1,n_1
end do
end do
end do
For each point it calculates the discretized
differentials by the following scheme with subroutine
set_box:
Scheme for
3D:
(i ,j
,k-1) (i ,j-1,k ) (i-1,j ,k )
(i ,j ,k )
u1 u2 u3
u1 u2 u3 u1 u2 u3
u1 u2 u3
--------------------------------------------------------------
1
2 3 4 5 6
7 8 9 10 11
12
(i+1,j ,k) (i ,j+1,k) (i ,j
,k+1)
u1 u2 u3 u1 u2 u3 u1 u2 u3
------------------------------------------
13 14 15 16 17 18 19 20 21
(i-1,j-1,k )
(i+1,j-1,k ) (i-1,j+1,k ) (i+1,j+1,k )
u1
u2 u3 u1 u2 u3
u1 u2 u3 u1 u2
u3
---------------------------------------------------------------
22 23 24 25 26 27 28
29 30 31 32 33
(i-1,j ,k-1) (i+1,j ,k-1) (i-1,j ,k+1) (i+1,j ,k+1)
u1
u2 u3 u1 u2 u3
u1 u2 u3 u1 u2
u3
---------------------------------------------------------------
34 35 36
37 38 39 40
41 42 43 44 45
(i ,j-1,k-1)
(i ,j+1,k-1) (i ,j-1,k+1) (i ,j+1,k+1)
u1
u2 u3 u1 u2 u3
u1 u2 u3 u1 u2
u3
---------------------------------------------------------------
46 47 48
49 50 51 52
53 54 55 56 57
Scheme for 2D:
(i,j-1) (i-1,j)
(i,j) (i+1,j) (i,j+1)
(i-1,j-1) (i+1,j-1)
u1 u2 u3
u1 u2 u3 u1 u2 u3 u1 u2 u3
u1 u2 u3 u1 u2 u3 u1 u2 u3
-------------------------------------------------------------------------------------------
1 2 3 4 5 6 7 8 9 10 11 12
13 14 15 16 17 18 19 20 21
(i-1,j+1)
(i+1,j+1)
u1 u2 u3 u1 u2 u3
--------------------------
22 23 24 25 26 27
Since we have a tensor, the
divergence has three components, so we have for each point
a computational molecule for the first,
second and third component.
A computational molecule contains the
discretized differential operators according to the above scheme.
So the differential operator
for each point can be described by a matrix
discM(1..3,1..57)
and the right hand side
rhsV(1,..,3).
After
discM has been
calculated, the nonzero matrix elements are written to the matrix:
The positions of the points
in uV are calculated via the internal subroutine
get_position. (ind_ver=1,..3,
ind_hor=1,..,57).
Note: For each
point the differential operator is set up for each box (numbox=1,...,8)
if the box is within the simulation region.
Discretization is done by box
integration, so the differential stars of all relevant boxes are simply added.
For Dirichlet points a constant
(here: 1d10)
is simply added to the diagonal elements which corresponds to the condition
uM(i,j,k,1:3) = 0d0.
Setup of differential operator with subroutine
set_box
Input parameters are (i,j,k),
boxnum (1..8 =
number of octant)
Output is
discretize_boxM(1..3,1..57)
and rhsV(1..3):
The program first finds
the corresponding point on the material grid (i_guent)
and gets the elasticity matrix on that point (elastM)
in the crystal system
(i.e. with crystal axes as basis)
with subroutine get_elasticity. Then the tensor is transformed to the (x,y,z)-system
(=calculation system)
in subroutine
transform_elasticity.
NOTE: -----
check if rot or transpose(rot) ---------
Now box integration over
specified octant is performed. See
Stefan Hackenbuchner's PhD thesis for mathematical background.
It is:
vec{t} = vec{n}.tensor{sigma}
and
integral [
div sigma dV ] = integral [ vec{t} dA ]
(Gauss)
vec{n}
is the unit vector normal to surface A and dA is the scalar surface element.
vec{v} = (1 0 0)
--> vec{t} = (sigma_11 , sigma_12 , sigma_13)vec{v} = (0 1 0)
--> vec{t} = (sigma_21 , sigma_22 , sigma_23)vec{v} = (0 0 1)
--> vec{t} = (sigma_31 , sigma_32 , sigma_33)
So integration over box1 gives
(octant i-,j-,k-;
normal vectors -e_x, -e_y, -e_z):
integral [vec{t} dA] = - (sigma_11, sigma_12, sigma_13) *
h_2V(j-1)*h_3V(k-1)/4 -
(sigma_21, sigma_22, sigma_23) * h_1V(i-1)*h_3V(k-1)/4 -
(sigma_31, sigma_32, sigma_33) * h_1V(i-1)*h_2V(j-1)/4 +
+ integration over inner surfaces
When the integration over all eight boxes are added, the fluxes through the inner
surfaces cancel;
SO THE INTEGRAL OVER INNER SURFACES IS SET TO 0!
Note that the box integration
provides an equation that consists of three components which are the three
components of the differential operator.
The differential operators
sigma_ij are calculated by subroutine
set_stress.
This subroutine returns a
matrix discretize_stressV(1..3,1..57)
which contains the discretized differentials.
Note: discretize_stressV
contains
only values different from zero in the jth row for sigma_ij.
The relation holds:
sigma_i = elasticity_matrix_ij epsilon_j (Voigt notation)
So in subroutine set_stress
the programm loops through the relevant matrix elements of the elasticity matrix
and if it is nonzero the according differential is set.
This is done via subroutine set_epsilon.
In subroutine
set_epsilon
relevant components of the strain tensor epsilon are calculated. The finite difference scheme
for that differential operator is provided by subroutine
dui_dxj.
- 2D: subroutine
dui_dxj:
!
! Creates differential du_ii/dx_jj for point (i,j)
! and box number boxnum (=1..4).
!
! discretizeV(1..27) contains corresponding coefficients.
!
! int_dir ... direction for integration
! 1 --> surface perpendicular to i-direction
! 2 --> surface perpendicular to j-direction
!
! Note: Numbering of grid points that enter into the
computational molecule.
! For each grid point we have three components (x,y,z), thus we
obtain
! a vector of length 3*9=27.
!
! 1 <--> (i ,j-1)
! 2 <--> (i-1,j )
! 3 <--> (i ,j )
! 4 <--> (i+1,j )
! 5 <--> (i ,j+1)
! 6 <--> (i-1,j-1)
! 7 <--> (i+1,j-1)
! 8 <--> (i-1,j+1)
! 9 <--> (i+1,j+1)
!
!
!
! i-1,j+1 i,j+1 i+1,j+1
: 8
5 9
!
:
3 4
! i-1,j i,j i+1,j
:
2 3 4
!
:
1 2
! i-1,j-1 i,j-1 i+1,j-1
: 6
1 7
!
!
The adjoining octants of point (i,j) are shown in
red colour.
!---------------------------------------------------------------
! - left point
! (pos2,pos4) (pos1,pos3)
!------------------
! octants: 3 |
4
! -----
!
1 | 2
!---------------------------------------------------------------
======================= octant 1
===============================
----------- surface perpendicular to
i-axis ---------------
-- d/dx F --
= F(i,j) - F(i-1,j)
= 3 - 2
= pos2 - pos1
-- d/dy F --
= 0.5*[[F(i,j)-F(i,j-1)] + [F(i-1,j)-F(i-1,j-1)]]
=
3 - 1 +
2 - 6
pos2 pos1
pos4 pos3
----------- surface perpendicular to j-axis
---------------
-- d/dx F --
= 0.5*[[F(i,j)-F(i-1,j)] + [F(i,j-1)-F(i-1,j-1)]]
=
3 - 2 +
1 - 6
pos2 pos1
pos4 pos3
-- d/dy F --
= F(i,j) - F(i,j-1)]
= 3 -
1
= pos2 - pos1
======================= octant 2
===============================
----------- surface perpendicular to
i-axis ---------------
-- d/dx F --
= F(i+1,j) - F(i,j)
= 4 -
3
-- d/dy F --
= 0.5*[[F(i,j)-F(i,j-1)] + [F(i+1,j)-F(i+1,j-1)]]
=
3 - 1 +
4 - 7
----------- surface perpendicular to j-axis
---------------
-- d/dx F --
= 0.5*[[F(i+1,j)-F(i,j)] + [F(i+1,j-1)-F(i,j-1)]]
=
4 - 3 +
7 - 1
-- d/dy F --
=
======================= octant 3
===============================
----------- surface perpendicular to
i-axis ---------------
-- d/dx F --
=
-- d/dy F --
= 0.5*[[F(i,j+1)-F(i,j)] + [F(i-1,j+1)-F(i-1,j)]]
=
5 - 3 +
8 - 2
----------- surface perpendicular to j-axis
---------------
-- d/dx F --
= 0.5*[[F(i,j)-F(i-1,j)] + [F(i,j+1)-F(i-1,j+1)]]
=
3 - 2 +
5 - 8
-- d/dy F --
= F(i,j+1) - F(i,j)]
= 5 -
3
======================= octant 4
===============================
----------- surface perpendicular to
i-axis ---------------
-- d/dx F --
=
-- d/dy F --
= 0.5*[[F(i,j+1)-F(i,j)] + [F(i+1,j+1)-F(i+1,j)]]
= 5
- 3 + 9
- 4
----------- surface perpendicular to j-axis
---------------
-- d/dx F --
= 0.5*[[F(i+1,j)-F(i,j)] + [F(i+1,j+1)-F(i,j+1)]]
= 4
- 3 + 9
- 5
-- d/dy F --
=
Note: The dislocations
u1,u2,u3 refer to the
reference frame, so epsilon must be corrected.
This is accomplished in the
following way:
- get actual lattice constants
a
and
c
a
is valid in [100] and [010]
direction and c
is valid in [001] direction in
the crystal coordinate system.
- The coordinates of these
basis vectors of the
crystal system are expressed in the (x',y',z')
calculation system:
cryst_xV, cryst_yV, cryst_zV
cryst_xV = rotM.(1 0 0 )
[ rotM is
rotation matrix of material grid point: ---------- check if rot or rotT
--------------]
- The projection to the direction
in which the differential points are calculated for each:
x_cos=dot_product(vecV,cryst_xV)
y_cos=dot_product(vecV,cryst_yV)
z_cos=dot_product(vecV,cryst_zV)
where
vecV=(100)
for d/dx, ...
- This correction is a constant,
so it enters the right hand side of the equation, so:
x_project = x_cos*x_cos*(a-a_ref)/a
y_project = y_cos*y_cos*(a-a_ref)/a
z_project = z_cos*z_cos*(c-c_ref)/c
sum
= x_cos*x_cos + y_cos*y_cos + z_cos*z_cos
rhs
=
(x_project+y_project+z_project)/sum
- Finally the
differentials are scaled accordingly:
x_project = x_cos*x_cos*a_ref/a
y_project = y_cos*y_cos*a_ref/a
z_project = z_cos*z_cos*c_ref/c
total_project = (x_project+y_project+z_project)/sum
disc1V = disc1V*total_project
disc2V = disc2V*total_project
(disc1V and disc2V
contain the discretized differential)
Note:
- The computational
molecule
is calculated in
subroutine
dui_dxj
for NEGATIVE differentials, i.e. - div ( )
is discretized.
- The correction of the strain
enters the right hand side, so it must have POSITIVE sign.
|
