As a mathematician, most people think I spend my time with writing mystical symbols on blackboards and solving problems nobody is interested in. Of course this is not true, I just try to find symplectic numerical integrators for Hamiltonian systems… But for some reason, telling people this doesn’t help to demystify math?!

In this post, I want to show you the beautiful world of geometry and some of the things we mathematicians do!

For mathematicians, physicists and other math-experts, this post might offer an illustrative introduction

to an interesting field, called geometric numerical integration, a mix between ideas from theoretical physics and applied math.

But primarily I want to address non-experts! Mathematics – especially geometry – is all about structures

which should help us humble humans to capture the very foundations of reality. I want to show how one particular mathematical

structure helps us to see more than numbers and how to use it!

*Notice: All content of this post is other people research and not my own work! I just love the topic and want to introduce it to other people as well.*

#### Mathematical structures matter…

Let me give you an example for mathematical structures from the real world.

The performance of students is often measured by grades, which are represented in Germany by numbers ranging from 1 to 6.

Now, numbers on its own have different mathematical structures assigned to it.

For example we can consider numbers together with addition and multiplication (N,+,·), which form a so-called algebraic *ring.*

In this mathematical structure we can calculate for example 2 + 3 = 5.

But considering grades together with addition and multiplication, i.e. (G,+,·) is nothing but bullshit. A calculation like 2 + 3 = 5 does have no meaning for grades.

(Maybe a weighted sum, like 2 “+” 3 = 2.5 would be reasonable for the space of grades.)

The message: Not everything which is represented by a number belongs to the same underlying mathematical space!

Even if I’ll talk about molecules and simulations, the key message of this post is to show,

how identifying mathematical structures can help to build better predictions of the reality.

Ignoring these structures will literally lead to (mathematical) explosions in our model.

### Molecules, the building blocks of life

My Bachelor’s thesis was about molecules and without any intention, I ended up doing geometry.

But more about that later.

Molecules are assembled out of many atoms (Hydrogen, Carbonate, Oxygen…). These atoms are bound

together since they love sharing some of their electrons.

How do they move? To answer this, we need to find the active forces.

#### Forces and their dynamics

If atoms are close to each other, complicated forces determine their movement.

For example a bond between two atoms pulls them together if they are too far apart

and pushes them apart if they are too close.

But there are many different forces, the image on the right hand side visualizes some of them.

If we know these forces, we arrive at the starting point for mathematical models

and can start to build simulations. For this we use the fundamental law

force = mass · acceleration

If the force acting on an atom are constant, we could easily compute its movement and positions at all times.

But when atoms change their positions, the acting forces change as well…

This has dramatic consequences: In each single moment, as short as it might be, the forces acting between the atoms might change at least a little bit. Which complicates the task of computing the movements of the atoms.

### (A bad) simulation of molecules (which does not respect the geometry)

The easiest way to build a simulation is always the assumption of small changes.

If we want to simulate the position of our atoms times

0, 2ps, 4ps, 6ps, … (ps = picoseconds = 0.000000001 seconds)

it is tempting to assume that forces between these time steps are constant.

This leads to the following simulation approach:

If our atom is initially at position 2.3nm (nanometer) and has velocity 1.1 nm/ps, then we can use the force model to compute the acting force at the current time.

Let the force acting on the atom be of strength 0.14 nN (piconewton).

The prediction for the next position and velocity of the atom is given by the two-step algorithm:

- The new position 2ps after the initial state is approximated by
2.3nm + 2ps · 1.1nm/ps = 4.5nm.

- The new velocity after 2ps is set to
1.1nm/ps + 2ps · 0.14nm/ps² = 1.38nm/ps.

We repeat these two steps each two picoseconds and with the new values as initial data.

This method is called *explicit Euler method, *named after the great physicist and mathematician Leonard Euler.

Back in the days, many people were employed just to calculate the two steps above over and over again.

In this way, the movement of planets, stars or ballistic objects were approximated. Which was important for

navigation on sea or war.

But it is important to notice, that we just assume the force to stay constant during these 2 picoseconds.

…at least, we update it every two picoseconds. Let’s find out if this is good enough!

Since I suck at doing calculations by myself, I implemented this method as a plugin for a software called Avogadro, hence the computer does all the work for me. Below you can see a simulation which works as described above!

If you have never seen an explosion in mathematics, that is what you see in the video above!

It is called explosion, since the small atoms will become arbitrarily fast. Guess what, that’s bad and has nothing to do with the real behavior of atoms.

#### Improving the algorithm without geometry

Of course, there are many non-geometric ways to get better results.

We could do smaller time steps. But if we compute an update after each picosecond,

we need twice as many calculations as before. In fact, even modern computers

are too slow, to make realistic predictions of molecules with the explicit Euler method from above,

since the time steps necessary for good accuracy are far too small.

(There are some other improvements as well, but we skip them,

since they all share a common disability.)

## The geometry of molecules

With tools from mathematics and physics, we can improve our simulation

such that the atoms stay together, as they do.

### An illustrative example

*We restrict our-self to the movement of one atom (the other atom is assumed to be fixed)!*

(That’s also typical for mathematicians, we ignore everything which is complicated just to find problems which are

simple enough for us to solve.)

The dynamics of this single atoms are boring. Our atom will just oscillate back and forth. Like shown in the image below.

Until now, we used two numbers, the position and the velocity, to describe the state of an atom.

For most scientist there is nothing more to find here, but as mathematicians we want to find the mathematical space, where these numbers belong to.

### The phase space

In our last sketch, the arrows indicated how fast our little atom is. But there should be two arrows in the central position, since the atom traversed this point in both directions.

Let us get rid of this ambiguity and move an arrow a bit upwards, if it is pointing to the right and we draw the arrow more on the bottom, if the velocity points to the right side.

The resulting picture shows the dynamics of the atom in the* (position-velocity) phase space.*

Which properties does this (blue) space have? We will just call it a *phase space *from now on.

(Usually, it is better to use a (position-momentum) phase space, but I will commit the crime of simplifying a lot of things!)

### Physical flow in phase space

One of the unifying ideas in math is to investigate spaces by means of their homomorphisms, i.e.

we look how maps transforms the space and it’s properties. A very interesting set of maps is given by physical flows in phase space.

A point in our phase space is a pair of a position and a velocity, i.e. two numbers.

For each initial position and velocity, we can associate a position and velocity after one time step (i.e. after 2 picoseconds).

We call this map the *flow* in phase space.

The get a feeling for the exact physical flow, we take a cloud of 18 points in the phase space (i.e 18 initial positions and velocities, show as purple dots in the picture below).

We look where these points end up after some times steps. In this way we can visualize how the physical flow transforms the phase space itself.

One interesting observation is, that the volume of the point cloud does not change after the application of the maps.

This is no incidence, but a true fact which holds for all physical systems! The statement is known as *Liouville’s Theorem:***
**Physical flows preserve volumes in phase space!

It seems that volume is kind of important?!

### The correct mathematical structure: (Phase space, Volume form)

Knowing about this law, allows us to assign a better mathematical structure to the phase space.

We call it a *volume form*, which is just a tool to compute the volume of (infinitely small) point clouds in the phase space.

Now each time we work with our phase space, we should also keep track how the *volume form* changes!

If the phase space is two-dimensional, the pair (phase space, volume form) is the same as a symplectic manifold.

Now, **Liouville’s Theorem (for 2D phase spaces)** can be restated as: *Physical flows preserve the symplectic structure of the phase space!*

For the curious reader: Symplectic manifolds are indeed a very proper structure for classical physics.

In higher dimensional spaces (i.e. for example one atom in 3D space, i.e. 6D phase space) the symplectic strucutre is a bit complicated.

We would have to consider the 2D phase spaces for each single space dimension individually (up-down, left-right and forward-backward).

and the symplectic form is given by the sum of volumes of these individual 2D phase spaces.

Hence understanding the case for an atom in 1D space is a good starting point for understanding the general theory.

Please notice that I used several brutal simplifications here! Exact definitions can be found in the reference or in my presentation about this topic.

### The flow of the (geometry ignoring) numerical simulation in phase space

Previously we just saw that the (bad) numerical simulation does lead to faster and faster atoms.

With the language developed above, we can understand more precisely, how the geometry of the phase space was violated!

We take the point cloud with 18 point, just as above, and compute the numerical flow, i.e. the new positions and velocities using the update rule.

One time step does look reasonable, but only since our eyes are not good enough to see the small violation of the geometry already.

Therefore we look where the point cloud will be after 100 time steps (i.e. 200 picoseconds).

The result is dramatic! Obviously the volume of our point cloud is not conserved!

We can also look at the background color gradient of the phase space, to see that the complete phase space is massively sketched.

In reality, the gradient should not change at all!

### A geometry preserving numerical simulation

It is possible to find an update rule, for which the flow preserves volume in phase space! (Even more, it preserves the symplectic structure).

- The new velocity after 2ps is computed by assuming constant velocity and computing the new velocity after 2ps
1.1nm/ps + 2ps · 0.14nm/ps² = 1.38nm/ps.

- The new position is approximated by adding the
*just computed new velocity*, not the initial velocity!`2.3nm + 2ps ·`

**1.38nm/ps**= 5.06nm.

We repeat these two steps each two picoseconds and with the new values as initial data. This method is called the *symplectic Euler method*.

Is this method now better? We preform the experiment for a pair of bonded atoms form above again!

Does it again lead to a (mathematical) explosion?

Now explosion this time! As we can see, the atoms start not to oscillate like crazy, which is the correct behavior.

### It works! But why?

The idea is as follows: For our approximations we always assume that the acceleration and velocity is constant for 2 picoseconds,

but we can choose between the initial velocity and acceleration (explicit Euler method) and the terminal velocity and acceleration (implicit Euler method).

Using the initial values will lead to an increase of volume in phase space.

Using the terminal velocity and acceleration leads to a decrease of volume in phase space.

Hence we just us a mix of both, which is of course also an approximation, but it preserves the volume in phase space as well, since it is a mix between

an overestimation and an underestimation!

Hence, we now have a simulation which does not only follow approximately the physical law of motion, but it also

respects the geometric structure of the phase space. Which leads to notably better results.

If we look at the point cloud under this numerical flow, we find that the volume will be strangely deformed by the flow, but it preserves

the volume of the point cloud itself.

Of course, this is not the end of the story, more just the beginning! The method has far more sweet properties

which can be found in the reference cited at the end of this post or in this presentation of mine about symplectic numerical integration.

## Applications in molecular dynamical simulations

Finally, we can use this technique to simulate more complex molecules.

In fact most scientists use variations of symplectic methods to simulate molecules. Some know about the underlying geometric reasons,

others just trust their predecessors.

As a final example we want to simulate a full atom, which is surrounded by small water molecules.

These water molecules are not visible here, but they hit the big molecules very often, which causes

the fluctuation of the atoms in the movie below.

If we look carefully at the long chain of atoms on the left, we see that it moves closer to the center.

(Which is interesting as a topic on its own and connected to the famous law that entropy always increases!)

As a result of the simulation we are often just interested in such characteristic values, like the length of a molecules.

Since these simple values might help us to understand how many molecules of one kind interact and arrange each other.

It is only reasonable, that if we want to measure something like a length of a molecule, we should respect the geometry of the physical space!

### Take home message

Not everything which can be represented by numbers belongs to the same mathematical space!

If the computations we do violated these structures, we might run into trouble.

Unfortunately, mathematical structures are invisible, but if we find them, we can use them to our advantage.

#### A personal note

If you are interested in more details, don’t hesitate to contact me, I love this topic. Or check out one of the references below or

my presentation. Since this is my first post in this blog, I also appreciate criticisms and suggestions how to improve further posts.

### References

Images and videos where create with the great molecular editor *Avogadro*

http://avogadro.cc/

For more on the scientific side of symplectic integration, we recommend the book:

*Geometric Numerical Integration, Structure-Preserving Algorithms for Ordinary Differential Equations,*

Ernst Hairer, Christian Lubich, Gerhard Wanner, Springer 2006