Empirical Potential Methods and Strain in Atomistic Simulations

Report
Outline
1.
2.
3.
4.
5.
6.
7.
8.
9.
Epitaxy in Semiconductor Crystal Growth
Elastic Description of Strain in Cubic Semiconductor Crystals
Atomistic Description of Strain
Molecular Statics and Force Fields
Keating's Valence Force Field
Stillinger-Weber Potential
Tersoff Potential
Simulation of Nanostructures
Piezoelectricity in Zincblende and Wurtzite crystals
TMCS III, Leeds 18th Jan 2012
Slide 1
Semiconductors
Group IV: Si, Ge, C
Group III-V: GaAs, InAs, AlAs, GaP, InP, AlP, GaN, InN, AlN, GaSb, InSb, AlSb
Group II-VI: CdSe, ZnSe, ZnS, CdS, MgSe, ZnTe
InN
TMCS III, Leeds 18th Jan 2012
Slide 2
Epitaxy
Epitaxy: (Greek; epi "above" and taxis "in ordered manner")
describes an ordered crystalline growth on a monocrystalline
substrate.
Homo-epitaxy (same layer and substrate material)
Hetero-epitaxy (different layer and substrate material).
In the Hetero-epitaxy case growth can be:
Lattice Matched: same, or very close, lattice constant of layer and
substrate e.g. GaSb/InAs or AlAs/GaAs
Lattice Mismatched: different lattice constant of layer and substrate
material e.g. InP/GaAs or InN/GaN.
TMCS III, Leeds 18th Jan 2012
Slide 3
Lattice Matched and Mismatched Epitaxy
Lattice Matched
TMCS III, Leeds 18th Jan 2012
Lattice Mismatched
Slide 4
Lattice Mismatched Epitaxy
In lattice mismatched heteroepitaxy the
layer material can be made to “adapt”
(can become smaller or larger) its in
plane lattice constant to match that of
the substrate (pseudomorphic growth).
Consequently volume conservation
(though volume is not perfectly
conserved) dictates that the lattice
constant in the growth direction needs
to become larger/smaller.
In this way the lattice periodicity is
maintained in the growth plane, but lost
in the growth direction.
Lattice Mismatched
TMCS III, Leeds 18th Jan 2012
Slide 5
Quantum Mechanics in action
Nanostructures: 2D Growth, 2D Growth + etching, 3D Growth
1D Quantum Well
2D Multi Quantum Wires
In0.52Al0.48As
AlAs
AlAs
In0.84Ga0.16As
In0.52Al0.48As
Taurino et al Mat Sci and Eng B, 67 (1999) 39
Growth direction
Scanning
Tunneling
Microscopy
TMCS III, Leeds 18th Jan 2012
Green: Free Carrier, Red: Confinement
Slide 6
Outline
1.
2.
3.
4.
5.
6.
7.
8.
9.
Epitaxy in Semiconductor Crystal Growth
Elastic Description of Strain in Cubic Semiconductor Crystals
Atomistic Description of Strain
Molecular Statics and Force Fields
Keating's Valence Force Field
Stillinger-Weber Potential
Tersoff Potential
Simulation of Nanostructures
Piezoelectricity in Zincblende and Wurtzite crystals
TMCS III, Leeds 18th Jan 2012
Slide 7
Elastic Strain
Semiconductors are produced by depositing liquid or gasses that
when coalesce and solidify follow the crystal structure of the “seed”,
usually a substrate of high crystalline quality.
During this deposition, often done in very small amounts (low growth
rate), as small as depositing one atomic layer at the time, if the layer
material has a bulk lattice constant larger than the substrate, then
the crystal will appear slightly deformed from its equilibrium state.
We refer to this material as “strained”.
z
zˆ
y
yˆ
xˆ
Unstrained
TMCS III, Leeds 18th Jan 2012
Strained
x
We chose the axes vectors
x,y,z arbitrarily, but need to
be linearly independent.
Note that while the axes
vectors are chosen to be
unitary (in units of the lattice
constant) in the unstrained
case, the strained axes are
not necessarily unitary.
Slide 8
zˆ
Elastic Strain
This picture is general
and valid for all types
of crystals, not just
simple cubic.
z
y
yˆ
xˆ
Unstrained
Strained
x
Unstrained and strained axis can be easily related:
x    1   xx  xˆ   xy yˆ   xz zˆ
y    yx xˆ  1   yy  yˆ   yz zˆ
z    zx xˆ   zy yˆ   1   zz  zˆ
TMCS III, Leeds 18th Jan 2012
The numerical coefficient εij
define the deformation of all the
atoms in the Crystal.
The diagonal terms εii control
the length of the axis, while the
off diagonal terms εij control the
angles between the axis.
Slide 9
Elastic Strain
This set of equations are in the
form of a mathematical entity
called Tensor.
x    1   xx  xˆ   xy yˆ   xz zˆ
The equations define the strained
position of any atom within the
crystal that upon strain moves
from R to R’.
z
zˆ
R
Unstrained
yˆ
R
y    yx xˆ  1   yy  yˆ   yz zˆ
z    zx xˆ   zy yˆ   1   zz  zˆ
Unstrained and strained
positions are written in terms of
the old and new axis:
y
R   xˆ   yˆ   zˆ
R    xˆ    yˆ    zˆ 
xˆ
TMCS III, Leeds 18th Jan 2012
Strained
x
Important: notice how the
coefficients α,β,γ are the same in the
unstrained and strained system
Slide 10
Elastic Strain
We now substitute the new axis with the expressions for the distortion:
R    xˆ    yˆ    zˆ 
x    1   xx  xˆ   xy yˆ   xz zˆ
y    yx xˆ  1   yy  yˆ   yz zˆ
z    zx xˆ   zy yˆ   1   zz  zˆ
After a little manipulation and taking into account the expression for R:
R  R 
 
xx
   yx   zx  xˆ     xy    yy   xz
 yˆ    
xz
   yz   zz  zˆ
Provided the original position and the distortion tensor are known, this
expression gives a practical way of calculating the position of any atom
inside a strained unit cell.
TMCS III, Leeds 18th Jan 2012
Slide 11
Strain Components
Often there is confusion between the terms strain and distortion. In
this lectures we follow the notation used in Jasprit Singh’s book, for
which the strain components eij are different from the distortion
components but related to them by:
e xx   xx
e yy   yy
e xy  x   y    yx   xy
e zz   zz
e yz  y   z    zy   yz
e xz  x   z    zx   xz
The final expressions for the off-diagonal terms eij are an
approximation in the limit of small strain.
Dilation: expresses how much the volume of the unit cell changes, and
in the limit of small strain is given by:
  e xx  e yy  e zz
Biaxial Strain: expresses how much the unit cell is strained in the z
direction compared to the x and y:
ebx  e zz 
1
2
e
xx
 e yy 
Uniaxial Strain: strain in one direction only, e.g. if eij= constant and
eii=0 then the strain is uniaxial in the [111]
TMCS III, Leeds 18th Jan 2012
Slide 12
Stress Components
Stress components: the force components (per unit area) that
causes the distortion of the unit cell.
There are 9 components:
Xx, Xy, Xz, Yx, Yy, Yz, Zx, Zy, Zz
Capital letters: direction of the force
Subscript: direction normal to the plane on which the stress is
applied (x is normal to yz, y is normal to xz, z is normal to xy, )
zˆ
xˆ
Xx
xy
TMCS III, Leeds 18th Jan 2012
yˆ
The number of independent
components reduces when we consider
that in cubic systems (like diamond or
zincblende) there is no torque on the
system (stress does not produce
angular acceleration).
Therefore Xy= Yx, Yz= Zy, Zx= Xz
And we are only left with 6:
Xx, Yy, Zz ; Yz, Zx, Xy
Slide 13
Elastic constants
The stress components are connected to the strain components via
the small strain elastic constants:
X x  c11 e xx  c12 e yy  c13 e zz  c14 e yz  c15 e zx  c16 e xy
Y y  c 21 e xx  c 22 e yy  c 23 e zz  c 24 e yz  c 25 e zx  c 26 e xy
Z z  c 31 e xx  c 32 e yy  c 33 e zz  c 34 e yz  c 35 e zx  c 36 e xy
Y z  c 41 e xx  c 42 e yy  c 43 e zz  c 44 e yz  c 45 e zx  c 46 e xy
Z x  c 51 e xx  c 52 e yy  c 53 e zz  c 54 e yz  c 55 e zx  c 56 e xy
X
y
 c 61 e xx  c 62 e yy  c 63 e zz  c 64 e yz  c 65 e zx  c 66 e xy
In practice we never have to deal with all 36 elastic constants.
First of all it is always the case that cij=cji which reduced the total to 21.
Second in real crystals, particularly cubic, the lattice symmetry reduces
the number even more.
Therefore in ZB we only have 3 independent constants: c11,c12,c44
In WZ there are 5: c11,c12,c13, c33, c44
TMCS III, Leeds 18th Jan 2012
Slide 14
Some more definitions
Elastic strain energy density for ZB:
U 

1
2
c11  e xx  e yy  e zz   c12  e xx e yy  e xx e zz  e zz e yy   c 44  e xy  e yz  e xz 
2
2
2
Bulk Modulus for ZB:
2
B 
2
 c11  2 c12 
3
Shear Constant for ZB:
C 
 c11  c12 
2
TMCS III, Leeds 18th Jan 2012
Slide 15
2

Properties of Semiconductors
ZB
Si
Ge
C
Ga-As
In-As
Al-As
Ga-P
In-P
Al-P
Ga-N
In-N
Al-N
Ga-Sb
In-Sb
Al-Sb
a
(Ǻ)
5.431
5.658
3.567
5.653
6.058
5.662
5.451
5.869
5.463
4.500
4.980
4.380
6.096
6.479
6.135
TMCS III, Leeds 18th Jan 2012
B
(Mbar)
0.980
0.713
0.442
0.757
0.617
0.747
0.921
0.736
0.886
2.060
1.476
2.030
0.567
0.476
0.855
C’
(Mbar)
0.502
0.410
0.478
0.364
0.229
0.288
0.440
0.269
0.329
0.825
0.424
0.698
0.270
0.183
0.414
c11
c12
c44
(Mbar) (Mbar) (Mbar)
1.660 0.640 0.796
1.260 0.440 0.677
10.79 1.24
5.78
1.242 0.514 0.634
0.922 0.465 0.444
1.131 0.555 0.547
1.507 0.628 0.763
1.095 0.556 0.526
1.325 0.667 0.627
3.159 1.510 1.976
2.040 1.190 1.141
2.961 1.565 2.004
0.927 0.378 0.462
0.720 0.354 0.341
1.407 0.579 0.399
Slide 16
Strain in Lattice Mismatched Epitaxy
Poisson ratio: is a measure of the tendency of materials to
stretch in one direction when compressed in another. This ratio
depends on the substrate orientation and the type of crystal. For
cubic crystals including ZB:
 
c11
 
[001]
2 c12
c11  2 c12  c 44
c11  c12  4 c 44
[111]
Strain: in pseudomorphic growth one can consider, independent of
the substrate orientation, strain to have only two components, one
parallel to the growth plane and one perpendicular.
e 
a Substrate
a Layer
1
e  
e

Important: in [001] growth: e = exx= eyy and e= ezz
TMCS III, Leeds 18th Jan 2012
Slide 17
Strain in [111] pseudomorphically
grown layers
Important: in [111] growth the combination of e and e results in a
strain tensor with exx = eyy = ezz and exy = exz = eyz
The distortions in this case are:
 xx   yy   zz
 2 1 2 c11  4 c12  4 c 44 
 
 e //
 3 3 c11  2 c12  4 c 44 
 xy   yz   xz
 1 1 2 c11  4 c12  4 c 44 
  
 e //
 3 3 c11  2 c12  4 c 44 
z
[111]
y
(1,1,1)
x
TMCS III, Leeds 18th Jan 2012
Important: the distortions are expressed
in the basis system where x, y and z are
aligned with the [100], [010] and [001]
directions.
Instead e and e are defined so that they
relate to strain in the (111) plane and the
[111] direction, respectively.
Slide 18
Outline
1.
2.
3.
4.
5.
6.
7.
8.
9.
Epitaxy in Semiconductor Crystal Growth
Elastic Description of Strain in Cubic Semiconductor Crystals
Atomistic Description of Strain
Molecular Statics and Force Fields
Keating's Valence Force Field
Stillinger-Weber Potential
Tersoff Potential
Simulation of Nanostructures
Piezoelectricity in Zincblende and Wurtzite crystals
TMCS III, Leeds 18th Jan 2012
Slide 19
Tetrahedral Bonding
In the Zincblende crystal, just like in the diamond one, atoms
bond together to form tetrahedrons.
Hence the individual atomic orbitals merge to form sp3 hybrid orbitals
TMCS III, Leeds 18th Jan 2012
Slide 20
Wurtzite
While Zinblende is the preferred crystal structure of III-As, III-P and
III-Sb, III-N tend to crystallize preferentially in hexagonal form. The
hexagonal crystal with a two atom basis consisting of cations and
anions is called Wurtzite.
View from the top
Perspective View
Zincblende
Wurtzite
Two adjacent tetrahedrons overlap in the z direction in WZ but not in ZB.
Hence second nearest neighbours in WZ are actually closer than in ZB at
equilibrium. The modified inter-atomic forces result in a slight reduction of
the interatomic distance between the first nearest neighbours.
TMCS III, Leeds 18th Jan 2012
Slide 21
The 7th elastic parameter
Is a description based on 6 strain components enough to describe
all deformations in a ZB or WZ crystals?
The distance that the atom is
displaced by is characterized by
the Kleinman parameter
z
a 
4
With a the lattice constant and γ
the shear strain.
This results in a crystal where the
atomic bonds are not all of the
same length.
Strain in the [111]
TMCS III, Leeds 18th Jan 2012
+
+
+
+
Only 3
identical
sp3
orbitals
Slide 22
Strain from atomic positions
Given the 5 coordinates of the atoms in a tetrahedron how do we
reverse engineer the strain?
 1   xx

  yx
 
zx

 1   xx

  yx
 
zx

 xy
 xz   x1 
1   yy
 yz
 zy
1   zz
 x1' 
   ' 
  y1    y1 
 z   z' 
 1   1 
 1   xx

  yx
 
zx

 xy
 xz   x 0 
1   yy
 yz
 zy
1   zz
 xy
 x 0' 
   ' 
  y0    y0 
 z   z' 
 0   0 
 xz   x 2 
1   yy
 yz
 zy
1   zz
 x 2' 
   ' 
  y2    y2 
 z   z' 
 2   2 
This become a simple system of linear equations easily solvable.
The solution gives the 6 components of the strain tensor.
However the deformation on the position of the yellow atom, dependant
on the Keinman parameter, is still undetermined and requires a separate
calculation.
TMCS III, Leeds 18th Jan 2012
Slide 23
The issue of local/global composition
Furthermore strain is a relative property (variation of e.g. bond
length compared to an initial state).
Microscopists refer to strain as difference in the bond lengths
compared to the host.
Theorists think of strain as deformation of a material from its bulk
state.
Everyone else does not usually know what they are talking about!!
If dealing with an alloy and if wanting to take the theorist approach,
one needs to know what the lattice constant of the alloy is.
But what does composition mean?
It makes sense for a large uniform block, not for non uniform.
We take the approach of counting atoms up to second nearest
neighbour form the centre of the tetrahedron
TMCS III, Leeds 18th Jan 2012
Slide 24
Outline
1.
2.
3.
4.
5.
6.
7.
8.
9.
Epitaxy in Semiconductor Crystal Growth
Elastic Description of Strain in Cubic Semiconductor Crystals
Atomistic Description of Strain
Molecular Statics and Force Fields
Keating's Valence Force Field
Stillinger-Weber Potential
Tersoff Potential
Simulation of Nanostructures
Piezoelectricity in Zincblende and Wurtzite crystals
TMCS III, Leeds 18th Jan 2012
Slide 25
Modelling Strain in Real Structures
Because of its impact on the electronic properties strain in
semiconductor nanostructures always needs to be evaluated with
the highest possible accuracy.
Measurements (usually involving electron microscopy analysis) are
not usually sufficiently accurate, so modelling is the only viable
alternative.
Simple elasticity formulas are acceptable when dealing with
standard cases where strains are uniform or approximating strains
as uniform is acceptable, e.g. a simple quantum well.
They become useless however in complex quantum well
structures, wires and dots where strains are non uniform.
In time several methods have been developed ranging from
continuum, finite element, analytic and atomistic.
Atomistic methods are now widely used for quantum dots while
continuum methods are the preferred methods for quantum wells.
TMCS III, Leeds 18th Jan 2012
Slide 26
Molecular Dynamics
• Molecular Dynamics is a computer simulation in which a starting set
of atoms or molecules is made to interact for a period of time
following the laws of Physics (e.g Newton’s Laws).
• In Semiconductor science one can build an atomistic model of a
strained crystal but if the strain is not known a priori then atoms are
not going to be in their equilibrium positions.
• Then their motion paths are dictated by the “force field” generated
by the potential of the solid.
TMCS III, Leeds 18th Jan 2012
Slide 27
Molecular Dynamics
Initial Position of the atoms r0i
Potential V (r0i)
Forces F=-grad V (r0i)
Velocities and acceleration
Evaluate the positions after Δt
Repeat till Forces are low
Often the simulation does not require very large atomic motion.
For instance for calculating strain one might only want to allow small
atomic displacements from the crystal structure, without atom switching.
When Energy minimisation is the fundamental criterion and forces are
used to direct the geometry optimisation rather than predicting the final
positions, we are using a “Molecular Statics” simulation.
TMCS III, Leeds 18th Jan 2012
Slide 28
Outline
1.
2.
3.
4.
5.
6.
7.
8.
9.
Epitaxy in Semiconductor Crystal Growth
Elastic Description of Strain in Cubic Semiconductor Crystals
Atomistic Description of Strain
Molecular Statics and Force Fields
Keating's Valence Force Field
Stillinger-Weber Potential
Tersoff Potential
Simulation of Nanostructures
Piezoelectricity in Zincblende and Wurtzite crystals
TMCS III, Leeds 18th Jan 2012
Slide 29
Valence Force Field
The “force field” that is generated by the potential of the atoms in
the solid can be represented as a 3 body potential.
In the Keating's Valence Force Field: The Potential is the sum of the
potential energy between the
pairs of atoms i and j (two
VV F F   V 2 ( R i  R j )   V 3 (ˆijk )
ij
ij
body), plus a term that depends
on the angle between i,j and a
nn
3 ij
2
1
third atom k (three body).
2
0 2


V2   
(
R

R
)

(
d
)
i
j
ij
0 2

Rj
2
8  (d ) 
i
V3 
j
nn
1

2
i
ij
3  i , jk
 ( R j  R i )  ( R k  R i )  ( d ) cos  0 


j ,k i 8  ( d )
0
ij
2
0
ij
2
2
jki
j
i
Ri
k
k
k
dij0 is the unstrained bond length of atoms i and j and 0 is the
unstrained bond angle (e.g. for zinc-blend cos0=-1/3), and ijk is the
angle between atoms i, j and k.
The local chemistry is contained in the parameters  and  , which are
fitted to the elastic constants
P.N. Keating, Phys. Rev. 145, 637 (1966)
TMCS III, Leeds 18th Jan 2012
Rk
Slide 30
Valence Force Field
The VFF is widely used for all types of nanostructures.
VFF is basically a parabolic approximation to the potential of solids
V(R)
Ω is the volume
occupied by one
atom
R0
2
Uniform: same
distortion in
x,y and z
B 
R
1 d E
 dv
2
E  Ecoh
1
3
R  v R
1
Binding Energy
The main limitation is that there
are only 2 parameters ( and  ) but
3 elastic constants even for
Zincblende!!!
Non Uniform: z
stretch, x,y compress
(by the same amount)
and viceversa
2
C '
1 d E
 d
c11 
a

 3

c12 
c 44 
TMCS III, Leeds 18th Jan 2012
1
a
4



1
R y  R y / (1   )
1
Rz  Rz

a   

E  Ecoh
R x  R x (1   )
1
1
2
Slide 31
Progress in Valence Force Field
Anharmonicity correction:
• Ability to reproduce anharmonic effects is linked to the quality of prediction
of the phonon spectrum.
• Some progress has been presented (e.g. Lazarenkova et al, Superlattices
and Microstructures 34, 553 (2003)).
• Not clear why phonon frequencies, elastic constants and mode Grüneisen
parameters are not correlated (Porter et al J. Appl. Phys. 81, 96 (1997)).
• For Ionicity in Zincblende to solve this problem check recent P. Han and G.
Bester, Phys. Rev. B 83 174304 (2011)
Distance between 2
ions, one of which is in
the central cell
U 
c
1
2
Z  Z e

l 
2
l 0 , 
0
rs
Ionicity and Wurtzite:
• Empirical potentials were historically developed for Si and Ge (pure covalent
bonds)
• III-V are mainly covalent, partially ionic. II-VI are both covalent and ionic
• Only for infinite crystals or systems were the charge is uniformly distributed
this it’s not a big deal.
• Important in III-N WZ (Grosse and Neugebauer, PRB 63, 085207 (2001)),
and can be incorporated following Ewald summation scheme (codes available).
• Also check Camacho et al (Physica E, Vol. 42, p. 1361 (2010)) “application
of Keating’s valence force field to non–ideal wurtzite materials”
TMCS III, Leeds 18th Jan 2012
Slide 32
Outline
1.
2.
3.
4.
5.
6.
7.
8.
9.
Epitaxy in Semiconductor Crystal Growth
Elastic Description of Strain in Cubic Semiconductor Crystals
Atomistic Description of Strain
Molecular Statics and Force Fields
Keating's Valence Force Field
Stillinger-Weber Potential
Tersoff Potential
Simulation of Nanostructures
Piezoelectricity in Zincblende and Wurtzite crystals
TMCS III, Leeds 18th Jan 2012
Slide 33
Stillinger-Weber
The “force field” that is generated by the potential of the atoms in
the solid can be represented as a 3 body potential.
In the Stillinger-Weber potential:
The Potential is the sum of the
2
nn
nn
potential energy between the
1

V SW   V 2 ( rij )   V 3 ( rij )   cos  ijk  
pairs of atoms i and j (two
3


ij
ijk
body), plus a term that depends
1
nn
on the angle between i,j and a
1
 rij  a 
p
q
V 2    A  B rij  rij   e
third atom k (three body).
2 i j
Rj
j
This in an adaptation of the well known Lennard-Jones
potential used for liquefied noble gasses.
jki
Ri
i
2
1
1 

nn
Rk
   rij  a     rik  a   
1
1



V3     e
  cos  ijk  
k
k
2 i j ,k  i
3

k
This potential works very well for Si in diamond structure where the bond
angle cos0=-1/3.
The local chemistry is contained in the parameters A, B , p, q , a, λ and γ
which are fitted to various material properties.
F. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985)
TMCS III, Leeds 18th Jan 2012
Slide 34
Stillinger-Weber
The SW is not as widely used as VFF, but it has his niche
(thermodynamics of Si mainly).
In a way it should perform much better than VFF as it is not a
parabolic approximation to the potential of solids.
Parameterisations take into account the crystal phase diagram and
check that diamond is the lowest energy structure
Works reasonably well for diamond-Si but not for other crystal structures.
TMCS III, Leeds 18th Jan 2012
Slide 35
Outline
1.
2.
3.
4.
5.
6.
7.
8.
9.
Epitaxy in Semiconductor Crystal Growth
Elastic Description of Strain in Cubic Semiconductor Crystals
Atomistic Description of Strain
Molecular Statics and Force Fields
Keating Valence Force Field
Stillinger-Weber Potential
Tersoff Potential
Simulation of Nanostructures
Piezoelectricity in Zincblende and Wurtzite crystals
TMCS III, Leeds 18th Jan 2012
Slide 36
Tersoff Potential
The “force field” that is generated by the potential of the atoms in
the solid can be represented as a 3 body potential.
In the Tersoff potential:
The Potential is the sum of the
potential energy between the
   ij ( rij  re ) 
   ij ( rij  re ) 
pairs of atoms i and j (two body),
V ij  Aij e
 bij B ij e
multiplied times a term (bij) that
depends on the angle between i,j


and a third atom k (three body).
bij  1    ij  ij 


Rj
 1
2 ni
ni
The expression for bij (known as bond order) is written
as to emulate the atomic coordination number Z.
Hence ζ is sometimes called the pseudo-coordination.
 ij 

jki
j
i
k
k i, j
g(θ) and ω describe the angular and radial forces dependence.
g ( ijk )  1 
 
ci
di
2
2

ci
d i  ( h i  cos  jki )
2
TMCS III, Leeds 18th Jan 2012
2
 ijk  e
Rk
k
k
f c ( rik ) g ( ijk )   ijk
Ri
  3 ( rij  rik ) 3 


Slide 37
Tersoff Potential
g ( ijk )  1 
jki
 
ci
di
2
2

ci
d i  ( h i  cos  jki )
2
 ijk  e
2
ω
g
j
  3 ( rij  rik ) 3 


i
k
k
k
θeq
θ
  3 ( rij  rik ) 3 


angular forces: resistance to bend
radial forces: resistance to stretch
• When fitting to Bulk Modulus g(θ) is always g(θeq) and ωijk==1
• When fitting to Shear Constant g(θ)≠ (θeq) but ωijk==1
• When fitting c44 then both g(θ) ≠ (θeq) and ωijk ≠ 1
• Hence the Kleinman parameter links angular and radial forces!!!
TMCS III, Leeds 18th Jan 2012
Slide 38
Tersoff Potential
This potential describes covalent bonding and works very
well for different crystal structures for group IV and despite
the partial ionicity of the bond, group III-V.
The local chemistry is contained in the parameters A, B , re,
α, β , γ, c, d, h, n and λ, which are fitted to various material
properties.
J. Tersoff, Phys Rev Lett 56, 632 (1986) & Phys Rev B 39, 5566 (1989)
Sayed et al, Nuclear Instruments and Methods in Physics Research 102, 232 (1995)
TMCS III, Leeds 18th Jan 2012
Slide 39
Tersoff Potential
The TP is not as widely used as VFF, but its use is rapidly increasing
as parameterizations are improved.
Again it should perform much better than VFF as it is not a parabolic
approximation to the potential of solids.
As there are many parameters, parameterisations can take into
account many things, including the crystal phase diagram, all the
cohesive and elastic properties and many more.
2
c 44 
1 d E coh

d
2
E  E coh
R x  R x   R y
R y  R y
R z  R z
Works rather well for zincblende and diamond group IV and III-V but it is
not yet optimized for thermodynamic and vibrational properties.
D. Powell, M.A. Migliorato and A.G. Cullis, Phys. Rev. B 75, 115202 (2007)
TMCS III, Leeds 18th Jan 2012
Slide 40
Progress in Tersoff
The Kleinman parameter
• The many parameters need putting to good use.
• Kleinman deformation is critical because expresses the balance between
radial and angular forces (Powell et al PRB 75, 115202 (2007))
hydrostatic distortion 
0.00
sublattice displacement 
0.9
0.02
0.04
0.06
0.08
0.10
circles: InAs
squares: GaAs
DFT
0.8
0.7
0.6
Tersoff
0.5
0.4
0.3
DFT
Range of physical shear strains
0.2
0.00
0.02
0.04
0.06
shear distortion
0.08
0.10
Tersoff
Ionicity and Phonons
• Ionicity, like VFF, is missing.
• Crystal growth only possible if ionic contribution is included (Nakamura et
al J. Cryst. Growth 209, 232 (2000)
• Phonons are still independent of elastic constants (Powell et al, Physica E
32, 270th(2006)
TMCS III, Leeds 18
Jan 2012
Slide 41
Beyond Tersoff: bond order potentials
Π versus σ –bonding
• Tersoff neglects Π–bonding. Is it of consequence?
• Tersoff can to some extent reproduce surface
reconstruction energies (Hammerschmidt, PhD thesis)
Beyond σ -bonding
 1


bij  1    ij  ij 


ni
g ( ijk )  1 
 
ci
di
2 ni
 ij 

f c ( rik ) g ( ijk )   ijk
k i, j
2
2

ci
d i  ( h i  cos  jki )
2
2
 ijk  e
  3 ( rij  rik ) 3 


• It is generally possible to rewrite the bij with expressions directly obtained
from tight binding. (D.G. Pettifor, “Many atom Interactions in Solids”,
Springer Proceedings in Physics 48, 1990, pag 64))
• In this way the “bond order” can be explicitly obtained analytically to any
order (Murdick et al, PRB 73, 045206 (2006)).
• The second moment approximation is essentially equivalent to Tersoff
(Conrad and Scheerschmidt, PRB 58, 4538 (1998))
TMCS III, Leeds 18th Jan 2012
Slide 42
Outline
1.
2.
3.
4.
5.
6.
7.
8.
9.
Epitaxy in Semiconductor Crystal Growth
Elastic Description of Strain in Cubic Semiconductor Crystals
Atomistic Description of Strain
Molecular Statics and Force Fields
Keating's Valence Force Field
Stillinger-Weber Potential
Tersoff Potential
Simulation of Nanostructures
Piezoelectricity in Zincblende and Wurtzite crystals
TMCS III, Leeds 18th Jan 2012
Slide 43
General Tips for MD
Building Models:
• If possible try and use existing software
• Try and guess final positions: it saves a lot of computational time
Empirical Potentials:
• Codes that use VFF, SW and Tersoff are usually freely available!
• IMD (Stuttgart), CPMD (IBM-Zurich) are parallel (for running on
clusters) and open source
• Nemo3 (Purdue) uses VFF
• Always check what version of the potentials are being used!!
Molecular Statics:
• Make sure that the parameters that control the length of time the
simulation is running for are set to reasonable values
• Build your simulation up in size to see what you can get away
with in terms of system sizes and check that results do not
depend on the size chosen
Strain:
• Good strain algorithms exist and are freely available
• If you write your own you need a nearest neighbour list. Usually
MD produces one
Gridding:
• Strain is first obtained onto the atomic grid. Then to use it often
it needs converting to an ordered grid. One can use various
methods like Gaussian smoothing or weighted average.
TMCS III, Leeds 18th Jan 2012
Slide 44
MD of QDs using Tersoff Potential
After MD
Before MD
300
300
-0.06500
-0.06690
250
-0.05500
250
-0.05690
-0.04500
-0.03500
-0.04690
200
200
-0.04500
-0.02500
-0.03690
-0.006900
150
-0.02690
[001] Å
[001] Å
0.01500
-0.01500
150
0.005000
-0.05500
-0.05500
-0.01690
100
100
-0.02500
-0.005000
-0.005000
-0.01500
-0.03500
0.005000
0.01500
-0.006900
0.02500
0.005000
50
εxx
0.03500
50
2.000E-4
0.04500
50
100
150
200
250
300
50
350
100
150
[010] Å
200
250
300
350
[001] Å
300
300
0.001500
-0.09850
-0.06690
-0.08850
250
250
-0.07850
-0.05690
[001] Å
-0.04690
-0.04690
-0.05690
-0.01690
-0.006900
150
-0.03850
-0.04850
-0.03690
-0.02690
-0.02690
-0.03690
100
-0.06850
-0.05850
200
[001] Å
Floating
200
0.01150
-0.01690
-0.04850
-0.03850
150
εyy
-0.02850
-0.008500
-0.01850
-0.02850
-0.01850
100
-0.008500
0.001500
-0.006900
0.01150
50
50
2.000E-4
0.02150
0.02900
0.001500
50
100
150
200
250
300
350
50
100
150
[010] Å
200
250
300
350
[010] Å
300
300
-0.09800
-0.008000
-0.02175
-0.08800
-0.01175
-0.07800
-0.001750
250
250
0.008250
PBC
-0.03800
0.02825
200
-0.02800
0.002000
0.01200
0.02200
0.06825
0.07825
[001] Å
[001] Å
0.04825
150
150
0.08825
0.05825
100
0.1082
100
-0.01800
-0.008000
0.002000
0.01200
0.02200
0.03200
0.03200
0.04200
0.09825
-0.03800
-0.02800
200
0.03825
0.05825
-0.05800
-0.04800
0.01825
0.05825
-0.06800
-0.01800
0.04200
-0.008000
-0.01800
0.07200
0.08200
0.1282
0.1382
Fixed
50
100
150
200
[010] Å
TMCS III, Leeds 18th Jan 2012
50
0.002000
0.002000
0.1482
0.1500
250
300
350
0.05200
0.06200
0.1182
50
0.04200
0.09200
0.1020
0.1120
0.1220
0.002000
50
100
0.1320
150
200
250
300
350
[010] Å
Slide 45
0.1420
0.1500
εz
z
Outline
1.
2.
3.
4.
5.
6.
7.
8.
9.
Epitaxy in Semiconductor Crystal Growth
Elastic Description of Strain in Cubic Semiconductor Crystals
Atomistic Description of Strain
Molecular Statics and Force Fields
Keating's Valence Force Field
Stillinger-Weber Potential
Tersoff Potential
Simulation of Nanostructures
Piezoelectricity in Zincblende and Wurtzite crystals
TMCS III, Leeds 18th Jan 2012
Slide 46
Kleinman Parameter
2
c 44 
1 d E coh

d
2
E  E coh
R x  R x   R y
R y  R y
R z  R z
The distance that the atom is displaced by is
characterized by the Kleinman parameter
z
a 
4
With a the lattice constant and γ the shear strain.
This results in a crystal where the atomic bonds
are not all of the same length.
TMCS III, Leeds 18th Jan 2012
Slide 47
Piezoelectricity
In the case of a uniaxial distortion the displacement is in the [111]
direction, and can still be characterized by the Kleinman parameter


R x  R x 
R y 
R z 

2

2
2
Ry 
2

Rx  R y 
Rx 
4
identical
sp3
orbitals

2
2
Rz
xy z 
Rz
4
R y  Rz
+
a 
Strain in the [111]
-
+
+
+
+
+
+
+
Only 3
identical
sp3
orbitals
The displacement of cations relative to anions in III-V semiconductors
results in the creation of electric dipoles in the polar direction which in
ZB is the direction that lacks inversion symmetry.
TMCS III, Leeds 18th Jan 2012
Slide 48
Piezoelectricity in Zincblende
The effect can be quantified by writing a general expression for the
polarization as a function of the so called “piezoelectric coefficients”
and the distortion components.
Px 

e xkl  kl
Convention is:
xx=1, yy=2, zz=3, yz=4, zx=5, xy=6

e ykl  kl
In ZB, for symmetry, the only non zero
coefficients are e14= e25= e36

e zkl  kl
k ,l  x , y , z
Py 
k ,l  x , y , z
Pz 
k ,l  x , y , z
P  e14    yz   zy   zx   xz   xy   yx 
In actual fact this picture is incomplete as only includes coefficients
linked to linear terms in the strain.
In the past 6 years the importance of including also coeffiecients linked
to quadratic terms in the strain (e.g. xx2 or xy xz ) has been
highlighted (so called non linear or second order Piezo effect).
• M.A. Migliorato et al, Phys. Rev. B 74, 245332 (2006)
• L. C. Lew Yan Voon and M. Willatzen, J. Appl. Phys. 109, 031101 (2011) REVIEW
• A. Beya-Wakata et al, Phys. Rev. B 84, 195207 (2011)
TMCS III, Leeds 18th Jan 2012
Slide 49
Piezoelectricity in Wurtzite
Zincblende
Wurtzite
Spontaneous
polarization
Strain induced polarization
Quadratic terms in the strain (e.g. xx2 or xy xz ) are also important.
There is still some controversy between the early accepted values of
mainly for the spontaneous polarization coefficients
• L. C. Lew Yan Voon and M. Willatzen, J. Appl. Phys. 109, 031101 (2011) REVIEW
• J. Pal et al, Phys. Rev. B 84, 085211 (2011)
TMCS III, Leeds 18th Jan 2012
Thank you!!!
Acknowledgments:
Joydeep Pal , Umberto Monteverde, Geoffrey Tse, Vesel Haxha, Raman Garg (University of Manchester)
Dave Powell (University of Sheffield)
GP Srivastava (University of Exeter)
TMCS III, Leeds 18th Jan 2012
Slide 51

similar documents