• Keine Ergebnisse gefunden

# The Frequency Shift of the Verlet Integrator

## Conclusions

### 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) =m2x(t)

∂t2 (A.3.7)

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

2x2. (A.3.8)

The well-known solution of this differential equation reads x(t) =C1cos

rk mt

+C2sin rk

mt

, (A.3.9)

where usually the substitutionω = qk

m is performed, yielding x(t) =C1cos

ωt

+C2sin

ωt

. (A.3.10)

The corresponding velocityv(t)then is obtained as v(t) =

∂tx(t) =−C1ωsin ωt

+C2ωcos ωt

. (A.3.11)

By lettingC1 = 1,C2 = 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 atransformedvelocityv0(t)is defined by v0(t) = v(t)

ω . (A.3.12)

It follows directly from this definition that in thex−v0 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) =−ω2x(t) (A.3.13)

A.3. The Frequency Shift of the Verlet Integrator

gives

x(t+t) =x(t) +tv(t)−t2ω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 v0 in the second equation:

v0(t+∆t) =v0(t)− ∆tω 2

2x(t) +∆tωv0(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

Tt(z) =Re(z) +tωIm(z)− t2ω2 2 Re(z) +i

"

Im(z)− 2

2Re(z) +tωIm(z)−t2ω2 2 Re(z)

# (A.3.17)

with

Re z(t) =x(t), (A.3.18) Im z(t) =v0(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−v0 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 =x2+v02= x2+ v

2

ω2 = 2

kEpot+ 2

mEkin, (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ωrsin(ϕ)−∆t

2ω2

2 rcos(ϕ) +i

"

rsin(ϕ)− ∆tω 2

2rcos(ϕ) +∆tωrsin(ϕ)− ∆t

2ω2

2 rcos(ϕ) #

(A.3.21) and

Tt(z)

2= rcos(ϕ) +∆tωrsin(ϕ)∆t

2ω2 2 rcos(ϕ)

!2

+ rsin(ϕ)∆tω 2

2rcos(ϕ) +∆tωrsin(ϕ)∆t

2ω2 2 rcos(ϕ)

!2

=r2 1+1

4∆t3ω3sin()1

4∆t4ω4cos()1

8∆t5ω5sin() + 1 32∆t6ω6

1+cos()

! ,

(A.3.22) thus

Tt(z)

=|z| r

1+1

4∆t3ω3sin()1

4∆t4ω4cos()1

8∆t5ω5sin() + 1 32∆t6ω6

1+cos().

(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 ∆t3ω3

(to see this, please note that√

1+x3 ≈1+ x23 +. . .). 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 ∆t6ω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−v0 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−v0 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−v0 space, and this angle only depends on∆tandω.

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

^(z1,z2) =arccos Re(z1)Re(z2) +Im(z1)Im(z2)

|z1| · |z2|

!

. (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

Tt(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) r2

!

. (A.3.26)

Inserting the terms from Equation A.3.21 for Re Tt(z) and Im Tt(z) and di-viding byr2leads to

^z,Tt(z)=arccos cos2(ϕ)− t

2ω2

2 cos2(ϕ) +sin2(ϕ)

∆t

2ω2

2 sin2(ϕ) + ∆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−x2

≈ x+O x3

, 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, asR

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,Tt(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−v0 space derived in Equation A.3.29, the cycle duration of the discrete system can be determined by

τcycle=

α ∆t (A.3.30)

= t

arccos 1−t22ω2. (A.3.31) Applyingω = τ

cycle, one then obtains

ωverlet= arccos 1t2ω22exact

∆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 cm1for a C–H vibration at 3000 cm1. At∆t =1.0 fs, the deviation even amounts to 40 cm1.

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

ωverlet

ωexact/ cm1 ∆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

Outline

ÄHNLICHE DOKUMENTE

Wähle deine Sprache

Die Website wird in die von dir ausgewählte Sprache übersetzt.

Vorgeschlagene Sprachen für dich:

Andere Sprachen: