** Conclusions**

**A. Appendix**

**A.3. The Frequency Shift of the Verlet Integrator**

This operator will be termedVerlet propagatorin the following.

The differential equation of a one-dimensional harmonic oscillator reads

− * ^{∂V x}*(t)

^{}

*∂x*(t) = −F x(t)^{} =m*∂*^{2}x(t)

*∂t*^{2} (A.3.7)

with a harmonic potentialV defined by
V(x) = ^{k}

2x^{2}. (A.3.8)

The well-known solution of this differential equation reads
x(t) =C_{1}cos

rk mt

+C_{2}sin
rk

mt

, (A.3.9)

where usually the substitution*ω* =
qk

m is performed, yielding
x(t) =C_{1}cos

*ωt*

+C_{2}sin

*ωt*

. (A.3.10)

The corresponding velocityv(t)then is obtained as
v(t) = ^{∂}

*∂t*x(t) =−C_{1}*ω*sin
*ωt*

+C_{2}*ω*cos
*ωt*

. (A.3.11)

By lettingC_{1} = 1,C_{2} = 0, it becomes clear that the trajectories of this system are
ellipses in thex−vphase space with the two radii 1 and*ω. Only in the caseω* =1
(i.e.,k=m), the trajectories are circles.

For the further derivation, it is desirable for all such trajectories to be circles.

Therefore atransformedvelocityv^{0}(t)is defined by
v^{0}(t) = ^{v}(t)

*ω* . (A.3.12)

It follows directly from this definition that in thex−v^{0} space all trajectories which
obey Equation A.3.7 are circles. In the following, only the transformed velocity will
be considered.

Applying the equations of the Verlet integrator from Equations A.3.1 and A.3.3 to the harmonic oscillator by substituting

a x(t)^{}= −^{1}
m

*∂V x*(t)^{}

*∂x*(t) =−^{k}

mx(t) =−*ω*^{2}x(t) (A.3.13)

A.3. The Frequency Shift of the Verlet Integrator

gives

x(t+_{∆}t) =x(t) +_{∆}tv(t)−^{∆}^{t}^{2}^{ω}^{2}

2 x(t), (A.3.14)

v(t+_{∆t}) =v(t)− ^{∆tω}

2

2

2x(t) +_{∆tv}(t)− ^{∆t}

2*ω*^{2}
2 x(t)

. (A.3.15)

with the velocity finally being substituted by the transformed velocity v^{0} in the
second equation:

v^{0}(t+_{∆t}) =v^{0}(t)− ^{∆tω}
2

2x(t) +_{∆tωv}^{0}(t)−^{∆t}

2*ω*^{2}
2 x(t)

. (A.3.16)

Combining Equations A.3.14 and A.3.16 into the form of a complex Verlet propa-gator such as in Equation A.3.6 yields

T_{∆}_{t}(z) =Re(z) +_{∆}tωIm(z)− ^{∆}^{t}^{2}^{ω}^{2}
2 Re(z)
+i

"

Im(z)−^{∆}^{tω}
2

2Re(z) +_{∆}tωIm(z)−^{∆}^{t}^{2}^{ω}^{2}
2 Re(z)

# (A.3.17)

with

Re z(t)^{} =_{x}(t), (A.3.18)
Im z(t)^{} =v^{0}(t). (A.3.19)
Please note that the parameters ∆tand*ω* always appear with equal exponents in
the above equation. The propagator is therefore invariant under any parameter
change which leaves the product ∆t*ω* invariant. This is in line with intuition, as
doubling the frequency*ω* and halving the time step∆t is just a re-parametrization
of time and should not change anything else.

As already stated above, exact solution trajectories of any harmonic oscillator are
circles in thex−v^{0} space, i.e. the absolute value|z|remains constant along these
trajectories.

**Claim 1:**The Verlet propagator from Equation A.3.17 keeps|z|approximately
con-stant, and the residuum vanishes quickly as∆tbecomes small.

Note: As

|z|^{2} =x^{2}+v^{0}^{2}= x^{2}+ ^{v}

2

*ω*^{2} = ^{2}

kE_{pot}+ ^{2}

mE_{kin}, (A.3.20)
Claim 1 is equivalent to the statement that the Verlet propagator approximately
conserves the total energy of the system.

**Proof:** Switching to the trigonometric representation of z by substituting z =
r cos(*ϕ*) +isin(*ϕ*)^{} on the right hand side of the following equation gives

T_{∆t} z

=rcos(*ϕ*) +_{∆tωr}sin(*ϕ*)−^{∆t}

2*ω*^{2}

2 rcos(*ϕ*)
+i

"

rsin(*ϕ*)− ^{∆tω}
2

2rcos(*ϕ*) +_{∆tωr}sin(*ϕ*)− ^{∆t}

2*ω*^{2}

2 rcos(*ϕ*)
#

(A.3.21) and

T_{∆}_{t}(z)^{}^{}

2= rcos(* _{ϕ}*) +

_{∆tωr}sin(

*)−*

_{ϕ}^{∆t}

2*ω*^{2}
2 rcos(* _{ϕ}*)

!2

+ rsin(*ϕ*)−^{∆tω}
2

2rcos(*ϕ*) +∆tωrsin(*ϕ*)−^{∆t}

2*ω*^{2}
2 rcos(*ϕ*)

!2

=r^{2} 1+^{1}

4∆t^{3}*ω*^{3}sin(2ϕ)−^{1}

4∆t^{4}*ω*^{4}cos(2ϕ)−^{1}

8∆t^{5}*ω*^{5}sin(2ϕ) + ^{1}
32∆t^{6}*ω*^{6}

1+cos(2ϕ)^{}

! ,

(A.3.22) thus

T_{∆}_{t}(z)^{}^{}

=|z| r

1+^{1}

4∆t^{3}*ω*^{3}sin(2ϕ)−^{1}

4∆t^{4}*ω*^{4}cos(2ϕ)−^{1}

8∆t^{5}*ω*^{5}sin(2ϕ) + ^{1}
32∆t^{6}*ω*^{6}

1+cos(2ϕ)^{}.

(A.3.23)
It can be seen that the Verlet propagator keeps the absolute value of the argument
approximately constant, and the residue is of the order ofO _{∆t}^{3}*ω*^{3}

(to see this, please note that√

1+x^{3} ≈1+ ^{x}_{2}^{3} +. . .). A time step∆t *ω* will lead to quickly
vanishing deviations. Apart from that, any real simulation will uniformly sample
the angle *ϕ, such that all terms which contain sin*(2*ϕ*)or cos(2*ϕ*)will cancel out
on average. Depending on the angle, some propagator steps will slightly enlarge

|z|, while others will slightly reduce it, causing a fluctuation of|z|around its exact
value. The only residual term in Equation A.3.23 which does not fluctuate around
zero is of the order of O _{∆t}^{6}*ω*^{6}

. This term is always positive, meaning that |z| (and therefore the total energy of the system) will always increase over the long term. However, a sufficiently small time step∆twill make this effect almost vanish.

As already discussed above, the exact solution trajectory of any harmonic
os-cillator constitutes a circle in the x−v^{0} space. Therefore it can be assumed that
the discrete Verlet propagator rotates the argument by a specific angle within the
complex plane in order to approximately resemble this exact trajectory.

A.3. The Frequency Shift of the Verlet Integrator

**Claim 2:** The angle by which the Verlet propagator rotates its argument around
the origin within the x−v^{0} space is approximately independent on the choice of
the argument itself.

This would mean that successive time steps of the Verlet propagator always
rotate by approximately the same angle in the x−v^{0} space, and this angle only
depends on∆tand*ω.*

**Proof:** The angle between two complex numbers can be determined by utilizing
the dot product:

^(z_{1},z_{2}) =arccos Re(z_{1})Re(z_{2}) +Im(z_{1})Im(z_{2})

|z_{1}| · |z_{2}|

!

. (A.3.24)

To obtain the angle by which the Verlet propagator rotates its argument, write

^^{}^{z,}^{T}∆t(z)^{}=arccos Re(z)Re T_{∆t}(z)^{}+Im(z)Im T_{∆t}(z)^{}

|z| ·^{}T_{∆t}(z)^{}

!

. (A.3.25)

Following from Claim 1, we can safely approximate

T_{∆}_{t}(z)^{}^{} ≈ |z|, and therefore

|z| ·^{}T_{∆t}(z)^{} ≈ |z|^{2}. Substitutingz =r cos(*ϕ*) +isin(*ϕ*)^{} yields

^^{}^{z,}^{T}∆t(z)^{}=arccos rcos(*ϕ*)Re T_{∆t}(z)^{}+rsin(*ϕ*)Im T_{∆t}(z)^{}
r^{2}

!

. (A.3.26)

Inserting the terms from Equation A.3.21 for Re T_{∆}_{t}(z)^{} and Im T_{∆}_{t}(z)^{} and
di-viding byr^{2}leads to

^^{}^{z,}^{T}∆t(z)^{}=_{arccos} _{cos}^{2}(*ϕ*)− ^{∆}^{t}

2*ω*^{2}

2 cos^{2}(*ϕ*) +sin^{2}(*ϕ*)

−^{∆t}

2*ω*^{2}

2 sin^{2}(*ϕ*) + ^{∆t}

3*ω*^{3}

4 sin(*ϕ*)cos(*ϕ*)

! (A.3.27)

=arccos 1− ^{∆t}

2*ω*^{2}
2 + ^{∆t}

3*ω*^{3}

4 sin(*ϕ*)cos(*ϕ*)

!

. (A.3.28)

Keeping in mind that arccos 1−x^{2}

≈ x+O x^{3}

, it is visible that the rotation
angle is approximately linearly dependent of the time step∆tfor fixed*ω, which is*
completely in line with the expectations, as the Verlet propagator need to traverse
larger pieces of the trajectory circle with larger time steps. The term which depends
on *ϕ*is of higher order with respect to∆tω, and therefore vanishes quickly. Apart

from that, like already discussed in the proof of claim 1, the angles*ϕ*will be
uni-formly sampled in a real simulation, canceling out the residual terms on average
over long runs, asR2π

0 sin(*ϕ*)cos(*ϕ*)dϕ=0. Some of the time steps will possess a
rotation angle above average, others below.

The final approximation for the angle (which neglects the residual terms) there-fore reads

^^{}^{z,}^{T}∆t(z)^{}≈arccos
1−^{∆}^{t}

2*ω*^{2}
2

. (A.3.29)

**Derivation of the Frequency Shift**

Based on the angle*α* which one time step covers in the x−v^{0} space derived in
Equation A.3.29, the cycle duration of the discrete system can be determined by

*τ*_{cycle}= ^{2π}

*α* ∆t (A.3.30)

= ^{2π}^{∆}^{t}

arccos 1−^{∆}^{t}^{2}_{2}^{ω}^{2}^{}^{.} ^{(A.3.31)}
Applying*ω* = _{τ}^{2π}

cycle, one then obtains

*ω*_{verlet}= ^{arccos 1}− ^{∆}^{t}^{2}^{ω}_{2}^{2}^{exact}^{}

∆t . (A.3.32)

In practical applications, it might be more useful to compute the exact frequency based on the approximate one. Fortunately, the inverse function of the above is easily given by

*ω*_{exact}=

p2−2 cos(_{∆t}*ω*_{verlet})

∆t . (A.3.33)

A.3. The Frequency Shift of the Verlet Integrator

**Results**

In the following, some results with parameters from typical simulations are shown.

For a typically used time step of 0.5 fs, the deviation of the frequency is found to
be around 10 cm^{−}^{1}for a C–H vibration at 3000 cm^{−}^{1}. At∆t =1.0 fs, the deviation
even amounts to 40 cm^{−}^{1}.

**Table A.1.:**Frequency shift of Verlet integrator for typical vibrations and time steps∆t.

*ω*_{verlet}

*ω*_{exact}/ cm^{−}^{1} ∆t=0.1 fs ∆t =0.5 fs ∆t=1.0 fs

10 10.000 000 014 8 10.000 000 369 10.000 001 48

100 100.000 014 8 100.000 369 100.001 48

1000 1000.0148 1000.370 1001.483

2000 2000.118 2002.967 2012.013

3000 3000.399 3010.064 3041.396

4000 4000.946 4024.025 4101.159

10 000 10 014.834 10 411.945 13 033.591