等圧等温(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 |
セル角度は保持するが、圧力の異方性は考慮可能 |
表中の2番目と3番目の手法はBerendsen barostatで本質的に同じものです。(ちなみに3番目の手法であるInhomogeneousBerendsenはASEのマニュアルで記載はないのですが、NPTBerendsenとともにASEのなかでクラスが定義されています。)したがって、ASEのフレームワーク内で使える機能としてはParrinello-Rahman法とBerendsen法の2種類のみになります。熱浴自体は6-2節で解説した手法が用いられているので、これらの熱浴法の特徴をよく考慮する必要があります。
それではまず、系の自由度と汎用性が比較的高いParrinnello-Rahman法を見てみましょう。
Parrinello-Rahman法の運動方程式¶
Parrinello-Rahman法はいわゆる拡張系の計算手法で、Nosé–Hoover熱浴法の時のように、計算対象の系が仮想的に外部の一定温度、一定圧力の系と接続されていることを仮定します。この場合の運動方程式は下記のように記述されます。(導出の詳細については参考文献[1-3]をご確認ください。)
ここで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のシミュレーションで十分に平衡状態に達します。
[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
結果を可視化してみると以下のようになります。
ここではread
関数に index="::100"
を指定することでTrajectoryファイルから間隔を間引いた結果に対して可視化しています。 Cuは固体の構造を保っており、大きく動いているわけではありませんが、振動の動きによりセルサイズに変化があることがわかります。
[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)
熱浴の時定数\(\tau_t\)は20 fs、圧力制御パラメーターであるpfactorは2e6 GPa\(\cdot\)fs\(^2\)としています。上記の計算条件で熱平衡状態にするため300 Kで20 psのNPT-MDシミュレーションを行った際のセル体積の時間変化の様子を確認してみましょう。以下のようなコードで特定温度におけるMDシミュレーションの結果であるtrajファイルを解析出来ます。
[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()
上記コードを実行すると、以下のような解析結果が得られます。
NPTアンサンブルで計算することでセル体積が約4%程の範囲で振動することが見て取れます。今回の計算対象は立方晶であるため、この体積の幾何平均から結晶の格子定数を計算出来ます。200 Kから1000 Kまで実行し、その各温度での規格化した格子定数の平均値をプロットすると以下の結果が得られます。(平衡状態に達したあとの格子定数を計算したいため、np.mean(vol[int(len(vol)/2):])**(1/3)
の部分では、Trajectoryの後半半分のみを使用して、体積から格子定数を計算しています。)
[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>
上記のコードを実行すると以下のようなプロットが得られます。
この結果にはさらにPFPを用いた計算結果、および参考までに実験値を比較対象としています。EMTおよびPFPのデータ共に非常に良い一致を示しており、小さな差ではありますがPFPはより実験値に近い傾向を示しています。
この規格化した格子定数の温度依存性から線熱膨張係数(CTE、\(\alpha\))を以下の関係式をもちいて算出します。
ここで\(\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です。
pfactorが小さい領域では高速でセル体積が振動しており、高周波と低周波の振動が混在していて挙動が不安定な領域もあるのであまり好ましくありません。pfactorが大きくなるにつれて振動の周期が徐々に長くなり、小さなpfactorの計算初期で起きるの大きなセル体積の変化も殆ど見られません。明確にpfactorのどの値から使うと良いという指標があるわけではないですが、10\(^6\) Ga\(\cdot\)fsec\(^2\)以上であれば小さな振動は確認できるものの中央値は相変わらず安定しているようです。上限についても同様で明確な基準はなく、pfactorが大きいと振動周期が長くなり扱いづらいので、上記の例に関していえば10\(^6\)から10\(^7\)ぐらいが妥当な領域かと考えます。
[Advanced] Berendsen barostatのパラメーター依存性¶
Berendsen barostatを用いた計算方法について説明します。Berendsen barostatは以下の方程式に従って圧力の時間発展が計算されます。(導出の詳細は参考文献[5]をご覧ください。)
上式から明らかなように、Berendsen圧力制御法では指数関数的に系の各瞬間における圧力を外圧\(\mathbf{P}_o\)に近づけていきます。その速度は時定数\(\tau_P\)により制御されます。
各MDステップ毎に、各原子の座標とセルベクトルは以下の式で表される係数でスケーリングされます。
熱浴のコントロールの際に\(\tau_T\)が時定数だったように、圧力制御法でも適切な時定数\({\tau_P}\)を設定してやる必要があります。それでは実際の計算事例を見てみましょう。
ASEの機能であるNPTBerendsenクラスを用いるとdynamicsを定義するオブジェクトは以下の形で記述されます。
[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
もほぼ同様の設定です。(クラス名をNPTBerendsen
からInhomogeneous_NPTBerendsen
とすることで利用可能です。)1点異なるのは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\)の方からみてみます。
時定数が小さいほど振動の周期が不安定で大きく暴れ、逆にあまり大きくとりすぎると変化があまりに緩やかで平衡に至るまでの時間がかかります。この辺りは熱浴法の時定数と同じ考え方です。力学的にfcc-Cuと類似の系に関しては体積変化の安定性と収束性から\(\tau_P\)は10\(^2\) fsから10\(^3\) fs付近が適切なようです。
次に\(\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