\[\def\CC{\bf C} \def\QQ{\bf Q} \def\RR{\bf R} \def\ZZ{\bf Z} \def\NN{\bf N}\]

Day 6 Stability and Accuracy of Molecular Dynamics¶

through the tutorial there are jupyter cells that can allow you to run the tasks… feel free to add more if they are needed or use the console to run.

Setup (optional)¶

This step is needed only if you run on google colab. Instructions may work on other systems but are not tested.

[ ]:

! pip install data_tutorials
! apt install gfortran

Now bring the data needed for the exercises.

[ ]:

from data_tutorials.data import get_data
get_data(
    url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_6/Chaos/",
    filename=["plots.py","chaos.f90"],
    folder=".",
)

Introduction¶

MD algorithms are traditionally very simple, but they are not all equally accurate or stable. This exercise allows you to explore some of the aspects of stability and accuracy for a small selection of them. The algorithms we shall look at are:

Velocity Verlet
Leapfrog
Verlet’s original method
Rahman’s algorithm

All the Verlet variants are simple, explict algorithms and are widely used. Rahman’s algorithm is a predictor-corrector algorithm, which is now of historical interest only. It was used by Rahman in his pioneering work on Argon. We shall investigate the accuracy and stability of these algorithms, using velocity Verlet as the standard.

The files for this exercise are in the directory ~/WORKSHOP/Day_6/Chaos;

Formulae¶

It will be useful to have to hand the relevant formulae for the different integrators.

Velocity Verlet:

\[\vec{r}_{n+1} = \vec{r}_{n} + \vec{v}_n \Delta t + \vec{n}_{n} \frac{\Delta t^{2}}{2} + O(\Delta t^{4})\]

\[\vec{v}_{n+1} = \vec{v}_{n} + \frac{\Delta t}{2}(\vec{a}_{n} + \vec{a}_{n+1}) + O(\Delta t^{3})\]

or equivalently:

\[\vec{v}_{n+\frac{1}{2}} = \vec{v}_{n} + \frac{\Delta t}{2}\vec{a}_{n} + O(\Delta t^{2})\]

\[\vec{r}_{n+1} = \vec{r}_{n} + \vec{v}_{n+\frac{1}{2}} \Delta t + O(\Delta t^{4})\]

\[\vec{v}_{n+1} = \vec{v}_{n+\frac{1}{2}} + \frac{\Delta t}{2}\vec{a}_{n+1} + O(\Delta t^{3})\]

Leapfrog:

\[\vec{v}_{n+\frac{1}{2}} = \vec{v}_{n-\frac{1}{2}} + \vec{a}_{n} \Delta t + O(\Delta t^{3})\]

\[\vec{r}_{n+1} = \vec{r}_{n} + \vec{v}_{n+\frac{1}{2}} \Delta t + O(\Delta t^{4})\]

Verlet:

\[\vec{r}_{n+1} = 2\vec{r}_{n} + \vec{r}_{n-1} + \vec{a}_{n} \Delta t^{2} + O(\Delta t^{4})\]

\[\vec{v}_{n} = \frac{1}{2 \Delta t}(\vec{r}_{n+1} - \vec{r}_{n-1}) + O(\Delta t^{2})\]

Rahman:

\[\vec{r}_{n+1} = \vec{r}_{n-1} + 2 \vec{v}_{n} \Delta t + O(\Delta t^{3})\]

\[\vec{v}_{n+1} = \vec{v}_{n} + (\vec{a}_{n} + \vec{a}_{n+1}) \frac{\Delta t}{2} + O(\Delta t^{3})\]

\[\vec{r}_{n+1} = \vec{r}_{n} + (\vec{v}_{n} + \vec{v}_{n+1}) \frac{\Delta t}{2} + O(\Delta t^{3})\]

In all the above formulae \(r_{n}\), \(v_{n}\) and \(a_{n}\) represent position, velocity and acceleration respectively, at the n’th time step, while \(\Delta t\) represents the time step. Terms like \(O(\Delta t^{2})\) indicate the order of error. Note that only one of these methods is self-starting, in the sense that only the velocity Verlet can begin immediately from a knowledge of \(r_{n}\), \(v_{n}\) and \(a_{n}\). The others require some additional calculation.

The Program Chaos¶

The CHAOS program is designed to simulate a single particle moving in a confining quartic potential of the form:

\[V(x,y) = x^{2}y^{2}(x^{2} + y^{2})\]

The long time trajectory shows random-like behaviour and will (we hope!) rapidly show up differences in trajectories calculated by different methods. All of the above integration algorithms are coded in the program, though you will only look at two of them at any one time; the velocity Verlet compared with one of the others in sequence.

Use the CHAOS and its input file IN. Take a good look at the program source (chaos.f90) and make sure you understand what it is supposed to do. Ask if you don’t! Pay particular attention to how each of the integration algorithms is started. Only the velocity Verlet method can start without extra calculation.

You may compile the program using the command:

[ ]:

! gfortran -o chaos chaos.f90

input looks like this

[ ]:

%%writefile IN
        Experiment number
     Number of time steps
        Print control
01       Time step

[ ]:

! ./chaos < IN

Where IN is the input file. You will be able to control the chaos program by setting the numbers in the IN file. The contents of the IN file are as follows, each is a number written on a separate line:

key The experiment number (see below) [integer] nsteps The required number of timesteps [integer] nprnt Print control (number of lines printed =nsteps/nprnt) [integer] tstep The time step [real]

You are advised to use a fresh copy of the supplied IN file every time you start a new experiment. If you accidentally use an IN file modified from a previous experiment, the things you were meant to observe may not occur!

The output from the program will mainly appear in an output file named OUT, though sometimes the average energy and fluctuation will be printed to the screen. The contents may vary with each experiment, as will be described below. It is a good idea to rename and store the OUT files for comparison between the different experiments below.

[ ]:

%%bash

cp OUT OUT.1

[ ]:

from plots import *

[ ]:

plotFile("OUT",xcol=1,ycol=2,xlabel="time",ylabel="y")
plotFile("OUT",xcol=1,ycol=3,xlabel="time",ylabel="y")

Experiment 1. Velocity Verlet and Leapfrog¶

Take the input file IN as you downloaded it. Run the program (it will run for 1000 time steps and will take only a few seconds) then you may examine the OUT file with your preferred editor. The file contains the following data on each line:

the integration time;
the x and y coordinates of the particle obtained using velocity verlet;
the x and y coordinates of the particle obtained using leapfrog; and
the difference in position between the two methods.

You should do the following:

Compare the coordinates obtained by the two methods. Do they agree closely or not at all? Is the difference in position a constant or is it increasing?

It is a good idea to plot the difference in position as a function of time using one of the graph plotters available. It should be apparent that the two algorithms follow the same trajectory quite closely but a gradual divergence takes place.

The default setting for the time step in the IN file is 0.01. Try a few experiments using larger and smaller values than this, to see if this causes the trajectories to depart from each other at a different rate. Beware of increasing the time step too much, or the program may crash!

It is also interesting to plot the x coordinate versus the y coordinate using a graph plotting program. You will need to cut the relevant two columns of data from the OUT file to do this, but the result will be interesting. You can use this method to compare the trajectories obtained from the different algorithms. Also, try increasing the number of timesteps to (say) 50,000 and plot the trajectory again. You will get an idea of how the potential energy function in this system confines the particle.

[ ]:

%%bash

Experiment 2. Velocity Verlet and Original Verlet¶

Edit a copy of the original IN file and set the ‘experiment key’ on the first line to 2. This will select the original Verlet algorithm for comparison with velocity Verlet. Run the code. You should obtain a similar result to the previous case but with differences in the actual numbers.

Once again you can experiment with the magnitude of time step and see if this changes the rate of divergence.

It should be apparent from these experiments that the Verlet algorithms generate closely similar trajectories (at least over the first few hundred time steps), though they are clearly diverging and after a while will follow significantly different trajectories.

It should also be apparent that the original Verlet method diverges slightly faster from the trajectory given by the velocity Verlet than does the leapfrog algorithm. This may be attributed to rounding in the machine arithmetic; Verlet’s original method adds a term of order \(\Delta t^{2}\) to the position, while both velocity verlet and leapfrog add a term in \(\Delta t\). This difference in arithmetic is enough to cause the difference in trajectory seen!

[ ]:

%%bash

Experiment 3. Velocity Verlet and Rahman’s Algorithm¶

Edit a fresh copy of the original IN file, setting the ‘experiment key’ to 3. This will select the Rahman predictor-corrector algorithm for comparison with velocity Verlet. Now run the code.

Compare again the particle trajectories and the deviations obtained. The behaviour may seem to resemble the previous cases, but the magnitude in the trajectory difference will be considerably larger. It will be clear that Rahman’s and Verlet’s methods provide trajectories that deviate significantly from each other, though over the first 1000 time steps the trajectories still look broadly the same. Of course this experiment alone does not tell us which of the two algorithms is the more accurate. i.e. which is giving the \(best\) trajectory.

Once again you should experiment with the magnitude of the time step and see how this affects the trajectory divergence. One thing you should find is that Rahman’s algorithm is more likely to crash at larger time steps. This is a common failing in predictor-corrector algorithms!

Rahman’s algorithm, being a predictor-corrector, is iterative. The number of iterations is determined by the variable icyc in the source code of chaos.f90, and is set to 2. Locate the assignment of icyc near the top of the source code and set it to a larger number and see if the trajectory obtained is significantly different. Note that each iteration demands a new calculation of the force, which makes this a relatively expensive algorithm computationally, so icyc should ideally be kept small.

[ ]:

%%bash

Experiment 4: Energy Conservation¶

Divergent trajectories are worrying, but how significant is this from a thermodynamic point of view? To investigate this, we will calculate the total energy of the particle as given by velocity Verlet and Rahman and compare them.

Take an original copy of the IN file and set the ‘experiment key’ to 4. Then run the CHAOS program. At then end of the run, the program will print on screen the calculated fluctuation in energy of the two methods. It will be apparent that, over this range of time steps, the two algorithms give the same energy, with about the same fluctuation.

(Note: In comparing the two algorithms in this way, over a relatively short trajectory, we must bear in mind that a quantitative comparison is only meaningful if the two trajectories are reasonably close. If they have widely different trajectories, comparison is invalid because they will be exploring different parts of phase space. The previous experiments show that the two trajectories in this case are close enough.)

The OUT file from this experiment contains the following data on each line:

The integration time;
The energy calculated by velocity Verlet;
The energy obtained by Rahman.

Now plot the energies are against time. It will appear that the two algorithms have a different pattern of fluctuations: where the deviations from the average energy are largest, the velocity Verlet algorithm tends to oscillate about the mean, whereas Rahman’s algorithm shows a single large excursion from the mean. It is not too fanciful to imagine that this difference may have consequences for the long term stability of these algorithms.

To check this hypothesis, run the same simulation for 10,000 steps (ten times longer) and plot the result. The Rahman algorithm will show a marked drift from the mean, while Verlet’s remains quite stable.

It is worth repeating these two simulations, with the time step set to double the original value. The run of 1000 time steps will show that the fluctuation in the energy for the Rahman algorithm is much larger than for Verlet. (We are assuming that the trajectories are reasonably similar! – Check this!) Furthermore the plot of energy against time shows a distinct drift in the Rahman case. The run of 10,000 steps will not finish properly as the Rahman algorithm breaks down. Nevertheless it is useful to plot the energy against time (using a truncated OUT file!) to see how the two algorithms compare up to the point of breakdown.

[ ]:

%%bash

Experiment 5: Time Reversibility¶

The experiments so far suggest that the velocity Verlet algorithm has a built-in long term stability, while Rahman’s does not and (in common with many predictor-corrector algorithms) accumulates significant errors with time.

Why is the velocity Verlet algorithm so stable? One reason is that is possesses the property of time reversibility, which is implicit in Newton’s laws of motion. Thus changing the sign of time has no effect on the essential physics of the system. Rahman’s algorithm, with its iterative corrector, does not have this property, which implies a deficiency of some kind. For example, solving the equations of motion using the Verlet algorithms ensures the Boltzmann distribution function is constant in time, but using Rahman’s algorithm does not. If the distribution function is constant in time, the mechanical properties of the system, such as energy, will not drift.

You can test the time reversibility of the velocity Verlet and Rahman algorithms with the program CHAOS. Start with a copy of the original IN file, and edit the ‘experiment key’ to 5. Then run the program. CHAOS will integrate the equations of motion in a positive time direction for half the simulation (500 steps) and then reverse the time direction (by reversing the velocity at the half-way point) and integrate back to zero time.

Now examine the OUT file. You need only look at the first and last lines in the file. You will find that for the time step of 0.01, the velocity Verlet algorithm has returned the particle exactly to its starting position. However, the Rahman algorithm does not quite achieve this.

Next you should try increasing the time step in the IN file and repeating the experiment several times. It will surprise you how large the time step can be and yet preserve the time reversiblity.

Alternatively, if you make the time step successively smaller you will, of course, eventually cause time reversibility apparently to appear in the Rahman algorithm as well, but this is purely circumstantial, and not an intrinsic property of the algorithm itself.

[ ]:

%%bash

[ ]: