Canonical (NVT) ensemble

TL;DR

  • NVT-MD simulations assume constant temperature and volume. This ensemble should be chosen when one wants to ignore volume and temperature changes.

  • Applications include ion diffusion in solids, adsorption and reaction on slab-structured surfaces, and adsorption and reaction on clusters.

  • There are several variations of heat baths in NVT-MD simulation methods to control temperature, such as the Berendsen, Langevin, and Nosé-Hoover heat baths. They are also implemented in ASE.

  • The Berendsen heat bath method is simple and has good convergence, but it changes the velocity distribution of atoms in the entire system uniformly, which can result in unnatural phenomena.

  • The Langevin heat bath method controls temperature by applying random forces to individual atoms, which are appropriately sampled using a statistical method. Since each atom is controlled individually, even mixed phases can be handled properly. It is not suitable for precise trajectory analysis because it does not reproduce the same trajectories even after repeated calculations due to its statistical nature.

  • The Nosé-Hoover heat bath method is one of the most universally used temperature control methods. In principle, it reproduces the correct NVT ensemble, but there are exceptions for special systems.

In the previous section, we learned about MD simulation in the NVE ensemble, the simplest form of MD simulation.

We will now discuss MD simulations for one of the most commonly used state distributions, canonical ensemble (NVT ensemble). The canonical ensemble is a statistical mechanical state distribution where the number of particles (N), volume (V), and temperature (T) are constant. It is applied when the volume change is negligible for the target system. A good example would be a simulation in a relatively low temperature region where the volume expansion is negligible.

Let’s take a look at some typical methods through some specific examples.

Calculation example 1: Melting of fcc-Al

As a first example, let’s look at the simulation for the melting state of metallic aluminum treated in section 6-1.

In section 6-1, the simulation was performed at a high initial temperature of 1600 K to ensure that the metallic aluminum appears sufficiently molten. However, this time we will run the simulation with an initial temperature of 1000 K, and the NVE ensemble is the same as in the previous section, only the temperature is set to 1000 K. The results are shown below.

fcc-Al_NVE_1000Kstart

Fig6-2a. Fcc-Al in NVE ensemble starting at 1000 K

(File: ../input/ch6/6-2_fcc-Al_NVE_1000Kstart.traj)

Although the lattice is indeed much perturbed, the original FCC crystal structure is clearly visible. The melting point of real Al metal is 933 K, but an initial temperature of 1000 K seems insufficient to melt the system.

We will then run the NVT ensemble at the set temperature of 1000 K. The calculation method is as follows. The calculation time is almost the same as the first NVE ensemble calculation.

[1]:
import os
from asap3 import EMT
calculator = EMT()

import ase
from ase.build import bulk
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution,Stationary
from ase.md.nvtberendsen import NVTBerendsen
from ase.md import MDLogger
from ase import units
from time import perf_counter

# Set up a crystal
atoms = bulk("Al","fcc",a=4.3,cubic=True)
atoms.pbc = True
atoms *= 3
print("atoms = ",atoms)

# input parameters
time_step    = 1.0    # fsec
temperature  = 1000    # Kelvin
num_md_steps = 1000000
num_interval = 10000
taut         = 1.0    # fs

atoms.calc = calculator

print(f"taut = {taut:.3f}")

output_filename = "./output/ch6/fcc-Al_NVT-Berendsen_1000K"
print("output_filename = ",output_filename)

log_filename = output_filename + ".log"
print("log_filename = ",log_filename)
traj_filename = output_filename + ".traj"
print("traj_filename = ",traj_filename)


# Set the momenta corresponding to the given "temperature"
MaxwellBoltzmannDistribution(atoms, temperature_K=temperature,force_temp=True)
Stationary(atoms)  # Set zero total momentum to avoid drifting

# run MD
dyn = NVTBerendsen(atoms, time_step*units.fs, temperature_K = temperature, taut=taut*units.fs, loginterval=num_interval, trajectory=traj_filename)

# Print statements
def print_dyn():
    imd = dyn.get_number_of_steps()
    etot  = atoms.get_total_energy()
    temp_K = atoms.get_temperature()
    stress = atoms.get_stress(include_ideal_gas=True)/units.GPa
    stress_ave = (stress[0]+stress[1]+stress[2])/3.0
    elapsed_time = perf_counter() - start_time
    print(f"  {imd: >3}   {etot:.3f}    {temp_K:.2f}    {stress_ave:.2f}  {stress[0]:.2f}  {stress[1]:.2f}  {stress[2]:.2f}  {stress[3]:.2f}  {stress[4]:.2f}  {stress[5]:.2f}    {elapsed_time:.3f}")


dyn.attach(print_dyn, interval=num_interval)
dyn.attach(MDLogger(dyn, atoms, output_filename+".log", header=True, stress=True, peratom=True, mode="a"), interval=num_interval)

# Now run the dynamics
start_time = perf_counter()
print(f"    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)")
dyn.run(num_md_steps)
atoms =  Atoms(symbols='Al108', pbc=True, cell=[12.899999999999999, 12.899999999999999, 12.899999999999999])
taut = 1.000
output_filename =  ./output/ch6/fcc-Al_NVT-Berendsen_1000K
log_filename =  ./output/ch6/fcc-Al_NVT-Berendsen_1000K.log
traj_filename =  ./output/ch6/fcc-Al_NVT-Berendsen_1000K.traj
    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)
    0   23.764    1000.00    7.17  7.14  7.08  7.30  0.08  0.12  0.07    0.002
  10000   42.640    1001.31    0.83  0.80  0.86  0.81  -0.12  -0.59  -0.15    2.098
  20000   43.716    996.81    0.48  0.38  0.44  0.62  0.30  -0.02  -0.00    4.170
  30000   43.096    994.42    0.83  0.89  0.68  0.90  0.06  0.06  -0.07    6.610
  40000   42.908    995.85    0.72  0.51  0.72  0.93  -0.31  -0.12  -0.09    8.731
  50000   43.791    992.83    0.46  1.08  0.03  0.25  -0.24  0.00  -0.73    10.865
  60000   42.992    1002.44    0.56  1.09  0.54  0.06  -0.17  -0.04  0.32    12.912
  70000   42.523    1004.48    0.73  0.86  1.32  0.02  -0.02  -0.02  -0.52    15.045
  80000   43.748    1004.11    0.46  0.68  -0.02  0.72  -0.02  -0.14  -0.41    17.109
  90000   44.248    996.46    0.35  0.03  0.87  0.15  0.17  0.11  -0.12    19.143
  100000   42.705    1005.16    0.77  0.55  1.04  0.70  0.20  -0.27  -0.03    21.176
  110000   41.748    1002.37    1.13  0.70  0.95  1.74  0.16  -0.26  0.55    23.247
  120000   43.692    994.98    0.34  -0.15  0.63  0.52  0.06  -0.17  -0.39    25.306
  130000   42.717    997.73    0.74  0.78  0.61  0.84  -0.59  0.25  -0.19    27.341
  140000   43.516    1005.41    0.55  1.00  0.51  0.13  -0.03  0.26  0.01    29.368
  150000   42.965    1001.53    0.61  0.67  0.48  0.67  0.14  -0.11  0.09    31.393
  160000   42.345    1003.07    0.57  0.10  0.65  0.95  -0.32  0.20  -0.27    33.422
  170000   43.904    999.60    0.62  0.60  0.44  0.80  -0.01  -0.14  -0.16    35.493
  180000   41.930    999.55    1.13  1.03  1.19  1.18  -0.18  0.04  0.38    37.554
  190000   41.640    994.53    1.01  0.46  1.26  1.31  -0.06  0.18  -0.21    39.647
  200000   43.522    999.87    0.56  0.14  0.56  0.99  -0.05  0.26  -0.13    41.684
  210000   44.519    994.04    0.20  0.72  0.32  -0.44  0.25  -0.40  -0.40    43.720
  220000   43.242    1004.29    0.73  1.21  0.33  0.66  -0.29  0.40  -0.03    45.782
  230000   42.976    994.38    0.65  0.87  0.15  0.92  -0.06  0.43  -0.28    47.818
  240000   41.835    997.14    0.91  0.94  0.77  1.02  -0.01  0.02  -0.22    49.852
  250000   43.087    1001.53    0.61  0.23  0.62  0.98  0.30  0.01  0.09    51.887
  260000   42.897    997.49    0.80  0.47  0.75  1.19  -0.19  0.01  -0.05    53.992
  270000   44.549    1001.97    0.25  0.06  0.52  0.16  -0.01  -0.26  0.12    56.176
  280000   41.939    999.00    0.86  0.17  1.03  1.39  -0.09  0.14  -0.41    58.349
  290000   41.940    1004.48    0.99  0.62  1.13  1.21  -0.03  -0.48  0.28    60.440
  300000   42.026    991.52    0.87  0.70  0.71  1.18  -0.41  0.04  -0.28    62.510
  310000   42.958    994.24    0.62  1.07  0.08  0.71  -0.13  0.41  0.12    64.563
  320000   44.344    998.56    0.29  -0.40  0.51  0.75  -0.20  -0.18  0.03    66.705
  330000   43.142    1007.59    0.60  0.46  0.82  0.52  -0.06  1.10  -0.17    68.739
  340000   43.896    995.24    0.40  0.21  0.12  0.86  -0.31  -0.54  0.04    70.887
  350000   43.730    1006.50    0.51  -0.33  0.97  0.89  -0.01  0.30  -0.47    72.917
  360000   42.838    1004.45    0.87  0.98  0.20  1.41  -0.70  0.07  -0.43    74.964
  370000   43.373    999.10    0.67  0.48  0.75  0.78  -0.04  -0.34  0.07    77.000
  380000   43.482    994.15    0.56  0.10  0.38  1.21  0.26  -0.21  0.10    79.037
  390000   43.064    999.20    0.62  0.14  0.37  1.33  -0.35  0.42  0.15    81.074
  400000   42.217    999.70    0.75  0.40  1.04  0.82  -0.00  -0.00  -0.08    83.121
  410000   43.197    1005.39    0.71  0.34  1.35  0.43  0.56  0.25  -0.12    85.211
  420000   43.210    1006.11    0.45  0.94  -0.47  0.89  -0.05  -0.42  -0.09    87.576
  430000   42.243    1003.94    0.85  0.62  0.98  0.96  -0.70  -0.18  -0.10    89.769
  440000   44.286    995.46    0.13  -0.11  0.42  0.07  -0.17  -0.45  -0.49    91.840
  450000   43.378    996.32    0.58  0.13  1.31  0.30  0.25  -0.05  0.31    93.889
  460000   41.765    999.12    1.02  1.41  1.01  0.65  0.32  -0.19  -0.03    96.016
  470000   43.689    999.34    0.28  0.24  0.46  0.14  0.15  -0.32  -0.18    98.073
  480000   42.378    994.10    0.77  0.65  0.76  0.91  -0.44  0.54  0.64    100.110
  490000   42.174    1003.58    0.79  0.84  0.56  0.96  -0.10  0.10  -0.24    102.156
  500000   43.513    1001.70    0.44  0.19  0.53  0.60  0.49  0.02  -0.01    104.197
  510000   42.466    997.10    0.72  0.94  0.84  0.39  -0.23  -0.28  -0.09    106.284
  520000   45.174    996.93    0.32  0.31  0.78  -0.13  -0.68  0.13  -0.15    108.340
  530000   41.321    1003.59    1.06  1.20  0.97  1.00  -0.12  0.05  0.33    110.376
  540000   43.508    1002.91    0.51  0.49  0.77  0.26  -0.00  -0.09  -0.09    112.405
  550000   42.929    1002.61    0.74  1.34  -0.29  1.16  0.08  0.52  0.31    114.470
  560000   43.748    994.74    0.16  0.05  0.94  -0.50  -0.11  -0.02  0.03    116.571
  570000   42.713    999.79    0.87  0.43  1.28  0.89  -0.10  -0.23  0.33    118.608
  580000   41.789    999.73    0.95  1.02  0.64  1.21  -0.01  0.05  0.35    120.640
  590000   43.109    995.22    0.53  0.80  0.49  0.30  0.05  -0.07  0.15    122.669
  600000   43.090    997.19    0.54  0.88  0.83  -0.09  0.53  0.35  0.17    124.702
  610000   43.561    998.84    0.47  -0.04  0.83  0.60  -0.29  -0.16  0.08    126.847
  620000   43.111    994.91    0.61  0.66  0.33  0.83  -0.43  0.01  0.31    128.889
  630000   43.008    996.17    0.72  0.48  1.05  0.63  0.43  0.35  0.50    130.946
  640000   41.568    1000.35    0.99  0.52  1.36  1.09  -0.20  0.05  -0.35    132.975
  650000   44.186    995.09    0.24  -0.01  0.36  0.36  0.13  -0.02  -0.28    135.008
  660000   43.030    998.40    0.50  1.01  0.14  0.36  -0.18  -0.27  -0.59    137.039
  670000   44.271    999.31    0.35  0.75  0.02  0.29  -0.21  -0.07  0.47    139.066
  680000   43.616    1003.14    0.53  0.13  0.86  0.60  -0.24  -0.13  -0.35    141.092
  690000   42.047    994.11    0.86  0.62  0.97  0.99  0.34  -0.23  -0.20    143.125
  700000   43.206    997.86    0.65  0.42  1.22  0.32  -0.17  -0.03  0.11    145.171
  710000   42.089    999.57    1.01  0.82  1.36  0.85  -0.80  0.30  0.31    147.213
  720000   41.764    996.35    0.99  0.39  0.84  1.75  0.10  0.38  -0.59    149.238
  730000   42.024    1004.12    1.01  0.60  1.74  0.68  0.19  -0.10  0.28    151.268
  740000   42.816    1004.54    0.79  0.52  0.88  0.98  0.09  0.21  -0.83    153.799
  750000   45.395    999.93    0.07  -0.29  0.32  0.17  0.00  0.04  0.39    156.148
  760000   42.748    1003.24    0.69  1.14  0.73  0.20  -0.44  0.56  0.14    158.340
  770000   43.949    1003.85    0.59  0.85  0.30  0.61  -0.49  0.06  -0.11    160.415
  780000   42.946    1004.73    0.91  1.21  0.73  0.77  -0.26  -0.42  0.25    162.832
  790000   42.096    996.55    0.79  0.88  0.93  0.55  -0.25  0.04  0.40    164.925
  800000   41.883    999.76    0.94  1.20  1.12  0.51  -0.14  0.31  -0.06    167.022
  810000   44.263    1007.70    0.30  0.21  0.40  0.29  0.38  0.03  0.16    169.131
  820000   43.193    1000.32    0.59  0.18  0.67  0.91  0.25  -0.13  0.13    171.386
  830000   43.690    999.11    0.57  0.41  0.95  0.36  0.15  0.55  -0.42    173.632
  840000   41.244    996.13    1.15  1.04  1.01  1.41  0.10  0.13  -0.06    175.916
  850000   43.915    1002.81    0.32  0.58  0.65  -0.26  -0.44  -0.03  -0.19    178.193
  860000   42.515    1000.75    0.82  1.12  0.34  1.01  0.12  -0.68  0.08    180.435
  870000   43.699    1003.19    0.38  -0.18  0.60  0.71  0.29  0.12  -0.43    182.817
  880000   40.910    996.99    1.09  0.68  1.45  1.14  0.00  0.27  0.20    185.049
  890000   41.492    1004.17    1.09  1.14  0.84  1.28  0.17  0.02  -0.13    187.212
  900000   43.619    1001.87    0.48  0.90  -0.19  0.73  0.26  -0.02  0.38    189.255
  910000   41.999    1001.88    0.83  0.81  0.82  0.86  -0.20  0.37  0.09    191.432
  920000   43.409    999.04    0.37  0.46  0.32  0.32  -0.26  -0.12  0.26    193.663
  930000   41.530    1001.77    1.13  0.79  1.15  1.46  -0.19  -0.40  -0.17    195.861
  940000   43.632    998.60    0.55  -0.38  1.39  0.63  -0.43  -0.25  -0.29    197.989
  950000   40.596    995.63    1.25  0.76  1.76  1.23  -0.24  0.25  0.13    200.084
  960000   42.383    999.50    0.88  0.84  0.25  1.54  -0.10  0.35  -0.01    202.278
  970000   41.641    1003.08    0.98  1.15  1.29  0.51  0.12  0.19  -0.04    204.411
  980000   40.679    997.51    1.24  1.68  1.03  1.00  -0.49  -0.13  0.06    206.572
  990000   43.507    1002.15    0.46  0.78  0.06  0.53  -0.30  -0.35  0.29    208.710
  1000000   43.088    1001.48    0.67  1.09  0.52  0.39  -0.43  -0.12  0.42    210.838
[1]:
True

In this case, we are using a method called the Berendsen heat bath method. In this method, the temperature is adjusted by scaling the velocity of the atoms so that they approach the target temperature [1-3].

Specifically, as explained in section 6-1, the relationship between temperature and velocity is given by

\[K = \sum_{i=1}^{N} \frac{1}{2} m_i {\mathbf{v}}_i^2 = \frac{3}{2} k_B T\]

To control the temperature, the velocity of the atoms in the entire system is scaled and controlled.

The coefficient that scales the velocity is scaled is given by the following equation

\[\lambda = \sqrt{1+\frac{\Delta t}{\tau_T}\left(\frac{T_0}{T(t)}-1\right)}\]

where \(\Delta t\) is the time step and \(\tau_T\) is the time constant of the heat bath. This method can reach the target temperature more gradually by controlling the time constant \(\tau\). When \(\tau_T=\Delta t\), it is identical to the simplest temperature control method, the so-called rate scaling method. (i.e. \(\lambda = \sqrt{T_0/T(t)}\)) In this case, the temperature is always a given temperature and the temperature constraint is considered to be maximized. Conversely, as \(\tau\) increases, the temperature control becomes weaker and it takes longer to reach the target temperature. Considering the actual state of matter, a certain amount of relaxation time should be necessary, and for practical use, it is common to set the relaxation time to be sufficiently short compared to the MD calculation time.

In this method, the temperature of the system is corrected to decay exponentially as follows.

\[\frac{dT}{dt} = \frac{T_0-T}{\tau_T}\]

Now let’s look at the results of the NVT-MD simulation above.

[2]:
from ase.io import Trajectory
from pfcc_extras.visualize.povray import traj_to_apng
from IPython.display import Image


traj = Trajectory(traj_filename)
traj_to_apng(traj, f"output/ch6/Fig6-2_fcc-Al_NVT-Berendsen_1000K_unwrap.png", rotation="0x,0y,0z", clean=True, n_jobs=16)

Image(url="output/ch6/Fig6-2_fcc-Al_NVT-Berendsen_1000K_unwrap.png")
[Parallel(n_jobs=16)]: Using backend ThreadingBackend with 16 concurrent workers.
[Parallel(n_jobs=16)]: Done  18 tasks      | elapsed:    4.5s
[Parallel(n_jobs=16)]: Done 101 out of 101 | elapsed:   10.2s finished
[2]:
[3]:
# Please remove comment out if you want to show trajectory in nglviewer.

from pfcc_extras.visualize.view import view_ngl
# view_ngl(traj)

fcc-Al_NVE_1000Kstart

Fig6-2b. Fcc-Al in NVT ensemble using Berendsen thermostat at 1000 K

(File: ../input/ch6/6-2_fcc-Al_NVT-Berendsen_1000K.traj”)

In contrast to the previous NVE calculation, the structure is very dynamic and its behavior is complex. The structure is in a randomly arranged melt state. This is because the temperature is always controlled and kept around 1000 K in NVT ensemble calculations, which is always above the melting point. (Note: The actual ASE results show that the atoms diffuse violently and spread outward from the cell. For display purposes, the above figure uses the wrap method to draw the atoms back inside the cell.)

The following plots show the time evolution of the total energy (Tot.E), potential energy (PE), kinetic energy (KE), and temperature (Temp.) of the entire system when calculated with the NVE and NVT ensembles. By definition, the total energy remains constant for the NVE ensemble, and the temperature is nealy fixed for the NVT ensemble. In the Berendsen heat bath method, which is a king of NVT methods, the total energy fluctuates, but the kinetic energy and temperature stay nearly constant. (Actually, there are small fluctuations, but they are too small to be seen due to the scale of the graph.)

66a541abbc0d44b2805ffcd7a0f1b9da 50f88b6047e0401cb56ee3a79ae0d4e7

Fig.6-2c. Comparison of energies and temperatures between NVE (top) and NVT (bottom) calculations.

The above plot regarding NVT can be prepared as follows. See previous section regarding NVE.

[11]:
import pandas as pd

df = pd.read_csv(
    log_filename,
    delim_whitespace=True,
    names=["Time[ps]", "Etot/N[eV]", "Epot/N[eV]", "Ekin/N[eV]", "T[K]",
           "stressxx", "stressyy", "stresszz", "stressyz", "stressxz", "stressxy"],
    skiprows=1,
    header=None,
)
df
[11]:
Time[ps] Etot/N[eV] Epot/N[eV] Ekin/N[eV] T[K] stressxx stressyy stresszz stressyz stressxz stressxy
0 0.1 -3.4368 -3.4670 0.0302 233.3 1.586 1.311 0.904 -0.118 -0.025 -0.255
1 0.2 -3.5662 -3.5964 0.0301 233.2 1.318 1.447 1.152 0.081 -0.109 -0.138
[12]:
import numpy as np
import matplotlib.pyplot as plt


fig = plt.figure(figsize=(10, 5))

ax1 = fig.add_subplot(4, 1, 1)
ax1.set_xticklabels([])  # x axis label
ax1.set_ylabel('Tot E (eV)')  # y axis label
ax1.plot(df["Time[ps]"], df["Etot/N[eV]"], color="blue",alpha=0.5)

ax2 = fig.add_subplot(4, 1, 2)
ax2.set_xticklabels([])  # x axis label
ax2.set_ylabel('P.E. (eV)')  # y axis label
ax2.plot(df["Time[ps]"], df["Epot/N[eV]"], color="green",alpha=0.5)

ax3 = fig.add_subplot(4, 1, 3)
ax3.set_xticklabels([])  # x axis label
ax3.set_ylabel('K.E. (eV)')  # y axis label
ax3.set_ylim([0.0, 0.2])
ax3.plot(df["Time[ps]"], df["Ekin/N[eV]"], color="orange",alpha=0.5)

ax4 = fig.add_subplot(4, 1, 4)
ax4.set_xlabel('time (ps)')  # x axis label
ax4.set_ylabel('Temp. (K)')  # y axis label
ax4.plot(df["Time[ps]"], df["T[K]"], color="red",alpha=0.5)
ax4.set_ylim([500., 1700])

fig.suptitle("Time evolution of total, potential, and kinetic energies, and temperature in NVT.", y=0)

#plt.savefig("Fig6-2_fcc-Al_NVT-Berendsen_1000K_E_vs_t.png")  # <- Use if saving to an image file is desired
plt.show()
_images/6_2_md-nvt_13_0.png

A note on calculations using the Berendsen heat bath method.

  1. Setting the appropriate time constant. The Berendsen heat bath method has a time constant \(\tau\) as a temperature control parameter. This \(\tau\) is a parameter that determines the relaxation time to reach the target temperature. The smaller \(\tau\) is, the faster the convergence speed is. If \(\tau\) is too small, the calculation becomes unstable and the temperature value shows unnatural behavior. In general, \(\tau\) should be several to five times of the time step of the MD simulation. The following graph shows the behavior of the overall system temperature when the time constant \(\tau\) is varied for the above fcc-Al system, see Fig. 6-2d. As mentioned above, when \(\tau\) is 0.2 fs, the system exhibits very unstable behavior, with large temperature transitions up and down, and the system is unable to maintain equilibrium. For \(\tau\) of 1 or 2 fs, the temperature is quite tightly controlled, and for \(\tau\) of 5 to 10 fs, the temperature is loosely controlled, and the fluctuations are relatively large but in a more realistic range.

b8a2b458651c4fff86e6da58005ea625

Fig.6-2d. Temperature evolution with various time constant \(\tau\).

  1. The case when the system converges to unrealistic state. The Berendsen dynamic is a velocity scaling method, in which the velocity of all atoms in the calculation is changed uniformly. By moving the entire velocity distribution, the temperature of the entire system is brought into alignment with the target temperature, but this operation itself is very artificial and differs from a thermodynamically correct canonical (NVT) ensemble. The major problem of this method is that it creates unrealistic states when applied to systems in which several substances with very different behaviors coexist. For example, if one wants to use NVT to calculate the interface between a solid and a liquid, Berendsen heat bath method might result in an unnatural situation where one will be at a high temperature and the other at a low temperature, and the total temperature will be at the target temperature. (A specific example is explained at the end of this section.)

If you are careful about the above points and do not require strictly correct thermodynamic distributions, the Berendsen heat bath method is a very convenient calculation method because it converges to an equilibrium state quickly and is simple and easy to understand in principle. For example, if you simply want to create a state with thermal fluctuation, this is a very effective calculation method.

Calculation example 2: Ion diffusion in solids

Let us consider ion diffusion in solid electrolytes as the next calculation example of the NVT ensemble.

Li\(_{10}\)GeP\(_{2}\)S\(_{12}\) (LGPS) is one of the materials that have attracted significant attention in recent years as a solid electrolyte for all-solid-state batteries. It is characterized by very high ionic conductivity (12 mS/cm @ 298 K) comparable to that of liquid electrolytes. The crystal structure is as follows (The reference structure is from Li10Ge(PS6)2_mp-696138_computed in Materials Project. The structure shown below is based on that structure and optimized with PFP).

Input cif file is modified from
A. Jain, S.P. Ong, G. Hautier, W. Chen, W.D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, K.A. Persson (*=equal contributions)
The Materials Project: A materials genome approach to accelerating materials innovation APL Materials, 2013, 1(1), 011002.
Licensed under CC BY 4.0
[13]:
atoms = ase.io.read("../input/ch6/LGPS_2x2x1_pfp-opt.cif")

view_ngl(atoms, representations=["ball+stick"])

247cca7eac1949ae8db632b433688b65

Fig.6-2e. Crystal structure of Li\(_{10}\)GeP\(_2\)S\(_{12}\).

For convenience of calculation, unit cell is expanded to 2x2x1. This is because we expect that increasing the number of Li ions will reduce the statistical error. Although the calculation time may increase slightly, there is no significant problem in practice, because the extended structure has only about 200 atoms in total. We would like to determine the Li ion diffusion coefficient within this material using MD simulations. Since the composition of this material is complex and there is no general classical force field available, we will use the PFP implemented in Matlantis. (It is possible that a general-purpose force field could exist if we do a more thorough search of academic papers, but since it was not immediately found and even if it were found, the difficulty of implementation is unknown, the time, effort, and accuracy involved in these processes make the use of a general-purpose machine-learning potential a major advantage.)

The script for running this NVT-MD simulation is shown below. The following is a test calculation with 10 MD steps. In the actual calculation for computing diffusivity, 10^6 steps calculation is performed at 1 fs/step and 1 ns MD simulation is performed.

[14]:
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
estimator = Estimator(model_version="v2.0.0",calc_mode=EstimatorCalcMode.CRYSTAL_U0)
calculator = ASECalculator(estimator)

import ase
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.langevin import Langevin
from ase.md import MDLogger
from ase import units
from time import perf_counter

# Set up a crystal
atoms = ase.io.read("../input/ch6/LGPS_2x2x1_pfp-opt.cif")
atoms.pbc=True
print("atoms = ",atoms)
atoms.calc = calculator

# input parameters
time_step    = 1.0    # fsec
temperature  = 773    # Kelvin
num_md_steps = 10     # Total MD step, for testing.
# num_md_steps = 1e6    # Total MD step, for actual run.
num_interval = 1      # Output printing interval
friction_coeff = 0.002
print(f"friction_coeff = {friction_coeff:.3f}")

# Define output filenames
output_filename = "./output/ch6/Li10Ge(PS6)2_mp-696138_Langevin_773K"
print("output_filename = ",output_filename)
log_filename = output_filename + ".log"
print("log_filename = ",log_filename)
traj_filename = output_filename + ".traj"
print("traj_filename = ",traj_filename)

# run MD
dyn = Langevin(atoms, time_step*units.fs, friction=friction_coeff, temperature_K=temperature, loginterval=num_interval, trajectory=traj_filename)

# Print statements
def print_dyn():
    imd = dyn.get_number_of_steps()
    etot  = atoms.get_total_energy()
    temp_K = atoms.get_temperature()
    stress = atoms.get_stress(include_ideal_gas=True)/units.GPa
    stress_ave = (stress[0]+stress[1]+stress[2])/3.0
    elapsed_time = perf_counter() - start_time
    print(f"  {imd: >3}   {etot:.3f}    {temp_K:.2f}    {stress_ave:.2f}  {stress[0]:.2f}  {stress[1]:.2f}  {stress[2]:.2f}  {stress[3]:.2f}  {stress[4]:.2f}  {stress[5]:.2f}    {elapsed_time:.3f}")


dyn.attach(print_dyn, interval=num_interval)
dyn.attach(MDLogger(dyn, atoms, log_filename, header=True, stress=True, peratom=True, mode="w"), interval=num_interval)

# Now run the dynamics
start_time = perf_counter()
print(f"    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)")
dyn.run(num_md_steps)
atoms =  Atoms(symbols='Ge8Li80P16S96', pbc=True, cell=[[17.2189, 0.0, 0.0], [-0.07190443131522067, 17.742454297891182, 0.0], [-0.13259731631828015, -0.4464045597143071, 13.005365377826626]], spacegroup_kinds=...)
friction_coeff = 0.002
output_filename =  ./output/ch6/Li10Ge(PS6)2_mp-696138_Langevin_773K
log_filename =  ./output/ch6/Li10Ge(PS6)2_mp-696138_Langevin_773K.log
traj_filename =  ./output/ch6/Li10Ge(PS6)2_mp-696138_Langevin_773K.traj
    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)
    0   -699.887    0.00    0.00  0.00  0.00  0.00  -0.00  -0.00  0.00    0.090
    1   -699.880    0.28    0.00  0.00  0.00  0.00  0.00  -0.00  -0.00    0.201
    2   -699.871    0.63    0.00  0.00  0.00  0.00  0.00  0.00  -0.00    0.276
    3   -699.864    0.93    0.00  0.00  0.00  0.00  0.00  0.00  -0.00    0.602
    4   -699.857    1.21    0.00  -0.00  0.00  0.00  0.00  0.00  -0.00    0.716
    5   -699.849    1.49    0.00  -0.00  0.00  -0.00  0.00  0.00  -0.00    0.963
    6   -699.840    1.85    -0.00  -0.00  0.00  -0.00  0.00  0.00  -0.00    1.216
    7   -699.833    2.11    -0.00  -0.00  -0.00  -0.00  0.00  0.00  -0.00    1.434
    8   -699.825    2.43    -0.01  -0.01  -0.00  -0.01  -0.00  0.00  -0.00    1.690
    9   -699.816    2.74    -0.01  -0.01  -0.01  -0.01  -0.00  0.00  -0.00    1.761
   10   -699.808    3.00    -0.01  -0.01  -0.01  -0.01  -0.00  0.00  -0.00    2.090
[14]:
True

When MD is performed in the time of actual (non-test) calculation, the following figure is obtained with trajectory visualization. You can see that the purple Li ions are moving (diffusing) significantly. Although it is shown as if it is going to an empty space, Ge, P, and S actually exist outside the unit cell due to the periodic boundary condition.

476a9615f8194a0d942c68e94448ed25

Li\(_{10}\)GeP\(_{2}\)S\(_{12}\) crystal at 773 K.

Here, we have performed NVT-MD simulations using the Langevin heat bath method. This method is a statistical approach that controls the temperature by applying a force on each atom that follows a certain probability distribution, often called a random force [1,4].

Specifically, the equation of motion for Langevin heat bath method can be written in the following form

\[m_i\ddot{\mathbf{x}_i}(t)=\mathbf{F}(\mathbf{x}_i(t)) -\gamma m_i\dot{\mathbf{x}_i}(t)+\sqrt{2k_B T\gamma m_i}\mathbf{\eta}_i(t)\]

where \(m_i\) is the mass of each atom, \(\mathbf{x}_i(t)\) is the coordinate vector of atom \(i\) at time \(t\), \(\mathbf{F}(\mathbf{x}_i(t))\) is the force on atom \(i\) and \(\gamma\) is the friction coefficient. The last term on the right-hand side is the random force, and the function \(\mathbf{\eta}_i(t)\) must satisfy the relation \(<\mathbf{\eta}_i(0)\mathbf{\eta}_j(t)>=\delta_{ij}\delta(t)\). This allows us to reproduce the canonical distribution.

As mentioned in the problem with the Berendsen heat bath method in example 1, if the motions of the fast diffusing Li and other anionic structures are extremely different, there is a risk that the Berendsen heat bath method will produce unrealistic results due to the divergence of these behaviors. The Langevin method controls the forces on individual atoms and thus avoids this problem.

The Langevin heat bath method requires a friction coefficient \(\gamma\) to be set, the value between 1e-2 and 1e-4 is generally considered to be in the appropriate range although it depends on the target of calculation. The larger the value, the tighter the convergence conditions and the faster the target temperature is reached. Generally, there is no problem with the above range of values, but we recommend checking the sensitivities for these friction coefficient values as a preliminary study when considering a new calculation target.

In this example, the simulation time is set to 1 ns at 773 K, which is a relatively high temperature state. The reason for performing MD at high temperatures is to guarantee sufficient Li ion diffusion, which may not be a realistic temperature range. While there is no guarantee that the actual material will be stable in this temperature range, the density is kept constant in the NVT ensemble, so the structure remains relatively stable and the diffusion coefficient can be calculated. Calculations at such high temperatures are a commonly used technique and are suitable for the purposes of this tutorial because they tend to produce relatively clean data.

We will calculate the diffusion coefficient from the above MD Trajectory. First, we calculate the mean-squared dispacement (MSD) from the initial position of the atom of interest [5].

\[MSD = \frac{\Sigma_i(\mathbf{x}_i(t)-\mathbf{x}_i(0))^2}{N}\]

where \(\mathbf{x}_i(t)\) is the coordinate of atom \(i\) at time \(t\) and \(N\) is the number of atoms. The following graph shows this MSD with time on the horizontal axis.

01f291c1a3424f8498a9b70f96bed7d3

Fig.6-2f. Mean-square displacement of Li-ion in Li\(_{10}\)GeP\(_{2}\)S\(_{12}\) crystal at 773 K.

From the MSD obtained here, the diffusion coefficient can be calculated using the following equation

\[D = \frac{MSD}{6\Delta t}\]

It is important to note that this equation assumes the system has fully reached equilibrium. If the calculation results include data that are not in equilibrium, it is necessary to remove that range of data. Also, if the MSD plot contains noise when the temperature is low, diffusion is unlikely to occur, or the number of atoms of interest is small, a more stable result can be obtained by calculating the slope using linear regression or other methods.

The result is a diffusion coefficient of 2.7x10\(^{-5}\) cm\(^2\)/s, which is generally in good agreement with the previously reported value [6]. However, we have not optimized the calculation conditions for the sake of simplicity. If you actually want to perform accurate calculations efficiently, you should examine the dependence on these parameters such as the number of calculation steps, time step, temperature, and size of the system to be calculated in advance and examine them carefully.

Thus, the NVT ensemble can be used to evaluate the diffusion of a substance at a specific temperature. For example, if one tries to study this using NVE, it is difficult to correctly evaluate the temperature because there is a difference between the initial temperature and the temperature at which an equilibrium state is reached. Therefore, the NVT ensemble is an effective method to observe the thermal equilibrium state at a specific temperature.

For more information on this calculation example, please also see the following

Calculation example 3: Oxidation reaction on a metal surface

As a final example, we consider the case of oxidation reaction on a metal surface. Consider a clean fcc-Al (111) surface, exposed to an oxygen atmosphere. Since Al readily oxidizes at room temperature, we will attempt to reproduce such a situation in our calculations.

The model to be calculated is a slab model of fcc-Al (111) with a 10 Å vacuum layer and randomly placed oxygen molecules.

[15]:
atoms = ase.io.read("../input/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K_wrapped.cif")

view_ngl(atoms, representations=["ball+stick"])

bf96cdbb68d14152b28aebe5bbb987bd

Fig.6-2g. Starting structure of fcc-Al(111) slab with 20 O\(_{2}\) molecules.

In this case, we will use PFP implemented in Matlantis since we need to deal with the reaction between the metal interface and gas-phase molecules. Generally, a calculation method such as DFT is necessary to study the reaction with high accuracy, but the large number of atoms and large cell size make it very computationally prohibiting. On the other hand, if PFP is used, the reaction process can be reproduced with accuracy comparable to DFT, and since sufficient number of atoms is can be calculated, it is possible to perform a trial study in more realistic condition at a very high speed. Let us now run an MD simulation of the NVT ensemble at 300 K. The calculation script is as follows (The num_md_steps is set small for the test calculation. Please increase it for a production run.)

[16]:
import os
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
estimator = Estimator(model_version="v2.0.0",calc_mode=EstimatorCalcMode.CRYSTAL_U0_PLUS_D3)
calculator = ASECalculator(estimator)

from ase.io import read,write
from ase.build import bulk
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution,Stationary
from ase.md.npt import NPT
from ase.md import MDLogger
from ase import units
from time import perf_counter

# Set up a crystal
atoms = read("../input/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K_wrapped.cif")
atoms.pbc = True
atoms.calc = calculator
print("atoms = ",atoms)

# input parameters
time_step    = 1.0    # fsec
temperature  = 300    # Kelvin
num_md_steps = 200    # For test
# num_md_steps = 10000  # For actual run
num_interval = 100
ttime        = 20.0   # tau_T thermostat time constant in fsec
print(f"ttime = {ttime:.3f}")

output_filename = "./output/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K"
print("output_filename = ",output_filename)

log_filename = output_filename + ".log"
print("log_filename = ",log_filename)
traj_filename = output_filename + ".traj"
print("traj_filename = ",traj_filename)


# Set the momenta corresponding to the given "temperature"
MaxwellBoltzmannDistribution(atoms, temperature_K=temperature,force_temp=True)
Stationary(atoms)  # Set zero total momentum to avoid drifting

# run MD
dyn = NPT(atoms,
          time_step*units.fs,
          temperature_K = temperature,
          externalstress = 0.1e-6 * units.GPa,  # Ignored in NVT
          ttime = ttime*units.fs,
          pfactor = None,   # None for NVT
          loginterval=num_interval,
          trajectory=traj_filename
          )

# Print statements
def print_dyn():
    imd = dyn.get_number_of_steps()
    etot  = atoms.get_total_energy()
    temp_K = atoms.get_temperature()
    stress = atoms.get_stress(include_ideal_gas=True)/units.GPa
    stress_ave = (stress[0]+stress[1]+stress[2])/3.0
    elapsed_time = perf_counter() - start_time
    print(f"  {imd: >3}   {etot:.3f}    {temp_K:.2f}    {stress_ave:.2f}  {stress[0]:.2f}  {stress[1]:.2f}  {stress[2]:.2f}  {stress[3]:.2f}  {stress[4]:.2f}  {stress[5]:.2f}    {elapsed_time:.3f}")


dyn.attach(print_dyn, interval=num_interval)
dyn.attach(MDLogger(dyn, atoms, log_filename, header=True, stress=True, peratom=True, mode="w"), interval=num_interval)

# Now run the dynamics
start_time = perf_counter()
print(f"    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)")
dyn.run(num_md_steps)
atoms =  Atoms(symbols='Al72O40', pbc=True, cell=[8.59135, 9.92043, 31.6913], spacegroup_kinds=..., calculator=ASECalculator(...))
ttime = 20.000
output_filename =  ./output/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K
log_filename =  ./output/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K.log
traj_filename =  ./output/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K.traj
    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)
  100   -384.135    330.91    1.12  1.05  1.33  0.99  -0.25  0.18  0.22    17.075
  200   -391.182    375.24    0.88  1.38  0.76  0.50  0.39  0.23  -0.03    34.796

This time, we use the Nosé-Hoover heat bath method. In the above code, the NPT module is specified, but you can simulate NVT using the Nosé-Hoover heat bath by specifying pfactor=None.

In this method, two new independent variables (\(\xi\) and \(s\)) are introduced into the equations of motion to create the contact with a virtual heat bath and control the temperature [1,3-5]. The equation of motion has the following form

\[\dot{\mathbf{x}_i} = \mathbf{p}_i/m_i\]
\[\dot{\mathbf{p}_i} = \mathbf{F_i}-\zeta \mathbf{p}_i\]
\[\dot{\zeta} = \frac{1}{\tau_T^2}\left(\frac{T(t)}{T_o}-1\right)\]
\[\dot{s} = 3N s\zeta\]

The control variable is a new time constant (\(\tau_T\)) which the user must set to an appropriate value. According to the ASE manual, a value around 25 fs is considered appropriate. Too small value will make the calculation unstable, while too large value will result in extremely poor convergence and take too long to reach the target temperature.

Let’s take a look at the resulting trajectory. (The following is an APNG file only, but the name of the trajectory file is given in the caption, so you can actually visualize the trajectory file on Matlantis to see how it behaves.)

[17]:
from pfcc_extras.visualize.view import view_ngl
from ase.io import Trajectory
traj = Trajectory("./output/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K_wrapped.traj")

v = view_ngl(traj, representations=["ball+stick"], replace_structure=True)
display(v)

3652890ec1c942baaef043ec7950ed36

Fig.6-2h. Starting structure of fcc-Al(111) slab with 20 O\(_{2}\) molecules.

As you can see, oxygen molecules, which were initially randomly placed in the vacuum, adsorb and break into individual oxygen atoms on the Al surface in a short period of time to form aluminum oxide. Aluminum forms an oxide phase even at room temperature under low oxygen partial pressure (\(pO_2\)=10\(^{-18}\) atm). So in a condition like 20 oxygen molecules are packed into a space of only 9x10x10 Å\(^3\), an aluminum oxide film is formed even in only 10 ps time (If you want to know more about the oxidation equilibrium of aluminum, see Ellignham diagram. This will give you an idea of the specific temperature and oxygen partial pressure regions under equilibrium condition.)

As can be seen from the above data, MD simulations allow us to observe individual reaction processes in detail. For example, in Fig.6-2e, the \(O_2\) molecule circled in green from the left is adsorbed on the Al bridge site on the surface, bonded to two Al atoms, and then the \(O_2\) is dissociated and incorporated into the surface structure. And the Al atoms coordinated to the three O atoms adsorb to the next \(O_2\) molecule.

d50309ce300c4cb68e183f77f28372a6

Fig.6-2i. Reaction process of an O\(_2\) molecule on fcc(111) Al.

As we have just seen above, NVT-MD simulations allow detailed observation of the dynamics of the adsorption process on gases and solid surfaces. If this were, a more complex molecule or surface structure, the analysis might be more difficult, but the basic process and concept would remain the same.

Here is the final remark on the Nosé-Hoover heat bath. This method is one of the most commonly used heat bath methods because it is known to reproduce statistically correct canonical ensembles. A drawback of this method is that it cannot reproduce the correct canonical ensemble in some special cases (e.g. harmonic oscillators) when applied to systems that satisfy unintended conservation laws. This means that the possible physical states cannot be sampled correctly. For practical use, however, this problem can be avoided by fixing the center of mass of the system to be computed, making it highly versatile. (Generally speaking, there are more advanced methods (such as the Nosé-Hoover chain method) which overcome above complication, so those interested should check the references [4,5].)

Summary

Finally, we summarize the simulation methods currently available for ASE’s NVT ensemble.

Class

Ensemble

Parameter

Description

Pros

Cons

NVTBerensen

NVT

time constant (\(\tau_T\))

A type of velocity scaling method. Deterministic.

Simple and general.

Non-canonical distribution. Controls only the temperature of the entire system.

Langevin

NVT

friction parameter (\(\gamma\))

Temperature control by random forces and virtual friction terms. Stochastic.

Reproduce canonical distribution. Temperature control of individual atoms.

Statistical method.

NPT(pfactor=None)

NVT

time constant (\(\tau_T\))

Nosé-Hoover heat bath. Deterministic.

Reproduces canonical distribution for many systems.

Noncanonically distributed systems also exist. Artificial modes of motion can be introduced. Only the temperature of the entire system is controlled.

Andersen

NVT

frequency factor (\(\nu\))

Modifies the velocity of a randomly selected group of atoms. Stochastic.

Simple. Temperature control of individual atoms.

Non-canonical distribution. Statistical Method.

If you are unsure of which method to apply when performing a NVT ensemble calculation, the Nosé-Hoover heat bath method is often a reasonable choice. The theoretical foundation of this methos is sound despite few drawbacks in some special cases. In addition, its implemented in a wide range of molecular simulation tools and is readily available. It is a deterministic method and is also suitable in terms of reproducibility when one wants to study the trajectories (i.e. coordinates, velocity trajectories) of each atom in phase space.

The Langevin heat bath method is often used to calculate liquid systems or to perform Brownian motion-like calculations. In principle, this method can reproduce the correct NVT ensemble and control temperature efficiently. However, it is a stochastic method and cannot be expected to reproduce precise atomic trajectories each time the simulation is run. The ASE manual states that this method is not suitable for controling the temperature when applied to systems with endothermic reactions. However, we have not been able to confirm such a fact, at least we have no data to support it.

The Berendsen heat bath method (NVTBerendsen) is a simple, generic method and easy to implement. This may be a good starting point method for test calculations. However, it cannot be said to reproduce the NVT ensemble strictly, and the thermal fluctuations are artificial and lack precision. It is useful for situations where you want to create a thermally fluctuating structure, create a non-equilibrium state, etc.

Finally, there is a method implemented in ASE called the Andersen heat bath method, but it has no particular advantage over the other methods and is not as widely used as other methods mentioned above. So we decided to save ourselves from giving the discussion on this method here.

[Advanced] Parameter dependency of the Nosé-Hoover heat bath method

Let us examine the dependence of the temperature of a system on the time constant \(\tau\) of the Nosé-Hoover heat bath. A fcc-Al crystal structure (3x3x3 unit cells) was equilibrated for different \(\tau_t\) as follows. PFP(v2.0.0) was used in the calculation, and no initial temperature was set, i.e. the starting temperature was 0 K.

8cd385e24c8e4453935119e1e7c9d265

Fig.6-2j. Time evolution of temperature in Nose-Hoover thermostat

Let’s look at three separate areas because the behavior varies greatly with the value of \(\tau\).

  • The smaller the value, the faster the speed of approaching the target temperature. But at 5 fs, even if the average value reaches the target temperature, the amplitude of the temperature oscillation is large and fluctuates from 750 K to 100 K. This makes it difficult to assume a “constant temperature” even if the average value is stable.

  • In the 10-40 range, the convergence to the target temperature is relatively fast and the convergence of temperatures is stable (bottom of the above figure). However, at 40 fs, there appears to be an artificial periodicity in the period of oscillation.

  • When \(\tau\) is more than 100 fs, the temperature rises slowly and oscillates in very large waves. Although the temperature may eventually converge to the target temperature, the temperature reaches several thousand K in some places in the large wave-like behavior, which may have a significant impact depending on the target under consideration.

In the above example, \(\tau\) seems to be reasonable around 10 or 20 fs. However, there is a point in the beginning of the calculation where the temperature exceeds 1000 K, momentarily. Is it possible to avoid this problem? Let’s try to initialize the temperature using the Maxwell-Boltzmann distribution as a conditioning method, or create an equilibrium state using a more rapidly converging method such as the Berendsen thermostat. The Berendsen heat bath was run for only 1 ps after the start of the calculation, and then switched to the Nosé-Hoover heat bath. How effective is such preconditioning in practice?

de488335121c42be8d1821c335963962

Initialized by the Maxwell-Boltzmann Distribution

9a07eee124614ff595df350b3e200b53

Initialized by Berendsen thermostat (τ=1.0fs)

12cacbe99e6c4617b45301c642bcf110

Initialized by Berendsen thermostat (τ=10.0fs)

9fb423871d924535b120032af6e266aa

Initialized by Berendsen thermostat (τ=30.0fs)

Fig.6-2k. Time evolution of temperature with various initialization methods:

Log-Log scale (left), Linear scale (right)

When the Nosé-Hoover method was simply applied to a structurally optimized system (i.e., 0 K), the temperature instantaneously rose to about 1000-2000 K. In contrast, the temperature increase in the initial stage of the calculation can be limited to about 500-600 K using either of these methods. The problem was not completely eliminated, but it seems to have been mitigated to a certain extent.

[Advanced] NVT ensemble in 2-phase coexisting state in Berendsen heat bath

Here, we present a typical example the problems with the Berendsen heat bath method. The target system is an Al atom cluster in water. The system consists of 165 molecules of water molecules and an Al atom cluster with 135 atoms. The water molecules are designed to be roughly 1 g/cc in density. The details of the setup are not very important, but we have set a relatively realistic range of densities and temperatures so that readers can see the relevance in practice. Below is the image of the system under consideration.

a4d5d7f530bb4ebc9be388c39fca1639

Fig.6-2l. Starting structure of fcc-Al(111) slab with 20 O\(_{2}\) molecules.

(File: ../input/ch6/165h2o_and_135Al_20x20x20_opt.cif)

Let us apply the Berendsen and Nosé-Hoover heat bath methods to this system in the MD simulation to bring it to thermal equilibrium at 100°C. The MD time step is 1 fs and the computation time is 10 ps. Below are plots of temperatures for water molecules, Al clusters, and the system as a whole.

b8091bbba8cb493db7bfa35c2d356e9d 124fee519b8646f9a68e492bcd3b234e

Fig.6-2m. Comparison of the temperature evolution in Berendsen and Nosé–Hoover thermostats. Each line represents the temperature of the entire system (blue), Al cluster (gray), and water (orange), respectively.

The lines represent the temperatures of the entire system (blue), Al cluster (gray), and water (orange), respectively in each graph. In the case of the Berendsen heat bath, the temperature of the Al cluster is about 420 K on average, whereas the temperature of H\(_2\)O remains around 360 K. This happens because the Berendsen bath method controls the velocities of all atoms with a single scaling parameter. Although the overall temperature of the system is good and stable, the atoms behave very differently from the Al clusters and the H\(_2\)O molecules, so their temperatures are decoupled. It is like a hot Al sphere floating in relatively warm water. This is a good example of a phenomenon that cannot happen in reality and does not reproduce the correct thermal equilibrium state.

In contrast, when the Nosé-Hoover heat bath method is used, the temperatures of the entire system and each of its components are consistent and physically reasonable.

Reference

[1] M.P. Allen and D.J. Tildesley, “Computer simulaiton of liquid”, 2nd Ed., Oxford University Press (2017) ISBN 978-0-19-880319-5. DOI:10.1093/oso/9780198803195.001.0001

[2] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak, “Molecular dynamics with coupling to an external bath”, J. Chem. Phys. (1984) 81 3684, DOI:10.1063/1.448118

[3] “GROMACS 2019-rc1 Reference Manual”, https://manual.gromacs.org/documentation/2019-rc1/reference-manual/index.html

[4] M.E. Tuckerman, “Statistical mechanics: Theory and molecular simulation”, Oxford University Press (2010) ISBN 978-0-19-852526-4. https://global.oup.com/academic/product/statistical-mechanics-9780198525264?q=Statistical%20mechanics:%20Theory%20and%20molecular%20simulation&cc=gb&lang=en#

[5] D. Frenkel and B. Smit “Understanding molecular simulation - from algorithms to applications”, 2nd Ed., Academic Press (2002) ISBN 978-0-12-267351-1. DOI:10.1016/B978-0-12-267351-1.X5000-7

[6] Mo et al. “First Principles Study of the Li10GeP2S12 Lithium Super Ionic Conductor Material”, Chem.Mater. (2012) 24, 15-17 https://pubs.acs.org/doi/10.1021/cm203303y