カノニカル（NVT）アンサンブル¶
TL;DR¶
NVT-MDシミュレーションは温度および体積を一定とするシミュレーションで、体積変化や温度変化を無視したいときに使える計算方法。
用途としては固体中のイオン拡散、スラブ構造表面での吸着や反応、クラスター上の吸着や反応など。
NVT-MDシミュレーションの方法は温度をコントロールする熱浴法によりバリエーションがいくつか存在し、Berendsen熱浴、Langevin熱浴、NoseHoover熱浴などが一般的でASEにも実装されている。
Berendsen熱浴法は簡便で収束性も良いが、系全体の原子の速度分布を一律に変更するため不自然な現象が起こり得る。
Langevin熱浴法は統計的手法で適切にサンプリングされたランダム力を個々の原子に与えることで温度を制御する。各原子を個別に制御するため混相のようなものでも適切に扱える。計算を繰り返しても同じ軌道を再現しないため精密に軌道を解析するのには不向き。
Nosé–Hoover熱浴法は最も汎用的に利用される温度制御法のひとつ。原理的に正しいNVTアンサンブルを再現するが特殊な系で例外が存在する。
前節ではMDシミュレーションの最もシンプルな形であるNVEアンサンブルのMDシミュレーションついて学びました。
つづいて、よく使われる状態分布のひとつであるカノニカルアンサンブル（正準状態分布とも呼びます）におけるMDシミュレーションについて説明します。カノニカルアンサンブルは粒子数（N）、体積（V）、温度（T）が一定の統計力学的な状態分布です。検討したい対象で体積変化が無視できるような場合に適用されます。たとえば比較的低温領域で体積膨張が無視できるようなシミュレーションが良い例になります。
いくつかの実践的な具体的な例を通して、代表的な手法についてみていきましょう。
計算事例１： fcc-Alの溶融¶
一つ目の事例として6-1節で扱った金属アルミニウムの溶融状態の計算方法についてみてみましょう。
6-1節では金属アルミニウムが十分に溶融して見えるように1600 Kという高い初期温度でシミュレーションを行いました。しかし、今回は初期温度を1000 Kで計算してみます。NVEアンサンブルでの計算方法は前節と一緒で温度を1000 Kに設定するだけです。これを実行した結果が以下のようになります。
Fig6-2a. Fcc-Al in NVE ensemble starting at 1000 K
(File: ../input/ch6/6-2_fcc-Al_NVE_1000Kstart.traj)
確かに格子はだいぶ揺らいではいますが、元のFCCの結晶構造が明らかに見て取れると思います。Alの融点は933 Kですが、初期温度1000 Kでは不十分なようです。
続いてNVTアンサンブルで設定温度1000 Kで実行してみます。計算方法は以下のとおりです。計算時間は最初のNVEアンサンブルの計算と殆ど変わらないと思います。
[1]:
import os
from asap3 import EMT
calculator = EMT()
import ase
from ase.build import bulk
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution,Stationary
from ase.md.nvtberendsen import NVTBerendsen
from ase.md import MDLogger
from ase import units
from time import perf_counter
# Set up a crystal
atoms = bulk("Al","fcc",a=4.3,cubic=True)
atoms.pbc = True
atoms *= 3
print("atoms = ",atoms)
# input parameters
time_step = 1.0 # fsec
temperature = 1000 # Kelvin
num_md_steps = 1000000
num_interval = 10000
taut = 1.0 # fs
atoms.calc = calculator
print(f"taut = {taut:.3f}")
output_filename = "./output/ch6/fcc-Al_NVT-Berendsen_1000K"
print("output_filename = ",output_filename)
log_filename = output_filename + ".log"
print("log_filename = ",log_filename)
traj_filename = output_filename + ".traj"
print("traj_filename = ",traj_filename)
# Set the momenta corresponding to the given "temperature"
MaxwellBoltzmannDistribution(atoms, temperature_K=temperature,force_temp=True)
Stationary(atoms) # Set zero total momentum to avoid drifting
# run MD
dyn = NVTBerendsen(atoms, time_step*units.fs, temperature_K = temperature, taut=taut*units.fs, loginterval=num_interval, trajectory=traj_filename)
# 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}")
dyn.attach(print_dyn, interval=num_interval)
dyn.attach(MDLogger(dyn, atoms, output_filename+".log", 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 = Atoms(symbols='Al108', pbc=True, cell=[12.899999999999999, 12.899999999999999, 12.899999999999999])
taut = 1.000
output_filename = ./output/ch6/fcc-Al_NVT-Berendsen_1000K
log_filename = ./output/ch6/fcc-Al_NVT-Berendsen_1000K.log
traj_filename = ./output/ch6/fcc-Al_NVT-Berendsen_1000K.traj
imd Etot(eV) T(K) stress(mean,xx,yy,zz,yz,xz,xy)(GPa) elapsed_time(sec)
0 23.764 1000.00 7.17 7.19 7.18 7.14 -0.03 0.10 -0.04 0.002
10000 45.661 997.10 0.06 -0.03 0.23 -0.00 0.49 0.11 -0.68 6.141
20000 40.488 1001.20 1.28 0.95 1.39 1.49 0.07 0.31 -0.22 12.240
30000 42.446 999.51 0.91 0.35 1.33 1.06 -0.07 0.11 0.36 19.200
40000 44.621 1000.50 0.39 0.66 -0.03 0.55 0.69 0.27 0.51 28.047
50000 41.777 1002.44 1.03 1.36 1.23 0.49 0.02 -0.12 0.16 37.952
60000 43.185 1001.37 0.70 0.61 0.86 0.62 -0.45 0.17 0.24 46.973
70000 42.578 1005.13 0.82 0.74 0.88 0.83 -0.08 0.35 0.48 56.703
80000 44.024 998.34 0.28 -0.25 0.38 0.71 -0.21 -0.69 0.41 64.885
90000 43.371 999.74 0.55 0.46 0.54 0.64 0.30 -0.07 -0.64 73.537
100000 41.107 1002.71 1.25 1.90 0.63 1.22 0.00 -0.09 0.13 81.986
110000 42.462 999.09 0.77 0.79 0.85 0.68 0.18 -0.04 -0.13 90.240
120000 40.450 1000.00 1.06 1.76 0.66 0.75 -0.20 0.19 -0.05 99.628
130000 43.466 997.90 0.43 0.24 0.58 0.46 -0.48 0.59 -0.04 108.143
140000 42.179 998.05 0.92 0.51 1.05 1.21 -0.15 0.55 0.14 116.789
150000 42.484 1001.88 0.80 0.86 0.96 0.59 -0.16 0.38 -0.31 125.978
160000 41.827 996.92 0.81 0.73 0.83 0.89 0.37 0.22 0.41 134.414
170000 44.053 998.86 0.61 0.58 0.46 0.78 0.02 -0.22 0.16 142.578
180000 43.217 1002.70 0.61 1.13 0.75 -0.06 0.19 -0.56 -0.02 151.821
190000 44.495 1007.19 0.09 -0.33 0.42 0.17 -0.33 -0.06 0.12 160.257
200000 44.564 1006.14 0.37 0.84 0.23 0.05 -0.24 -0.03 0.07 168.267
210000 43.669 999.60 0.67 0.57 0.43 1.02 -0.14 0.25 -0.38 176.887
220000 42.006 1004.46 1.07 0.60 1.05 1.55 -0.23 -0.08 0.20 185.297
230000 43.512 1002.65 0.60 0.25 1.34 0.21 -0.11 0.18 -0.24 193.943
240000 42.741 1000.95 0.78 0.67 0.73 0.95 0.35 0.08 -0.13 202.552
250000 43.862 996.97 0.25 0.15 0.47 0.12 0.12 0.24 0.17 211.166
260000 43.631 1001.84 0.41 0.47 0.21 0.54 -0.07 0.19 0.37 220.092
270000 42.637 1005.87 0.61 0.68 0.99 0.15 0.48 0.24 0.56 229.247
280000 42.994 1002.63 0.83 0.29 1.52 0.69 -0.04 0.53 0.26 237.561
290000 42.622 996.87 0.72 0.19 1.38 0.58 -0.18 -0.29 -0.46 245.825
300000 42.120 1000.96 0.97 1.40 0.48 1.03 0.10 0.44 0.35 253.356
310000 41.948 999.98 0.91 1.08 0.95 0.72 0.31 0.07 0.05 261.935
320000 42.457 1003.75 0.83 0.66 1.23 0.59 -0.27 -0.01 -0.36 331.649
330000 43.615 1003.10 0.44 0.73 0.12 0.47 0.07 -0.15 0.47 363.026
340000 42.523 1001.14 0.64 0.59 0.62 0.69 -0.03 -0.09 -0.03 370.787
350000 42.058 999.63 0.89 0.86 1.11 0.71 -0.28 0.54 -0.31 379.067
360000 43.424 997.89 0.59 0.68 0.27 0.81 0.09 0.11 -0.25 387.963
370000 42.529 999.25 0.73 0.66 0.79 0.73 0.50 0.10 -0.42 394.554
380000 41.670 999.99 1.05 1.11 0.43 1.60 0.27 0.13 0.07 400.607
390000 42.868 994.82 0.61 -0.16 1.12 0.87 -0.03 0.28 -0.01 406.777
400000 42.262 998.90 0.96 0.88 0.61 1.38 -0.29 -0.25 -0.56 412.805
410000 42.061 998.11 0.80 0.94 1.11 0.35 0.08 -0.03 0.33 418.863
420000 42.227 1001.73 0.75 1.20 1.13 -0.09 0.23 -0.38 -0.10 425.068
430000 43.208 1006.62 0.57 0.54 0.70 0.47 -0.44 -0.59 0.03 431.148
440000 41.026 999.03 1.25 1.19 1.05 1.52 -0.05 -0.42 -0.11 437.426
450000 40.831 999.36 1.20 0.89 1.80 0.90 0.27 -0.10 0.21 443.452
460000 42.171 997.59 0.75 0.98 0.19 1.08 -0.09 -0.37 0.10 449.604
470000 42.703 1004.33 0.78 0.77 1.00 0.56 -0.00 0.51 0.10 455.938
480000 42.322 1005.63 0.90 0.50 1.15 1.07 -0.16 -0.15 0.30 461.979
490000 44.869 1006.86 0.33 0.54 0.10 0.34 0.07 0.05 0.76 468.042
500000 43.368 999.61 0.56 0.78 0.32 0.57 -0.08 -0.53 0.10 474.149
510000 42.557 1001.21 0.77 1.26 0.57 0.48 0.04 0.07 0.22 480.161
520000 43.193 1000.85 0.40 -0.39 0.64 0.94 -0.21 -0.38 0.30 486.252
530000 43.417 999.61 0.69 0.92 0.60 0.54 -0.25 0.27 0.11 492.244
540000 41.011 1001.73 1.10 0.95 1.22 1.13 -0.10 -0.16 0.27 498.538
550000 42.295 999.79 0.95 1.03 1.10 0.71 0.28 0.27 0.12 504.475
560000 42.773 1002.67 0.61 1.17 0.57 0.10 0.32 0.07 -0.28 510.951
570000 42.253 1000.13 0.85 0.80 1.03 0.73 0.13 0.37 0.18 517.089
580000 43.353 996.70 0.56 0.90 0.33 0.47 0.24 -0.01 0.16 523.234
590000 42.717 1001.11 0.68 0.63 0.96 0.44 -0.12 0.19 0.23 529.451
600000 41.756 996.23 1.04 0.50 1.59 1.04 0.33 -0.20 -0.09 535.695
610000 42.011 998.61 0.89 0.77 1.21 0.70 0.25 -0.38 -0.03 541.673
620000 43.435 1003.28 0.58 0.48 1.07 0.20 0.15 -0.06 0.23 547.782
630000 44.957 996.35 0.30 0.45 0.35 0.11 0.05 0.06 -0.35 553.840
640000 42.937 996.20 0.65 0.43 0.75 0.76 0.23 -0.56 -0.36 559.810
650000 42.769 1001.59 0.67 0.13 1.44 0.45 -0.42 -0.27 0.38 565.347
660000 43.322 1004.52 0.59 0.56 0.05 1.16 -0.09 0.16 -0.03 582.402
670000 41.388 997.71 1.05 0.72 1.17 1.26 -0.12 0.03 0.03 587.981
680000 41.888 1001.65 1.02 0.73 1.00 1.33 0.02 -0.13 -0.31 593.655
690000 42.636 1000.20 0.91 0.97 0.98 0.76 0.32 -0.13 -0.10 599.724
700000 43.108 992.74 0.67 0.99 0.49 0.51 -0.15 0.36 0.45 604.735
710000 43.644 1010.39 0.62 1.07 0.64 0.15 -0.37 -0.22 -0.21 607.843
720000 44.271 1004.63 0.27 0.03 0.26 0.51 -0.02 0.12 -0.41 611.116
730000 42.490 1003.67 0.97 0.28 1.53 1.10 -0.17 0.41 0.79 614.577
740000 43.006 1001.16 0.46 0.61 0.24 0.52 -0.04 0.09 -0.17 618.229
750000 42.048 994.53 0.89 0.70 0.87 1.11 -0.01 -0.18 0.23 621.761
760000 42.784 1004.37 0.82 1.08 0.46 0.93 -0.24 -0.15 0.10 625.236
770000 45.002 993.18 -0.04 0.44 -0.66 0.09 -0.09 0.54 0.17 629.138
780000 41.006 999.13 1.04 1.17 1.20 0.76 0.56 -0.23 0.36 632.664
790000 44.627 998.88 0.29 0.11 0.84 -0.08 0.42 -0.49 0.28 635.900
800000 43.433 1004.60 0.83 0.64 0.98 0.86 -0.02 0.23 0.01 639.517
810000 42.961 996.03 0.83 0.63 0.69 1.16 -0.08 0.44 0.24 642.676
820000 43.261 1008.76 0.47 0.20 0.56 0.64 -0.17 0.01 0.29 646.588
830000 44.443 998.68 0.19 0.47 -0.28 0.39 0.03 0.14 0.10 650.348
840000 43.358 1000.20 0.49 0.61 0.17 0.70 0.00 0.12 0.18 654.203
850000 40.909 999.72 1.21 1.41 1.07 1.16 0.45 -0.10 0.03 657.643
860000 41.880 1001.53 0.87 1.55 -0.20 1.26 -0.19 0.30 -0.33 661.020
870000 42.494 1001.03 0.64 1.11 -0.10 0.90 -0.38 -0.14 0.43 664.592
880000 42.393 999.43 0.71 0.99 0.85 0.28 -0.21 0.06 -0.37 668.199
890000 43.094 999.00 0.55 0.19 0.72 0.75 -0.19 0.15 0.08 671.932
900000 42.943 994.53 0.79 0.87 1.23 0.25 -0.20 -0.23 0.14 675.086
910000 44.272 999.36 0.31 -0.25 1.01 0.17 -0.01 0.33 -0.22 678.699
920000 42.750 994.83 0.76 0.56 0.97 0.76 0.10 -0.10 -0.11 681.836
930000 43.651 995.37 0.45 0.95 1.06 -0.65 -0.31 0.29 -0.05 685.068
940000 41.926 999.87 1.09 1.37 1.19 0.72 -0.30 -0.62 0.15 688.486
950000 42.630 1000.63 0.80 0.42 0.94 1.05 0.08 0.08 0.41 691.777
960000 44.108 999.31 0.40 0.54 0.40 0.26 0.22 -0.30 -0.54 695.243
970000 42.049 995.80 0.95 1.31 1.07 0.46 -0.21 0.07 0.47 698.910
980000 40.368 997.00 1.30 1.65 0.95 1.31 0.23 -0.49 0.37 702.153
990000 43.381 993.22 0.50 0.72 0.89 -0.10 -0.54 -0.12 -0.11 705.736
1000000 42.183 1002.23 0.88 0.89 1.06 0.70 -0.15 -0.01 -0.21 709.365
[1]:
True
今回はBerendsen熱浴法という方法を使っています。この計算手法では原子の速度を目標温度に近づくようにスケーリングして温度調整する方法です。（参考文献[1-3]）
具体的には、6-1節でも解説したように、温度と速度の関係は
で与えられるため、温度をコントロールするために計算対象である系全体の原子の速度をスケールして制御しています。実際に速度をスケールする係数は以下の式で記述されます。
ここで\(\Delta t\)は時間ステップサイズ、\(\tau_T\)は熱浴の時定数です。この方法では時定数\(\tau\)をコントロールすることでより緩やかに目標温度に到達することが可能です。\(\tau_T=\Delta t\)のとき、いわゆる速度スケーリング法と呼ばれる最も単純な温度制御法と同一になります。（つまり \(\lambda = \sqrt{T_0/T(t)}\) ）この場合、温度は常に所定の温度一定になり、温度の拘束条件が最大になっていると考えられます。逆に\(\tau\)が大きくなるにつれて温度制御は緩くなり、目標温度到達するのに時間がかかることになります。現実の物質の状態を考えるとある程度の緩和時間はあってしかるべきで、実用上はMDの計算時間に対して十分短くて済むような設定しするのが一般的です。ちなみに、この方法では系の温度は以下のように指数関数的に減衰するように補正されます。
それでは上記のNVT-MDシミュレーションの結果を見てみましょう。
[2]:
from ase.io import Trajectory
from pfcc_extras.visualize.povray import traj_to_apng
from IPython.display import Image
traj = Trajectory(traj_filename)
traj_to_apng(traj, f"output/ch6/Fig6-2_fcc-Al_NVT-Berendsen_1000K_unwrap.png", rotation="0x,0y,0z", clean=True, n_jobs=16)
Image(url="output/ch6/Fig6-2_fcc-Al_NVT-Berendsen_1000K_unwrap.png")
[Parallel(n_jobs=16)]: Using backend ThreadingBackend with 16 concurrent workers.
[Parallel(n_jobs=16)]: Done 18 tasks | elapsed: 6.1s
[Parallel(n_jobs=16)]: Done 101 out of 101 | elapsed: 12.5s finished
[2]:
[3]:
# Please remove comment out if you want to show trajectory in nglviewer.
from pfcc_extras.visualize.view import view_ngl
# view_ngl(traj)
Fig6-2b. Fcc-Al in NVT ensemble using Berendsen thermostat at 1000 K
(File: ../input/ch6/6-2_fcc-Al_NVT-Berendsen_1000K.traj」)
先ほどのNVE計算と打って変わって構造は非常にダイナミックに動いており、挙動は複雑でほぼランダムに配置された溶融状態であることが見て取れます。これはNVTアンサンブルでの計算で常に温度が1000 Kになるようにコントロールされているので常に融点を超えた状態になっているためです。（注釈：実際ASEで実行した結果を見ると原子が激しく拡散してセルの外に広がっていくことが分かります。表示の便宜上、上図はwrap関数を用いて原子をセルの内側に戻して描画しています。）
以下のプロットがNVEとNVTアンサンブルで計算した場合の全エネルギー（Tot.E）、ポテンシャルエネルギー（PE）、運動エネルギー（KE）、系全体の温度（Temp.）の時間発展になります。当然ではありますが、NVEアンサンブルでは全エネルギー、NVTアンサンブルでは温度が一定になっていることが見て取れます。それに対してBerendsen 熱浴法では全エネルギーは揺らぎますが、運動エネルギーと温度はほぼ一定に保たれています。（実際は小さな揺らぎらありますが、グラフのスケールの問題で揺らぎが小さくて見えていません。）
Fig.6-2c. Comparison of energies and temperatures between NVE (top) and NVT (bottom) calculations.
上記のNVTに関するプロットは以下のように行うことができます。NVEに関しては前節をご参照ください。
[4]:
import pandas as pd
df = pd.read_csv(
log_filename,
delim_whitespace=True,
names=["Time[ps]", "Etot/N[eV]", "Epot/N[eV]", "Ekin/N[eV]", "T[K]",
"stressxx", "stressyy", "stresszz", "stressyz", "stressxz", "stressxy"],
skiprows=1,
header=None,
)
df
/tmp/ipykernel_49628/2151791511.py:3: FutureWarning: The 'delim_whitespace' keyword in pd.read_csv is deprecated and will be removed in a future version. Use ``sep='\s+'`` instead
df = pd.read_csv(
[4]:
| Time[ps] | Etot/N[eV] | Epot/N[eV] | Ekin/N[eV] | T[K] | stressxx | stressyy | stresszz | stressyz | stressxz | stressxy | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.0 | 0.2200 | 0.0908 | 0.1293 | 1000.0 | 7.194 | 7.179 | 7.141 | -0.029 | 0.101 | -0.042 |
| 1 | 10.0 | 0.4228 | 0.2939 | 0.1289 | 997.1 | -0.032 | 0.225 | -0.002 | 0.492 | 0.109 | -0.681 |
| 2 | 20.0 | 0.3749 | 0.2455 | 0.1294 | 1001.2 | 0.949 | 1.388 | 1.489 | 0.070 | 0.313 | -0.224 |
| 3 | 30.0 | 0.3930 | 0.2638 | 0.1292 | 999.5 | 0.347 | 1.332 | 1.060 | -0.073 | 0.114 | 0.362 |
| 4 | 40.0 | 0.4132 | 0.2838 | 0.1293 | 1000.5 | 0.657 | -0.031 | 0.554 | 0.694 | 0.271 | 0.509 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 96 | 960.0 | 0.4084 | 0.2792 | 0.1292 | 999.3 | 0.537 | 0.401 | 0.261 | 0.221 | -0.302 | -0.541 |
| 97 | 970.0 | 0.3893 | 0.2606 | 0.1287 | 995.8 | 1.312 | 1.070 | 0.456 | -0.212 | 0.069 | 0.466 |
| 98 | 980.0 | 0.3738 | 0.2449 | 0.1289 | 997.0 | 1.651 | 0.947 | 1.314 | 0.232 | -0.487 | 0.365 |
| 99 | 990.0 | 0.4017 | 0.2733 | 0.1284 | 993.2 | 0.720 | 0.889 | -0.099 | -0.540 | -0.122 | -0.110 |
| 100 | 1000.0 | 0.3906 | 0.2610 | 0.1295 | 1002.2 | 0.886 | 1.061 | 0.699 | -0.155 | -0.012 | -0.212 |
101 rows × 11 columns
[5]:
import numpy as np
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(10, 5))
ax1 = fig.add_subplot(4, 1, 1)
ax1.set_xticklabels([]) # x axis label
ax1.set_ylabel('Tot E (eV)') # y axis label
ax1.plot(df["Time[ps]"], df["Etot/N[eV]"], color="blue",alpha=0.5)
ax2 = fig.add_subplot(4, 1, 2)
ax2.set_xticklabels([]) # x axis label
ax2.set_ylabel('P.E. (eV)') # y axis label
ax2.plot(df["Time[ps]"], df["Epot/N[eV]"], color="green",alpha=0.5)
ax3 = fig.add_subplot(4, 1, 3)
ax3.set_xticklabels([]) # x axis label
ax3.set_ylabel('K.E. (eV)') # y axis label
ax3.set_ylim([0.0, 0.2])
ax3.plot(df["Time[ps]"], df["Ekin/N[eV]"], color="orange",alpha=0.5)
ax4 = fig.add_subplot(4, 1, 4)
ax4.set_xlabel('time (ps)') # x axis label
ax4.set_ylabel('Temp. (K)') # y axis label
ax4.plot(df["Time[ps]"], df["T[K]"], color="red",alpha=0.5)
ax4.set_ylim([500., 1700])
fig.suptitle("Time evolution of total, potential, and kinetic energies, and temperature in NVT.", y=0)
#plt.savefig("Fig6-2_fcc-Al_NVT-Berendsen_1000K_E_vs_t.png") # <- Use if saving to an image file is desired
plt.show()
Berendsen熱浴法を用いた計算に関する注意点について。
適切な時定数の設定。 Berendsen熱浴法には温度制御パラメーターに時定数\(\tau\)があります。この\(\tau\)は目的温度に到達するための緩和時間を決めるパラメーターであり、\(\tau\)が小さいと収束の速度が速まります。\(\tau\)が余りに小さいと計算が不安定になり温度の値が不自然な挙動を示します。一般に\(\tau\)はMDシミュレーションのタイムステップの数倍から5倍程度が適切です。上記のfcc-Alの系に対して時定数\(\tau\)を変化した際の系全体の温度の挙動が以下のグラフで示しています。(Fig.6-2d参照）前記のとおり、\(\tau\)が0.2 fsの場合非常に不安定な挙動を示しており、大きく温度が上下に遷移して平衡状態を保てなくなっています。\(\tau\)が1もしくは2 fsではかなりタイトに温度が制御されており、5～10 fsで温度制御が緩くなり、比較的揺らぎは大きいですがより現実的な領域の揺らぎになっていると考えられます。
Fig.6-2d. Temperature evolution with various time constant \(\tau\).
非現実的な状態に陥る可能性。Berendsen熱浴法はいわゆる速度スケーリング法と呼ばれる部類の手法で、計算対象の全原子の速度を一律に変更します。速度分布全体を移動することで系全体の温度を目標温度に合わせに行きますが、この操作自体が非常に人工的で、熱力学的に正しいカノニカル（NVT）アンサンブルとは異なります。この手法の大きな問題は複数の挙動が大きく異なる物質が共存している系に対して適用すると非現実的な状態を作ってしまうことです。例えば固体と液体の界面の計算をNVTで行いたい場合、安易にBerendsen熱浴法を適用すると片方は高温に、もう一方は低温になり、トータルでは目標温度になるという不自然なことが起こり得ます。(本節末尾で具体例を解説しています。)
上記の点を気を付けて熱力学的な厳密性を求めるのでなければ、Berendsen熱浴法は平衡状態への収束が早く、原理的にもシンプルで分かりやすいので非常に便利な計算手法です。例えば単純に「熱による揺らぎを与えた状態を作りたい」などの場合には非常に有効な計算手法です。
計算事例２： 固体中のイオン拡散¶
さて、次のNVTアンサンブルの計算事例として固体電解質中におけるイオン拡散を考えてみます。
Li\(_{10}\)GeP\(_{2}\)S\(_{12}\)（LGPS）は近年注目を浴びている全固体電池の固体電解質として注目を集めている物質のひとつで、液系電解質に匹敵する非常に高いイオン伝導度（12 mS/cm @ 298 K）が特徴です。結晶構造は以下のとおりです。（基準となる構造はMaterials ProjectのLi10Ge(PS6)2_mp-696138_computedを用いています。以下に示す構造は、その構造をベースにPFPで最適化したものを使用しています。）
[6]:
atoms = ase.io.read("../input/ch6/LGPS_2x2x1_pfp-opt.cif")
view_ngl(atoms, representations=["ball+stick"])
Fig.6-2e. Crystal structure of Li\(_{10}\)GeP\(_2\)S\(_{12}\).
計算の便宜上、Unit Cellを2x2x1に拡張しています。これはLiイオンの数を増やすことで統計的に誤差が小さくなることを期待しているからです。計算時間は若干増えるかもしれませんが、実用上、拡張した構造でも全原子含め200原子程度なので特に問題はありません。
この材料内でのLiイオン拡散係数をMDシミュレーションを用いて求めてみたいと思います。この材料は組成も複雑で汎用の古典力場が存在しないのでMatlantisに実装されているPFPを使います。（学術論文を詳しく調べれば存在する可能性はありますが、すぐに見つからなかった事、見つかったとしても実装の難易度が不明なのでこれらのプロセスに係る時間や労力、精度の問題を考慮しても汎用的な機械学習ポテンシャルであるPFPを使えることは大きな利点です。）
今回のNVT-MDシミュレーションの実行スクリプトは以下のとおりです。 以下はテスト計算のためMDステップを10ステップとしています。 本計算では、1 fs/step で 10^6 steps計算を行い、1 nsのＭＤシミュレーションを行っています。
[7]:
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)
import ase
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.langevin import Langevin
from ase.md import MDLogger
from ase import units
from time import perf_counter
# Set up a crystal
atoms = ase.io.read("../input/ch6/LGPS_2x2x1_pfp-opt.cif")
atoms.pbc=True
print("atoms = ",atoms)
atoms.calc = calculator
# input parameters
time_step = 1.0 # fsec
temperature = 773 # Kelvin
num_md_steps = 10 # Total MD step, for testing.
# num_md_steps = 1e6 # Total MD step, for actual run.
num_interval = 1 # Output printing interval
friction_coeff = 0.002
print(f"friction_coeff = {friction_coeff:.3f}")
# Define output filenames
output_filename = "./output/ch6/Li10Ge(PS6)2_mp-696138_Langevin_773K"
print("output_filename = ",output_filename)
log_filename = output_filename + ".log"
print("log_filename = ",log_filename)
traj_filename = output_filename + ".traj"
print("traj_filename = ",traj_filename)
# run MD
dyn = Langevin(atoms, time_step*units.fs, friction=friction_coeff, temperature_K=temperature, loginterval=num_interval, trajectory=traj_filename)
# 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}")
dyn.attach(print_dyn, interval=num_interval)
dyn.attach(MDLogger(dyn, atoms, log_filename, header=True, stress=True, peratom=True, mode="w"), 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 = Atoms(symbols='Ge8Li80P16S96', pbc=True, cell=[[17.2189, 0.0, 0.0], [-0.07190443131522067, 17.742454297891182, 0.0], [-0.13259731631828015, -0.4464045597143071, 13.005365377826626]], spacegroup_kinds=...)
friction_coeff = 0.002
output_filename = ./output/ch6/Li10Ge(PS6)2_mp-696138_Langevin_773K
log_filename = ./output/ch6/Li10Ge(PS6)2_mp-696138_Langevin_773K.log
traj_filename = ./output/ch6/Li10Ge(PS6)2_mp-696138_Langevin_773K.traj
imd Etot(eV) T(K) stress(mean,xx,yy,zz,yz,xz,xy)(GPa) elapsed_time(sec)
0 -699.887 0.00 0.00 0.00 0.00 0.00 -0.00 -0.00 0.00 0.045
1 -699.879 0.32 0.00 0.00 0.00 0.00 -0.00 -0.00 0.00 0.105
2 -699.871 0.65 0.00 0.00 0.00 0.00 -0.00 -0.00 -0.00 0.164
3 -699.862 0.98 -0.00 -0.00 0.00 -0.00 -0.00 -0.00 -0.00 0.229
4 -699.854 1.31 -0.00 -0.00 -0.00 -0.00 -0.00 -0.00 -0.00 0.313
5 -699.842 1.76 -0.00 -0.00 -0.00 -0.01 -0.00 0.00 -0.00 0.368
6 -699.834 2.07 -0.00 -0.01 0.00 -0.01 -0.01 0.00 -0.00 0.439
7 -699.826 2.36 -0.01 -0.01 -0.00 -0.01 -0.01 0.00 -0.00 0.492
8 -699.819 2.62 -0.01 -0.01 -0.00 -0.01 -0.01 0.00 -0.00 0.567
9 -699.811 2.86 -0.01 -0.01 0.00 -0.01 -0.01 0.00 -0.00 0.619
10 -699.805 3.04 -0.01 -0.01 -0.00 -0.01 -0.01 0.01 -0.00 0.672
[7]:
True
本計算の時間でMDを実行した場合、Trajectoryを可視化すると以下の図が得られました。 紫色のLi イオンが大きく動いている(拡散している) ことが見てとれます。 なにもない場所に行っているように表示されていますが、周期境界条件のため、実際にはUnit cell外にもGe, P, Sが存在しています。
Li\(_{10}\)GeP\(_{2}\)S\(_{12}\) crystal at 773 K.
今回はLangevin熱浴法を用いてNVT-MDシミュレーションを行います。本手法では個々の原子に係る力をランダムフォースと呼ばれるある確率分布に従う力を作用させることで温度を制御する統計的なアプローチです。（参考文献[1,4]）
具体的には、Langevin動力学の運動方程式は以下のような形で書けます。
ここで\(m_i\)は対象原子の質量、\(\mathbf{x}_i(t)\)は時間\(t\)における原子\(i\)の座標ベクトル、\(\mathbf{F}(\mathbf{x}_i(t))\)が原子\(i\)にかかる力、\(\gamma\)は摩擦係数です。右辺最後の項がランダムフォースで、関数\(\mathbf{\eta}_i(t)\)が\(<\mathbf{\eta}_i(0)\mathbf{\eta}_j(t)>=\delta_{ij}\delta(t)\)の関係を満たす必要があります。これによりカノニカル分布を再現することが出来ます。
計算事例１で出てきたBerendsen熱浴法の問題点のところで言及したように、高速に拡散するLiとそれ以外のアニオン構造の運動が極端に異なる場合、Berendsen熱浴法を用いるとこれらの挙動が乖離して非現実的な結果をもたらす危険があります。Langevin熱浴法では個々の原子にかかる力を制御するためこのような問題が起きません。
Langevin熱浴法では摩擦係数\(\gamma\)を設定する必要があり、計算対象にもよりますが、一般的に1e-4から1e-2が概ね適切な範囲とされています。（値が大きいほど収束条件がタイトになり、目標温度に到達する速度が早まります。）概ね上記の値の範囲で問題ないことが多いですが、新規の計算対象を検討する際には事前検討として、この摩擦係数の値に対するセンシティビティを確認することをオススメします。
今回の計算事例は比較的高温状態の773 Kで、シミュレーション時間を1 nsに設定しています。高温でMDを行う理由は十分なLiイオンの拡散を保証するためで、あまり現実的ではない温度領域かもしれません。実際の物質がこの温度領域で安定な保証はありませんが、NVTアンサンブルでは密度が一定に保たれるため、構造は比較的安定に保たれるので拡散係数を計算することが可能です。このような高温状態での計算は一般的によく用いられる手法であり、比較的きれいなデータになりやすいこともあって本Tutorialの目的に適しています。
それでは上記のMDのTrajectoryから拡散係数を求めます。まずは対象の原子の初期位置からの平均二乗距離（mean-squared dispacement, MSD）を計算します。（参考文献[5]）
ここで\(\mathbf{x}_i(t)\)は時間\(t\)における原子\(i\)の座標,\(N\)は対象の原子数です。このMSDを横軸を時間にしてプロットしたものが以下のグラフになります。
Fig.6-2f. Mean-square displacement of Li-ion in Li\(_{10}\)GeP\(_{2}\)S\(_{12}\) crystal at 773 K.
ここで得られたMSDから拡散係数は以下の式で計算出来ます。
注意点としては、この式は系が十分に平衡状態に達した状態を仮定していて、計算結果に平衡状態でないデータが含まれる場合はその範囲のデータを除去しておく必要があることです。また、温度が低く拡散が起こりにくい、または対象の原子数が少ないなどの場合にMSDのプロットにノイズが含まれる場合、線形回帰などで傾きを計算することでより安定した結果を得ることが出来ます。
ここでの計算結果は拡散係数が2.7x10\(^{-5}\) cm\(^2\)/sとなり、概ね既報の計算値とよく一致しています。(参考文献[6]) ただし今回は簡単のため計算条件を最適化していません。実際に精度良い計算を効率よく行いたい場合には、あらかじめ計算ステップ数、計算ステップサイズ、温度、計算対象のサイズなど事前検討でこれらのパラメーターに対する依存性を調べ、良く吟味しておく必要があります。
このようにNVTアンサンブルを使うと特定温度での物質の拡散状態を評価できます。これを例えばNVEで検討しようとすると初期温度と定常状態に達した時の温度で差があるので正しく評価することが難しくなります。したがって特定温度での熱平衡状態を観察するのにNVTアンサンブルは有効な手法になります。
本計算事例に関して詳しくは、以下もご覧ください。
計算事例３： 金属表面での酸化反応¶
最後の事例として、金属表面における酸化反応の事例を検討します。キレイなfcc-Alの(111)表面があり、酸素雰囲気下に暴露されている状態を考えます。Alは室温で酸化するので、計算上そのような状態が再現出来るのか試行してみます。
計算対象のモデルとしては、fcc-Alの(111)のスラブモデルに10Åの真空層を作り、そこにランダムに酸素分子を配置した状態を初期状態とします。
[8]:
atoms = ase.io.read("../input/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K_wrapped.cif")
view_ngl(atoms, representations=["ball+stick"])
Fig.6-2g. Starting structure of fcc-Al(111) slab with 20 O\(_{2}\) molecules.
この事例においては金属界面と気相分子の反応を扱う必要があるのでMatlantisに実装されているPFPを使います。一般的には精度よく検討するにはDFTのような計算手法が必要になりますが、原子数も多くセルサイズも大きいため計算負荷が非常に高くなります。その点、PFPをもちいればDFTに匹敵する精度で反応プロセスも再現でき、原子数も十分確保できるのでより現実的な条件を非常に高速に試行検討することが可能です。
それではこれを300 KでNVTアンサンブルのMDシミュレーションを行ってみます。計算スクリプトは以下のとおりです。 (テスト計算のため``num_md_steps``を小さく設定しています。実際の計算時は大きくして実行してください。)
[9]:
import os
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_PLUS_D3)
calculator = ASECalculator(estimator)
from ase.io import read,write
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
# Set up a crystal
atoms = read("../input/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K_wrapped.cif")
atoms.pbc = True
atoms.calc = calculator
print("atoms = ",atoms)
# input parameters
time_step = 1.0 # fsec
temperature = 300 # Kelvin
num_md_steps = 200 # For test
# num_md_steps = 10000 # For actual run
num_interval = 100
ttime = 20.0 # tau_T thermostat time constant in fsec
print(f"ttime = {ttime:.3f}")
output_filename = "./output/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K"
print("output_filename = ",output_filename)
log_filename = output_filename + ".log"
print("log_filename = ",log_filename)
traj_filename = output_filename + ".traj"
print("traj_filename = ",traj_filename)
# Set the momenta corresponding to the given "temperature"
MaxwellBoltzmannDistribution(atoms, temperature_K=temperature,force_temp=True)
Stationary(atoms) # Set zero total momentum to avoid drifting
# run MD
dyn = NPT(atoms,
time_step*units.fs,
temperature_K = temperature,
externalstress = 0.1e-6 * units.GPa, # Ignored in NVT
ttime = ttime*units.fs,
pfactor = None, # None for NVT
loginterval=num_interval,
trajectory=traj_filename
)
# 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}")
dyn.attach(print_dyn, interval=num_interval)
dyn.attach(MDLogger(dyn, atoms, log_filename, header=True, stress=True, peratom=True, mode="w"), 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 = Atoms(symbols='Al72O40', pbc=True, cell=[8.59135, 9.92043, 31.6913], spacegroup_kinds=..., calculator=ASECalculator(...))
ttime = 20.000
output_filename = ./output/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K
log_filename = ./output/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K.log
traj_filename = ./output/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K.traj
imd Etot(eV) T(K) stress(mean,xx,yy,zz,yz,xz,xy)(GPa) elapsed_time(sec)
100 -390.145 269.61 0.82 1.16 0.86 0.43 0.12 -0.26 -0.09 6.812
200 -398.080 195.53 1.12 1.41 1.34 0.61 0.07 0.21 0.19 13.612
今回、熱浴法としてはNosé–Hoover熱浴を使用しています。 上記コードではNPT moduleが指定されていますが、pfactor=None と指定することでNosé–Hoover熱浴を用いたNVTのシミュレーションを行う事ができます。
この方法では2つの新たな独立変数（\(\xi\)と\(s\)）を運動方程式に導入することで、仮想の熱浴と接している状態を再現し温度を制御します。（参考文献[1,3-5]）解かなければいけない運動方程式は以下のような形になります。
制御変数は新たな時定数（\(\tau_T\)）でユーザーはこれを適切な値に設定する必要があります。ASEのマニュアルに従うと、25 fs付近が適切な値とされています。あまり小さい値は計算を不安定にし、大きすぎる値は収束が極端に悪くなり目標温度に到達するまでの時間がかかりすぎることになります。
それでは計算結果のtrajectoryを見てみましょう。(以下はAPNGファイルのみですが、キャプションにtrajファイル名が記載されているので、Matlantis上で実際にtrajectoryファイルを可視化して挙動を確認してみてください。
[10]:
from pfcc_extras.visualize.view import view_ngl
from ase.io import Trajectory
traj = Trajectory("./output/ch6/fcc111_Al_3x4x6_vac10A_20O2_NVT-NoseHoover_300K_wrapped.traj")
v = view_ngl(traj, representations=["ball+stick"], replace_structure=True)
display(v)
Fig.6-2h. Starting structure of fcc-Al(111) slab with 20 O\(_{2}\) molecules.
ご覧のとおり、最初はランダムに空隙中に配置された酸素分子が短時間のうちにAl表面に吸着、乖離して酸化アルミニウムを形成していく様子が見てとれます。アルミニウムは室温でも低酸素分圧下（\(pO_2\)=10\(^{-18}\) atm程度）で酸化物相を形成するので、今回のような僅か9x10x10Å\(^3\)程の空間に20もの酸素分子が詰め込まれている状態では、たった10 psの時間にも関わらず酸化アルミ被膜が形成されていきます。（アルミニウムの酸化平衡について知りたい方はEllignham diagramというキーワードで調べてみてください。具体的な温度と酸素分圧領域のイメージが分かるかとおもいます。）
上記のデータから確認できるとおり、MDシミュレーションを行うことによって個々の反応プロセスをつぶさに観察することが可能です。例えば以下のような状態がシミュレーションの初期段階で確認出来ます。Fig.6-2eでは左から緑の丸で囲まれた\(O_2\)分子が表面のAlのbridgeサイトに吸着し2つのAl原子に結合、その後\(O_2\)結合が開裂して表面構造に取り込まれます。その後３つのO原子に配位されたAlが次の\(O_2\)分子と吸着する様子が見て取れます。
Fig.6-2i. Reaction process of an O\(_2\) molecule on fcc(111) Al.
このようにNVT－MDシミュレーションを使うとガスと固体表面の吸着プロセスのダイナミクスを詳細に観察できます。これが例えばより複雑な分子や表面構造になっても、解析は難しくなるかもしれませんが基本的なプロセスと考え方は変わりません。
最後にNosé–Hoover熱浴について補足します。Nosé–Hoover熱浴は統計力学的に正しいカノニカルアンサンブルを再現することが知られており、非常によく用いられる熱浴法のひとつです。この手法の欠点は特殊なケース（例えば調和振動子等）で、意図しない保存則を満たすような系に対して適用すると正しいカノニカルアンサンブルを再現出来なくなります。これが何を意味するかというと、物理的に起こり得る状態を正しくサンプリング出来ないということになります。ただし実用的には、計算対象の系の重心を固定するなどでこのような問題を回避することが出来きるので非常に汎用的な手法です。（一般には、Nosé–Hoover法以外にも、Nosé–Hoover chain法のようなより発展的な手法も存在するので興味のある方は参考文献をご確認ください。[4,5]）
まとめ¶
最後に本節のまとめとして、現時点で利用可能なASEのNVTアンサンブルでのシミュレーション手法を整理してみます。
クラス名 |
ensemble |
パラメーター |
解説 |
メリット |
デメリット |
|---|---|---|---|---|---|
NVTBerensen |
NVT |
時定数(\(\tau_T\)) |
速度スケーリング法の一種。決定論的手法。 |
シンプルで汎用的 |
非カノニカル分布。系全体の温度のみ制御。 |
Langevin |
NVT |
摩擦係数(\(\gamma\)) |
確率的手法。ランダムな力と仮想的な摩擦項によって温度制御。 | カノニカル分布を再現。個別の原子の温度制御。 |
統計的手法 |
|
NPT(pfactor=None) |
NVT |
時定数(\(\tau_T\)) |
Nosé–Hoover熱浴。決定論的手法。 |
多くの系でカノニカル分布を再現。 |
非カノニカル分布な系も存在。人工的な運動モードが入り得る。系全体の温度のみ制御。 |
Andersen |
NVT |
頻度因子(\(\nu\)) |
確率的手法。ランダムに選択された原子群の速度を変更。 | シンプル。個別の原子の温度制御。 |
非カノニカル分布。統計的手法 |
|
実際にNVTアンサンブルの計算を行う際に、どの手法を適用するか迷ったらNosé–Hoover熱浴法を用いるのが無難な選択肢であると思います。原理的に正しく、欠点が少ないため多くの分子シミュレーションプログラムに実装されており、ユーザーがすぐに利用できるコトが多いかと思います。本手法は決定論的な手法であり、各原子の位相空間上の軌道（i.e.座標、速度の軌跡）を調べたいときにも再現性の観点で適しています。
Langevin熱浴法は液系の計算やブラウン運動のような計算を行いたいときによく用いられます。原理的にもNVTアンサンブルを再現し、効率よく温度制御出来るため良く用いられる手法です。ただし確率論的手法であり、精密な原子の軌道の再現を期待できるものではありません。ASEのマニュアルには吸熱反応を伴う系に適用すると温度制御がうまくできないというような記述もありますが、少なくとも我々の調査の中ではそのような事実示すデータが確認できていません。
Berendsen熱浴法（NVTBerendsen）はシンプルで汎用的手法であり、実装も簡単です。検討初期のテスト計算などを行いたいときに始めるには良いでしょう。ただし、厳密な意味でNVTアンサンブルを再現できるとは言えず、熱的な揺らぎも人工的なので精度を欠きます。熱的に揺らいだ構造を作成したい、非平衡状態を作り出したいなどのような条件の時に有用な手法と言えるでしょう。
最後にASEに実装されている手法でAndersen熱浴法というのがありますが、特に他の手法と比較して特筆すべき利点もなく、上の3つの計算方法があれば事足りることが殆どで、活用される事例が少ないため今回は説明は割愛しています。
[Advanced] Nosé–Hoover熱浴法のパラメーター依存性¶
温度制御に用いられるNosé–Hoover熱浴の時定数\(\tau\)に対する系全体の温度の依存性をみてみます。fcc-Al（3x3x3 unit cells）を\(\tau_t\)を変えて計算した結果が以下の通りです。 計算にはPFP(v2.0.0)を用いています。初期温度は特に設定していません。
Fig.6-2j. Time evolution of temperature in Nose-Hoover thermostat
\(\tau\)の値で挙動が大きく異なるので３つの領域に分けて見てみましょう。
小さいほど目的温度に近づく速度は早くなりますが、5 fsでは平均値が目的温度に到達しても温度振動の振幅が大きく上は750 Kから下は100 Kまで変動しています。これでは平均値は安定していても「定温」と仮定するのは難しいと思われます。
10-40の範囲では比較的目標温度への収束も早く、温度の収束性も安定しています。(上図下部）ただし、40 fsの場合振動の周期に人工的な周期性があるように見えます。
\(\tau\)が100 fs以上の場合、温度上昇も遅く非常に大きく波打つように振動しています。いずれは目標温度に収束するのかもしれませんが、大きく波打つ挙動の中で温度が数千Kに達しているところもみられ、検討している対象によっては大きな影響を受ける可能性があります。
上記の例では\(\tau\)が10か20 fs辺りが妥当な時定数なようです。ただし、計算の初期で瞬間的にですが1000 Kを超える温度になるところがあります。この問題を回避することが出来るのでしょうか？ためしに事前にコンディショニングとしてMB分布を利用した温度の初期化、もしくはBerendsen熱浴法等の方法による平衡状態を作って計算してみましょう。Berendsen熱浴は計算開始から1 psだけ実行したのちにNosé–Hoover熱浴に切り替えます。実際このようなプレコンディショニングはどの程度有効でしょうか？
Initialized by the Maxwell-Boltzmann Distribution
Initialized by Berendsen thermostat (τ=1.0fs)
Initialized by Berendsen thermostat (τ=10.0fs)
Initialized by Berendsen thermostat (τ=30.0fs)
Fig.6-2k. Time evolution of temperature with various initialization methods:
Log-Log scale (left), Linear scale (right)
構造最適化された系（i.e. 0 K）に対して、単純にNosé–Hoover法を適用した時は、瞬間的にですが、約1000-2000Kまで温度が上昇していました。これに対してこれらいずれの方法を用いても計算初期の温度上昇が約500-600 Kまでに抑えられていることが確認できます。問題は完全に解消されませんでしたが、それなりに緩和されたと考えられそうです。
[Advanced] Berendsen熱浴法の問題 - ２相共存状態でのNVTアンサンブル¶
Berendsen熱浴法の問題点を示すための典型的な事例として、水中に沈んだAl原子クラスターの例を見てみます。計算対象は水分子が165分子、Al原子クラスター（135原子）です。水分子はおおよそ1g/ccとなるように設計しています。 細かい設定はあまり重要ではないのですが、比較的現実的な密度や温度の範囲の設定にしています。以下が計算対象の構造を可視化したものです。
Fig.6-2l. Starting structure of fcc-Al(111) slab with 20 O\(_{2}\) molecules.
(File: ../input/ch6/165h2o_and_135Al_20x20x20_opt.cif)
さて、この系をMDシミュレーションでBerendsen熱浴法とNosé–Hoover熱浴法を適用して100℃で熱平衡状態にしてみます。MDの時間ステップは1 fs, 計算時間は10 psとします。以下に水分子、Alクラスター、そして系全体の温度のプロットを示します。
Fig.6-2m. Comparison of the temperature evolution in Berendsen and Nosé–Hoover thermostats. Each line represents the temperature of the entire system (blue), Al cluster (gray), and water (orange), respectively.
それぞれのグラフで青色の線が系全体の温度、グレーとオレンジがそれぞれAlクラスターと水の温度です。Berendsen熱浴の場合、Alクラスターの温度が 平均420Kぐらいなのに対し、H\(_2\)Oは360Kぐらいを推移しています。これはBerendsen熱浴法で全ての原子を一つの速度スケーリングパラメーターで制御しているためです。 系全体の温度は非常に良く安定しているのですが、AlクラスターとH\(_2\)O分子とでは原子の挙動が大きく異なるので、それぞれの温度が分離（decouple）してしまっています。 比較的ぬるいお湯の中に熱いAlの球が浮いているようなものです。現実的には起こり得ない現象で正しい熱平衡状態が再現されていない良い例です。
これに対して、Nosé–Hoover熱浴法を用いた方は、系全体の温度と構成要素のそれぞれの温度が概ねよく一致していて問題ないことが確認出来ます。
参考文献¶
[1] M.P. Allen and D.J. Tildesley, “Computer simulaiton of liquid”, 2nd Ed., Oxford University Press (2017) ISBN 978-0-19-880319-5. DOI:10.1093/oso/9780198803195.001.0001
[2] 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, DOI:10.1063/1.448118
[3] “GROMACS 2019-rc1 Reference Manual”, https://manual.gromacs.org/documentation/2019-rc1/reference-manual/index.html
[4] 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#
[5] D. Frenkel and B. Smit “Understanding molecular simulation - from algorithms to applications”, 2nd Ed., Academic Press (2002) ISBN 978-0-12-267351-1. DOI:10.1016/B978-0-12-267351-1.X5000-7
[6] Mo et al. “First Principles Study of the Li10GeP2S12 Lithium Super Ionic Conductor Material”, Chem.Mater. (2012) 24, 15-17 https://pubs.acs.org/doi/10.1021/cm203303y