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]).
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()
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.
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>
The following plot is obtained by running the above code.
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
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.
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]).
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
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\).
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.
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