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   4.409    175.46    0.51  0.51  0.22  0.79  0.18  0.14  -0.15    3.081
  2000   4.554    168.67    0.81  0.62  0.49  1.32  -0.05  0.11  -0.37    6.234
  3000   5.053    205.40    0.66  0.42  0.53  1.02  -0.17  -0.15  -0.19    8.929
  4000   4.447    185.29    0.45  0.43  0.59  0.34  0.04  -0.31  -0.11    12.269
  5000   4.475    177.11    -0.04  0.11  0.17  -0.41  0.19  0.05  -0.15    15.029
  6000   5.026    162.44    -0.18  -0.10  -0.10  -0.34  0.19  0.29  0.04    17.755
  7000   6.189    216.28    0.06  0.15  -0.17  0.19  -0.10  -0.10  0.12    21.232
  8000   5.431    183.85    0.44  0.41  0.37  0.54  -0.47  0.04  0.11    24.111
  9000   5.293    206.53    0.23  0.03  0.23  0.44  -0.50  0.22  0.24    27.106
  10000   5.408    198.72    -0.42  -0.52  -0.08  -0.64  -0.18  0.23  -0.15    29.854
  11000   5.527    199.20    -0.56  -0.44  -0.48  -0.76  -0.27  0.17  -0.23    32.392
  12000   5.505    215.48    0.17  0.68  -0.16  -0.01  -0.07  0.12  -0.34    35.132
  13000   4.902    186.62    0.82  1.07  0.47  0.92  -0.10  0.15  0.00    37.888
  14000   5.960    222.68    0.19  -0.03  0.20  0.41  0.08  0.07  0.02    41.139
  15000   5.288    178.66    -0.46  -0.94  -0.03  -0.41  0.36  -0.26  -0.04    44.053
  16000   5.017    206.58    -0.37  -0.52  -0.07  -0.51  0.28  -0.22  0.13    47.279
  17000   4.912    179.07    0.04  0.30  -0.09  -0.10  0.67  0.03  -0.11    50.617
  18000   5.376    200.65    0.68  1.21  0.29  0.52  0.60  0.07  -0.04    53.710
  19000   5.138    183.70    0.84  1.47  0.27  0.77  0.65  -0.03  0.04    57.105
  20000   5.343    190.34    0.01  0.07  -0.15  0.10  0.71  -0.19  0.04    60.639
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   6.948    290.55    0.14  -0.21  0.44  0.18  0.24  0.44  0.02    2.960
  2000   7.892    318.00    1.45  1.21  1.42  1.74  0.26  0.30  -0.46    5.941
  3000   9.439    371.59    1.60  1.55  1.01  2.23  0.37  0.55  -0.62    8.924
  4000   9.264    311.59    0.57  0.71  0.17  0.84  0.40  0.21  -0.62    12.344
  5000   9.115    288.89    -1.17  -1.01  -0.79  -1.71  0.25  0.47  -0.65    15.636
  6000   7.865    266.53    -1.53  -1.73  -0.65  -2.20  0.17  0.22  -0.69    18.228
  7000   8.026    318.54    -0.80  -0.96  -0.51  -0.93  0.32  0.47  -0.11    21.306
  8000   7.964    294.51    0.61  0.79  0.02  1.01  0.46  0.49  0.01    24.509
  9000   8.997    318.42    1.30  1.42  0.76  1.70  0.24  0.46  -0.07    27.221
  10000   8.712    305.76    1.00  1.10  0.70  1.19  0.08  0.62  -0.54    30.628
  11000   7.166    236.95    0.01  -0.03  0.49  -0.43  -0.75  0.12  -0.32    33.836
  12000   8.156    287.28    -1.48  -1.65  -1.30  -1.51  -0.63  0.06  0.03    36.733
  13000   8.132    284.67    -1.49  -1.58  -1.55  -1.32  -0.59  0.24  0.03    39.784
  14000   7.562    280.99    0.29  0.34  0.05  0.49  -0.49  0.15  -0.19    42.367
  15000   9.115    282.56    0.95  1.25  0.79  0.82  0.09  0.55  -0.20    45.349
  16000   8.767    297.23    1.42  1.49  1.58  1.20  -0.18  0.25  -0.08    48.828
  17000   8.571    324.21    0.38  0.40  0.64  0.09  -0.26  0.24  0.00    51.943
  18000   7.356    273.88    -0.93  -1.05  -1.05  -0.69  -0.65  0.02  -0.14    54.991
  19000   7.944    288.90    -1.59  -1.55  -2.02  -1.21  -0.40  0.12  -0.02    58.106
  20000   7.669    278.09    -0.19  -0.20  -0.46  0.07  -0.70  0.60  -0.02    61.495
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   9.901    380.54    -0.73  -0.88  -0.25  -1.06  0.41  -0.12  -0.01    3.310
  2000   10.953    424.68    0.90  0.70  0.91  1.07  0.34  -0.26  -0.04    5.950
  3000   11.837    412.66    1.90  2.16  1.42  2.14  -0.47  -0.03  0.37    9.051
  4000   12.274    439.53    1.86  2.02  1.84  1.73  -0.41  -0.42  0.27    12.127
  5000   10.531    380.62    1.04  1.05  1.03  1.03  -0.64  -0.40  0.08    15.048
  6000   10.598    379.05    -0.80  -1.25  -0.61  -0.54  -1.23  -0.08  0.07    17.199
  7000   10.303    414.59    -1.51  -1.22  -1.94  -1.38  -0.94  0.04  -0.28    19.426
  8000   11.592    415.85    -1.66  -1.27  -1.77  -1.94  -0.85  0.08  -0.10    22.226
  9000   10.502    380.82    0.41  0.71  0.53  -0.01  -0.62  -0.20  -0.17    25.079
  10000   11.888    399.96    1.09  1.04  1.50  0.73  -0.15  0.35  -0.04    27.839
  11000   10.849    362.57    1.82  1.51  1.90  2.03  0.13  -0.22  0.41    30.935
  12000   11.039    372.20    0.78  0.24  0.99  1.12  -0.11  -0.13  0.24    33.910
  13000   10.181    343.80    -0.48  -0.55  -0.34  -0.55  -0.11  -0.15  0.13    36.696
  14000   11.158    385.41    -1.91  -1.46  -1.82  -2.44  -0.58  0.04  0.37    39.738
  15000   11.076    378.35    -1.73  -1.31  -1.99  -1.90  0.02  0.27  0.45    42.240
  16000   11.301    417.99    -0.36  -0.40  -0.44  -0.24  -0.12  0.49  0.29    45.421
  17000   12.074    430.58    0.88  0.97  0.66  1.00  0.47  0.66  0.46    48.240
  18000   11.617    391.27    1.58  1.77  1.37  1.60  0.45  0.33  0.10    50.898
  19000   13.698    552.14    0.86  0.81  1.02  0.75  0.03  0.55  0.42    53.981
  20000   11.356    403.27    -0.28  -0.42  -0.11  -0.29  0.21  0.62  -0.19    57.308
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   12.937    509.57    -1.61  -1.42  -1.57  -1.85  -0.50  0.18  0.00    2.853
  2000   12.492    450.87    -0.12  -0.56  -0.23  0.42  0.07  0.10  -0.45    5.125
  3000   13.562    474.38    1.96  1.67  1.87  2.34  0.22  0.42  -0.35    7.974
  4000   15.870    481.72    2.16  1.91  2.00  2.58  0.28  0.06  -0.50    10.835
  5000   15.115    497.90    3.16  3.41  3.06  3.00  0.61  0.43  0.21    13.974
  6000   14.462    499.74    2.67  3.02  2.73  2.27  0.59  0.32  0.31    16.821
  7000   13.970    496.22    0.74  1.07  0.68  0.49  0.19  -0.10  0.26    19.647
  8000   14.336    480.33    -1.85  -2.14  -2.31  -1.11  -0.34  -0.81  -0.48    22.709
  9000   13.104    484.38    -2.65  -3.46  -2.59  -1.91  -0.07  -1.67  -0.26    25.913
  10000   15.766    567.79    -3.54  -3.66  -3.10  -3.85  0.68  -1.39  -1.18    29.078
  11000   13.868    512.80    -1.40  -0.86  -1.01  -2.34  0.42  -1.00  -0.03    32.009
  12000   14.059    500.51    0.57  0.88  1.02  -0.18  0.46  -0.81  -0.48    34.781
  13000   16.236    537.63    1.53  1.18  1.14  2.26  0.27  -0.79  0.48    37.879
  14000   15.715    531.19    2.89  2.46  2.30  3.90  0.50  -0.35  0.46    41.016
  15000   15.919    492.59    2.12  2.05  2.26  2.07  -0.17  -0.78  0.20    43.721
  16000   14.143    532.73    1.64  2.26  1.72  0.95  0.24  -0.50  0.39    47.002
  17000   14.338    487.14    -1.13  -1.23  -1.40  -0.77  0.36  -1.19  0.52    49.995
  18000   13.942    491.84    -2.39  -2.94  -2.41  -1.81  -0.43  -1.45  0.71    53.389
  19000   12.268    448.83    -2.47  -2.78  -2.07  -2.56  -0.15  -1.24  0.99    56.305
  20000   13.945    535.41    -1.69  -1.27  -1.58  -2.23  -0.11  -0.69  0.36    59.510
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   16.087    641.72    -3.71  -3.55  -4.01  -3.56  -0.94  -0.10  -0.10    3.070
  2000   16.481    603.98    -2.17  -1.77  -2.60  -2.14  -1.18  0.60  -1.15    5.762
  3000   17.627    604.66    -0.56  -0.95  -0.97  0.23  -0.82  0.45  -0.64    30.637
  4000   16.322    512.69    2.23  2.02  2.14  2.54  -0.14  0.55  -1.16    69.998
  5000   17.990    530.96    2.99  2.52  3.70  2.76  -0.04  0.58  -0.93    100.589
  6000   18.332    593.33    4.17  4.60  4.02  3.90  0.34  0.44  -0.64    103.463
  7000   20.593    703.42    3.37  3.25  3.51  3.36  0.46  0.38  -0.51    106.795
  8000   17.471    588.41    2.53  2.50  2.40  2.68  0.92  1.32  -1.43    110.182
  9000   16.230    550.43    0.87  0.38  0.95  1.28  1.41  1.33  -0.76    113.385
  10000   17.951    648.30    -1.83  -2.22  -1.43  -1.84  2.11  1.18  0.14    116.471
  11000   16.131    574.53    -2.30  -2.31  -2.72  -1.88  1.25  1.15  0.04    119.182
  12000   16.993    577.14    -3.42  -3.49  -3.64  -3.13  0.94  0.39  -0.22    121.976
  13000   15.844    565.26    -2.62  -2.79  -2.66  -2.40  -0.06  -0.97  0.75    124.385
  14000   16.340    600.13    -1.73  -2.17  -1.36  -1.65  -0.14  -1.12  -0.20    127.670
  15000   17.722    646.40    -0.89  -0.93  -0.71  -1.03  -0.90  -0.41  -0.19    130.387
  16000   17.251    559.00    0.32  0.59  -0.29  0.67  -0.35  -0.74  0.14    133.107
  17000   17.865    586.11    0.83  1.11  0.70  0.69  0.71  -0.96  0.04    135.098
  18000   18.496    591.65    1.10  0.43  1.65  1.23  0.26  -0.85  0.24    137.047
  19000   17.045    557.54    1.45  1.30  1.72  1.34  0.03  -0.24  0.35    138.984
  20000   17.930    649.16    1.01  1.35  1.06  0.62  -0.04  -0.85  -0.54    140.904
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   21.130    717.63    -8.10  -7.76  -9.00  -7.53  -0.50  0.03  -0.46    1.926
  2000   18.584    654.84    -5.86  -4.20  -7.07  -6.31  -0.36  -1.02  -0.08    3.894
  3000   18.151    672.98    -5.27  -5.15  -5.65  -4.99  -1.52  0.27  -0.06    5.894
  4000   18.098    681.02    -4.13  -5.21  -3.67  -3.51  -1.42  0.19  0.12    7.915
  5000   19.217    726.64    -3.66  -3.42  -3.52  -4.02  -2.07  1.15  -0.36    9.953
  6000   20.184    782.32    -2.69  -1.73  -3.07  -3.26  -1.38  0.01  -0.22    11.895
  7000   20.636    711.38    -2.05  -1.98  -2.44  -1.74  -1.91  -0.04  0.19    13.915
  8000   21.412    794.20    -0.88  -1.46  -0.99  -0.20  -1.21  1.14  0.07    15.915
  9000   21.412    692.04    0.45  -0.09  0.83  0.60  -1.29  1.36  1.36    17.828
  10000   20.974    715.60    2.27  1.49  3.24  2.07  -1.39  0.71  1.52    19.732
  11000   20.436    680.67    3.30  3.45  3.41  3.05  -0.80  1.51  0.68    21.754
  12000   21.224    615.08    2.82  2.44  2.54  3.49  -0.69  1.41  0.76    23.649
  13000   21.853    649.51    2.97  2.96  2.82  3.12  -0.77  1.32  0.18    25.699
  14000   22.078    691.09    3.29  2.91  3.57  3.40  -0.67  1.50  0.19    27.687
  15000   22.062    609.07    2.38  2.46  2.19  2.49  -1.17  0.73  -0.21    29.632
  16000   22.444    694.37    2.59  2.18  2.40  3.21  -1.01  1.74  -1.03    31.583
  17000   20.422    607.63    3.06  3.09  2.91  3.17  -0.59  1.32  -0.61    33.578
  18000   21.280    693.12    2.60  2.65  3.27  1.87  -0.72  1.05  -0.64    35.613
  19000   21.836    685.83    1.20  1.59  1.50  0.51  -0.16  -0.00  -0.87    37.611
  20000   20.453    706.17    1.69  1.60  0.97  2.50  -0.54  -0.04  -0.31    39.638
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   21.992    746.13    -8.04  -6.10  -8.79  -9.24  -0.01  -0.46  -1.04    2.031
  2000   19.312    683.29    -5.32  -5.73  -5.06  -5.17  0.87  -1.82  -1.32    3.942
  3000   20.427    723.90    -4.14  -5.04  -3.40  -3.98  1.59  -1.85  -2.34    6.007
  4000   24.735    918.51    -4.42  -5.42  -3.57  -4.27  1.64  -1.24  -2.25    8.020
  5000   23.588    809.15    -4.55  -3.35  -5.15  -5.14  1.16  -0.60  -1.53    9.936
  6000   21.399    741.90    -3.83  -3.59  -4.46  -3.43  0.61  0.01  -0.95    11.921
  7000   21.007    730.71    -3.13  -3.64  -2.55  -3.21  0.82  1.35  -0.48    13.907
  8000   21.089    696.61    -2.25  -3.55  -0.74  -2.47  1.49  1.80  0.13    15.887
  9000   23.073    852.01    -2.45  -2.75  -2.36  -2.24  0.51  1.46  0.72    17.806
  10000   21.750    731.26    -2.54  -1.78  -3.12  -2.71  1.41  0.52  1.06    19.723
  11000   23.415    868.18    -1.70  -0.61  -2.82  -1.67  2.38  0.37  0.76    21.641
  12000   22.295    728.36    -2.01  -2.36  -1.08  -2.59  1.93  -0.67  0.35    23.638
  13000   23.454    801.19    -2.08  -2.84  -1.60  -1.82  2.03  0.38  0.04    25.699
  14000   20.723    690.10    -0.98  -0.33  -1.35  -1.27  0.27  -0.96  -1.26    27.780
  15000   22.425    732.43    -1.06  -0.71  -1.44  -1.04  -0.77  0.08  -0.88    29.822
  16000   24.276    787.47    -0.19  -1.47  1.10  -0.19  -0.63  -0.06  0.24    31.709
  17000   24.079    761.50    -0.04  0.46  -0.80  0.23  -1.64  0.28  -0.68    33.677
  18000   23.593    717.72    0.63  2.66  -0.44  -0.33  -0.85  0.44  -0.09    35.814
  19000   23.733    717.47    1.07  0.35  2.40  0.46  -0.11  -0.16  -0.30    37.802
  20000   23.627    782.00    2.00  1.25  3.11  1.66  0.56  0.73  -0.04    39.794
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   24.053    856.20    -6.40  -6.13  -5.54  -7.54  -1.60  -1.35  0.64    1.963
  2000   29.027    833.89    0.01  -0.23  -0.49  0.77  -1.99  -0.56  0.11    3.886
  3000   29.311    832.15    1.97  2.12  0.89  2.90  -2.61  -0.93  0.76    5.972
  4000   30.075    981.67    2.49  2.61  2.97  1.90  -2.07  -1.12  0.30    8.011
  5000   29.670    989.71    1.63  0.89  2.92  1.08  -1.05  -0.49  -0.11    9.968
  6000   27.201    904.73    2.18  2.36  1.78  2.39  -0.28  -0.20  0.11    11.899
  7000   27.076    917.39    3.43  4.50  2.78  3.00  0.80  0.24  0.18    13.895
  8000   29.889    995.00    2.59  2.89  2.43  2.46  0.71  0.28  0.24    15.882
  9000   30.146    974.61    1.56  1.04  1.64  2.01  -0.48  0.92  -0.78    17.822
  10000   25.306    806.60    1.53  1.33  0.97  2.28  0.27  0.35  -1.43    19.831
  11000   26.427    813.28    0.39  1.61  -0.28  -0.16  0.69  -0.37  -1.40    21.923
  12000   25.481    950.54    0.28  -0.68  1.36  0.16  0.13  -0.80  -2.72    23.897
  13000   25.082    875.69    -1.90  -2.96  -0.29  -2.45  -0.31  -0.15  -1.23    25.888
  14000   25.217    926.64    -2.53  -1.26  -3.28  -3.04  -1.71  -0.83  -0.87    27.920
  15000   22.603    817.02    -0.83  -1.01  -0.65  -0.82  -2.14  -1.40  -0.49    30.011
  16000   28.101    918.84    -1.67  -2.51  -0.99  -1.51  -1.79  -2.30  -0.37    32.024
  17000   27.314    1031.26    -0.01  -0.25  0.51  -0.29  -1.49  -3.04  0.07    34.046
  18000   24.763    713.64    -0.86  -0.12  -2.05  -0.41  -1.61  0.18  0.10    36.032
  19000   27.750    904.59    -1.58  -1.35  -1.78  -1.61  -0.33  1.07  -0.01    37.931
  20000   27.337    893.66    -0.31  -1.50  0.23  0.36  1.00  -0.40  1.17    39.895
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   28.119    1011.01    -2.14  -1.45  -2.38  -2.59  -1.14  0.95  -0.38    2.036
  2000   33.702    1004.21    2.73  2.65  1.75  3.79  -0.86  0.01  -1.29    3.995
  3000   35.155    958.17    2.76  2.73  2.89  2.65  -0.50  0.46  -0.13    6.055
  4000   32.087    968.53    3.13  3.16  3.04  3.17  0.08  0.16  -0.13    8.055
  5000   29.413    931.13    1.13  0.94  0.64  1.80  0.52  -1.28  1.53    9.984
  6000   26.729    1009.75    0.28  0.99  -0.17  0.01  -0.41  -1.27  0.95    12.082
  7000   30.541    1027.19    -2.95  -2.34  -3.19  -3.31  1.24  -0.82  0.47    14.091
  8000   28.584    955.48    -1.39  -2.16  -1.24  -0.78  1.24  -1.43  0.37    16.028
  9000   30.511    1133.23    -1.47  -1.90  -1.53  -0.98  0.84  -0.83  -0.61    17.962
  10000   30.152    1092.36    0.25  0.39  0.35  0.03  -0.05  -1.58  0.10    19.996
  11000   28.550    969.00    1.32  1.63  1.39  0.94  0.48  -0.76  0.47    22.044
  12000   31.246    872.10    0.64  0.87  1.79  -0.75  -0.00  -0.60  0.44    24.042
  13000   31.001    1081.19    0.86  0.98  0.13  1.48  0.57  1.33  0.70    26.003
  14000   28.189    1040.45    -0.33  0.35  -1.06  -0.28  0.50  0.56  0.83    28.044
  15000   25.827    831.57    -0.21  0.01  0.98  -1.63  1.15  -1.57  0.27    29.965
  16000   31.611    1140.17    -0.48  -0.54  -1.45  0.54  0.90  -1.73  0.20    32.034
  17000   30.605    1003.21    0.89  0.19  1.46  1.02  -0.53  -0.41  0.07    34.082
  18000   30.230    1017.39    1.31  2.48  1.71  -0.28  -0.74  -0.64  -0.74    36.084
  19000   32.260    1140.83    -0.35  0.10  -1.10  -0.07  -1.00  -0.62  -0.84    38.118
  20000   25.343    859.61    -0.01  0.37  -0.64  0.25  -1.96  -0.99  0.62    40.049

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:   10.7s remaining:   10.7s
[Parallel(n_jobs=16)]: Done  20 out of  20 | elapsed:   13.7s 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.

1b1be1220daf44c9afc3a23a9ff51a2a

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.82147984000271, 10.842812387233078, 10.86569914139755, 10.88904498199313, 10.911654633645538, 10.939929859456134, 10.962750012570796, 10.99208239153427, 11.017357433409266]
norm_lat_a =  [0.99803256 1.         1.00211078 1.00426389 1.00634911 1.00895685
 1.01106149 1.01376672 1.01609777]

[5]:

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

The following plot is obtained by running the above code.

19d4895926c1495bba81bc01e1ea66fd

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.

3f11aae0d75c464b8f64cb6cf015bc8c

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\).

1498baf50df841718d01230e21491a12

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.

bf51276c298646e08bcbe72148dc5062

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