Molecular dynamics thesis

r(t
+ ∆
t)
=
r(t)
+
v(t)
∆
t
+
1
2
a(t)(
∆
t)
2
,
(23)
where a is determined from the positions of the atoms using eqn. (2), and from these new
positions the new accelerations a(t + 
∆
t) are calculated. The velocities at t + 
∆
t are
calculated in two steps:
v(t
+
1
2
∆
t)
=
v(t)
+
1
2
∆
ta(t),
(24)
v(t
+ ∆
t)
=
v(t
+
1
2
∆
t)
+
1
2
∆
ta(t
+ ∆
t).
(25)
When positions, velocities, and accelerations at t + 
∆
t are known, the process can be
repeated to calculate the next configuration.
The timestep 
∆
t  is not a constant in a simulation run. Timestep control is discussed
in the next section.
2.3.3 Timestep control
The timestep 
∆
t should be chosen in such a way that no atom moves more than a
small distance 
∆
r  over which the forces on the atoms do not change significantly. In this
way, the numerical result after 
∆
t will be a close approximation of the analytical result. An
accurate criterion for the magnitude of 
∆
r  would involve the computation of countless
forces with atoms at various trial positions. As this would lead to excessive computer time,
we have followed the usual approach of using a fixed distance limit 
∆
r  
≤
R. R is chosen in
such a way that even in the situation in which dF /dr  has its highest value, the forces on
the atoms do not change significantly over a distance R. In the simulations used to calculate
the results in this thesis the value R  = 0.020 Å was chosen.
During IBAD simulations, velocities of atoms vary over a very broad range.
Therefore a constant timestep is impractical, because the timestep would always have to be
so small that the fastest atom of the entire simulation does not move more than 
∆
r. Because
of this, a variable timestep is chosen, allowing the fastest particle of one configuration to
move a distance R. That way, a much larger timestep can be used for most configurations,
because only few configurations contain very fast atoms. Another way to improve
efficiency is to exclude atoms from limiting the timestep as long as they are in free flight.
The algorithm to determine the timestep in this way is somewhat complicated. The
general idea for timestep control with variable timestep and free flight algorithm is as
follows:

---

19
- for the first atom in a loop over all atoms experiencing a force (a 
≠
0), the distance r  that
it will move away from its starting point during an initial trial timestep
*
∆
t is calculated 
using eqn. (26),
**
∆
r
( )
2
=
v
2
∆
t
( )
2
+
(a
•
v)
∆
t
( )
3
+
1
4
a
2
∆
t
( )
4
.
(26)
The acceleration is included in eqn. (26) because for atoms with a small velocity, for 
instance close to ultimate collision impact, it can contribute significantly to 
∆
r  and 
neglecting it could result in a very inaccurate approximation for systems with small 
numbers of atoms (neglecting (a•v)(
∆
t)
3 
+ 1/4*a
2
*(
∆
t)
4
will result in a timestep that is too 
long. In systems of thousands of atoms this is not important because there will always be 
other atoms that limit the time step. However, theoretically, in systems with few atoms, 
such as two-particle systems used in accuracy tests, situations can arise in which any time 
step could be selected).
Depending on whether this distance is smaller or larger than 
R
, ∆
t is increased or 
decreased and 
∆
r is calculated again for the atom under consideration. This process is 
repeated until 
∆
r corresponds to R.
***
- for the other atoms experiencing a force, the travel distance r in time 
∆
t  of the previous 
atom is calculated and if it exceeds R, 
∆
t is further reduced. This establishes a first 
estimate of t .
- free atoms (a =0) are allowed to travel an unlimited distance, but only for as long as 
they do not come near other atoms during the trial timestep. If they do, the trial time is 
further reduced to allow them to come just inside the influence of another atom.
- atoms for which the cosine of the angle between a  and v   is smaller than 
−
8 9
are 
treated in a special way. The criterion cos(v ,a ) < 
−
8 9
selects those atoms whose 
velocity and acceleration point almost exactly in opposite directions. These atoms are about
to ‘make a turn’. This requires a special approach because the distance they travel away 
from their starting points is not a monotonously increasing function of 
∆
t.. This can be 
seen as follows. From eqn. (26) we get
d
∆
r
( )
2
( )
d(
∆
t)
=
2v
2
∆
t
+
3(a
•
v)
∆
t
( )
2
+
a
2
∆
t
( )
3
.
(27)
( r )
2
  is not a monotonous function of 
∆
t if there are real, non-zero roots for the 
equation
d
∆
r
( )
2
( )
d(
∆
t)
=
0
⇒
a
2
∆
t
( )
2
+
3(a
•
v)
∆
t
+
2v
2
=
0 .
(28)
In general, the roots 
∆
t
1
and 
∆
t
2
of eqn. (28) are
∆
t
1
,
∆
t
2
= −
3(a
•