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

Day 8 Thermostats¶

These exercises use DL_POLY to demonstrate the use of temperature scaling and extended system methods. The working directory for these exercises is ~/WORKSHOP/Day_8/MD_Ensembles. The DL_POLY v5 manual can be found at this link for reference whilst the executable, DLPOLY.Z.

You shall get same answers using any molecular dynamics code, just be sure you actually read the manual since subtle differences may exist.

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 weas_widget ase
! 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_8/MD_Ensembles/Exercise_1/",
    filename=["CONTROL","FIELD","CONFIG"],
    folder="Exercise_1/",
)
get_data(
    url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_8/MD_Ensembles/Exercise_1/high_T/",
    filename=["CONTROL","FIELD","CONFIG"],
    folder="Exercise_1/high_T/",
)
get_data(
    url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_8/MD_Ensembles/Exercise_1/no_scale/",
    filename=["CONTROL","FIELD","CONFIG"],
    folder="Exercise_1/no_scale/",
)
get_data(
    url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_8/MD_Ensembles/Exercise_2/",
    filename=["CONTROL","FIELD","CONFIG"],
    folder="Exercise_2/",
)
get_data(
    url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_8/MD_Ensembles/Exercise_3/",
    filename=["CONTROL","FIELD","CONFIG"],
    folder="Exercise_3/",
)

get_data(
    url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_8/MD_Ensembles/",
    filename=["rdf.py","statis.py"],
    folder=".",
)

get_data(
    url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/.raw/",
    filename=["DLPOLY.Z"],
    folder="/usr/local/bin",
)
! chmod 755 /usr/local/bin/DLPOLY.Z

[ ]:

from statis import *
from rdf import *
import numpy as np
from weas_widget import WeasWidget
from ase.io import read

Exercise 1: Velocity scaling¶

This exercise demonstrates the use of velocity scaling to set the initial temperature of an MD run. The test case is comprised of rigid SPC water molecules. The FIELD, CONTROL and CONFIG files are in the directory Exercise_1 and there is a copy of all the input files with the scale command removed(rescale_frequency) in the sub-directory no_scale. Study the files and try to understand what is being requested of the simulation. The DL_POLY keyword rescale_frequency is used to define the number of MD steps between velocity scaling events. In this example we have rescale_frequency 10 steps, so that every 10 time steps during the equilibration phase of the simulation run the velocities will be scaled to bring the temperature into line with the set temperature (295 K).

Run DLPOLY.Z for the system with and without scaling in the respective directories.

start by visualising your data

[ ]:

water = read("Exercise_1/CONFIG")
water.center()

v=WeasWidget()
v.from_ase(water)
v.avr.model_style = 1
v.avr.show_hydrogen_bonds = True
v

[ ]:

%%bash

cd Exercise_1
rm -f STATIS
DLPOLY.Z

you can view now the trajectory created

[ ]:

traj = read("Exercise_1/HISTORY",index=":")
_ = [f.center() for f in traj]

v=WeasWidget()
v.from_ase(traj)
v.avr.model_style = 1
v.avr.show_hydrogen_bonds = True
v

[ ]:

%%bash

cd Exercise_1/no_scale
rm -f STATIS
DLPOLY.Z

We will print the instantaneous temperature from STATIS files. We use a simple python script to extract data and plot… if you are curious you can look into it and extend it to plot more quantities.

[ ]:

s1 = Statis("Exercise_1/STATIS")
s2 = Statis("Exercise_1/no_scale/STATIS")

s1.plot("temperature")
s2.plot("temperature")

rescale_frequencyFor the RDF plot use the another small python script (see below) This will print the g(r) of all species in your system, O and H in this case.For this example the oxygen atoms of the water are labelled OW.

The start point for these simulations was generated using a high temperature run (see Exercise_1/high_T). How has this affected the runs with and without velocity scaling?

[ ]:

%%bash

cd Exercise_1/high_T
rm -f STATIS
DLPOLY.Z

[ ]:

r1 = rdf("Exercise_1/RDFDAT")
r2 = rdf("Exercise_1/no_scale/RDFDAT")

r1.plot()
r2.plot()

Exercise 2: Thermostats¶

This exercise looks at the choice of coupling time constant for a simulation using a thermostat to obtain the NVT ensemble using Hoover thermostat as modified by Melchiona et al. The system is the same as in Exercise 1, an SPC model of water which has been initially run at 500 K. We will test coupling time constants of 2.0 ps, 0.5 ps and 0.02 ps. Go into the Exercise_2 directory and make some appropriately named sub-directories, e.g.: ‘’mkdir 0p002’’ Copy the CONFIG, CONTROL and FIELD files into your sub-directories, edit the CONTROL file to contain the right value (look for ensemble_thermostat_coupling line) then save and run the DLPOLY.Z command.

[ ]:

%%bash

cd Exercise_2
mkdir -p 0p002
cd 0p002
cp ../CONFIG .
cp ../FIELD .
cp ../CONTROL .
# do not forget to change the CONTROL file

[ ]:

%%bash

cd Exercise_2/0p002
DLPOLY.Z

Add cells as needed to do all the simulations

Remember the starting structure in this case was at high temperature so the first thing to consider is the way that the thermostat brings the temperature into line with the target temperature and once this is achieved the fluctuations observed.

Using the same scrtipt as above plot the temperature vs time plots for each of your simulations as we did in Exercise 1. What is the effect of the coupling time constant on the rate at which the target temperature is achieved and the fluctuations in the temperature in general?

[ ]:

s1 = Statis("Exercise_2/0p002/STATIS")

s1.plot("temperature")

Fluctuations contain information. Remember that the aim is not to eliminate fluctuations altogether but to obtain a simulation that correctly samples the NVT ensemble. In this part of the exercise we will look at the fluctuations in the total energy for the NVT simulation using the different coupling constants. We will use the restart file in DL_POLY to obtain a structure for each of the coupling constants. For example, in the 0p002 directory type.

[ ]:

%%bash

cd Exercise_2
mkdir -p equilib_run
cp CONTROL FIELD equilib_run
cp 0p002/REVCON equilib_run/CONFIG

cd equilib_run
DLPOLY.Z

This makes a copy of the files but with the final structure of the first run (REVCON) as a new start point. Now run the three equilib_run jobs and use the statis.py to obtain the time evolution of the total energy. The fluctuation in the total energy is given by:

\[\langle \delta E^{2} \rangle = \langle E^{2} \rangle - \langle E \rangle ^{2}\]

and this quantity is related to the heat capacity via:

\[C_{V} = \frac{N_{m}}{kT^{2}} \langle \delta E^{2} \rangle _{NVT}\]

where \(N_{m}\) is the number of particles in the simulation.

Use the python script to extract the Etot timeseries so that the required averages can be calculated. The result is a numpy array so things shall be straight forward. If you are not comfortable with python, the timeseries is saved in etot.dat. Use that with your favourite tool.

How does the time constant affect the observed fluctuations? Is this consistent with what you saw in the first part of the exercise?

[ ]:

s = Statis("Exercise_2/equilib_run/STATIS")

s.plot("total_energy")
Etot = s.get_total_energy()
print(Etot)
np.savetxt('etot.dat', Etot, delimiter='\n')

If you have time try to estimate the specific heat capacity of SPC water using your results, it is difficult to get accurate values with short runs (you could create better data once you return to your lab), have a look at the variation of your values when using different sections of the data.

Caveat¶

Finally, here is another reminder that some “thermostats”, even some widely-used ones, do not generate the canonical \(NVT\) ensemble. They may not generate any equilibrium ensemble at all! Even if they do, in principle, correspond to a well-defined ensemble, they may, in practice, produce unphysical velocity distributions (sometimes called the “flying ice cube effect”). This can be dangerous: there are published papers containing incorrect results which have been traced back to a poor choice of thermostat. You may like to read Braun et al, J Chem Theor Comput, 14, 5262 (2018) for some guidance in this area, and possibly also Petersen and Searles, Phys Chem Chem Phys, 24, 6383 (2022) for an interesting theoretical result. In any case, be aware of the issue.

[ ]: