Advanced Physical and Theoretical Chemistry: Molecular Dynamics

Whether you want to see molecules move on an atomistic level or want to investigate gas particles interacting with each other, or even follow stars or galaxies as they brush past each other, then the simulation of particles is the way to go. The mathematical and algorithmic foundations are surprisingly similar in all these cases.

In these two modules (one on fundamentals and one more in-depth on applications), we will write a particle simulation in two dimensions.

Whether it's planets or atoms you are considering, Newton's laws apply reasonably well, most importantly

\[\mathbf{F}(\mathbf{x})=m\mathbf{a}(\mathbf{x})\]

which relates the force acting on a particle \(\mathbf{F}(\mathbf{x})\) to the acceleration \(\mathbf{a}(\mathbf{x})\) that this particle of mass \(m\) feels. Since we usually can obtain the force acting on a particle quite easily, we now also know about the acceleration which allows us to obtain the trajectory (i.e. the position over time \(\mathbf{x}(t)\)) by numerical integration. This means solving an ordinary differential equation (ODE) for which there is a wealth of possible techniques. We will focus on the most simple ones: the Euler approach and the (velocity) verlet integrator.

The core idea of the Euler Integrator is to take the current (\(t=t_0\)) position of a particle \(\mathbf{x}(t_0)\) move it in the direction of the position derivative, i.e. the velocity \(\mathbf{v}(t_0)\). If this is done in tiny steps, chances are that the velocity direction has not changed (much) and the new position still lies on (or very close to) the trajectory of that particle. Technically, this can be seen as updating velocity and position in very small time steps \(\Delta t\):

\[\mathbf{x}(t+\Delta t) = \mathbf{x}(t) + \mathbf{v}(t)\Delta t\]

\[\mathbf{v}(t+\Delta t) = \mathbf{v}(t) + \mathbf{a}(t)\Delta t\]

Repeating this over and over with very small \(\Delta t\) yields the full trajectory in discrete positions. That is a bit as if you watch somebody walk who is illuminated by a stroboscope: you cannot see the full motion, but you know what the motion looked like at equidistant points in time. This integrator is simple, easy to implement and highly inaccurate. This can easily be seen if you consider a particle moving horizontally while a force is pushing it down, creating a curved path. At the end of the time step \(\Delta t\), the velocity vector is less steep than it should be, i.e. the resulting position \(\mathbf{x}(t+\Delta t)\) is always higher than the true position where it should have been, which yields a systematic error at every step. Decreasing a time step \(\Delta t\) makes any numerical integrator more accurate but at the same time also increases the cost as now more steps (\(T/\Delta t\) many) are necessary to describe the same total time \(T\).

A better alternative is the Verlet integrator. The idea is the same, but a slighly different formalism improves its accuracy over the Euler integrator. It is derived from considering two Taylor expansions of the trajectory forward and backwards in time:

\[\mathbf{x}(t+\Delta t) = \mathbf{x}(t) + \mathbf{v}(t)\Delta t + \frac{1}{2}\mathbf{a}(t)\Delta t^2 + \frac{1}{6}\mathbf{j}(t)\Delta t^3 + \mathcal{O}(\Delta t^4)\]

\[\mathbf{x}(t-\Delta t) = \mathbf{x}(t) - \mathbf{v}(t)\Delta t + \frac{1}{2}\mathbf{a}(t)\Delta t^2 - \frac{1}{6}\mathbf{j}(t)\Delta t^3 + \mathcal{O}(\Delta t^4)\]

Here \(\mathbf{j}\) is the jerk (or jolt), i.e. the derivative of the acceleration w.r.t. time. \(\mathcal{O}(\Delta t^4)\) groups all higher order terms in the expansion. Adding the two expansions yields

\[\mathbf{x}(t+\Delta t) = 2\mathbf{x}(t)-\mathbf{x}(t-\Delta t)+\mathbf{a}(t)\Delta t^2 \]

which is an update that does not require any velocity information.

Finally, the Velocity Verlet integrator yields explicit positions and velocities by following the Verlet integrator idea for the velocities rather than positions:

\[\mathbf{x}(t+\Delta t) = \mathbf{x}(t) + \mathbf{v}(t)\Delta t + \frac{1}{2}\mathbf{a}(t)\Delta t^2\]

\[\mathbf{v}(t+\Delta t) = \mathbf{v}(t) + \frac{1}{2}\left[\mathbf{a}(t+\Delta t) + \mathbf{a}(t)\right]\]

Here the trick is to obtain \(\mathbf{a}(t+\Delta t)\) by evaluating the force acting on the particle at the new position, therefore \(\mathbf{x}(t+\Delta t)\) needs to be calculated before \(\mathbf{v}(t+\Delta t)\) can be obtained.

Using this background, you will now implement Python code to calculate the trajectory of \(n\) charged particles in an immovable box. All particles bounce off the box boundary without energy transfer, i.e. their velocity component orthogonal to the box surface is inverted.

You can develop and run code online without installation. Please submit your code with your report. It is recommended to read all tasks first before starting to develop code. Code examples given here will deliberately use beginner-level Python features to remain accessible, even at the cost of performance. Feel free to deviate from the suggested function naming and function signatures if you are familiar with Python. Implementations using NumPy/JAX/numba are accepted as well, provided they solve the same problem.

Part 1: Fundamentals

First, please familiarise yourself with the basics of Python, e.g. via the self-learning course in moodle or via the material for the block course in Python which contains all you need. Knowing NumPy or JAX is advantageous but not necessary to solve the tasks.
Let's create a single particle. Write a function init_universe(bounds: tuple[float, float]) -> tuple[tuple[float, float], tuple[float, float]] which randomly places a particle in the box defined by the tuple bounds and returns its position and a random velocity. The box is rectangular and goes from the origin to the point bounds.
Now the particle shall move. Write a function get_trajectory(total_time: float, time_step: float, bounds) -> list[tuple[float, float]] which uses the Euler integrator and init_universe to calculate the trajectory of a single particle by iteration over small steps \(\Delta t\) (given as time_step) until the total time \(T\) (total_time) is elapsed. This function returns a list of coordinates where the particle has been. You can visualise your trajectory using
```
bounds = (10, 10)
trajectory = get_trajectory(10, 0.01, bounds)

import matplotlib.pyplot as plt
plt.plot([_[0] for _ in trajectory], [_[1] for _ in trajectory])
plt.plot((0, 0, bounds[0], bounds[0], 0), (0, bounds[1], bounds[1], 0, 0), color="red")
```
So far the particle starts in the box and leaves it quickly. That should not happen. Now constrain your particle to the box by writing a new function enforce_boundaries(position: tuple[float, float], velocity: tuple[float, float], bounds: tuple[float, float]) -> tuple[tuple[float, float], tuple[float, float]] which makes sure the particle bounces off the box walls. You should assume that only one box boundary is crossed between time steps. Add this function to get_trajectory and verify that the particle stays in the box.
What happens if you increase the time step? Explore a wide range of time steps and document your findings by showing the trajectories. Explain your findings and discuss what they mean for simulating single particles in a box.
Describe what would be needed to be done if we had an irregular boundary, i.e. not a simple rectangle? Think of a circle or a polygon. No implementation is needed, just a description of how to deal with the reflection in that case.
Which sources of deviations from the true trajectory do you expect in this particular case? Is the Euler integrator sufficient for this problem?

Part 2: Application

Now generalize your code to multiple particles using NumPy. We will consider charged particles avoiding each other in the box, i.e. all charges are positive.

First write a function to calculate the Coulomb force acting on a set of charged particles: coulomb_repulsion_force(charges: np.ndarray, positions: np.ndarray) -> np.ndarray. You may use atomic units, i.e. the prefactor in the force expression is 1. The force for the \(i\)-th particle is

\[\mathbf{F}_i = q_i\sum_j q_j\frac{(\mathbf{x}_i - \mathbf{x}_j)}{|\mathbf{x}_i - \mathbf{x}_j|^3}\]
Write a function coulomb_energy(charges, positions) which calculates the Coulomb energy from the pairwise distances \(d_{ij}\):

\[E_\textrm{pot} = \frac{1}{2}\sum_i \sum_{j\neq i} \frac{q_iq_j}{d_{ij}}\]
Update the universe creation to generate \(n\) particles of random positions, velocities, masses and charges (which are the return values): init_charged_universe(n: int, bounds: tuple[float, float]) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray].
Implement the Euler integrator for this problem, writing get_euler_trajectory(n_particles: int, total_time: int, time_step: float, bounds: tuple[float, float]) -> tuple[np.ndarray, np.ndarray] which returns the trajectory and the total energy. This includes the kinetic energy:

\[E_\textrm{kin} = \frac{1}{2}\sum_i m_i \mathbf{v}_i^2\]
Write get_velocity_verlet_trajectory with the same function arguments as get_euler_trajectory to use velocity verlet instead.
Choose a number of particles between two and five and run both integrators such that the particles interacted with each other/the box at least a few times. Running the integrators should take roughly 10 seconds. Suggest a metric for the integration error which should be derived from the total energy data. Make a table with the integration error as a function of time step for the Euler integrator and the Velocity Verlet integrator and discuss your findings. (If the metric does not improve with smaller time steps, then either the integrator is faulty or the metric is flawed.)
Suggest a way how to implement the direction reversal upon touching the box in the Verlet algorithm where no velocities are available. No actual implementation needed.