等圧等温（NPT）アンサンブル

TL;DR

• NPT-MDシミュレーションでは圧力と温度が一定に制御される手法で、圧力制御（barostat）に関してASEではParrinello-Rahman法、Berendsen法の2種類の一般的な手法が存在。

• 用途としては固体や流体の熱膨張、相転移、加圧状態のシミュレーション等

• Parrinello-Rahman法ではシミュレーションのセルの自由度が全て可変。pfactorを適切に設定する必要があり。

• Berendsen barostatではBerendsen thermostatと同様、非常に効率的に収束性良く制御可能。Berendsen barostatはセルの角度は固定して各セル長を独立に可変、もしくは各セル長の比率を固定する2つのモードで計算可能。compressibilityを入力パラメーターとして適切に設定が必要。

本節では圧力と温度の両方が一定となるような平衡状態を作る計算手法について説明します。一般に圧力制御の仕組みはbarostatと呼ばれ、一般に6-2節で出てきた熱浴法（thermostat）と同時に用いることで等温等圧アンサンブル（isothermal-isobaric ensemble、またはNPT）と呼ばれる状態分布を生成することを目的としてます。NPT-MDシミュレーションで（原理的に）検討可能な現象には

• 固体の熱膨張率

• 融点の予測

• 固体の相転移

• 流体（ガス、液体）の密度予測

etc.

が考えられます。原理的にという枕詞がついている理由は、これらの現象の再現性は計算に用いる力場の精度に大きく依存するためで、特に小さなエネルギー差に依存する分子間力や流体の状態予測は非常に難しいことが知られているからです。本チュートリアルでは比較的精度が高いと考えられる固体の事例を通してNPT-MDについて学んでいきます。

まずは本チュートリアルで利用するASEに実装されているNPT-MDの手法について確認します。2022年6月現在、ASEで標準的に利用可能な実装は以下の3種類が存在します。

クラス名	ensemble	パラメーター	熱浴	圧力制御	コメント
NPT	NPT	時定数(\(\tau_t\)),圧力因子(pfactor)	Nosé–Hoover	Parrinello-Rahman	セルの全自由度が可変、制御可能
NPTBerendsen	NPT	\(\tau_t\),\(\tau_P\),\(\beta_T\)	Berendsen	Berendsen	セル形状は維持し体積変化のみ
InhomogeneousBerendsen	NPT	\(\tau_t\),\(\tau_P\),\(\beta_T\)	Berendsen	Berendsen	セル角度は保持するが、圧力の異方性は考慮可能

表中の２番目と３番目の手法はBerendsen barostatで本質的に同じものです。（ちなみに3番目の手法であるInhomogeneousBerendsenはASEのマニュアルで記載はないのですが、NPTBerendsenとともにASEのなかでクラスが定義されています。）したがって、ASEのフレームワーク内で使える機能としてはParrinello-Rahman法とBerendsen法の２種類のみになります。熱浴自体は6-2節で解説した手法が用いられているので、これらの熱浴法の特徴をよく考慮する必要があります。

それではまず、系の自由度と汎用性が比較的高いParrinnello-Rahman法を見てみましょう。

Parrinello-Rahman法の運動方程式

Parrinello-Rahman法はいわゆる拡張系の計算手法で、Nosé–Hoover熱浴法の時のように、計算対象の系が仮想的に外部の一定温度、一定圧力の系と接続されていることを仮定します。この場合の運動方程式は下記のように記述されます。（導出の詳細については参考文献[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}\]

ここでNosé–Hoover熱浴以外の項に限定すると、圧力制御の時定数\(\tau_P\)、系の重心\(R_o\), 目標の外圧\(P_o\),シミュレーションセル体積\(V\)があります。\(\eta\)が圧力制御自由度を表す変数になっています。\(\mathbf{h} = (\mathbf{a}, \mathbf{b}, \mathbf{c})\)で\(\mathbf{a},\mathbf{b},\mathbf{c}\)はそれぞれシミュレーションセルの各辺を定義するセルベクトルです。

上記の方程式で温度\(T_o\)や圧力\(P_o\)以外にユーザーが設定しなければいけない値は\(\tau_T\)と\(\tau_P\)の2つがあります。それではまず、単純に\(\tau_T\)と\(\tau_P\)にを任意の値に設定した場合を以下に示します。ここで\(\tau_T\)を20 fsecとしています。\(\tau_P\)を直接指定しませんが、pfactorと呼ばれる値が\(\tau_P^2B\)となります。\(B\)はbulk modulusを指しており、この値をあらかじめ計算して指定する必要があります。しかし、\(\tau_P\)自体の正確な値が事前にわからないのでpfactorを特定する術がありません。また\(B\)の値も異方性を示す構造や異なるタイプの物質が混在する際は計算しようがありません。したがって、pfactorそのものの概ねの値を設定してbarostatの挙動を調べます。以下の例でも言及しますが、少なくとも金属の結晶系では約10\(^6\) GPa\(\cdot\)fs\(^2\) から10\(^7\) GPa\(\cdot\)fs\(^2\) のオーダーの値を用いると計算が安定して収束性も良いようです。それ以外の材料となる場合は、念のため事前検討として、pfactorの値を振ってみて体積変化の様子を確認しておくことをオススメします。

計算事例：熱膨張係数の計算

それではNPTアンサンブルの計算事例としてNosé–Hoover thermostatとParrinello-Rahman barostat（ASEのNPTクラス）を使って固体の熱膨張係数を計算してみます。温度を変化させながら熱平衡状態を作り出し、その時の格子定数の平均値から各温度における熱膨張係数を算出します。

今回は簡単のためfcc-Cuを用いて計算を行います。計算に用いられるスクリプトは以下のとおりです。温度は200 Kから1000 Kまで100 K刻み、外圧は1 barとします。対象の構造はfcc-Cuを3x3x3 unit cellsに拡張して108原子で計算を行います。以下では高速性のためASAP3-EMT力場を用いていますが、PFPで計算する場合も全く同様です。時間ステップは1 fsで20 psのシミュレーションで十分に平衡状態に達します。

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

結果を可視化してみると以下のようになります。

ここではread関数に index="::100" を指定することでTrajectoryファイルから間隔を間引いた結果に対して可視化しています。 Cuは固体の構造を保っており、大きく動いているわけではありませんが、振動の動きによりセルサイズに変化があることがわかります。

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

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)

熱浴の時定数\(\tau_t\)は20 fs、圧力制御パラメーターであるpfactorは2e6 GPa\(\cdot\)fs\(^2\)としています。上記の計算条件で熱平衡状態にするため300 Kで20 psのNPT-MDシミュレーションを行った際のセル体積の時間変化の様子を確認してみましょう。以下のようなコードで特定温度におけるMDシミュレーションの結果であるtrajファイルを解析出来ます。

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()

上記コードを実行すると、以下のような解析結果が得られます。

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

NPTアンサンブルで計算することでセル体積が約4%程の範囲で振動することが見て取れます。今回の計算対象は立方晶であるため、この体積の幾何平均から結晶の格子定数を計算出来ます。200 Kから1000 Kまで実行し、その各温度での規格化した格子定数の平均値をプロットすると以下の結果が得られます。（平衡状態に達したあとの格子定数を計算したいため、np.mean(vol[int(len(vol)/2):])**(1/3) の部分では、Trajectoryの後半半分のみを使用して、体積から格子定数を計算しています。）

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]

<matplotlib.legend.Legend at 0x7ff29dee4dc0>

上記のコードを実行すると以下のようなプロットが得られます。

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

Experimental data is taken from Reference [4].

この結果にはさらにPFPを用いた計算結果、および参考までに実験値を比較対象としています。EMTおよびPFPのデータ共に非常に良い一致を示しており、小さな差ではありますがPFPはより実験値に近い傾向を示しています。

この規格化した格子定数の温度依存性から線熱膨張係数（CTE、\(\alpha\)）を以下の関係式をもちいて算出します。

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

ここで\(\alpha\)が線熱膨張係数で、\(a(T)\)と\(a_{RT}\)は温度\(T\)と室温における格子定数です。以下がASAP3-EMT、PFP、実験値[4]のサマリーです。

	\(\alpha\) ( \(10^{-5}\) /K)
emt	2.23
pfp	2.13
exp	1.74

計算誤差を考えると、いずれの計算結果も実験値とリーズナブルな範囲で一致していると考えてよいと思います。

このような手法でNPTアンサンブルのMDシミュレーションを用いることによって熱膨張係数を算出することが可能です。

最後に補足として言及しますと、「液体や気体の熱膨張も同様の手法で再現できるのでは？」と考えるかと思います。原理的には確かにその通りなのですが、現時点では再現性の良いモデルというのは限定的と思われます。古典力場で特定の液体や気体のみに特化したモデルは存在するのかもしれませんが、第一原理計算を用いた場合、液体や気体では固体に比べて非常に小さい分子間相互作用の精度が重要になり、現在ではそのような精度の量子化学計算自体が非常に困難です。特にNNPを作成する際には膨大な計算データが必要になり、高精度な量子化学計算を多量に行うのは現時点であまり現実的ではありません。従って今後この領域での精度が高いモデルの開発が期待されるところです。

[Advanced] Parrinello-Rahman barostatのパラメーター依存性

NPTアンサンブルでのMDシミュレーションを行う際にpfactorというパラメーターを設定する必要があり、適切な値の範囲が 10\(^6\) GPa\(\cdot\)fs\(^2\)から10\(^7\) GPa\(\cdot\)fs\(^2\)付近という説明がありました。ご参考までにpfactorの値を変化し際の結果を以下に示します。計算対象は上記と同じfcc-Cuの3x3x3 unit cellsで、温度は300 Kです。

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

pfactorが小さい領域では高速でセル体積が振動しており、高周波と低周波の振動が混在していて挙動が不安定な領域もあるのであまり好ましくありません。pfactorが大きくなるにつれて振動の周期が徐々に長くなり、小さなpfactorの計算初期で起きるの大きなセル体積の変化も殆ど見られません。明確にpfactorのどの値から使うと良いという指標があるわけではないですが、10\(^6\) Ga\(\cdot\)fsec\(^2\)以上であれば小さな振動は確認できるものの中央値は相変わらず安定しているようです。上限についても同様で明確な基準はなく、pfactorが大きいと振動周期が長くなり扱いづらいので、上記の例に関していえば10\(^6\)から10\(^7\)ぐらいが妥当な領域かと考えます。

[Advanced] Berendsen barostatのパラメーター依存性

Berendsen barostatを用いた計算方法について説明します。Berendsen barostatは以下の方程式に従って圧力の時間発展が計算されます。（導出の詳細は参考文献[5]をご覧ください。）

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

上式から明らかなように、Berendsen圧力制御法では指数関数的に系の各瞬間における圧力を外圧\(\mathbf{P}_o\)に近づけていきます。その速度は時定数\(\tau_P\)により制御されます。

各MDステップ毎に、各原子の座標とセルベクトルは以下の式で表される係数でスケーリングされます。

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

熱浴のコントロールの際に\(\tau_T\)が時定数だったように、圧力制御法でも適切な時定数\({\tau_P}\)を設定してやる必要があります。それでは実際の計算事例を見てみましょう。

ASEの機能であるNPTBerendsenクラスを用いるとdynamicsを定義するオブジェクトは以下の形で記述されます。

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もほぼ同様の設定です。（クラス名をNPTBerendsenからInhomogeneous_NPTBerendsenとすることで利用可能です。）１点異なるのはmaskが設定でき、mask=(1, 1, 1)とするとa,b,c方向全て独立して変化可能です。このtupleの要素を0と設定することでその方向でcellが固定されます。

いずれのBerendsen barostatでも注意が必要なのは温度、圧力のような計算条件以外に、圧力制御時定数(taup、\(\tau_P\)）や圧縮率（compressibility、\(\beta_T\)）と呼ばれるパラメーターを設定する必要があることです。それぞれの値に対するセル体積の時間変化の依存性を300 Kにおけるfcc-Cuの例を用いて確認します。まずは\(\tau_P\)の方からみてみます。

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

時定数が小さいほど振動の周期が不安定で大きく暴れ、逆にあまり大きくとりすぎると変化があまりに緩やかで平衡に至るまでの時間がかかります。この辺りは熱浴法の時定数と同じ考え方です。力学的にfcc-Cuと類似の系に関しては体積変化の安定性と収束性から\(\tau_P\)は10\(^2\) fsから10\(^3\) fs付近が適切なようです。

次に\(\beta_T\)に対する依存性です。結果は以下のとおりです。

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

\(\beta_T\)が小さいほど収束が非常に遅く、高い値にあるほど不安定になっていく様子がうかがえます。グラフの傾向から\(\beta_T\)は10\(^{-7}\)から10\(^{-6}\)fs位が良さそうです。

\(\tau_P\)と\(\beta_T\)に関してはいずれも厳密に正しい値というものは無く、おおよそ桁で値を変えた時にこの程度傾向が変化するということを理解しておけば概ね問題はなさそうです。

最後に補足すると、これらの数値はあくまでfcc-Cuのような金属で密な構造を持つ物質について適用出来ますが、もし全く異なる系統の物質（例えばポリマー、液体、気体等）をNPTで扱いたいとなれば事前検討でこれらの適切な値の領域を確認しておく必要があります。手間を惜しんで事前検討を省くと意図しない結果が得られて余計に時間がかかるなどになりかねないので、通常、新しい材料系に取り組む際は注意されることをお勧めします。

参考文献

[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