Constraints¶
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 h5py data_tutorials weas_widget ase
! apt install gfortran libhdf5-dev
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/Constraints/",
filename=["hdf5_module.py","dat_to_ase.py"],
folder=".",
)
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_6/Constraints/",
filename=["Makefile","config_io_module.f90","force_module.f90","hdf5_module.f90","maths_module.f90","md_constraints.f90","md_module.f90","md_springs.f90","md_springs_mts.f90"],
folder=".",
)
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_6/Constraints/",
filename=["md_springs.dat", "md_constraints.dat"],
folder=".",
)
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 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 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. 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.
[ ]:
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import numpy as np
import shutil
from glob import glob
from hdf5_module import read_file
from dat_to_ase import argon
from weas_widget import WeasWidget
plt.style.use(['seaborn-v0_8-talk','seaborn-v0_8-darkgrid','seaborn-v0_8-colorblind'])
plt.rc('image',cmap='viridis')
plt.rc('legend',frameon=True,framealpha=1.0)
The model¶
In this exercise we shall consider a simple system: a flexible chain molecule of \(N=64\) atoms in periodic boundary conditions. The nonbonded interaction between all pairs of atoms (except for the ones that are bonded to each other) is the WCA repulsive Lennard-Jones potential, \begin{equation*} u_{\mathrm{WCA}}(r) = \begin{cases} u_{\mathrm{LJ}}(r) + \varepsilon & r \leq r_{\mathrm{min}} \\ 0 & r> r_{\mathrm{min}} \end{cases} \end{equation*} where \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. \(r_{\mathrm{min}}=2^{1/6}\sigma\) is the position of the minimum of the Lennard-Jones potential. 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\).
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, and \(\kappa\) is the spring constant, set here to a fairly large value \(\kappa=10000\) (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*} 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*} An additional simplification is that 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 cell will build all the codes to be used here: md_constraints,md_springs, and md_springs_mts.
[ ]:
!make
The programs¶
The main programs are in the files 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 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, md_springs.dat and md_constraints.dat, one for each model. The format of these files is
n
xbox ybox zbox
x1 y1 z1 px1 py1 pz1
x2 y2 z2 px2 py2 pz2
x3 y3 z3 px3 py3 pz3
: : : : : :
xn yn yn pxn pyn pzn
where the first line gives the number of atoms, the second line gives the box dimensions (in this case we use cubic boxes) and the subsequent lines give the coordinates and momenta of each atom. If you wish to view these configurations, first use the Python script dat_to_xyz.py (outside this notebook) to convert the files to XYZ format, and use one of the supplied molecular graphics programs. In both cases, the systems have been equilibrated at a temperature \(T\approx1.0\) and a density
\(\rho\approx 0.7\). They should also have zero total momentum \(\sum_i \mathbf{p}_i=\mathbf{0}\), and the programs print this at the start and the end of the run, for you to check.
[ ]:
atoms=argon('md_springs.dat')
v=WeasWidget()
v.from_ase(atoms)
v
The programs all 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 /' | ./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}=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 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.
[ ]:
params, data = read_file('md_springs.hdf5')
[ ]:
print(params['Title'].astype(str))
N = params['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 \(3\) from the total to account for the three conserved (\(=0\)) components of total momentum.
[ ]:
Kavg = data['K'].mean()
Nfree = (3*N-3)
Tavg = 2*Kavg/Nfree
print(f'Average temperature T = {Tavg:10.4f}')
Hopefully the average temperature is very close to 1.
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()
estd = e.std()
print(f'Average energy/atom e = {eavg:10.5f}')
print(f'RMS deviation e_RMS = {estd: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. Once more, we use default parameters. Note that the default timestep dt is \(10\) times as long as for md_springs.
[ ]:
!echo '&nml /' | ./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.
[ ]:
params, data = read_file('md_constraints.hdf5')
[ ]:
print(params['Title'].astype(str))
N = params['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()
estd = e.std()
print(f'Average energy/atom e = {eavg:10.4f}')
print(f'RMS deviation e_RMS = {estd: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 very close to 1.0.
[ ]:
Kavg = data['K'].mean()
Nfree = (3*N-3) # This expression is INCORRECT, please insert corrected version
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, at \(100\) reduced time units. After each run the log and HDF5 output files are renamed, to save them for further analysis. These runs will just take a few minutes: wait until they are all complete before proceeding.
[ ]:
nsteps = [20000,50000,100000,200000]
t_run = 100.0
print(' nstep dt')
for run, nstep in enumerate(nsteps):
dt = t_run/nstep
! echo "&nml nstep=$nstep, dt=$dt /" | ./md_springs >& md_springs.log
shutil.move('md_springs.hdf5','md_springs_'+str(run)+'.hdf5')
shutil.move('md_springs.log', 'md_springs_'+str(run)+'.log')
print(f'{nstep:10d}{dt:10.4f}')
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 which we will call e_RMS. 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 ( pattern ):
"""Gathers desired data from all HDF5 files matching the supplied pattern.
Argument
--------
pattern : string
Matching pattern for names of HDF5 files to read
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(glob(pattern)) # Sorted list of files matching supplied pattern
d = {'dt':[],'CPU':[],'e_RMS':[]} # Initial dictionary of empty lists
for file in files: # Loop over files
params,data=read_file(file) # Read from the HDF5 file
d['dt'].append(params['dt']) # Append timestep to corresponding list
d['CPU'].append(params['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/params['N'] # Add energy/atom contribution into total
d['e_RMS'].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 md_springs.
[ ]:
md_springs = gather_data('md_springs_?.hdf5')
The next cell plots md_springs['e_RMS'] vs md_springs['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 \(y \propto x^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(figsize=(8,5))
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.plot(md_springs['dt'],md_springs['e_RMS'],'o',label='springs')
# Draw a crudely fitted line with slope 2 as a guide
ax.set_prop_cycle(None)
c=np.average(md_springs['e_RMS']/md_springs['dt']**2)
dt=np.array(ax.get_xlim())
ax.plot(dt,c*dt**2,'--')
ax.legend()
plt.tight_layout()
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. As before, this will not take long, you should wait until all the runs are complete before proceeding. The HDF5 files are renamed md_rattle_0.hdf5 etc, to reflect the constraint algorithm being used.
[ ]:
nsteps = [20000,50000,100000,200000]
t_run = 100.0
print(' nstep dt')
for run, nstep in enumerate(nsteps):
dt = t_run/nstep
! echo "&nml nstep=$nstep, dt=$dt /" | ./md_constraints >& md_constraints.log
shutil.move('md_constraints.hdf5','md_rattle_'+str(run)+'.hdf5')
shutil.move('md_constraints.log', 'md_rattle_'+str(run)+'.log')
print(f'{nstep:10d}{dt:10.4f}')
print('All done')
When the runs are all done, collect the dt, CPU and e_RMS data from these log files, as before.
[ ]:
md_rattle = gather_data('md_rattle_?.hdf5')
[ ]:
fig, ax = plt.subplots(figsize=(8,5))
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 for both cases
ax.plot(md_springs['dt'],md_springs['e_RMS'],'o',label='springs')
ax.plot(md_rattle['dt'],md_rattle['e_RMS'],'s',label='RATTLE')
# Draw crudely fitted lines with slope 2 as a guide
ax.set_prop_cycle(None)
dt=np.array(ax.get_xlim())
c=np.average(md_springs['e_RMS']/md_springs['dt']**2)
ax.plot(dt,c*dt**2,'--')
c=np.average(md_rattle['e_RMS']/md_rattle['dt']**2)
ax.plot(dt,c*dt**2,'--')
ax.legend()
plt.tight_layout()
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.
[ ]:
fig, ax = plt.subplots(figsize=(8,5))
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,4,6)))
ax.xaxis.set_minor_formatter(ticker.ScalarFormatter())
# Plot energy conservation vs CPU time for both cases
ax.plot(md_springs['CPU'],md_springs['e_RMS'],'o',label='springs')
ax.plot(md_rattle['CPU'],md_rattle['e_RMS'],'s',label='RATTLE')
# Draw crudely fitted lines with slope -2 as a guide
ax.set_prop_cycle(None)
cpu=np.array(ax.get_xlim())
c=np.average(md_springs['e_RMS']*md_springs['CPU']**2)
ax.plot(cpu,c/cpu**2,'--')
c=np.average(md_rattle['e_RMS']*md_rattle['CPU']**2)
ax.plot(cpu,c/cpu**2,'--')
ax.legend()
plt.tight_layout()
The RATTLE approach is still better, but the difference is not so dramatic. This is because the 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 Bailey et al. J Comput Phys, 227, 8949 (2008) and 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.
[ ]:
!echo '&nml milcshake=.true. /' | ./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. The HDF5 files are renamed md_milcshake_0.hdf5 etc, to reflect the constraint algorithm.
[ ]:
nsteps = [20000,50000,100000,200000]
t_run = 100.0
print(' nstep dt')
for run, nstep in enumerate(nsteps):
dt = t_run/nstep
! echo "&nml nstep=$nstep, dt=$dt, milcshake=.true. /" | ./md_constraints >& md_constraints.log
shutil.move('md_constraints.hdf5','md_milcshake_'+str(run)+'.hdf5')
shutil.move('md_constraints.log', 'md_milcshake_'+str(run)+'.log')
print(f'{nstep:10d}{dt:10.4f}')
print('All done')
[ ]:
md_milcshake = gather_data('md_milcshake_?.hdf5')
[ ]:
fig, ax = plt.subplots(figsize=(8,5))
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,4,6)))
ax.xaxis.set_minor_formatter(ticker.ScalarFormatter())
# Plot energy conservation vs CPU time for all 3 cases
ax.plot(md_springs['CPU'],md_springs['e_RMS'],'o',label='springs')
ax.plot(md_rattle['CPU'],md_rattle['e_RMS'],'s',label='RATTLE')
ax.plot(md_milcshake['CPU'],md_milcshake['e_RMS'],'^',label='MILC SHAKE')
# Draw crudely fitted lines with slope -2 as a guide
ax.set_prop_cycle(None)
cpu=np.array(ax.get_xlim())
c=np.average(md_springs['e_RMS']*md_springs['CPU']**2)
ax.plot(cpu,c/cpu**2,'--')
c=np.average(md_rattle['e_RMS']*md_rattle['CPU']**2)
ax.plot(cpu,c/cpu**2,'--')
c=np.average(md_milcshake['e_RMS']*md_milcshake['CPU']**2)
ax.plot(cpu,c/cpu**2,'--')
ax.legend()
plt.tight_layout()
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 you have time, make a comparison with the multiple time step (MTS) approach to the model with spring bonds. The main point of difference between md_springs.f90 and md_springs_mts.f90 lies in the advancement of coordinates and momenta within the main loop. 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=1
you should be able to see that this is just the standard velocity Verlet algorithm, as in the md_springs.f90 program. However, in this section, you are going to set n_mts=10.
The following cell carries out several runs, as before, renaming the HDF5 output file md_springs_mts_0.hdf5 etc. As always, wait until all the runs have finished before proceeding.
[ ]:
nsteps = [20000,50000,100000,200000]
t_run = 100.0
n_mts = 10
print(' nstep n_mts*dt n_mts*nstep dt')
for run, nstep in enumerate(nsteps):
dt = t_run/(nstep*n_mts)
! echo "&nml n_mts=$n_mts, nstep=$nstep, dt=$dt /" | ./md_springs_mts >& md_springs_mts.log
shutil.move('md_springs_mts.hdf5','md_springs_mts_'+str(run)+'.hdf5')
shutil.move('md_springs_mts.log', 'md_springs_mts_'+str(run)+'.log')
print(f'{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. We call the dictionary md_springs_mts after the method just used.
[ ]:
md_springs_mts = gather_data ('md_springs_mts_?.hdf5')
Try plotting the ['e_RMS'] data vs ['CPU'], again on log-log axes, for all four methods. You should find that RATTLE and MTS (with n_mts=10) give fairly similar performance to each other, as measured this way; both are better than the non-MTS 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.
[ ]:
fig, ax = plt.subplots(figsize=(8,5))
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,4,6)))
ax.xaxis.set_minor_formatter(ticker.ScalarFormatter())
# Plot energy conservation vs CPU time for all 4 cases
ax.plot(md_springs['CPU'],md_springs['e_RMS'],'o',label='springs')
ax.plot(md_rattle['CPU'],md_rattle['e_RMS'],'s',label='RATTLE')
ax.plot(md_milcshake['CPU'],md_milcshake['e_RMS'],'^',label='MILC SHAKE')
ax.plot(md_springs_mts['CPU'],md_springs_mts['e_RMS'],'D',label='springs MTS')
# Draw crudely fitted lines with slope -2 as a guide
ax.set_prop_cycle(None)
cpu=np.array(ax.get_xlim())
c=np.average(md_springs['e_RMS']*md_springs['CPU']**2)
ax.plot(cpu,c/cpu**2,'--')
c=np.average(md_rattle['e_RMS']*md_rattle['CPU']**2)
ax.plot(cpu,c/cpu**2,'--')
c=np.average(md_milcshake['e_RMS']*md_milcshake['CPU']**2)
ax.plot(cpu,c/cpu**2,'--')
c=np.average(md_springs_mts['e_RMS']*md_springs_mts['CPU']**2)
ax.plot(cpu,c/cpu**2,'--')
ax.legend()
plt.tight_layout()
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 Andersen, J Chem Phys, 72, 2384 (1980) in which atomic momenta are reselected at intervals from the Maxwell-Boltzmann distribution. The routine for this, andersen, is given in md_module.f90. In the md_springs.f90 program,
this is called (if required) at the start of each timestep. (Note, in many applications it would be more appropriate to call the routine much less frequently, and perhaps apply it to individual atom momenta rather than all at once).
The working of this thermostat is checked in the following cell.
[ ]:
!echo '&nml nvt=.true., temperature=1.5 /' | ./md_springs
Take a look at the average kinetic temperature for this run.
[ ]:
params, data = read_file('md_springs.hdf5')
[ ]:
print(params['Title'].astype(str))
N = params['N']
T = params['T']
print(f'Number of atoms N = {N:10d}')
print(f'Set temperature T = {T:10.4f}')
K = data['K']
Nfree = (3*N-3)
T = 2*K/Nfree
Tavg = T[2000:].mean()
print(f'Average temp T = {Tavg:10.4f}')
We (somewhat arbitrarily) discarded the first 2000 steps in calculating the mean, so as to allow for equilibration from \(T=1.0\) to \(T=1.5\).
Hopefully, you know that the kinetic energy, and hence the instantaneous kinetic temperature, are not fixed in the canonical ensemble: they fluctuate about the mean value. If you wish, you can plot T against step number (time), to see these fluctuations, or plot a histogram of values of T to quantify them.
[ ]:
# Do the plots here, if you wish
Now take a look at the md_constraints.f90 file. Just inside the loop over steps, you will see a comment about this thermostatting. Clearly, just calling the andersen routine is wrong: a random, isotropic, selection of atomic momenta would violate the time derivatives of the constraints. Fortunately, it was shown by Ryckaert and Ciccotti, Mol Phys, 58, 1125 (1986) that a correct set of momenta will be obtained if the constraint
conditions on the relative velocities of bonded atoms are imposed immediately after this selection. Remember that the adjustment of momenta is performed by the constraints_b routine (either rattle_b or milcshake_b). So, insert an extra call to constraints_b, so that section of code looks like
if ( nvt ) then
call andersen ( temperature, p )
call constraints_b ( box, r, p, iter )
end if
Recompile the program, and run, in the next cells.
[ ]:
!make
[ ]:
!echo '&nml nvt=.true., temperature=1.5 /' | ./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.
[ ]:
params, data = read_file('md_constraints.hdf5')
[ ]:
print(params['Title'].astype(str))
N = params['N']
T = params['T']
print(f'Number of atoms N = {N:10d}')
print(f'Set temperature T = {T:10.4f}')
Nfree = (3*N-3) # Expression INCORRECT, please insert corrected version as you did before
K = data['K']
T = 2*K/Nfree
Tavg = T[200:].mean()
print(f'Kinetic 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.
Optional exercise¶
There is a danger lurking for the unwary, in thermostatting the constrained system this way. If you have time, you can look into this.
Instead of reselecting all the atomic momenta in the Andersen method, quite often the reselection is done for just a fraction of them, chosen at random. Inside the andersen routine in the md_module.f90 file is a parameter fraction which allows this to be done. Change this from its current value of \(1.0\) to a lower value, say \(0.4\) or \(0.2\), and then recompile both the programs.
Re-run the md_springs program, with nvt=.true., temperature=1.5; you should see that all is well, i.e. the kinetic temperature still agrees with the thermostat temperature.
However, if you re-run md_constraints with nvt=.true., temperature=1.5, the run-average temperature should be noticeably different from the specified temperature. This may seem a surprising outcome, for such an apparently innocuous change!
This particular issue, the reasons for it, and a simple cure for the problem, have been described by Peters et al, J Chem Theor Comput, 10, 4208 (2014). For now, we will just leave it as a warning to take care when combining algorithms: things are not always as simple as they appear. If you intend to use the md_constraints.f90 program again, you should re-edit the md_module.f90 file, making sure that fraction=1.0 in the andersen routine.
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) and 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.
[ ]: