等圧等温（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   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

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

ここでは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:    8.9s remaining:    8.9s
[Parallel(n_jobs=16)]: Done  20 out of  20 | elapsed:   10.9s 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.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]

<matplotlib.legend.Legend at 0x7f4475de0e10>

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

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