NPT ensemble

TL;DR

  • In NPT-MD simulations, pressure and temperature are controlled and remain constant once the system reaches equilibrium. There are two common pressure control (barostat) methods available in ASE: the Parrinello-Rahman method and the Berendsen method.

  • Typical applications include simulations of thermal expansion, phase transitions, and pressurization of solids and fluids.

  • In the Parrinello-Rahman method, all degrees of freedom of the simulation cell are variable, and the control parameter (pfactor) must be set appropriately.

  • Berendsen barostat, like Berendsen thermostat, can control pressure efficiently for convergence; Berendsen barostat can be calculated in two modes: with fixed cell angles and independently variable cell lengths, or with fixed cell length ratios. compressibility must be properly set as an input parameter.

In this section, we will discuss a method to create an equilibrium state in which both pressure and temperature are constant. The pressure control mechanism is generally referred to as barostat, and is used in conjunction with the thermostat method described in section 6-2 to generate a state distribution called an isothermal-isobaric ensemble (also known as NPT).

Phenomena that can be studied, in principle, in NPT-MD simulations include

  • Coefficient of thermal expansion of solids

  • Prediction of melting point

  • Phase transitions in solids

  • Density prediction of fluids (gases and liquids)

etc.

The reason why the word “in principle” is used is that the reproducibility of these phenomena depends greatly on the accuracy of the force field used in the calculation, and in particular, it is known to be very difficult to predict intermolecular forces and fluid states that depend on small energy differences. In this tutorial, we will learn about NPT-MD through the case study of solids, which are considered to have relatively high accuracy.

First, we will review the NPT-MD methods implemented in ASE that we will use in this tutorial. There are three types of implementations available for ASE

Class

Ensemble

Parameter

thermostat

barostat

description

NPT

NPT

time constant (\(\tau_t\)), pressure factor (pfactor)

Nosé–Hoover

Parrinello-Rahman

セルの全All degrees of freedom of the cell are variable and controllable

NPTBerendsen

NPT

\(\tau_t\),\(\tau_P\),\(\beta_T\)

Berendsen

Berendsen

Cell shape is maintained and only volume changes

InhomogeneousBerendsen

NPT

\(\tau_t\),\(\tau_P\),\(\beta_T\)

Berendsen

Berendsen

Cell angles are preserved, but pressure anisotropy can be taken into account

The second and third methods in the table are essentially the same Berendsen barostat. (The third method, InhomogeneousBerendsen, is not mentioned in the ASE manual, but is defined as a class within ASE along with NPTBerendsen.) Thus, there are only two methods available within the ASE framework: the Parrinello-Rahman method and the Berendsen method. The heat bath itself uses the methods described in section 6-2, so the characteristics of these heat bath methods should be carefully considered.

Let us first look at the Parrinello-Rahman method, which has relatively high system flexibility and versatility.

Equations of motion for the Parrinello-Rahman method

The Parrinello-Rahman method is a so-called extended system calculation method, which assumes that the system under consideration is hypothetically connected to an external system of constant temperature and pressure, as in the Nosé-Hoover heat bath method. In this case, the equations of motion are written as follows (For details of the derivation, please refer to references [1-3]).

\[\dot{\mathbf{x}_i} = \mathbf{p}_i/m_i + \eta(\mathbf{x}_i-\mathbf{R}_o)\]
\[\dot{\mathbf{p}_i} = \mathbf{F_i}-(\eta+\zeta)\mathbf{p}_i\]
\[\dot{\zeta} = \frac{1}{\tau_T^2}\left(\frac{T(t)}{T_o}-1\right)-3\eta\zeta\]
\[\dot{s} = 3(N-1) s\zeta\]
\[\dot{\mathbf{\eta}} = \frac{V}{\tau_P^2Nk_BT_o}\left(\mathbf{P(t)}-P_o\mathbf{I}\right)+3\frac{\tau_T^2}{\tau_P^2}\zeta^2\mathbf{I}\]
\[\dot{\mathbf{h}}=\mathbf{\eta}\mathbf{h}\]

Besides the terms for the Nosé-Hoover heat bath, we have the pressure control time constant \(\tau_P\), the system center of mass \(R_o\), the target external pressure \(P_o\), and the simulation cell volume \(V\). The \(\eta\) is the variable for the pressure control degrees of freedom. \(\mathbf{h} = (\mathbf{a}, \mathbf{b}, \mathbf{c})\) and \(\mathbf{a},\mathbf{b},\mathbf{c}\) are cell vectors defining each edge of the simulation cell, respectively.

In addition to temperature \(T_o\) and pressure \(P_o\), there are two other values that the user must set in the above equation: \(\tau_T\) and \(\tau_P\). First, the case where \(\tau_T\) and \(\tau_P\) are simply set to random values is shown below. Here \(\tau_T\) is set to 20 fsec. Although \(\tau_P\) is not specified directly, the value called pfactor is \(\tau_P^2B\). The \(B\) refers to the bulk modulus, and this value must be calculated and specified in advance. However, since the exact value of \(\tau_P\) itself is not known in advance, there is no way to specify pfactor. Also, the value of \(B\) cannot be calculated for anisotropic structures or when different types of materials coexist in the system. Therefore, we will set an approximate value of pfactor to examine the behavior of the barostat. As mentioned in the following example, at least for crystalline metal systems, a value on the order of about 10\(^6\) GPa\(\cdot\)fs\(^2\) to 10\(^7\) GPa\(\cdot\)fs\(^2\) seems to give good convergence and stability in the calculation. For other materials, it is recommended to play with the value of pfactor and check the behavior of the volume change as a preliminary study.

Calculation example: Coefficient of thermal expansion

Now, we will use the Nosé-Hoover thermostat and the Parrinello-Rahman barostat (ASE’s NPT class) to compute the coefficient of thermal expansion of a solid as an example. The system is equilibrated at temperatures of given increment, and the thermal expansion coefficient is calculated from the average value of the lattice constant at each temperature.

For simplicity, we will use fcc-Cu for this calculation. A sample script is shown below. The temperature is varied from 200 K to 1000 K in 100 K interval and the external pressure is set to 1 bar. The structure is fcc-Cu extended to 3x3x3 unit cells with 108 atoms. In the example, the ASAP3-EMT force field is used for speed, but the same scheme can be used with PFP. The time step size is set to 1 fs, and the 20 ps simulation is found to be sufficient to reach equilibrium.

[1]:

import ase
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

calc_type = "EMT"
# calc_type = "PFP"

if calc_type == "EMT":
    # ASAP3-EMT calculator
    from asap3 import EMT
    calculator = EMT()
elif calc_type == "PFP":
    # PFP calculator
    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)
else:
    raise ValueError(f"Wrong calc_type = {calc_type}!")


# Set up a crystal
atoms_in = bulk("Cu",cubic=True)
atoms_in *= 3
atoms_in.pbc = True
print("atoms_in = ",atoms_in)

# input parameters
time_step    = 1.0    # fsec
#temperature = 300    # Kelvin
num_md_steps = 20000
num_interval = 10

sigma   = 1.0     # External pressure in bar
ttime   = 20.0    # Time constant in fs
pfactor = 2e6     # Barostat parameter in GPa
temperature_list = [200,300,400,500,600,700,800,900,1000]

# 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}")


# run MD
for i,temperature in enumerate(temperature_list):
    print("i,temperature = ",i,temperature)

    print(f"sigma = {sigma:.1e} bar")
    print(f"ttime = {ttime:.3f} fs")
    print(f"pfactor = {pfactor:.3f} GPa*fs^2")

    temperature_str = str(int(temperature)).zfill(4)
    output_filename = f"./output/ch6/fcc-Cu_3x3x3_NPT_{calc_type}_{temperature_str}K"
    log_filename = output_filename + ".log"
    traj_filename = output_filename + ".traj"
    print("log_filename = ",log_filename)
    print("traj_filename = ",traj_filename)

    atoms = atoms_in.copy()
    atoms.calc = calculator

    # Set the momenta corresponding to T=300K
    MaxwellBoltzmannDistribution(atoms, temperature_K=temperature,force_temp=True)
    Stationary(atoms)

    dyn = NPT(atoms,
          time_step*units.fs,
          temperature_K = temperature,
          externalstress = sigma*units.bar,
          ttime = ttime*units.fs,
          pfactor = pfactor*units.GPa*(units.fs**2),
          logfile = log_filename,
          trajectory = traj_filename,
          loginterval=num_interval
          )

    print_interval = 1000 if calc_type == "EMT" else num_interval
    dyn.attach(print_dyn, interval=print_interval)
    dyn.attach(MDLogger(dyn, atoms, log_filename, 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_in =  Atoms(symbols='Cu108', pbc=True, cell=[10.83, 10.83, 10.83])
i,temperature =  0 200
sigma = 1.0e+00 bar
ttime = 20.000 fs
pfactor = 2000000.000 GPa*fs^2
log_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0200K.log
traj_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0200K.traj
    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)
  1000   3.978    166.23    0.80  0.81  0.76  0.83  -0.01  0.27  0.04    1.133
  2000   4.698    164.97    0.67  0.61  0.63  0.77  -0.09  0.19  0.04    2.063
  3000   5.725    201.32    0.13  -0.12  0.25  0.28  0.08  -0.07  -0.02    3.005
  4000   6.272    209.89    -0.60  -0.85  -0.19  -0.74  0.09  -0.14  -0.14    3.898
  5000   4.987    194.84    -0.32  -0.15  -0.25  -0.55  0.07  0.02  0.06    4.764
  6000   5.138    194.18    -0.03  0.16  -0.35  0.08  0.18  -0.01  0.43    5.630
  7000   6.007    192.54    0.12  0.35  -0.45  0.45  0.54  -0.11  0.41    6.510
  8000   5.714    213.34    0.58  0.50  0.45  0.78  0.49  -0.06  0.25    7.411
  9000   4.851    187.68    0.19  0.12  0.35  0.09  0.11  -0.44  0.06    8.289
  10000   5.044    185.88    -0.53  -0.64  -0.18  -0.77  0.11  -0.61  -0.23    9.174
  11000   6.226    212.41    -0.88  -0.97  -0.86  -0.82  -0.10  -0.16  -0.12    10.012
  12000   5.503    187.05    0.19  0.22  0.03  0.33  0.16  -0.49  0.11    10.836
  13000   5.462    217.03    0.78  0.72  0.77  0.84  0.09  -0.43  0.15    11.658
  14000   5.579    195.67    0.18  -0.09  0.43  0.20  0.19  -0.51  -0.08    12.468
  15000   5.650    214.89    -0.43  -0.54  -0.08  -0.66  0.22  -0.47  0.17    13.278
  16000   5.339    170.60    -0.73  -0.83  -0.65  -0.70  0.28  -0.49  -0.01    14.089
  17000   4.925    187.46    0.20  0.18  -0.04  0.45  -0.15  -0.30  -0.06    14.897
  18000   4.952    175.70    0.72  0.63  0.45  1.06  -0.17  0.01  -0.04    15.708
  19000   4.997    176.19    0.58  0.54  0.63  0.57  -0.14  -0.38  -0.04    16.517
  20000   4.938    185.59    0.02  0.01  0.18  -0.14  -0.38  -0.44  0.34    17.355
i,temperature =  1 300
sigma = 1.0e+00 bar
ttime = 20.000 fs
pfactor = 2000000.000 GPa*fs^2
log_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0300K.log
traj_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0300K.traj
    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)
  1000   7.654    299.88    -0.07  0.33  -0.11  -0.43  -0.16  0.12  -0.51    0.832
  2000   9.322    332.30    0.81  1.18  0.62  0.64  -0.09  0.02  0.14    1.652
  3000   7.319    249.97    2.00  2.24  1.63  2.14  -0.04  -0.13  -0.12    2.472
  4000   7.429    277.97    1.04  1.16  0.74  1.23  -0.05  -0.13  0.04    3.299
  5000   7.835    316.40    -0.59  -0.60  -0.77  -0.42  0.44  0.01  0.48    4.118
  6000   7.766    289.53    -1.38  -1.40  -1.20  -1.54  0.39  -0.21  0.73    4.932
  7000   9.249    312.54    -1.00  -0.94  -0.79  -1.26  0.11  0.33  0.73    5.749
  8000   9.578    325.69    0.45  0.21  0.60  0.54  0.14  -0.09  0.54    6.571
  9000   9.292    293.75    1.21  0.79  1.50  1.33  -0.13  -0.71  0.76    7.387
  10000   8.428    301.90    1.00  0.84  1.15  1.01  0.03  -0.07  0.75    8.232
  11000   7.516    282.29    -0.30  0.02  -0.47  -0.45  -0.16  0.12  0.45    9.052
  12000   8.321    304.88    -1.49  -1.19  -1.47  -1.80  -0.05  -0.50  0.26    9.872
  13000   8.785    295.99    -1.28  -1.41  -1.31  -1.12  0.06  -0.26  0.18    10.691
  14000   7.991    272.72    0.55  0.01  0.72  0.91  -0.10  -0.33  0.21    11.508
  15000   8.030    281.03    1.53  1.22  1.58  1.78  -0.20  -0.10  0.09    12.327
  16000   8.779    318.99    0.96  0.98  0.94  0.94  -0.61  -0.36  0.04    13.153
  17000   7.856    280.24    -0.14  0.29  -0.33  -0.39  -0.77  -0.05  -0.35    14.034
  18000   8.726    356.50    -1.18  -1.16  -1.17  -1.21  -0.49  0.14  0.16    14.862
  19000   9.173    328.31    -1.50  -1.84  -1.40  -1.27  -0.60  0.01  0.14    15.749
  20000   7.327    267.16    0.47  -0.02  0.60  0.82  -0.23  -0.34  -0.13    16.770
i,temperature =  2 400
sigma = 1.0e+00 bar
ttime = 20.000 fs
pfactor = 2000000.000 GPa*fs^2
log_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0400K.log
traj_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0400K.traj
    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)
  1000   10.098    397.37    -1.39  -1.40  -1.59  -1.17  0.40  0.01  0.18    1.230
  2000   10.254    405.68    1.08  0.95  1.01  1.28  0.13  -0.22  0.04    2.132
  3000   10.624    380.96    3.05  2.86  3.07  3.20  0.40  0.22  0.31    3.247
  4000   11.342    354.85    3.15  3.24  3.27  2.94  0.19  0.50  0.19    4.358
  5000   10.732    357.70    1.49  1.68  1.32  1.48  0.29  0.08  -0.17    5.202
  6000   11.075    429.36    -1.11  -0.78  -1.23  -1.32  0.67  -0.12  0.19    6.050
  7000   10.846    382.45    -3.32  -3.71  -2.82  -3.44  0.59  0.49  0.21    6.884
  8000   10.633    392.35    -2.98  -3.51  -2.75  -2.67  0.54  0.10  0.52    7.717
  9000   10.706    436.02    -0.44  -0.67  -0.28  -0.36  0.26  0.01  0.35    8.542
  10000   10.366    363.55    2.09  2.27  1.74  2.26  0.08  -0.18  0.76    9.368
  11000   11.086    365.36    3.19  3.37  3.03  3.18  0.21  -0.38  0.18    10.199
  12000   11.349    384.52    2.52  2.57  2.77  2.22  0.17  -0.34  -0.12    11.028
  13000   11.061    369.82    -0.09  -0.59  0.15  0.18  -0.30  -0.66  0.42    11.850
  14000   11.056    423.24    -2.49  -2.83  -2.21  -2.45  0.13  -0.56  0.62    12.674
  15000   11.261    383.06    -3.35  -3.01  -3.48  -3.56  0.38  -0.49  0.33    13.494
  16000   10.551    376.12    -1.44  -1.12  -1.75  -1.45  0.32  -0.16  0.65    14.317
  17000   11.942    434.60    0.54  0.55  0.57  0.52  -0.12  -0.28  -0.09    15.137
  18000   12.471    435.71    2.31  2.00  2.75  2.19  0.45  -0.53  -0.36    15.985
  19000   12.968    409.56    1.98  1.67  2.15  2.14  -0.44  -0.24  -0.01    16.822
  20000   10.564    354.54    1.47  1.72  1.24  1.46  -0.39  -0.45  0.02    17.654
i,temperature =  3 500
sigma = 1.0e+00 bar
ttime = 20.000 fs
pfactor = 2000000.000 GPa*fs^2
log_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0500K.log
traj_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0500K.traj
    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)
  1000   14.673    520.26    -3.31  -3.38  -3.53  -3.03  -0.11  0.31  0.39    0.824
  2000   14.639    536.98    -1.21  -1.35  -1.42  -0.86  -0.04  -0.14  0.01    1.753
  3000   13.517    458.72    1.10  0.64  1.27  1.39  0.10  -0.17  0.04    2.687
  4000   14.045    478.80    2.65  2.58  2.85  2.52  0.41  0.05  0.00    3.598
  5000   15.437    484.76    2.34  2.45  2.13  2.45  -0.24  0.16  -0.48    4.469
  6000   15.525    516.64    1.79  2.26  1.21  1.91  0.39  -0.09  0.41    5.300
  7000   13.855    471.34    0.77  1.16  0.31  0.84  0.72  0.16  -0.58    6.141
  8000   13.553    484.43    -1.06  -1.28  -1.09  -0.81  0.86  0.70  -0.07    6.970
  9000   14.460    530.01    -2.49  -2.74  -2.32  -2.42  0.69  0.89  -0.44    7.797
  10000   14.742    534.49    -2.72  -2.71  -2.65  -2.81  1.05  0.24  -0.85    8.624
  11000   13.574    465.51    -1.36  -1.33  -1.50  -1.26  0.35  0.33  -0.30    9.451
  12000   14.380    504.15    0.38  0.17  0.60  0.36  0.46  0.49  -0.72    10.281
  13000   14.959    464.34    1.59  1.41  1.75  1.60  0.00  0.42  -0.23    11.106
  14000   13.694    435.18    2.74  3.18  2.46  2.56  -0.14  0.67  -0.29    11.948
  15000   14.828    468.54    1.60  1.77  1.33  1.72  -0.38  0.57  -0.43    12.771
  16000   14.646    498.14    0.30  0.43  0.42  0.05  0.32  0.14  -0.47    13.649
  17000   13.559    495.21    -0.86  -0.78  -0.58  -1.20  0.33  0.75  -0.63    14.475
  18000   14.272    518.76    -2.41  -2.63  -2.26  -2.35  -0.06  0.67  -0.53    15.338
  19000   12.809    504.62    -1.59  -1.63  -1.70  -1.45  -0.58  1.07  -0.71    16.204
  20000   14.694    505.54    -1.32  -1.44  -1.41  -1.10  -0.14  1.15  0.53    17.087
i,temperature =  4 600
sigma = 1.0e+00 bar
ttime = 20.000 fs
pfactor = 2000000.000 GPa*fs^2
log_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0600K.log
traj_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0600K.traj
    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)
  1000   14.309    488.11    -3.44  -3.10  -4.03  -3.19  -0.34  -0.01  -0.08    0.927
  2000   15.377    533.39    -3.38  -3.23  -3.69  -3.22  -0.11  0.38  -0.01    1.764
  3000   14.439    503.41    -1.09  -1.79  -0.48  -1.02  0.18  0.90  0.26    2.596
  4000   17.430    671.39    0.15  0.32  0.15  -0.02  -0.36  -0.20  0.44    3.425
  5000   15.922    538.48    2.68  2.79  2.95  2.30  0.18  0.39  -0.05    4.251
  6000   16.884    523.87    3.28  3.55  3.29  2.99  0.15  0.45  -0.19    5.078
  7000   17.547    560.94    3.75  3.28  3.94  4.03  0.35  0.51  0.09    5.905
  8000   19.650    613.92    2.91  2.38  3.60  2.73  0.24  0.97  -0.31    6.727
  9000   17.656    558.54    2.59  2.58  2.40  2.78  1.23  1.02  -0.19    7.577
  10000   17.574    631.69    1.78  2.24  0.86  2.22  0.38  1.09  0.27    8.410
  11000   17.762    636.11    -0.50  -0.71  -0.65  -0.15  0.30  1.19  0.27    9.253
  12000   13.458    461.03    -1.11  -1.54  -0.12  -1.67  0.20  1.09  0.47    10.191
  13000   14.940    555.56    -1.93  -1.88  -1.31  -2.59  -0.80  1.23  0.49    11.024
  14000   15.085    534.58    -2.66  -2.55  -3.09  -2.34  -0.57  1.17  0.58    11.866
  15000   16.737    591.48    -3.02  -2.98  -3.28  -2.80  -0.62  1.10  0.62    12.738
  16000   16.577    595.75    -1.98  -1.66  -1.85  -2.42  -0.14  0.34  0.62    13.611
  17000   17.616    640.28    -0.82  -0.38  -0.33  -1.74  0.49  0.66  0.44    14.444
  18000   15.836    530.52    0.90  1.28  0.55  0.86  -0.09  0.30  -0.37    15.278
  19000   16.765    574.35    2.02  1.37  1.98  2.72  -0.35  0.62  0.32    16.114
  20000   17.810    631.17    2.75  2.65  2.65  2.94  0.13  0.36  -0.19    16.981
i,temperature =  5 700
sigma = 1.0e+00 bar
ttime = 20.000 fs
pfactor = 2000000.000 GPa*fs^2
log_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0700K.log
traj_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0700K.traj
    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)
  1000   20.594    818.34    -7.78  -7.91  -7.12  -8.30  -0.17  -0.88  1.03    0.869
  2000   18.389    640.07    -7.07  -7.01  -7.39  -6.80  -0.19  -1.83  1.51    1.709
  3000   18.585    717.27    -6.12  -6.24  -6.68  -5.43  -0.54  -1.84  0.99    2.546
  4000   18.007    598.71    -6.10  -6.54  -5.32  -6.44  -0.23  -1.68  1.43    3.380
  5000   19.528    697.28    -5.77  -6.63  -5.48  -5.19  -0.24  -1.36  0.96    4.221
  6000   19.696    727.41    -5.21  -5.57  -5.03  -5.02  -0.49  -1.55  0.33    5.055
  7000   21.996    754.22    -5.51  -5.01  -5.78  -5.75  0.26  -1.33  0.52    6.081
  8000   20.226    752.49    -3.93  -3.80  -4.07  -3.92  -0.78  -0.37  -0.32    6.965
  9000   18.830    660.58    -2.78  -2.97  -2.73  -2.63  -0.77  0.39  0.26    7.896
  10000   20.346    803.03    -0.99  -0.97  -1.24  -0.77  -0.93  0.02  0.18    8.749
  11000   19.890    763.27    -1.02  -0.08  -1.38  -1.60  -1.22  -0.42  1.06    9.595
  12000   18.734    664.62    0.07  0.88  0.04  -0.73  -1.97  -0.02  0.78    10.430
  13000   19.332    621.22    0.40  0.20  1.26  -0.27  -1.37  -0.82  0.14    11.315
  14000   20.502    657.99    1.07  0.34  0.92  1.96  -1.80  -0.38  -0.47    12.160
  15000   20.546    633.89    1.32  1.56  1.12  1.28  -1.03  -1.55  -0.50    13.000
  16000   21.013    661.66    2.10  2.82  1.65  1.82  -1.78  -0.52  -0.29    13.839
  17000   22.299    711.83    2.18  2.66  2.39  1.51  -1.58  -1.63  0.06    14.751
  18000   23.074    760.55    1.93  1.80  2.15  1.84  -0.76  -1.42  0.95    15.653
  19000   22.322    730.21    1.86  2.19  1.40  2.00  -0.66  -0.96  -0.12    16.504
  20000   21.319    685.32    1.59  1.71  1.81  1.26  -0.85  -1.59  -0.00    17.343
i,temperature =  6 800
sigma = 1.0e+00 bar
ttime = 20.000 fs
pfactor = 2000000.000 GPa*fs^2
log_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0800K.log
traj_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0800K.traj
    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)
  1000   19.615    699.02    -7.17  -7.17  -7.03  -7.30  0.45  0.60  1.00    0.838
  2000   22.314    737.39    -7.19  -6.79  -7.49  -7.29  -0.19  0.42  1.92    1.670
  3000   22.217    905.76    -5.00  -5.64  -4.91  -4.44  -0.83  1.62  2.00    2.499
  4000   21.086    766.84    -3.02  -2.85  -2.47  -3.76  -0.50  0.63  2.62    3.329
  5000   19.333    691.99    -1.10  -0.48  -1.21  -1.60  -1.02  1.20  2.26    4.160
  6000   22.632    769.16    -0.51  -1.08  -0.34  -0.10  -0.29  1.13  1.81    5.005
  7000   22.952    749.84    0.28  0.31  -0.18  0.71  -0.33  0.67  1.51    5.835
  8000   22.408    679.48    0.24  -0.03  0.45  0.29  -0.44  -0.01  1.50    6.693
  9000   25.025    887.54    -0.93  -1.40  -0.56  -0.82  -0.13  0.31  1.59    7.587
  10000   23.422    848.35    0.64  0.89  0.68  0.36  0.86  0.94  0.33    8.528
  11000   26.432    840.75    0.59  0.19  0.09  1.50  0.92  0.73  -1.06    9.393
  12000   23.960    769.43    1.35  0.90  1.32  1.84  1.39  0.82  -0.92    10.238
  13000   23.571    869.00    1.94  1.60  2.73  1.50  1.45  0.25  0.01    11.074
  14000   22.455    755.61    1.09  1.35  0.76  1.17  0.82  -0.78  -0.08    11.911
  15000   24.966    841.11    0.45  0.46  0.02  0.88  -0.05  -0.92  -0.30    12.793
  16000   23.094    776.38    1.10  0.93  1.01  1.36  0.61  -0.83  -0.08    13.680
  17000   25.367    860.67    0.45  1.32  -0.10  0.13  -0.55  -0.23  0.19    14.569
  18000   24.627    821.41    1.38  1.82  1.07  1.24  -1.23  -0.52  0.40    15.541
  19000   23.767    761.80    1.47  1.25  2.39  0.77  0.13  -0.34  1.33    16.403
  20000   24.721    787.23    1.39  0.82  1.37  1.99  0.09  0.48  0.51    17.277
i,temperature =  7 900
sigma = 1.0e+00 bar
ttime = 20.000 fs
pfactor = 2000000.000 GPa*fs^2
log_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0900K.log
traj_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0900K.traj
    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)
  1000   23.834    788.37    -7.03  -6.71  -7.68  -6.70  -1.12  -1.30  0.58    0.857
  2000   26.979    841.52    -0.89  -0.63  -0.05  -1.98  -2.60  -1.04  0.28    1.708
  3000   31.115    1010.60    0.63  -0.16  1.30  0.76  -2.05  -1.11  0.44    2.555
  4000   30.643    915.34    0.85  0.83  0.33  1.39  -1.24  -1.05  -0.06    3.400
  5000   29.591    901.35    1.41  1.80  1.77  0.66  -0.50  -0.96  -0.77    4.329
  6000   27.527    1016.01    2.52  2.98  2.27  2.32  -0.35  -0.64  -0.24    5.236
  7000   28.071    781.00    1.04  1.12  -0.68  2.68  -0.18  0.70  0.68    6.089
  8000   28.884    959.52    2.75  2.16  2.90  3.18  1.21  1.19  0.80    6.966
  9000   28.956    965.19    1.74  2.36  2.11  0.73  0.42  0.93  0.73    7.818
  10000   26.964    861.06    1.83  2.50  1.01  1.99  0.33  1.58  0.91    8.665
  11000   28.701    938.36    0.05  -0.23  -0.47  0.84  0.40  2.28  -0.96    9.514
  12000   25.423    881.77    -0.39  -1.31  1.33  -1.19  -0.86  0.08  -0.66    10.365
  13000   26.335    872.98    -1.40  -0.45  -1.51  -2.25  -0.75  -0.35  -0.18    11.210
  14000   26.152    860.34    -0.32  -0.27  -1.10  0.40  -1.12  -0.53  -1.35    12.049
  15000   25.269    843.19    -0.70  -0.51  -0.51  -1.08  -0.71  -1.39  -0.26    12.887
  16000   25.793    881.28    -2.06  -1.22  -2.18  -2.79  0.16  -0.10  -0.66    13.724
  17000   24.458    860.41    0.58  0.25  0.37  1.11  0.85  0.41  0.13    14.565
  18000   27.555    958.01    0.51  0.37  0.55  0.62  0.07  0.03  0.36    15.405
  19000   26.487    922.41    0.65  1.34  0.28  0.32  0.63  0.11  0.67    16.242
  20000   28.488    960.88    1.55  2.10  1.33  1.22  -0.17  -0.58  0.39    17.079
i,temperature =  8 1000
sigma = 1.0e+00 bar
ttime = 20.000 fs
pfactor = 2000000.000 GPa*fs^2
log_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_1000K.log
traj_filename =  ./output/ch6/fcc-Cu_3x3x3_NPT_EMT_1000K.traj
    imd     Etot(eV)    T(K)    stress(mean,xx,yy,zz,yz,xz,xy)(GPa)  elapsed_time(sec)
  1000   29.100    1014.78    -2.97  -3.23  -2.36  -3.33  -0.16  -1.28  1.33    0.834
  2000   31.246    909.61    2.51  2.59  2.62  2.30  -1.70  -1.09  0.62    1.666
  3000   34.402    1067.05    3.76  4.51  3.11  3.66  -0.66  -0.25  0.16    2.510
  4000   30.820    979.53    2.29  1.83  2.25  2.80  0.12  0.32  -1.23    3.348
  5000   29.635    930.94    -0.98  -1.81  -0.37  -0.78  0.71  1.13  -0.77    4.195
  6000   28.940    932.28    -3.34  -3.73  -3.81  -2.47  1.75  2.05  -0.25    5.026
  7000   28.848    1049.14    -3.11  -2.61  -3.32  -3.40  1.30  0.51  -0.26    5.861
  8000   26.307    1015.33    -0.34  -0.90  0.64  -0.76  0.31  0.33  -0.50    6.697
  9000   29.778    1009.50    -0.18  0.07  -1.00  0.40  -0.01  -0.30  0.50    7.530
  10000   28.630    941.93    1.81  2.49  1.19  1.75  0.57  -0.08  1.06    8.362
  11000   29.936    1045.13    2.83  2.71  2.02  3.75  0.45  -0.36  0.79    9.224
  12000   28.866    958.14    2.02  1.38  2.35  2.32  0.77  0.58  1.12    10.061
  13000   28.790    993.03    -0.65  -0.76  -0.44  -0.75  -0.24  0.65  -0.15    10.895
  14000   30.622    1030.82    -0.46  1.57  -1.53  -1.42  -0.07  0.71  0.29    11.730
  15000   29.641    1020.87    -0.54  -1.40  0.12  -0.34  -0.80  -0.08  -0.12    12.572
  16000   30.512    999.98    0.61  -0.51  0.64  1.71  -0.16  -1.61  -0.52    13.431
  17000   31.112    994.90    0.81  1.41  0.32  0.71  0.49  -1.65  -0.79    14.266
  18000   30.148    973.65    0.01  1.02  0.15  -1.14  2.06  -1.40  0.54    15.118
  19000   31.246    1096.26    -1.39  -3.01  -1.26  0.10  0.68  -0.28  2.50    16.034
  20000   29.739    1111.59    -0.95  -1.34  -1.12  -0.37  2.73  1.02  1.12    16.936

A visualization of the results is shown below.

Here, the read method has been set to index="::100" to visualize the results of interval thinning from the Trajectory file. Although Cu retains its solid structure and is not moving significantly, we can see that there is a change in the cell size due to the vibrational movement.

[2]:

from ase.io import Trajectory, read
from pfcc_extras.visualize.povray import traj_to_apng
from IPython.display import Image


traj = read("./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0900K.traj", index="::100")
traj_to_apng(traj, f"output/ch6/Fig6-3_fcc-Cu_3x3x3_NPT_EMT_0900K.png", rotation="10x,-10y,0z", clean=True, n_jobs=16)

Image(url="output/ch6/Fig6-3_fcc-Cu_3x3x3_NPT_EMT_0900K.png")

[Parallel(n_jobs=16)]: Using backend ThreadingBackend with 16 concurrent workers.
[Parallel(n_jobs=16)]: Done  10 out of  20 | elapsed:    8.9s remaining:    8.9s
[Parallel(n_jobs=16)]: Done  20 out of  20 | elapsed:   10.9s finished

[2]:

[3]:

from ase.io import Trajectory, read
from pfcc_extras.visualize.view import view_ngl


traj = read("./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0900K.traj", index="::100")
view_ngl(traj, replace_structure=True)

The time constant \(\tau_t\) of the heat bath is 20 fs and the barostat parameter, pfactor, is 2e6 GPa\(\cdot\)fs\(^2\). Let’s check how the cell volume changes with time during a 20 ps NPT-MD simulation at 300 K to achieve thermal equilibrium under the above calculation conditions. You can analyze the traj file, which is the result of MD simulation at a specific temperature, with the following code.

[4]:

import matplotlib.pyplot as plt
from pathlib import Path
from ase.io import read,Trajectory

time_step = 0.01  # Time step size in ps between each snapshots recorded in traj
traj = Trajectory("./output/ch6/fcc-Cu_3x3x3_NPT_EMT_0300K.traj")

time = [ i*time_step for i in range(len(traj)) ]
volume = [ atoms.get_volume() for atoms in traj ]

# Create graph
fig = plt.figure(figsize=(8,3))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('time (ps)')  # x axis label
ax.set_ylabel('Cell volume (Å^3)')  # y axis label
ax.plot(time,volume, alpha=0.5)
ax.set_ylim([1100,1400])
#plt.savefig("filename.png")  # Set filename to be saved
plt.show()
_images/6_3_md-npt_8_0.png

The time evolution of the cell volume for a 20 ps NPT-MD simulation at 300 K to achieve thermal equilibrium under these calculation conditions is shown below.

6a8fa35e8f2d4d1ab17186671bfa9bc9

Fig.6-3a. Time evolution of cell volume. (fcc-Cu_3x3x3, @300 K and 1 bar)

It is shown that the cell volume oscillates within a range of about 4% with the NPT ensemble. Since we are dealing with cubic crystals, we can calculate the lattice constant of the crystal from the geometric mean of this volume, and by plotting the average normalized lattice constant at each temperature from 200 K to 1000 K, we obtain the following results. (Since we want to calculate the lattice constant after equilibrium is reached, we use only the second half of the Trajectory in the np.mean(vol[int(len(vol)/2):])**(1/3) part to calculate the lattice constant from the volume.)

[5]:

import matplotlib.pyplot as plt
import numpy as np
from pathlib import Path
from ase.io import read,Trajectory

time_step = 10.0  # Time step size between each snapshots recorded in traj
paths = Path("./output/ch6/").glob(f"**/fcc-Cu_3x3x3_NPT_{calc_type}_*K.traj")
path_list = sorted([ p for p in paths ])

# Temperature list extracted from the filename
temperature = [ float(p.stem.split("_")[-1].replace("K","")) for p in path_list ]
print("temperature = ",temperature)

# Compute lattice parameter
lat_a = []
for path in path_list:
    print(f"path = {path}")
    traj = Trajectory(path)
    vol = [ atoms.get_volume() for atoms in traj ]
    lat_a.append(np.mean(vol[int(len(vol)/2):])**(1/3))

print("lat_a = ",lat_a)

# Normalize relative to the value at 300 K
norm_lat_a = lat_a/lat_a[1]
print("norm_lat_a = ",norm_lat_a)

# Plot
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('Temperature (K)')  # x axis label
ax.set_ylabel('Normalized lattice parameter')  # y axis label
ax.scatter(temperature[:len(norm_lat_a)],norm_lat_a, alpha=0.5,label=calc_type.lower())
ax.legend(loc="upper left")

temperature =  [200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1000.0]
path = output/ch6/fcc-Cu_3x3x3_NPT_EMT_0200K.traj
path = output/ch6/fcc-Cu_3x3x3_NPT_EMT_0300K.traj
path = output/ch6/fcc-Cu_3x3x3_NPT_EMT_0400K.traj
path = output/ch6/fcc-Cu_3x3x3_NPT_EMT_0500K.traj
path = output/ch6/fcc-Cu_3x3x3_NPT_EMT_0600K.traj
path = output/ch6/fcc-Cu_3x3x3_NPT_EMT_0700K.traj
path = output/ch6/fcc-Cu_3x3x3_NPT_EMT_0800K.traj
path = output/ch6/fcc-Cu_3x3x3_NPT_EMT_0900K.traj
path = output/ch6/fcc-Cu_3x3x3_NPT_EMT_1000K.traj
lat_a =  [10.821572832083033, 10.842619005347451, 10.867506249620094, 10.88752168215781, 10.910928969034376, 10.938882785804418, 10.96019451392851, 10.990325464402156, 11.01840205198354]
norm_lat_a =  [0.99805894 1.         1.00229532 1.00414131 1.00630014 1.00887828
 1.01084383 1.01362277 1.01621223]

[5]:

<matplotlib.legend.Legend at 0x7f4475de0e10>
_images/6_3_md-npt_12_2.png

The following plot is obtained by running the above code.

7f79d08df15542c3ae500f690d4ff1b2

Fig.6-3b. Normalized lattice parameter as a function of temperature.

Experimental data is taken from Reference [4].

The results are further compared with the PFP calculations and, with the experimental data for reference. There is very good agreement for both the EMT and PFP data, while the PFP data appears to be closer to the experimental data with a small margin.

The coefficient of linear thermal expansion (CTE, \(\alpha\)) is calculated from the temperature dependence of the normalized lattice parameter using the following equation

\[\frac{a(T)}{a_{RT}} = 1 + \alpha \cdot T\]

where \(\alpha\) is the coefficient of linear thermal expansion and \(a(T)\) and \(a_{RT}\) are the lattice constants at temperature \(T\) and room temperature. The following is a summary of ASAP3-EMT, PFP, and experimental values [4].

\(\alpha\) ( \(10^{-5}\) /K)

emt

2.23

pfp

2.13

exp

1.74

Considering the calculation error, both calculation results are in reasonable agreement with the experimental values.

The coefficient of thermal expansion can be calculated by using MD simulations of the NPT ensemble in this way.

You may be thinking, “Can’t we reproduce thermal expansion of liquids and gases in the same way?” In principle, this is true, but at present there are only a limited number of models that can reproduce the thermal expansion of liquids and gases. There may be models that are specialized for specific liquids and gases in the classical force field, but when using first-principles calculations, the accuracy of intermolecular interactions in liquids and gases, which are much smaller than those in solids, is critical, and quantum chemical calculations with such accuracy are currently very difficult to perform. In particular, a large amount of computational data is required to create NNPs, and it is currently not very realistic to perform a large number of high-precision quantum chemical calculations. Therefore, the development of models with high accuracy in this area is expected in the future.

[Advanced] Parrinello-Rahman barostat parameter dependency

When performing MD simulations in the NPT ensemble, a parameter called pfactor needs to be chosen. It was explained that the appropriate value range is around 10\(^6\) GPa\(\cdot\)fs\(^2\) to 10\(^7\) GPa\(\cdot\)fs\(^2\). For your reference, the results when changing the value of pfactor are shown below. The calculation targets are the same fcc-Cu 3x3x3 unit cells as above, and the temperature is 300 K.

11ce7669c83f4b44a58da3b8082923fa

Fig.6-3c. Time evolution of cell volume as a function of pfactor.

When pfactor is small, the cell volume is oscillating at high speeds which is not desirable, and some regions have a mixture of high and low frequency oscillations making the system unstable. As pfactor increases, the period of oscillation gradually increases, and the large volume changes at the beginning of the calculation are hardly noticeable. Although there is no clear cut indicator of which value of pfactor is correct, the median value seems to be as stable as ever, although small oscillations can be observed for values of 10\(^6\) Ga\(\cdot\)fsec\(^2\) or higher. Similarly, there is no well-defined reference for the upper limit, and since the larger the pfactor, the longer the oscillation period becomes and the more difficult it is to handle. A range between 10\(^6\) and 10\(^7\) seems appropriate for the above example.

[Advanced] Berendsen barostat parameter dependency

Here, we explain the calculation method using Berendsen barostat. The time evolution of pressure is calculated according to the following equation. (For details of the derivation, see reference [5]).

\[\frac{d\mathbf{P}}{dt}=\frac{\mathbf{P}_o-\mathbf{P}}{\tau_P}\]

It is clear from the above equation that the Berendsen pressure control method exponentially increases the pressure at each time in the system closer to the external pressure \(\mathbf{P}_o\). Its speed is controlled by the time constant \(\tau_P\).

At each MD step, the coordinates and cell vectors of each atom are scaled by a factor expressed by the following equation

\[\mu_{ij}=\delta_{ij}-\frac{\Delta t}{3\tau_P}\beta_{ij}\big\{ P_{oij} - P_{ij}(t) \big\}\]

Just as \(\tau_T\) was a time constant when controlling a heat bath, an appropriate time constant \({\tau_P}\) must be set for the pressure control method. Let’s look at an actual calculation example.

Using the NPTBerendsen class, the object defining the dynamics is written in the following form

[6]:

from ase.md.nptberendsen import NPTBerendsen

dyn = NPTBerendsen(
    atoms,
    time_step*units.fs,
    temperature_K = temperature,
    pressure_au = 1.0 * units.bar,
    taut = 5.0 * units.fs,
    taup = 500.0 * units.fs,
    compressibility_au = 5e-7 / units.bar,
    logfile = log_filename,
    trajectory = traj_filename,
    loginterval=num_interval
)

Inhomogeneous_NPTBerendsen, which handles anisotropic pressure control, can be used with the same setting. (Simply relacing the class name from NPTBerendsen to Inhomogeneous_NPTBerendsen should work.) One difference between the two classes is that Inhomogeneous_NPTBerendsen accept mask argument can be set, i.e. setting mask=(1, 1, 1) allows the cell paramters a, b, and c can vary independently. If the element of this tuple is set to 0, the cell is fixed in that direction.

In addition to the calculation conditions such as temperature and pressure, it is necessary to set parameters called barostat time constant (taup, \(\tau_P\)) and compression ratio (compressibility, \(\beta_T\)). The dependence of the time variation of the cell volume on each of these values will be checked using the example of fcc-Cu at 300 K. We will start with \(\tau_P\).

c6e39e3724b64fa29b23350675245575

Fig.6-3d. Time evolution of cell volume as a function of \(\tau_P\).

The smaller the time constant, the more unstable and violent the period of oscillation is. On the other hand, if the time constant is too large, the change is too gradual and it takes time to reach equilibrium. This is the same concept as the time constant of the heat bath method. For a system mechanically similar to fcc-Cu, the appropriate value for \(\tau_P\) is around 10\(^2\) fs to 10\(^3\) fs due to the stability and convergence of the volume change.

Next is the dependence on \(\beta_T\). The results are as follows.

85d7190aecb3475b8527a97b3cdab355

Fig.6-3e. Time evolution of cell volume as a function of \(\beta_T\).

The smaller \(\beta_T\) is, the slower the convergence is, and the higher the value, the more unstable it becomes. From the trend of the graph, 10\(^{-7}\) to 10\(^{-6}\)fs seems to be a good range for \(\beta_T\).

There is no strict guideline for choosing the values of \(\tau_P\) and \(\beta_T\). But it is better to have the sense of the trend when these values are changed by an order of magnitude.

Finally, it should be noted that these values are only applicable to dense metallic materials such as fcc-Cu. If one wants to treat a completely different kind of materials (e.g., polymers, liquids, gases, etc.) with NPT ensemble, the appropriate range of values should be studied in advance. In general, it is recommended that great care be taken when working with new material systems, as omitting this preliminary study may lead to unintended results and extra time.

Reference

[1] 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#

[2] Melchionna S. (2000) “Constrained systems and statistical distribution”, Physical Review E 61 (6) 6165 https://journals.aps.org/pre/abstract/10.1103/PhysRevE.61.6165

[3] S. Melchionna, G. Ciccotti, B.L. Holian, “Hoover NPT dynamics for systems varying in shape and size”, Molecular Physics, (1993) 78 (3) 533 https://doi.org/10.1080/00268979300100371

[4] F.C. Nix, D. MacNair, “NIST:The Thermal Expansion of Pure Metals: Copper, Gold, Aluminum, Nickel, and Iron” https://materialsdata.nist.gov/handle/11256/32

[5] 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 https://aip.scitation.org/doi/10.1063/1.448118