6 Introduction

V

T 2

V (r) = D 1 ’ e’a(r’b)

D

0 Er

0 b

Figure 1.1 Morse curve with a = 2/b (solid curve). Dotted curve: parabola with

same curvature as Morse curve at r = b: V = Da2 (r ’ b)2 .

implied. This book is meant to provide the necessary scienti¬c background

and to promote awareness for the limitations and inaccuracies of simulating

the “real world”.

Example: An oscillating bond

In this example we use brute-force simulation to attack a problem that could be

approached analytically, albeit with great di¬culty. Consider the classical bond

length oscillation of a simple diatomic molecule, using the molecule hydrogen ¬‚uo-

ride (HF) as an example. In the simplest approximation the potential function is

a parabola:

V (r) = 1 k(r ’ b)2 , (1.1)

2

with r the H“F distance, k the force constant and b the equilibrium distance. A

better description of the potential function is the Morse function (see Fig. 1.1)

2

V (r) = D 1 ’ e’a(r’b) , (1.2)

where D is the dissociation energy and a is a constant related to the steepness of

the potential. The Morse curve is approximated near the minimum at r = b by a

parabola with force constant k = 2Da2 .

The Morse curve (Morse, 1929) is only a convenient analytical expression that has

some essential features of a diatomic potential, including a fairly good agreement

with vibration spectra of diatomic molecules, but there is no theoretical justi¬cation

for this particular form. In many occasions we may not even have an analytical form

for the potential, but know the potential at a number of discrete points, e.g., from

quantum-chemical calculations. In that case the best way to proceed is to construct

the potential function from cubic spline interpolation of the computed points. Be-

1.1 What is this book about? 7

Table 1.1 Data for hydrogen ¬‚uoride

mass H mH 1.0079 u

mass F mF 18.9984 u

dissocation constant D 569.87 kJ/mol

equilibrium bond length b 0.09169 nm

kJ mol’1 nm’2

5.82 — 105

force constant k

cause cubic splines (see Chapter 19) have continuous second derivatives, the forces

will behave smoothly as they will have continuous ¬rst derivatives everywhere.

A little elementary mechanics shows that we can split o¬ the translational motion

of the molecule as a whole, and that “ in the absence of rotational motion “ the

bond will vibrate according to the equation of motion:

dV

μ¨ = ’

r , (1.3)

dr

where μ = mH mF /(mH + mF ) is the reduced mass of the two particles. When we

start at time t = 0 with a displacement ”r and zero velocity, the solution for the

harmonic oscillator is

r(t) = b + ”r cos ωt, (1.4)

with ω = k/μ. So the analytical solution is simple, and we do not need any nu-

merical simulation to derive the frequency of the oscillator. For the Morse oscillator

the solution is not as straightforward, although we can predict that it should look

much like the harmonic oscillator with k = 2Da2 for small-amplitude vibrations.

But we may expect anharmonic behavior for larger amplitudes. Now numerical sim-

ulation is the easiest way to derive the dynamics of the oscillator. For a spline-¬tted

potential we must resort to numerical solutions. The extension to more complex

problems, like the vibrations of a molecule consisting of several interconnected har-

monic oscillators, is quite straightforward in a simulation program, while analytical

solutions require sophisticated mathematical techniques.

The reader is invited to write a simple molecular dynamics program that uses

the following very general routine mdstep to perform one dynamics step with the

velocity-Verlet algorithm (see Chapter 6, (6.83) on page 191). De¬ne a function

force(r) that provides an array of forces F , as well as the total potential energy

V , given the coordinates r, both for the harmonic and the Morse potential. You

may start with a one-dimensional version. Try out a few initial conditions and

time steps and look for energy conservation and stability in long runs. As a rule

of thumb: start with a time step such that the fastest oscillation period contains

50 steps (¬rst compute what the oscillation period will be). You may generate

curves like those in Fig. 1.2. See what happens if you give the molecule a rotational

velocity! In this case you of course need a two- or three-dimensional version. Keep

to “molecular units”: mass: u, length: nm, time: ps, energy: kJ/mol. Use the data

for hydrogen ¬‚uoride from Table 1.1.

The following function performs one velocity-Verlet time step of MD on a system

of n particles, in m (one or more) dimensions. Given initial positions r, velocities v

and forces F (at position r), each as arrays of shape (n, m), it returns r, v, F and

8 Introduction

0.13

H-F distance (nm)

0.12

0.11

0.1

0.09

0.08

5 10 15 20 25 30 35 40

time (fs)

Figure 1.2 Oscillation of the HF bond length, simulated with the harmonic oscil-

lator (solid curve) and the Morse curve (long dash), both with initial deviation

from the equilibrium bond length of 0.01 nm, Dotted curve: Morse oscillator with

initial deviation of 0.03 nm, showing increased anharmonic behavior. Note that the

frequency of the Morse oscillator is lower than that of the harmonic oscillator. A

time step of 0.2 fs was used; the harmonic oscillator simulation is indistinguishable

from the analytical solution.

the potential energy V one time step later. For convenience in programming, the

inverse mass should be given as an array of the same shape (n, m) with repeats of

the same mass for all m dimensions. In Python this n — m array invmass is easily

generated from a one-dimensional array mass of arbitrary length n:

invmass=reshape(repeat(1./mass,m),(alen(mass),m)),

or equivalently

invmass=reshape((1./mass).repeat(m),(alen(mass),m))

An external function force(r) must be provided that returns [F, V ], given r. V is

not actually used in the time step; it may contain any property for further analysis,

even as a list.

python program 1.1 mdstep(invmass,r,v,F,force,delt)

General velocity-Verlet Molecular Dynamics time step

01 def mdstep(invmass,r,v,F,force,delt):

02 # invmass: inverse masses [array (n,m)] repeated over spatial dim. m

03 # r,v,F: initial coordinates, velocities, forces [array (n,m)]

04 # force(r): external routine returning [F,V]

05 # delt: timestep

1.2 A modeling hierarchy 9