Molecular dynamics thesis

i
=
f
j
(r
ij
)
j
=
1: j
≠
i
N
∑
,
(4)
where f
j 
is the electron density contributed by atom j  at the position of atom i.
Based on a comparison between different EAM models made by Morishita and De la
Rubia [4] for bcc vanadium the Johnson-Oh EAM [5] potential was chosen. Johnson and
Oh have proposed the following general expressions for the embedding function, the
electron density distribution, and the pair potential for a number of bcc metals:
F( )
= −
E
C
−
E
1V
UF
(
)
1
−
ln
e
 
 
  
 
 
n
 
 
 
 
 
 
e
 
 
  
 
 
n
,
(5)
f (r)
=
f
e
r
1e
r
 
 
 
 
,
(6)
(r)
=
K
3
r
r
1e
−
1
 
 
 
 
 
 
3
+
K
2
r
r
1e
−
1
 
 
 
 
 
 
2
+
K
1
r
r
1e
−
1
 
 
 
 
 
 
1
+
K
0
,
(7)
in which
n
=
9
Ω
B
−
15
Ω
G
2
E
C
−
E
1V
UF
(
)
,
(8)
K
3
= −
15
Ω
G
3A
+
2
3
4
A
−
1
2
S
 
 
 
 
,
(9)
K
2
= −
15
Ω
G
3A
+
2
15
16
A
−
3
4
S
−
9
8
AS
+
7
8
 
 
 
 
,
(10)
K
1
= −
15
Ω
G
7
8
−
3
4
S
 
 
 
 
,
(11)
K
0
= −
15
Ω
G
3A
+
2
201
112
A
−
27
28
S
−
87
56
AS
+
187
168
 
 
 
 −
1
7
E
1V
UF
,
(12)
in which
S
=
r
2e
r
1e
,
(13)
A
=
2C
44
C
11
−
C
12
,
(14)
G
=
3C
44
+
C
11
−
C
12
.
(15)

---

12
In these equations E
C 
is the cohesive energy, E
1V
UF
is the unrelaxed monovacancy formation
energy, 
e
is the (dimensionless) equilibrium electron density, f
e 
is a dimensionless factor
that cancels out for single-element potentials, r
1e
and r
2e
are the first and second equilibrium
neighbour distances, 
is a constant (of little influence, a value of 6 was chosen by Johnson
and Oh), is the atomic volume, B  is the bulk modulus, G is the Voigt average shear
modulus, and C
11
, C
12
, and C
44
are the cubic elastic constants.
The cohesive energy, equilibrium lattice constant a , the three elastic constants, and
the unrelaxed monovacancy formation energy are the input parameters for the potential. The
values used for constructing a bcc Mo-Mo EAM potential are given in table I [5]. Note that
for the lattice constant a more accurate value is used that differs slightly from the value of
3.150 Å used by Johnson and Oh.
Table I. Values used for constructing the EAM potential for bcc molybdenum. The values are the lattice
parameter a, atomic volume 
, the bulk modulus B, the Voigt average shearmodulus G, the anisotropy
ratio A, the cohesive energy E
C
, and the unrelaxed monovacancy formation energy 
E
1V
UF
.
quantity
a (Å)
3.1472
Ω
B (eV)
25.68
Ω
G (eV)
12.28
A
0.78
E
C
(eV)
6.81
E
1V
UF
(eV)
3.1
The embedding function of eqn. (5) is very steep near 
ρ
=0, the derivative tends
to - . As a result of discretization errors (see also section 2.3.4) this singularity can produce
unreliable results. To alleviate this problem, the embedding function has been multiplied by
the function
h( )
=
1
−
e
−
1
2
s
 
 
 
 
 
 
 
 
2
,
(16)
in which 
ρ
s
was given the value 0.5. This fixes the problem with the singularity. There is no
theoretical basis for doing this, but only the low-density part of the embedding function is
significantly affected. The embedding energy for low densities is small, consequently
reducing it further or increasing it by a small fraction has no significant impact. The Mo
embedding function is shown in fig. 1.
The electron density function has also been altered slightly from eqn. (6) to limit the
radial range and thereby the number of interactions. This is done by substituting for large
distances r
s
 < r < r
c
the density function f(r) by the function f
c
(r) given by
f
c
(r)
=
−
2 f(r
s
)
+
(r
s
−
r
c
) f'( r
s
)
(r
s
−
r
c
)
3
(r
−
r
c
)
3
+
3 f (r
s
)
−
(r
s
−
r
c
) f '(r
s
)
(r
s
−
r
c
)
2
(r
−
r
c
)
2
.
(17)
This function has the properties f
c
(r
s
) = f(r
s
), f
c
’(r
s
)  = f’(r
s
), f
c
(r
c
) = f
c
’(r
c
) = 0 . Therefore
the corrected function smoothly connects to the original function at r
s 
and then gradually

---

13
Figure 1. Embedding function F  versus electron density for Mo.
The solid curve represents the embedding function used in simulations.
The dashed curve is the unmodified embedding function (eqn. (5)).
goes to zero at r
c
. For r
s
and r
c
the values r
2e
+ 0.1*(r
3e
- r
2e 
) and r
2e
+ 0.5*(r
3e
- r
2e 
) were
chosen, in which r
3e
is the third equilibrium neighbour distance. This limits the number of
interactions to an average of 14, the sum of the number of nearest and next-nearest
neighbours. The pairpotential has been treated in a similar manner. A second modification
of the electron density function, at short distances, will be discussed later.
Since the Johnson-Oh EAM is derived for low energy situations, different potentials
are used for short distances. For distances smaller than 2.70 Å the Johnson-Oh pair
potential is stiffened. Instead of (r ) one should use 
S
(r ), given by
S
r
( )
=
r
( )
+
4.5 1
+
4
A
−
0.1
 
 
 
 
r
( )
−
r
1e
( )
(
)
r
r
1e
−
1
 
 
 
 
 
 
2
.
(18)
For distances smaller than 1.59 Å the Screened Coulomb pair potential with Molière weight
factors and Firsov screening length [6] is used. This purely repulsive potential has the form
(r)
=
Z
1
Z
2
e
2
4
0
r
r
a
 
 
 
 
,
(19)
in which
r
a
 
 
 
 =
c
i
e
−
d
i
r
a
1
3
∑
,
c
i
=
(0)
=
1
1
3
∑
,
(20)
and

---

14
a
=
9
2
128
 
 
 
 
 
1
3
a
B
Z
1
1
2
+
Z
2
1
2
 
 
 
 
 
 
−
2
3
,
(21)
in which Z
1 
and Z
2
are the atom numbers of the interacting particles at distance r, e is the
elementary charge, 
ε
0
is the vacuum permittivity, is the screening function, a  is the
screening length, c
i
and d
i
are constants, and a
B
is the Bohr radius. The Molière values for
c
i
and d
i
are given in table II.
Table II: Molière constants in the screening function.
i
c
i
d
i
1
0.35
0.3
2
0.55
1.2
3
0.10
6.0
Since the EAM potential is based on a combination of a pair potential and an electron
density and the Screened Coulomb potential is a pure pair potential, a method has to be
devised to smoothly go from one potential to the other. In other words, the electron density
function should fall off to zero when the distance from the nucleus decreases towards the
range where the Screened Coulomb pair potential takes effect. This is done by multiplying
f(r) with a Fermi-Dirac-like function,
g(r)
=
1
e
r
z
−
r
∆
r
z
+
1
,
(22)
in which r
z 
is a constant that determines where the function is equal to 1/2 and 
∆
r
z
is a
constant that determines how steeply the function rises from 0 to 1. For r
z 
and 
∆
r
z 
the
values of 1.8 Å and 0.025 Å were chosen. In this way the electron density and hence the
embedding energy are practically zero near the atom core. The stiffened Johnson-Oh and
Screened Coulomb pair potentials have been smoothly connected in the r -range where their
values are almost the same (1.59 to 2.08 Å). At the transition radius from EAM interaction
to pure pair interaction the electron density interaction is not as repulsive as the original
EAM potential as a result of this modification of the electron density. This causes a small
shoulder in the potential. This error has to be accepted in going from the EAM potential to
the pair potential. The Mo electron density distribution and Mo-Mo pair potential for low
electron densities and energies have been plotted in figs. 2 and 3. For an isolated system of
two molybdenum atoms the potential well shape for an isolated pair of atoms is depicted in
figure 4. The potential well only has this shape in the absence of other atoms. In the
presence of, for example, a third atom at a constant distance from the other two, the well
will not only be shifted up or down, but it will also change shape due to the extra electron
density contribution of the third atom.
All results in this work have been calculated using the modified Johnson-Oh EAM
potential in combination with the Firsov-Molière Screened Coulomb potential, unless
mentioned otherwise. A few simulations have been calculated using a different, Finnis-
Sinclair EAM potential without modifications. This potential is described in [7].
For noble gas interactions the Firsov-Molière Screened Coulomb pair potential has
been used with the values in table II. No attractive potential was added to this purely
repulsive potential. This is justified because the noble gas atoms are practically always in
high-energy situations where neglecting a few meV of binding energy has no influence.

---

15
Figure 2. Electron density f versus distance r for Mo. The solid
curve represents the electron density used in simulations. The dashed
curve is the unmodified density (eqn. (6)).
Figure 3. Mo-Mo pairpotential versus distance r for low energies.

---

16
Figure 4. Potential energy U  of an isolated pair of Mo atoms at a distance
r . Note the shoulder in the transition region around r = 1.8 Å.