Constraints¶
Colab Setup (optional)¶
The next cells should only be needed if you use Google Colab, so we test the environment for this. Commands may work in other environments too but are not tested. If you need to try, just replace the following by SETUP = True.
[ ]:
import sys
SETUP = "google.colab" in sys.modules
if SETUP:
print('Fetching files for Google Colab')
[ ]:
if SETUP:
# Install additional packages
! pip install h5py data_tutorials weas_widget ase meson ninja
! apt install gfortran libhdf5-dev
[ ]:
if SETUP:
from data_tutorials.data import get_data
# Get the python modules and top-level meson.build needed for this tutorial
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_6/Constraints/",
filename=["hdf5_module.py","meson.build"],
folder=".",
)
# Get the starting configuration, in the DATA subdirectory
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_6/Constraints/DATA/",
filename=["md_springs.xyz", "md_constraints.xyz"],
folder="DATA/",
)
# Get the Fortran source files, and meson.build file, in the src subdirectory
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_6/Constraints/src/",
filename=["meson.build","config_io_module.f90","force_module.f90","hdf5_module.f90","maths_module.f90","md_constraints.f90","md_module.f90","md_springs.f90"],
folder="src/",
)
Introduction¶
How do we simulate molecules which have strong internal bonds, as well as nonbonded interactions (such as Lennard-Jones) between the atoms? One approach is to allow the bond lengths etc to evolve under the influence of appropriate terms (for example, strong harmonic springs) in the potential energy. However, this requires a correspondingly small timestep, which can make the program very inefficient.
An alternative approach is to fix the bond lengths by applying constraints. The introduction of a practical way of doing this is one of the most important developments in molecular dynamics. It enables the simulation of molecules of great complexity with a high degree of stability and reliability. The original SHAKE algorithm was derived by JP Ryckaert et al., J Comput Phys, 23, 327 (1977), based on the Verlet algorithm: at each step, an iterative method is used to ensure that the bond lengths, angles, or other geometrical quantities, satisfy the imposed constraints. Somewhat later, a version based on the velocity Verlet algorithm was derived by HC Andersen, J Comput Phys, 52, 24 (1983), called RATTLE. This is essentially equivalent to SHAKE for the positions, but also ensures that the atomic velocities are consistent with the constraint conditions. Subsequently, several other algorithms have been developed to tackle constraints, with better performance in certain circumstances (e.g. LINCS, SETTLE). To compare with RATTLE, here we look at MILC SHAKE, which applies to linear chain topologies.
If we still wish to simulate the bond vibrations explicitly, instead of using constraints, one idea is to use multiple timesteps (MTS). This permits the expensive intermolecular interactions to be evaluated much less frequently than the rapidly-changing intramolecular spring forces, making the program more efficient; see e.g. M Tuckerman et al., J Chem Phys, 97, 1990 (1992). However, there are still dangers and limitations on the timesteps that can be used. This method is illustrated at the end of the notebook.
One should bear in mind that these two models (springs and constraints) are not equivalent to each other, even in the limit of very strong springs, and that they are both just approximations to the real system. Bond vibrations in most molecules are best described by quantum mechanics, not classical mechanics.
Preliminaries¶
Start by importing some useful Python modules and functions. Here are links to NumPy, Matplotlib, pathlib, and shutil, in case you need any details. We shall use ase, and the weas-widget tool to visualize atomic configurations. The remaining modules are supplied as part of this workshop. We also apply some Matplotlib settings.
[ ]:
# Standard and external modules
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from pathlib import Path
from shutil import copy
from itertools import cycle, islice
from ase.atoms import Atoms
from ase.io import read
from ase.data import vdw_radii, atomic_numbers
from weas_widget import WeasWidget
# Modules supplied in this workshop
from hdf5_module import read_hdf5_file
# Settings for Matplotlib plots
plt.style.use(['seaborn-v0_8-talk','seaborn-v0_8-darkgrid','seaborn-v0_8-colorblind'])
plt.rc('image',cmap='viridis') # Default colour map
plt.rc('figure',figsize=(8,5)) # Default figure dimensions in inches
plt.rc('figure.constrained_layout',use=True) # Default figure layout option
plt.rc('legend',frameon=True,framealpha=1.0) # Default appearance of figure legend
plt.rc('hist',bins=100) # Default number of bins to use in histograms
# Create a marker cycle for plotting, similar to prop_cycle
prop_cycle = plt.rcParams['axes.prop_cycle']
colors = prop_cycle.by_key()['color']
markers = ['o','s','^','D','v']
markers = list(islice(cycle(markers),len(colors))) # Match length of colors, cycle if necessary
marker_cycle = plt.cycler(color=colors,marker=markers)
# Setting for weas-widget viewer
weas_shape = {"width":"800px","height":"600px"} # Default shape for weas widget
The model¶
In this exercise we shall consider a simple system: a flexible chain molecule of \(N=64\) atoms in isolation. The nonbonded interaction between all pairs of atoms (excepting bonded pairs) is the full Lennard-Jones potential, \begin{equation*} u_{\mathrm{LJ}}(r) = 4\varepsilon \left[\left( \frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^6 \right]. \end{equation*} Here \(r\) is the distance between the atoms, \(\varepsilon\) is an energy characterizing the strength of the interaction, and \(\sigma\) is a length scale that characterizes the size of the atoms. As usual, reduced units are employed, so \(\varepsilon=1\) and \(\sigma=1\), and in addition Boltzmann’s constant is taken to be unity, \(k_{\mathrm{B}}=1\). Note that there are NO periodic boundary conditions in this system, and there is NO potential cutoff!
In the harmonic spring model, there are \(N-1\) bonds between successive atoms \(i\) and \(j=i+1\) in the chain, each having a potential energy \(u_{\mathrm{spring}}(r_{ij})\) where \(r_{ij}=|\mathbf{r}_{ij}|\), \(\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}\), and \begin{equation*} u_{\mathrm{spring}}(r) = \tfrac{1}{2} \kappa ( r - d )^2 . \end{equation*} In this potential \(d\) is the bond length, set to \(d=0.85\). and \(\kappa\) is the spring constant, set here to a fairly large value \(\kappa=10000\) (all in Lennard-Jones reduced units).
In the constrained-bond model, there are \(N-1\) constraints of the form \(r_{ij}=d\) between successive atoms \(i\) and \(j=i+1\). It is convenient to write these as \begin{equation*} r_{ij}^2 = \mathbf{r}_{ij}\cdot\mathbf{r}_{ij} = d^2 . \end{equation*} Once again, \(d=0.85\). By time-differentiating this equation, you can see that there is a corresponding constraint on the relative velocities \(\mathbf{v}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}\): \begin{equation*} \mathbf{r}_{ij}\cdot\mathbf{v}_{ij} = 0. \end{equation*}
In this very basic model, we shall take all the atoms to have the same mass, \(m\). Bear in mind that this will simplify several of the equations, particularly in the constraint algorithms. In the program we actually set \(m=1\).
Getting started¶
The next two cells will build, in the build directory, both the codes to be used here: md_constraints and md_springs.
[ ]:
!meson setup build
[ ]:
!meson compile -C build
The programs¶
The program source files are in the subdirectory src. The main programs are md_springs.f90 and md_constraints.f90. The force and constraint routines appear in the files force_module.f90 and md_module.f90 respectively. You may like to take a quick look at these files, to see how the programs operate.
An additional module config_io_module.f90 handles configuration input, and hdf5_module.f90 deals with the HDF5 output. There is no particular need to study those files.
The RATTLE algorithm¶
Here we give a bit more detail on the algorithm in the rattle_a and rattle_b routines within the md_module.f90 file. Since this workshop is quite short, you might like to skip this section, and return to it later, at your convenience.
We follow the description of HC Andersen, J Comput Phys, 52, 24 (1983). Referring to the lecture, after implementing the first kick \(\mathbf{p}_i(t)\rightarrow\mathbf{p}_i'(t+\tfrac{1}{2}\Delta t)\), and the drift \(\mathbf{r}_i(t)\rightarrow\mathbf{r}_i'(t+\Delta t)\), without including the effects of bond-length constraints, we arrive at the following equations for the remaining constraint corrections: \begin{align*} \mathbf{p}_i(t+\tfrac{1}{2}\Delta t) &= \mathbf{p}_i'(t+\tfrac{1}{2}\Delta t) + \frac{m}{\Delta t}\sum_j\lambda_{ij}\mathbf{r}_{ij}(t), \\ \mathbf{r}_i(t+\Delta t) &= \mathbf{r}_i'(t+\Delta t) + \sum_j\lambda_{ij} \mathbf{r}_{ij}(t). \end{align*} Each constraint acts along a bond vector \(\mathbf{r}_{ij}(t)=\mathbf{r}_i(t)-\mathbf{r}_j(t)\), and we just need to determine the multipliers \(\lambda_{ij}\). For our chain, the sum over \(j\) is limited to, at most, two terms \(j=i\pm 1\), i.e. all the atoms, except the ones at the end, are involved in two bonds. There are \(N-1\) constraint equations, and \(N-1\) unknowns \(\lambda_{ij}\), but the equations are coupled together.
However, in rattle_a, the constraints are tackled one at a time in an iterative process. At each iteration, there is a loop over all the constraints. For a given constraint \(ij\) where \(j=i+1\), set \(\lambda=\lambda_{ij}=\lambda_{ji}\), and \(\mathbf{q}_{ij}=\mathbf{r}_{ij}(t)=-\mathbf{r}_{ji}(t)\). Ignoring the effects of other constraints, the above equations lead to an updating scheme \begin{align*}
\mathbf{p}_i'' &=
\mathbf{p}_i' + \frac{m}{\Delta t}\lambda\,\mathbf{q}_{ij}, \quad
\mathbf{r}_i'' = \mathbf{r}_i' + \lambda\,\mathbf{q}_{ij},
\\
\mathbf{p}_j'' &=
\mathbf{p}_j' - \frac{m}{\Delta t}\lambda\,\mathbf{q}_{ij}, \quad
\mathbf{r}_j'' = \mathbf{r}_j' - \lambda\,\mathbf{q}_{ij} .
\end{align*} where \('\) and \(''\) simply denote successive estimates of momenta and positions.
It follows that \(\mathbf{r}_{ij}'' = \mathbf{r}_{ij}' + 2\lambda\, \mathbf{q}_{ij}\). In our simple example, all bond lengths are \(d\), and so the equation for \(\lambda\) is \begin{equation*}
\mathbf{r}_{ij}'' \cdot \mathbf{r}_{ij}'' =
\mathbf{r}_{ij}' \cdot \mathbf{r}_{ij}' + 4\mathbf{r}_{ij}' \cdot \mathbf{q}_{ij}\,\lambda
+ 4\,d^2\,\lambda^2 = d^2 .
\end{equation*} Dropping the (small) quadratic term, this becomes \begin{equation*}
A\lambda = \sigma
\quad\text{where}\quad A=4\mathbf{r}_{ij}' \cdot \mathbf{q}_{ij},
\quad\text{and}\quad \sigma=d^2 - \mathbf{r}_{ij}' \cdot \mathbf{r}_{ij}' .
\end{equation*} The corresponding equation \(\lambda=\sigma/A\) appears in the rattle_a routine of md_module.f90. The corrections are used to update the positions and momenta of the two atoms. Then, the program moves on to the next constraint. Of course, moving \(i\) to satisfy the \(i,i+1\) constraint may cause the \(i-1,i\) constraint to be violated once more. After all constraints have been considered, the loop begins the next iteration, during which all pairs will be
checked and possibly updated again. The loop concludes when all constraints are satisfied to within a prescribed tolerance.
The second constraint correction is made after implementing the second kick \(\mathbf{p}_i(t+\frac{1}{2}\Delta t)\rightarrow\mathbf{p}_i'(t+\Delta t)\) using the freshly calculated nonbonded forces at time \(t+\Delta t\), but without including the effects of constraints. These constraints act along the new bond vectors \(\mathbf{r}_{ij}(t+\Delta t)\): \begin{equation*}
\mathbf{p}_i(t+\Delta t) = \mathbf{p}_i'(t+\Delta t) + \sum_j \mu_{ij} \mathbf{r}_{ij}(t+\Delta t) .
\end{equation*} Once more, an iterative scheme is used. For bond \(ij\), setting \(\mu=\mu_{ij}=\mu_{ji}\), \(\mathbf{r}_{ij}=\mathbf{r}_{ij}(t+\Delta t)=-\mathbf{r}_{ji}(t+\Delta t)\), and ignoring the effects of the other constraints, we get an equation for the velocity difference \(\mathbf{v}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}\) \begin{equation*}
\mathbf{p}_i'' = \mathbf{p}_i' + \mu\, \mathbf{r}_{ij}, \quad
\mathbf{p}_j'' = \mathbf{p}_j' - \mu\, \mathbf{r}_{ij}, \quad
\Rightarrow\quad m\mathbf{v}_{ij}'' = m\mathbf{v}_{ij}' + 2\mu\, \mathbf{r}_{ij} .
\end{equation*} If all the masses were not the same, as in our simple example here, this equation would be more complicated. This gives \begin{equation*}
m\mathbf{v}_{ij}'' \cdot \mathbf{r}_{ij} =m\mathbf{v}_{ij}'\cdot \mathbf{r}_{ij} + 2\mu\, \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} = 0
\end{equation*} or \begin{equation*}
B\mu = \tau,
\quad\text{where}\quad B = 2\,\mathbf{r}_{ij} \cdot \mathbf{r}_{ij} = 2d^2,
\quad\text{and}\quad \tau= -m\,\mathbf{v}_{ij}'\cdot \mathbf{r}_{ij}.
\end{equation*} The corresponding equation \(\mu=\tau/B\) is in the rattle_b routine of md_module.f90. Once the correction \(\pm \mu \mathbf{r}_{ij}\) has been applied to the momenta, the rest of the constraints are considered. Subsequent iterations reconsider all the bonds again, and the iterative loop concludes when all the relative velocities satisfy their constraints to within a prescribed tolerance.
Input data¶
Two configuration files are supplied, DATA/md_springs.xyz and DATA/md_constraints.xyz, one for each model. These files are in extended XYZ format, and look like this:
n
Properties=species:S:1:pos:R:3:velo:R:3
X r(1,1) r(2,1) r(3,1) p(1,1) p(2,1) p(3,1)
X r(1,2) r(2,2) r(3,2) p(1,2) p(2,2) p(3,2)
X r(1,3) r(2,3) r(3,3) p(1,3) p(2,3) p(3,3)
: : : : : :
X r(1,n) r(2,n) r(3,n) p(1,n) p(2,n) p(3,n)
The first line gives the number of atoms, the second line gives details of the columns in the rest of the file, and the subsequent lines give the species, (x,y,z)-coordinates, and (x,y,z)-momenta, of each atom. Again, a reminder, there is no box. Everything is in LJ reduced units. We use X as the atomic species for each atom: in realistic simulations, a chemical element symbol would appear instead, and lengths would most commonly be given in Angstroms.
In both cases, the systems have been equilibrated at a temperature \(T\approx1.0\). They should also have zero total linear momentum \(\sum_i \mathbf{p}_i=\mathbf{0}\), and zero total angular momentum \(\sum_i \mathbf{r}_i\times\mathbf{p}_i=\mathbf{0}\) about the centre of mass, which is located at the origin. The programs print these values at the start and the end of the run, for you to check. Note that, in the absence of boundary conditions, both linear and angular momentum are conserved.
The following cell should allow you to visualize DATA/md_springs.xyz, converted into units appropriate for argon atoms, using weas-widget. The viewer size may be changed with weas_shape (defined in Preliminaries above) or with the widget’s ⛶ button (may not work in Colab). For clarity, ball-and-stick format (model_style=1) is selected, but in the GUI you may switch to Atoms/Model Style/Ball and Atoms/Radius Type/VDW to
give a space-filling picture.
[ ]:
this_dir = Path.cwd() # Current working directory
data_dir = this_dir / 'DATA' # Data directory
old_conf = data_dir / 'md_springs.xyz' # Supplied configuration file
new_conf = copy(old_conf,this_dir) # Copy supplied file into current directory
atoms=read(new_conf)
atoms.numbers[:] = atomic_numbers['Ar'] # Visualize as argon atoms
sigma = vdw_radii[atomic_numbers['Ar']]*2 # Argon diameter in Angstroms
atoms.positions *= sigma # Scale all positions
viewer=WeasWidget(from_ase=atoms,viewerStyle=weas_shape)
viewer.avr.bond.add_bond_pair('Ar','Ar',max=0.89*sigma)
viewer.avr.model_style=1
viewer
The programs both accept user input of run parameters via a namelist, which allows them to be run easily with default values, as well as allowing the user to specify selected values through keywords. In the following cells, these programs are run with default parameters. Start with the md_springs program.
[ ]:
!echo '&nml /' | build/md_springs
The program reports (at the start and end of the run) the worst deviation of bond length from the specified value: this should be quite small, due to the strong springs (we expect something of the order of \(\sqrt{k_{\text{B}}T/\kappa}\approx 0.01\)). Notice, however, that the worst (i.e. largest) time derivative of the bond length is not correspondingly small. Is this what you expect? Give this a moment’s thought.
Next, we analyze the run output. We are interested in the energy conservation, and the measured temperature (calculated from the kinetic energy). The relevant simulation attributes (parameters), and step-by-step datasets, may be extracted from the HDF5 file written by the simulation program, and this is done in next cells.
[ ]:
attr, data = read_hdf5_file('md_springs.hdf5')
[ ]:
print(attr['Title'].astype(str))
N = attr['N']
print(f'Number of atoms N = {N:10d}')
We expect the average kinetic energy for each quadratic degree of freedom to be \(\frac{1}{2}k_{\text{B}}T\). So for this system we can use the average total kinetic energy to estimate the temperature. Let \(N_\text{free}\) be the number of degrees of freedom: this is \(3\) for each atom, representing translation in each coordinate direction, but remember that we must subtract \(6\) from the total to account for the conserved total linear momentum and, in this case, angular momentum (feel free to ask if you need more explanation!).
[ ]:
Kavg = data['K'].mean()
Nfree = (3*N-6)
Tavg = 2*Kavg/Nfree
print(f'Average temperature T = {Tavg:10.4f}')
Hopefully the average temperature is not too far from \(1\). Actually, in this type of sytem, equilibration between the different types of energy in the system can be quite slow, so the agreement may not be precise.
The potential energy has two parts: \(U\) (the non-bonded potential) and \(V\) (the spring bonds). They both contribute to the total energy, along with the kinetic energy \(K\). We calculate the mean energy per atom \(e=E/N\), and the root-mean-square deviation of this same quantity.
[ ]:
e = ( data['K'] + data['U'] + data['V'] ) / N
eavg = e.mean()
eRMS = e.std()
print(f'Average energy/atom e = {eavg:10.5f}')
print(f'RMS deviation eRMS = {eRMS:10.1e}')
Hopefully the energy fluctuations are satisfactorily small. We shall look further into this, shortly.
Now we turn our attention to the constraints program. First, take a look at the supplied configuration.
[ ]:
old_conf = data_dir / 'md_constraints.xyz' # Supplied configuration file
new_conf = copy(old_conf,this_dir) # Copy supplied file into current directory
atoms=read(new_conf)
atoms.numbers[:] = atomic_numbers['Ar'] # Visualize as argon atoms
sigma = vdw_radii[atomic_numbers['Ar']]*2 # Argon diameter in Angstroms
atoms.positions *= sigma # Scale all positions
viewer=WeasWidget(from_ase=atoms,viewerStyle=weas_shape)
viewer.avr.bond.add_bond_pair('Ar','Ar',max=0.89*sigma)
viewer.avr.model_style=1
viewer
Once more, we use default parameters, defined within the program. Note that the default timestep dt is \(10\) times as long as for md_springs, and the default number of timesteps is a factor \(10\) smaller.
[ ]:
!echo '&nml /' | build/md_constraints
Now the program prints out the cumulative average number of iterations in both constraint stages, alongside the step number and CPU time, as a guide.
The worst bond length deviation figures are extremely small, and this time the worst bond length derivatives are also small. This is what we should expect if the constraints are being applied correctly. This highlights an important difference between the models. The momentum distribution for the model with springs is just the standard Maxwell-Boltzmann distribution, despite the fact that the bond lengths are staying very close to the desired values.
Again, we import the key simulation parameters and datasets, for analysis.
[ ]:
attr, data = read_hdf5_file('md_constraints.hdf5')
[ ]:
print(attr['Title'].astype(str))
N = attr['N']
print(f'Number of atoms N = {N:10d}')
The energy conservation may be checked as before. There is just one kind of potential energy, \(U\), the non-bonded potential.
[ ]:
e = ( data['K'] + data['U'] ) / N
eavg = e.mean()
eRMS = e.std()
print(f'Average energy/atom e = {eavg:10.4f}')
print(f'RMS deviation eRMS = {eRMS:10.1e}')
Temperature and degrees of freedom¶
Now it is your turn to do something: compute the kinetic temperature. This will be similar to the previous calculation, but there is a crucial difference: the number of degrees of freedom \(N_\text{free}\) is significantly different in this system, due to the bond constraints. Each constraint reduces the degrees of freedom by one.
Put the correct formula for \(N_\text{free}\) into the following cell. Your answer for T should be around \(1.0\).
[ ]:
Kavg = data['K'].mean()
Nfree = (3*N-6) # Correct this expression for Nfree
Tavg = 2*Kavg/Nfree
print(f'Kinetic temperature T = {Tavg:10.4f}')
Comparing performance¶
Now you are going to investigate energy conservation, and program speed, for the two approaches.
The following cell does several runs of md_springs, all from the same starting configuration. The time step dt is varied, and the number of steps nstep is chosen to keep the product nstep*dt constant, as specified by t_run. After each run the log and HDF5 output files are renamed according to run number, and moved into a subdirectory, to save them for further analysis. These runs will just take a few minutes: wait until they are all complete before proceeding.
[ ]:
result_dir = this_dir / 'SPRINGS/' # Subdirectory to contain results
result_dir.mkdir(exist_ok=True) # Create subdirectory
hdf5_file = this_dir / 'md_springs.hdf5' # HDF5 file produced by program
log_file = this_dir / 'md_springs.log' # Log file produced by program
old_conf = data_dir / 'md_springs.xyz' # Supplied configuration file
t_run, nsteps = 50.0, [10000,25000,50000,100000] # Use these values for all sets of runs
print('Run nstep dt')
for run, nstep in enumerate(nsteps):
dt = t_run/nstep
copy(old_conf,this_dir)
! echo "&nml nstep=$nstep, dt=$dt /" | build/md_springs >& md_springs.log
name = str(run)+'.hdf5'
hdf5_file.rename(result_dir/name) # Move HDF5 file into result directory, renaming it
name = str(run)+'.log'
log_file.rename(result_dir/name) # Move log file into result directory, renaming it
print(f'{run:3d}{nstep:10d}{dt:10.4f} done')
print('All done')
The following cell defines a function which gathers the important information from the HDF5 files: the timestep dt, CPU time taken CPU, and RMS energy per atom eRMS. Hopefully it is self-explanatory, but feel free to ask if anything is not clear. The function returns a dictionary of NumPy arrays, with the corresponding keys.
[ ]:
def gather_data ( dir ):
"""Gathers desired data from all HDF5 files within the specified directory.
Argument
--------
dir : Path
Pathlib Path to directory in which to search for HDF5 files
Returns
-------
dictionary of NumPy arrays
Keys in the dictionary correspond to timestep, CPU time, and RMS energy per atom
Successive array elements correspond to each simulation HDF5 file matching the pattern
"""
files = sorted(dir.glob('*.hdf5')) # Sorted list of HDF5 files within directory
d = {'dt':[],'CPU':[],'eRMS':[]} # Initial dictionary of empty lists
for file in files: # Loop over files
attr,data=read_hdf5_file(file) # Read from the HDF5 file
d['dt'].append(attr['dt']) # Append timestep to corresponding list
d['CPU'].append(attr['CPU']) # Append CPU time to corresponding list
e = np.zeros_like(data['K']) # Zero total energy/atom array
for energy in data.values(): # Loop over datasets (all types of energy)
e = e + energy # Add energy contribution into total
e = e / attr['N'] # Divide by number of atoms
d['eRMS'].append(e.std()) # Append RMS energy per atom to corresponding list
for key in d: # Loop over dictionary items
d[key] = np.array(d[key]) # Convert each list into a NumPy array
return d # Return dictionary of NumPy arrays
We gather the data from the above runs into a Python dictionary named results, and also store this as an item in a dictionary all_results, which will build up similar data for all the methods we are comparing.
[ ]:
results = gather_data(result_dir)
all_results = { 'springs': results }
The next cell plots results['eRMS'] vs results['dt'] on log-log axes. You should see that the points lie roughly on a straight line of slope \(\approx 2\). There is a line of the form \(e_{\text{RMS}} \propto \Delta t^2\) on the same graph, as a guide to the eye. This is expected for velocity Verlet: the root-mean-square energy fluctuations should be approximately proportional to the square of the time step.
[ ]:
fig, ax = plt.subplots()
ax.grid(True,which='both')
ax.set_xscale('log')
ax.set_xlabel(r'$\Delta t$')
ax.set_yscale('log')
ax.set_ylabel(r'$e_{\mathrm{RMS}}$')
ax.xaxis.set_major_locator(ticker.NullLocator())
ax.xaxis.set_minor_locator(ticker.LogLocator(subs='all'))
ax.xaxis.set_minor_formatter(ticker.LogFormatterSciNotation(minor_thresholds=(5,0.5)))
# Plot energy conservation vs timestep
ax.set_prop_cycle(marker_cycle)
ax.plot(results['dt'],results['eRMS'],ls='None',label='springs')
# Draw a crudely fitted line with slope 2 as a guide
ax.set_prop_cycle(None)
c=np.average(results['eRMS']/results['dt']**2)
dt=np.array(ax.get_xlim())
ax.plot(dt,c*dt**2,'--')
ax.legend();
Constraints¶
It will be interesting to carry out a similar series of runs, using constraints. If you have a look at the main loop in md_constraints.f90 you should see that the first step of velocity Verlet, in which the momenta are half-advanced and the positions are completely advanced (without constraints), is followed by a call to constraints_a, which points to the rattle_a routine in md_module.f90. This routine adjusts the positions along the bond vectors defined by the old (i.e.
un-advanced) coordinates, and makes some consistent adjustments to the momenta. The second half-kick of the momenta is followed by constraints_b, which points to the rattle_b routine, and makes momentum adjustments along bond vectors defined by the new positions.
The following cell carries out runs for different timesteps, with the same total run time, all from the same starting configuration. As before, this will not take long, you should wait until all the runs are complete before proceeding. The output files are moved into an appropriately-named subdirectory.
[ ]:
result_dir = this_dir / 'RATTLE/' # Subdirectory to contain results
result_dir.mkdir(exist_ok=True) # Create subdirectory
hdf5_file = this_dir / 'md_constraints.hdf5' # HDF5 file produced by program
log_file = this_dir / 'md_constraints.log' # Log file produced by program
old_conf = data_dir / 'md_constraints.xyz' # Supplied configuration file
print('Run nstep dt')
for run, nstep in enumerate(nsteps):
dt = t_run/nstep
copy(old_conf,this_dir)
! echo "&nml nstep=$nstep, dt=$dt /" | build/md_constraints >& md_constraints.log
name = str(run)+'.hdf5'
hdf5_file.rename(result_dir/name) # Move HDF5 file into result directory, renaming it
name = str(run)+'.log'
log_file.rename(result_dir/name) # Move log file into result directory, renaming it
print(f'{run:3d}{nstep:10d}{dt:10.4f} done')
print('All done')
When the runs are all done, collect the dt, CPU and eRMS data from the HDF5 files, as before.
[ ]:
results = gather_data(result_dir)
all_results['RATTLE'] = results
[ ]:
fig, ax = plt.subplots()
ax.grid(True,which='both')
ax.set_xscale('log')
ax.set_xlabel(r'$\Delta t$')
ax.set_yscale('log')
ax.set_ylabel(r'$e_{\mathrm{RMS}}$')
ax.xaxis.set_major_locator(ticker.NullLocator())
ax.xaxis.set_minor_locator(ticker.LogLocator(subs=(1,2,5)))
ax.xaxis.set_minor_formatter(ticker.ScalarFormatter())
# Plot energy conservation vs timestep for both cases
ax.set_prop_cycle(marker_cycle)
for method, results in all_results.items():
ax.plot(results['dt'],results['eRMS'],ls='None',label=method)
# Draw crudely fitted lines with slope 2 as a guide
ax.set_prop_cycle(None)
dt=np.array(ax.get_xlim())
for results in all_results.values():
c=np.average(results['eRMS']/results['dt']**2)
ax.plot(dt,c*dt**2,'--')
ax.legend();
Hopefully it will be clear that the energy conservation of RATTLE is more than an order of magnitude better, for each timestep. Really, the longer timesteps are not suitable for the model with spring bonds. In any case, it might be fairer to compare against CPU time. The following function will do this; we can use it again later.
[ ]:
def plot_all_results():
fig, ax = plt.subplots()
ax.grid(True,which='both')
ax.set_xscale('log')
ax.set_xlabel('CPU (seconds)')
ax.set_yscale('log')
ax.set_ylabel(r'$e_{\mathrm{RMS}}$')
ax.xaxis.set_major_locator(ticker.NullLocator())
ax.xaxis.set_minor_locator(ticker.LogLocator(subs=(1,2,5)))
ax.xaxis.set_minor_formatter(ticker.ScalarFormatter())
# Plot energy conservation vs CPU time for all cases
ax.set_prop_cycle(marker_cycle)
for method, results in all_results.items():
ax.plot(results['CPU'],results['eRMS'],ls='None',label=method)
# Draw crudely fitted lines with slope -2 as a guide
ax.set_prop_cycle(None)
cpu=np.array(ax.get_xlim())
for results in all_results.values():
c=np.average(results['eRMS']*results['CPU']**2)
ax.plot(cpu,c/cpu**2,'--')
ax.legend();
[ ]:
plot_all_results()
The RATTLE approach is still better, but the difference is not so dramatic. This is because the constraint iterations are more expensive than the spring force evaluation.
An alternative to RATTLE, for a simple linear chain such as this, is the MILC SHAKE algorithm, due to AG Bailey et al., J Comput Phys, 227, 8949 (2008); see also AG Bailey et al., Comput Phys Commun, 180, 594 (2009). If you are curious about this algorithm, more details are given at the end of this notebook. Suffice it to say that the routines milcshake_a and milcshake_b, within
md_module.f90, do exactly the same job as their rattle counterparts, but more efficiently. However, the algorithm only applies to linear chains.
This is how we run the program.
[ ]:
old_conf = data_dir / 'md_constraints.xyz' # Supplied configuration file
copy(old_conf,this_dir) # Copy supplied file into current directory
!echo '&nml milcshake=.true. /' | build/md_constraints
Notice that the average iteration counts per step are much reduced. In fact, part B of the algorithm only requires a single iteration.
Now we look at the performance of this algorithm with a sequence of runs, similar to before. Once again, the output files are moved into a subdirectory.
[ ]:
result_dir = this_dir / 'MILCSHAKE/' # Subdirectory to contain results
result_dir.mkdir(exist_ok=True) # Create subdirectory
hdf5_file = this_dir / 'md_constraints.hdf5' # HDF5 file produced by program
log_file = this_dir / 'md_constraints.log' # Log file produced by program
old_conf = data_dir / 'md_constraints.xyz' # Supplied configuration file
print('Run nstep dt')
for run, nstep in enumerate(nsteps):
dt = t_run/nstep
copy(old_conf,this_dir)
! echo "&nml nstep=$nstep, dt=$dt, milcshake=.true. /" | build/md_constraints >& md_constraints.log
name = str(run)+'.hdf5'
hdf5_file.rename(result_dir/name) # Move HDF5 file into result directory, renaming it
name = str(run)+'.log'
log_file.rename(result_dir/name) # Move log file into result directory, renaming it
print(f'{run:3d}{nstep:10d}{dt:10.4f} done')
print('All done')
[ ]:
results = gather_data(result_dir)
all_results['MILC SHAKE'] = results
[ ]:
plot_all_results()
The MILC SHAKE algorithm gives essentially the same energy conservation as RATTLE for the same timestep (provided the constraint tolerances are the same) but is faster.
Multiple timesteps¶
If time is short, you can skip this section.
If you have time, make a comparison with the multiple time step (MTS) approach to the model with spring bonds. To do this, we run the md_springs program with parameters such as nstep=10000, n_mts=10, instead of nstep=100000. The essential point of difference lies in the advancement of coordinates and momenta within the main loop of md_springs.f90. Take a look at this. If n_mts is not equal to 1, the kick half-steps, acting on the momenta, are applied separately for the
nonbonded forces f (which are recalculated at intervals n_mts*dt), and the spring bonds g (which are computed at intervals dt). If n_mts is equal to 1, you should be able to see that this is identical to the standard velocity Verlet algorithm, which we have explicitly (and redundantly) written out in the if ( n_mts == 1 ) then branch within the main loop.
The following cell carries out several runs, moving the output files into a subdirectory as before. As always, wait until all the runs have finished before proceeding.
[ ]:
result_dir = this_dir / 'MULTISTEP/' # Subdirectory to contain results
result_dir.mkdir(exist_ok=True) # Create subdirectory
hdf5_file = this_dir / 'md_springs.hdf5' # HDF5 file produced by program
log_file = this_dir / 'md_springs.log' # Log file produced by program
old_conf = data_dir / 'md_springs.xyz' # Supplied configuration file
n_mts = 10
print('Run nstep n_mts*dt n_mts*nstep dt')
for run, nstep in enumerate(nsteps):
dt = t_run/(nstep*n_mts)
copy(old_conf,this_dir)
! echo "&nml n_mts=$n_mts, nstep=$nstep, dt=$dt /" | build/md_springs >& md_springs.log
name = str(run)+'.hdf5'
hdf5_file.rename(result_dir/name) # Move HDF5 file into result directory, renaming it
name = str(run)+'.log'
log_file.rename(result_dir/name) # Move log file into result directory, renaming it
print(f'{run:3d}{nstep:15d}{dt*n_mts:10.4f}{n_mts*nstep:15d}{dt:10.5f}')
print('All done')
Collect together the relevant results, just as before.
[ ]:
results = gather_data (result_dir)
all_results['springs MTS'] = results
Try plotting the ['eRMS'] data vs ['CPU'], again on log-log axes, for all four methods. You should find that RATTLE and multiple timesteps (with n_mts=10) give fairly similar performance to each other, as measured this way; both are better than the single-step springs program, and worse than the special MILC SHAKE algorithm. Of course, more generally, this kind of comparison is dependent on the system studied, and on some of the parameters chosen.
[ ]:
plot_all_results()
Thermostatting with constraints¶
So far, all these simulations have been in the constant-\(NVE\) ensemble. There are a few subtleties associated with applying a thermostat to a system which includes constraints. One of the simpler thermostats is due to DM Heyes, Chem. Phys., 82, 285 (1983), corrected by M Hecht et al, Phys. Rev. E, 72, 011408 (2005). It is well explained in D Frenkel, B Smit, Understanding Molecular
Simulation, 3rd ed (2023) with an erratum. The momenta are scaled up or down by a common factor at each step. The factor is chosen by a Monte Carlo procedure which guarantees the correct sampling of the total kinetic energy. This is implemented in the routine thermostat, in the file md_module.f90. In the programs, this is called (if required) at the start of
each timestep. The thermostat due to G Bussi et al, J. Chem. Phys., 126, 014101 (2007) works in a similar way.
Provided the total linear and angular momenta, and the bond constraint derivatives, are zero, these thermostats do not cause them to become nonzero, since they simply scale all momenta up and down by a common factor.
Other methods, such as that of HC Andersen, J Chem Phys, 72, 2384 (1980), in which atomic momenta are reselected at intervals from the Maxwell-Boltzmann distribution, may be used. However, some care needs to be taken if it is wished to conserve total momentum, a correction described by JP Ryckaert and G Ciccotti, Mol Phys, 58, 1125 (1986) is needed to satisfy the bond constraints, and some pitfalls exist for the unwary, as explained by EAJF Peters et al, J Chem Theor Comput, 10, 4208 (2014).
You may like to consider these points in more detail, and perhaps jot down some of your thoughts in the cell below.
Enter your thoughts here.
The working of the thermostat is checked, for the constraints case, in the following cell. Starting from the supplied \(T=1.0\) configuration, the program applies a thermostat with \(T=2.0\).
[ ]:
old_conf = data_dir / 'md_constraints.xyz'
new_conf = copy(old_conf,this_dir)
!echo '&nml nvt=.true., temperature=2.0 /' | build/md_constraints
Once this has finished, compare the average kinetic temperature with the desired value. You’ll need to amend the expression for \(N_\text{free}\) below, recalling the exercise on this that you did earlier in the worksheet.
[ ]:
attr, data = read_hdf5_file('md_constraints.hdf5')
[ ]:
print(attr['Title'].astype(str))
N = attr['N']
T = attr['T']
print(f'Number of atoms N = {N:10d}')
print(f'Set temperature T = {T:10.4f}')
Nfree = (3*N-6) # Correct this expression again
U = data['U']
K = data['K']
T = 2*K/Nfree
discard=1000
Tavg = T[discard:].mean()
print(f'Kinetic temp T = {Tavg:10.4f}')
We (somewhat arbitrarily) discarded the first 1000 steps in calculating the mean, so as to allow for equilibration from \(T=1.0\) to \(T=2.0\). Hopefully the average is not too far from the specified value.
If you wish, you can plot T against step number (time), to see the fluctuations around its mean value, or plot a histogram of values of T to quantify them. These serve to emphasize that the kinetic temperature is not fixed in this ensemble!
[ ]:
# Insert your plotting code here
Although the kinetic energy is controlled quite well, the other energies typically respond more slowly to the sudden change in temperature. The polymer chain becomes more open, at the higher temperature. You can see this if you plot \(U\) against time, and if you view a snapshot of the final configuration new_conf in the current working directory, and compare with the initial one old_conf.
[ ]:
# Insert your plotting code here
[ ]:
atoms=read(new_conf)
atoms.numbers[:] = atomic_numbers['Ar'] # Visualize as argon atoms
sigma = vdw_radii[atomic_numbers['Ar']]*2 # Argon diameter in Angstroms
atoms.positions *= sigma # Scale all positions
viewer=WeasWidget(from_ase=atoms,viewerStyle=weas_shape)
viewer.avr.bond.add_bond_pair('Ar','Ar',max=0.89*sigma)
viewer.avr.model_style=1
viewer
Skip this last section if time is short.
Here, we just illustrate the operation of the thermostat for the model with springs.
[ ]:
old_conf = data_dir / 'md_springs.xyz'
new_conf = copy(old_conf,this_dir)
!echo '&nml nvt=.true., temperature=2.0 /' | build/md_springs
Take a look at the average kinetic temperature for this run.
[ ]:
attr, data = read_hdf5_file('md_springs.hdf5')
[ ]:
print(attr['Title'].astype(str))
N = attr['N']
T = attr['T']
print(f'Number of atoms N = {N:10d}')
print(f'Set temperature T = {T:10.4f}')
K = data['K']
U = data['U']
Nfree = (3*N-6)
T = 2*K/Nfree
discard=10000
Tavg = T[discard:].mean()
print(f'Average temp T = {Tavg:10.4f}')
Again, we discard the initial phase of the data, to allow for equilibration. Hopefully the result should agree with the specified value, showing that the thermostat works.
If you wish, insert some code here (similar to the code above, for the constraint model) to plot \(T\) vs step number, and to show a histogram of the distribution of \(T\).
[ ]:
# Insert your plotting code here
Finally, insert (if you wish) some code to plot the evolution of \(U\) with time, and take a look at the final configuration.
[ ]:
# Insert your plotting code here
[ ]:
atoms=read(new_conf)
atoms.numbers[:] = atomic_numbers['Ar'] # Visualize as argon atoms
sigma = vdw_radii[atomic_numbers['Ar']]*2 # Argon diameter in Angstroms
atoms.positions *= sigma # Scale all positions
viewer=WeasWidget(from_ase=atoms,viewerStyle=weas_shape)
viewer.avr.bond.add_bond_pair('Ar','Ar',max=0.89*sigma)
viewer.avr.model_style=1
viewer
This concludes the workshop. The remainder of the notebook gives a bit more background to the MILC SHAKE algorithm, but can be read later, at your convenience.
The MILC SHAKE algorithm¶
An advantage of SHAKE/RATTLE is that it is fairly general, applying to a wide range of bond topologies. Nonetheless, some molecular structures can be tackled more efficiently using a variant of the approach: the water molecule is an obvious case. Here, we have a linear chain molecule, and so we give a bit more detail regarding the MILC SHAKE algorithm, which applies to this situation. None of this is required reading, in order to complete the workshop.
The algorithm is well described by the original authors Bailey et al., J Comput Phys, 227, 8949 (2008); see also Bailey et al., Comput Phys Commun, 180, 594 (2009). Our implementation does not follow theirs precisely, but is equivalent. Here, we aim to highlight the resemblance to RATTLE, and we use the same notation as before.
For our chain molecule, we identify each constraint by successive atom indices \(i,i+1\). We do not tackle them one at a time. Instead, in milcshake_a, following the (unconstrained) first kick and drift, we aim to compute all \(\lambda_{i,i+1}\) together. The updating scheme for each atom \(i\), bonded to \(i+1\) and \(i-1\), is \begin{align*}
\mathbf{p}_i'' &=
\mathbf{p}_i'
+ \frac{m}{\Delta t}\bigl(\lambda_{i,i+1}\,\mathbf{q}_{i,i+1}
-\lambda_{i-1,i}\,\mathbf{q}_{i-1,i}\bigr),\quad i = 1\ldots N,
\\
\mathbf{r}_i'' &= \mathbf{r}_i'
+ \bigl(\lambda_{i,i+1}\,\mathbf{q}_{i,i+1}
-\lambda_{i-1,i}\,\mathbf{q}_{i-1,i}\bigr), \quad i = 1\ldots N ,
\end{align*} except that terms corresponding to non-existent bonds (at the ends of the chain) should be omitted. The \(C=N-1\) bond constraint equations are :nbsphinx-math:`begin{equation*}
bigl|mathbf{r}_{i,i+1}’’bigr|^2 =
Bigl|mathbf{r}_{i,i+1}’
+ bigl(-lambda_{i-1,i}mathbf{q}_{i-1,i} + 2lambda_{i,i+1}mathbf{q}_{i,i+1}
lambda_{i+1,i+2}mathbf{q}_{i+1,i+2}bigr)Bigr|^2
= d^2, quad i=1 ldots C. end{equation*}` Expanding, and dropping the quadratic terms, :nbsphinx-math:`begin{equation*} -2(mathbf{r}_{i,i+1}’cdotmathbf{q}_{i-1,i})lambda_{i-1,i} +4(mathbf{r}_{i,i+1}’cdotmathbf{q}_{i,i+1})lambda_{i,i+1} -2(mathbf{r}_{i,i+1}’cdotmathbf{q}_{i+1,i+2})lambda_{i+1,i+2}
= d^2 - mathbf{r}_{i,i+1}’cdotmathbf{r}_{i,i+1}’ .
end{equation*}` This has the form of a matrix equation \(\mathbb{A}\cdot\boldsymbol{\lambda} = \boldsymbol{\sigma}\), where \(\lambda_i\equiv\lambda_{i,i+1}\), \(\sigma_i=d^2 - \mathbf{r}_{i,i+1}'\cdot\mathbf{r}_{i,i+1}'\), \begin{align*}
A_{i,i-1} &= -2(\mathbf{r}_{i,i+1}' \cdot \mathbf{q}_{i-1,i}), \quad i=2 \ldots C ,
\\
A_{i,i} &= \quad 4(\mathbf{r}_{i,i+1}' \cdot \mathbf{q}_{i,i+1}), \quad i=1 \ldots C ,
\\
A_{i,i+1} &= -2(\mathbf{r}_{i,i+1}' \cdot \mathbf{q}_{i+1,i+2}), \quad i=1 \ldots C-1 ,
\end{align*} and all other \(A_{ij}=0\). If we ignored the off-diagonal elements, we would be left with \(\lambda_{i,i+1}=\sigma_i/A_{i,i}\), as we saw for rattle_a. Because \(\mathbb{A}\) is tridiagonal, the full matrix equation is easily solved for all the \(\lambda_{i,i+1}\) multipliers. These are used to update all the coordinates and momenta at once. The updated coordinates will satisfy all the bond constraint equations except for the neglected quadratic terms.
Because of this, iteration is still needed. The new coordinates are used to re-calculate the \(\sigma_i\), the matrix equation solved again, and the process repeated, until all constraints are satisfied to within the prescribed tolerance. This is similar to rattle_a; a slight difference is that the matrix \(\mathbb{A}\) is only calculated once, at the start, whereas the quantity \(A\) in rattle_a is recomputed at each iteration.
The milcshake_b routine, following the second kick, will implement the updating scheme \begin{equation*}
\mathbf{p}_i'' = \mathbf{p}_i'
+ (\mu_{i,i+1}\mathbf{r}_{i,i+1} - \mu_{i-1,i}\mathbf{r}_{i-1,i}),
\quad i=1 \ldots N,
\end{equation*} where, again, terms corresponding to nonexistent bonds at the ends should be omitted. The constraints on time derivatives of the bond length, \(m\mathbf{v}_{i,i+1}''\cdot\mathbf{r}_{i,i+1}=0\) in this case of equal masses \(m\), become \begin{equation*}
-(\mathbf{r}_{i-1,i}\cdot\mathbf{r}_{i,i+1})\mu_{i-1,i}
+ 2(\mathbf{r}_{i,i+1}\cdot\mathbf{r}_{i,i+1})\mu_{i,i+1}
-(\mathbf{r}_{i+1,i+2}\cdot\mathbf{r}_{i,i+1})\mu_{i+1,i+2}
= -m\mathbf{v}_{i,i+1}'\cdot\mathbf{r}_{i,i+1} .
\end{equation*} This is a matrix equation \(\mathbb{B}\cdot\boldsymbol{\mu}=\boldsymbol{\tau}\), with \(\mu_i\equiv\mu_{i,i+1}\), \(\tau_i=-m\mathbf{v}_{i,i+1}'\cdot\mathbf{r}_{i,i+1}\), \begin{align*}
B_{i,i-1} &= -(\mathbf{r}_{i-1,i}\cdot\mathbf{r}_{i,i+1}), \quad i=2 \ldots C,
\\
B_{i,i} &= \quad 2(\mathbf{r}_{i,i+1}\cdot\mathbf{r}_{i,i+1}), \quad i=1 \ldots C,
\\
B_{i,i+1} &= -(\mathbf{r}_{i+1,i+2}\cdot\mathbf{r}_{i,i+1}), \quad i=1 \ldots C-1,
\end{align*} and all other \(B_{ij}=0\). If we ignored the off-diagonal elements, we would have \(\mu_{i,i+1}=\tau_i/B_{i,i}\), as per rattle_b. Solving the triadiagonal matrix equation gives all the \(\mu_{i,i+1}\), and using these to update the momenta will exactly satisfy all the velocity constraints. No iterations are necessary, since there are no quadratic terms.
This concludes the notebook.
[ ]: