分子動力学シミュレーション¶
TL;DR¶
MDシミュレーションでは古典力学の運動方程式に則った計算方法で、全ての原子の位置と速度の時間発展を露わに観察できる。
時間スケールの制約があり一般的に数nsから数十nsぐらいが上限。その時間内に発現しない現象は別の計算手法を検討する必要がある。
MDシミュレーションには再現したい状態に応じて様々な状態（アンサンブル）があり、最も単純なアンサンブルはNVEアンサンブルという。
NVEアンサンブルでは全エネルギー保存則が成立し、計算精度および計算時間の観点から、時間ステップサイズを適切に設定する必要がある。
この章では系の時間発展をシミュレーションする、 分子動力学シミュレーション（Molecular dynamics, MD） について学びます。
MDシミュレーションでは個々の原子の軌道の時間発展を露わに取り扱うシミュレーションで、古典力学の運動方程式を積分することによって計算対象である原子の座標と速度を逐次計算していく手法です。この計算手法自体は原子間に働く力やエネルギーのモデルとは独立した理論で、古くから分子シミュレーションの分野で用いられてきました。従って、理論的な背景や事例等は多くの書籍や文献が多く存在するのでそれらを確認していただければと思います。（参考文献：[1-3]） 本チュートリアルではあくまでMatlantisを使って実践的にこれらの計算を実行するために必要な知識を習得するところが目的となっています。
早速ですが具体的な事例を通して、どんなことがMDシミュレーションで観察出来るかみてみましょう。
事前準備 - 必要なライブラリのインストール¶
本章ではところどころ、ASE上で簡易に利用できる古典力場であるASAP3のEMT力場を利用しています。古典力場であるため精度や用途はかなり限定的ですがシンプルかつ高速で実行可能なのでチュートリアルとして事例を示すには十分な機能を持っています。ちなみにASAP3-EMTで利用可能な元素はNi、Cu、Pd、Ag、Pt、Auとこれらの合金に限定されます。インストールの仕方は以下のとおりです。
[ ]:
!pip install --upgrade asap3
ASAP3-EMTに関する詳細はASAPのサイトをご確認ください。
MDシミュレーションでなにが出来るか¶
以下の事例では金属Al構造の溶融状態をMDシミュレーションで再現しています。ご存知の通りアルミニウムは我々がよく見かける物質で工業的に非常に重要な金属であり、その特性も良く知られています。融点は660.3 °C（933.45 K）で、常圧下では幅広い温度領域で面心立方格子（face-centered cubic, fcc）構造を持ちます。以下はそのfcc-Al構造を初期温度1600 Kに過熱して熔融させる過程を示したものです。計算時間は100 psecです。
Fig6-1a. Melting of fcc-Al in NVE ensemble starting at 1600 K
(File: ../input/ch6/6-1_fcc-Al_NVE_1600Kstart.traj)
最初はきれいにfccの結晶構造に配置されたAl原子ですが、計算が進むにつれ徐々に構造が乱れ、段々とシミュレーションセルの外側に拡散していく様子がうかがえます。このようにMDシミュレーションでは実際に原子が時間発展とともにどのような軌道を描くか詳細に観察することが可能です。
アルミニウムの溶融シミュレーション¶
それでは以下に、上記のfcc-Alの溶融過程を再現するのに用いたMDシミュレーションのサンプルコードを示します。今回は計算の高速化のためASAP3のEMT力場を用いています。（ユーザーのpython環境でこの力場を用いる際は始めにpip install asap3を実行してasap3のパッケージをインストールしてください。）以下の計算は100 psec実行しますが、数十秒ほどで計算が完了します。
[2]:
%%time
import os
from asap3 import EMT
calculator = EMT()
from ase.build import bulk
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution,Stationary
from ase.md.verlet import VelocityVerlet
from ase.md import MDLogger
from ase import units
from time import perf_counter
import numpy as np
# Set up a fcc-Al crystal
atoms = bulk("Al","fcc",a=4.3,cubic=True)
atoms.pbc = True
atoms *= 3
print("atoms = ",atoms)
# Set calculator (EMT in this case)
atoms.calc = calculator
# input parameters
time_step = 1.0 # MD step size in fsec
temperature = 1600 # Temperature in Kelvin
num_md_steps = 100000 # Total number of MD steps
num_interval = 1000 # Print out interval for .log and .traj
# 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
# Set output filenames
output_filename = "./output/ch6/liquid-Al_NVE_1.0fs_test"
log_filename = output_filename + ".log"
print("log_filename = ",log_filename)
traj_filename = output_filename + ".traj"
print("traj_filename = ",traj_filename)
# Remove old files if they exist
if os.path.exists(log_filename): os.remove(log_filename)
if os.path.exists(traj_filename): os.remove(traj_filename)
# Define the MD dynamics class object
dyn = VelocityVerlet(atoms,
time_step * units.fs,
trajectory = traj_filename,
loginterval=num_interval
)
# Print statements
def print_dyn():
imd = dyn.get_number_of_steps()
time_md = time_step*imd
etot = atoms.get_total_energy()
ekin = atoms.get_kinetic_energy()
epot = atoms.get_potential_energy()
temp_K = atoms.get_temperature()
print(f" {imd: >3} {etot:.9f} {ekin:.9f} {epot:.9f} {temp_K:.2f}")
dyn.attach(print_dyn, interval=num_interval)
# Set MD logger
dyn.attach(MDLogger(dyn, atoms, log_filename, header=True, stress=False,peratom=False, mode="w"), interval=num_interval)
# Now run MD simulation
print(f"\n imd Etot(eV) Ekin(eV) Epot(eV) T(K)")
dyn.run(num_md_steps)
print("\nNormal termination of the MD run!")
atoms = Atoms(symbols='Al108', pbc=True, cell=[12.899999999999999, 12.899999999999999, 12.899999999999999])
log_filename = ./output/ch6/liquid-Al_NVE_1.0fs_test.log
traj_filename = ./output/ch6/liquid-Al_NVE_1.0fs_test.traj
imd Etot(eV) Ekin(eV) Epot(eV) T(K)
0 32.139701294 22.336120234 9.803581060 1600.00
1000 32.144889111 9.851785137 22.293103974 705.71
2000 32.144634514 10.538523047 21.606111466 754.90
3000 32.144643608 11.205207875 20.939435733 802.66
4000 32.144360371 10.628376834 21.515983537 761.34
5000 32.144745056 10.642859695 21.501885362 762.38
6000 32.144964346 9.673146214 22.471818131 692.92
7000 32.145003412 10.329285833 21.815717579 739.92
8000 32.144790609 10.798458531 21.346332078 773.52
9000 32.145198824 9.717769786 22.427429038 696.11
10000 32.144979661 9.871641249 22.273338412 707.13
11000 32.144876202 10.343316303 21.801559899 740.92
12000 32.144993729 9.071014781 23.073978949 649.78
13000 32.144884399 9.615342610 22.529541788 688.77
14000 32.144833389 10.027070827 22.117762562 718.27
15000 32.144621196 10.539268847 21.605352348 754.96
16000 32.144626765 11.053743640 21.090883125 791.81
17000 32.144344782 10.994105364 21.150239418 787.54
18000 32.144140960 11.682574217 20.461566743 836.86
19000 32.144493219 11.282753092 20.861740128 808.22
20000 32.144620243 11.021678911 21.122941332 789.51
21000 32.144572004 10.528263721 21.616308283 754.17
22000 32.144853422 10.071815410 22.073038012 721.47
23000 32.144606778 11.575417579 20.569189199 829.18
24000 32.144877568 10.336973331 21.807904236 740.47
25000 32.144795922 11.360248814 20.784547108 813.77
26000 32.144841118 9.949809548 22.195031571 712.73
27000 32.144434541 10.326124400 21.818310141 739.69
28000 32.144740173 10.553984691 21.590755481 756.01
29000 32.144825049 10.503967331 21.640857718 752.43
30000 32.144716217 11.174413776 20.970302440 800.46
31000 32.145292268 9.117849480 23.027442787 653.14
32000 32.144760180 10.088349826 22.056410354 722.66
33000 32.144459801 11.068906755 21.075553046 792.90
34000 32.144621995 10.171031207 21.973590788 728.58
35000 32.144700888 11.363733751 20.780967137 814.02
36000 32.144550709 10.653176904 21.491373805 763.12
37000 32.144756902 10.381793060 21.762963842 743.68
38000 32.144730151 10.257004183 21.887725967 734.74
39000 32.144781327 10.380126523 21.764654804 743.56
40000 32.144900626 9.716120395 22.428780230 695.99
41000 32.144277311 11.816934913 20.327342398 846.48
42000 32.144732214 10.483903392 21.660828822 750.99
43000 32.144234909 12.011591863 20.132643045 860.42
44000 32.144975889 10.557127683 21.587848207 756.24
45000 32.144777993 11.268795111 20.875982882 807.22
46000 32.144417779 11.242316416 20.902101363 805.32
47000 32.144748539 11.097972071 21.046776468 794.98
48000 32.144715016 10.685419887 21.459295128 765.43
49000 32.145091005 9.685290825 22.459800180 693.79
50000 32.144660151 11.184341544 20.960318608 801.17
51000 32.144902738 10.747928937 21.396973802 769.90
52000 32.145163772 10.593081930 21.552081842 758.81
53000 32.144823393 11.091951675 21.052871718 794.55
54000 32.145110651 10.543491337 21.601619315 755.26
55000 32.145089909 10.704306098 21.440783811 766.78
56000 32.144971397 11.003250775 21.141720621 788.19
57000 32.144961685 11.822987403 20.321974282 846.91
58000 32.145234068 10.566794373 21.578439695 756.93
59000 32.145382006 10.000078187 22.145303819 716.33
60000 32.145418293 9.846825578 22.298592715 705.36
61000 32.145188405 10.333914594 21.811273811 740.25
62000 32.145072699 10.727550676 21.417522023 768.45
63000 32.145171536 11.116251652 21.028919884 796.29
64000 32.145034329 10.619028790 21.526005538 760.67
65000 32.144953516 10.635215744 21.509737772 761.83
66000 32.145025745 10.145599740 21.999426006 726.76
67000 32.145164126 9.595965458 22.549198668 687.39
68000 32.145200107 9.649000052 22.496200055 691.19
69000 32.145283749 9.871641792 22.273641957 707.13
70000 32.144927551 10.609717496 21.535210055 760.00
71000 32.144951658 10.544159100 21.600792558 755.31
72000 32.145100877 10.154726831 21.990374046 727.41
73000 32.144707777 10.546472670 21.598235107 755.47
74000 32.144862218 10.316668026 21.828194192 739.01
75000 32.144806403 10.065312180 22.079494223 721.01
76000 32.145314242 10.206742428 21.938571815 731.14
77000 32.144580203 10.995857023 21.148723180 787.66
78000 32.145020684 9.931523677 22.213497007 711.42
79000 32.145054354 9.985611439 22.159442915 715.30
80000 32.144343492 11.346276398 20.798067095 812.77
81000 32.144937808 10.538805485 21.606132323 754.92
82000 32.144803289 10.789040395 21.355762895 772.85
83000 32.144754963 10.774891056 21.369863907 771.84
84000 32.145431604 9.008525953 23.136905650 645.31
85000 32.144760194 10.057425630 22.087334564 720.44
86000 32.144698741 10.685036611 21.459662130 765.40
87000 32.144540251 10.640615691 21.503924560 762.22
88000 32.144090376 10.314197022 21.829893354 738.84
89000 32.144502722 11.396688453 20.747814270 816.38
90000 32.144564034 11.460520590 20.684043444 820.95
91000 32.144473663 10.635723413 21.508750250 761.87
92000 32.144908536 9.629389204 22.515519332 689.78
93000 32.144629395 9.741495751 22.403133644 697.81
94000 32.144635633 9.661071297 22.483564336 692.05
95000 32.144676379 10.231992764 21.912683615 732.95
96000 32.145286282 8.958930245 23.186356037 641.75
97000 32.144243927 10.629861599 21.514382329 761.45
98000 32.144669759 10.105523781 22.039145978 723.89
99000 32.144610378 9.758667508 22.385942870 699.04
100000 32.144478686 11.325155603 20.819323083 811.25
Normal termination of the MD run!
CPU times: user 15.3 s, sys: 102 ms, total: 15.4 s
Wall time: 16.9 s
プログラムの流れはスクリプト内のコメントを見ながら流れを追っていけばだいたい理解できるとは思いますが重要なポイントのみについて解説します。
（１）初速度分布
構造作成、計算パラメーターの設定が終わったのち、各原子の初速度は指定の温度に対応したMaxwell-Boltzmann分布に従う速度分布を与えます。上記コードでは、 MaxwellBoltzmannDistribution
を用いて行っています。この方法で与えた初速度には系全体の運動量の任意性があるため、このままでは系全体がドリフトする可能性があります。なのでMaxwell-BoltzmannDistributionで初速度を与えた後にStationaryという関数で系全体の運動量をゼロとし、質量中心の座標を固定します。NVEの事例だけでなく、この後出てくる温度制御や圧力制御を伴うシミュレーションでも重要になります。
（２）MDの実行
今回は、VelocityVerletというクラスを用いることで、速度ベルレ法によるMDを実行しています。 dyn.run(num_md_steps) の箇所で実際にMDの実行が行われています。 今回の設定では、time_step=1.0として、1stepあたり1fs = \(1 \times 10^{-15}\) secの時間幅で、num_md_steps=100000 step実行しています。
（3）計算結果の記録
このスクリプトにはprint_dynという標準出力（notebook上へのデータの出力）にエネルギーや温度の値を出力する関数があります。その他にもMDLoggerという関数があり、これは指定したLogファイルの中にエネルギー、温度、応力情報を出力します。独自に計算結果をファイル出力する機構を作っても良いですが、デフォルトで備わっているLoggerで大抵のアプリケーションに必要な情報は得られるので活用をお薦めします。
MD実行後、軌跡の可視化は以下のように行なえます。
[3]:
from ase.io import Trajectory
from pfcc_extras.visualize.view import view_ngl
traj = Trajectory(traj_filename)
view_ngl(traj)
[4]:
from pfcc_extras.visualize.povray import traj_to_apng
from IPython.display import Image
traj_to_apng(traj, f"output/Fig6-1_fcc-Al_NVE_1600Kstart.png", rotation="0x,0y,0z", clean=True, n_jobs=16)
# See Fig6-1a
#Image(url="output/Fig6-1_fcc-Al_NVE_1600Kstart.png")
[Parallel(n_jobs=16)]: Using backend ThreadingBackend with 16 concurrent workers.
[Parallel(n_jobs=16)]: Done 18 tasks | elapsed: 13.9s
[Parallel(n_jobs=16)]: Done 101 out of 101 | elapsed: 59.4s finished
NVEアンサンブルにおける物性値のヒストリー¶
各原子の速度（すなわち運動量）も分かるため各種エネルギーの履歴を追うことも出来ます。例えば、以下は上のfcc-Alの溶融の際の全エネルギー（Tot.E.）、ポテンシャルエネルギー（P.E.）、運動エネルギー（K.E.）、そして温度（Temp.）の時間プロファイルです。
[5]:
import pandas as pd
df = pd.read_csv(log_filename, delim_whitespace=True)
df
[5]:
| Time[ps] | Etot[eV] | Epot[eV] | Ekin[eV] | T[K] | |
|---|---|---|---|---|---|
| 0 | 0.0 | 32.140 | 9.804 | 22.336 | 1600.0 |
| 1 | 1.0 | 32.145 | 22.293 | 9.852 | 705.7 |
| 2 | 2.0 | 32.145 | 21.606 | 10.539 | 754.9 |
| 3 | 3.0 | 32.145 | 20.939 | 11.205 | 802.7 |
| 4 | 4.0 | 32.144 | 21.516 | 10.628 | 761.3 |
| ... | ... | ... | ... | ... | ... |
| 96 | 96.0 | 32.145 | 23.186 | 8.959 | 641.8 |
| 97 | 97.0 | 32.144 | 21.514 | 10.630 | 761.4 |
| 98 | 98.0 | 32.145 | 22.039 | 10.106 | 723.9 |
| 99 | 99.0 | 32.145 | 22.386 | 9.759 | 699.0 |
| 100 | 100.0 | 32.144 | 20.819 | 11.325 | 811.3 |
101 rows × 5 columns
[6]:
import numpy as np
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(10, 5))
#color = 'tab:grey'
ax1 = fig.add_subplot(4, 1, 1)
ax1.set_xticklabels([])
ax1.set_ylabel('Tot E (eV)')
ax1.set_ylim([31.,33.])
ax1.plot(df["Time[ps]"], df["Etot[eV]"], color="blue",alpha=0.5)
ax2 = fig.add_subplot(4, 1, 2)
ax2.set_xticklabels([])
ax2.set_ylabel('P.E. (eV)')
ax2.plot(df["Time[ps]"], df["Epot[eV]"], color="green",alpha=0.5)
ax3 = fig.add_subplot(4, 1, 3)
ax3.set_xticklabels([])
ax3.set_ylabel('K.E. (eV)')
ax3.plot(df["Time[ps]"], df["Ekin[eV]"], color="orange",alpha=0.5)
ax4 = fig.add_subplot(4, 1, 4)
ax4.set_xlabel('time (ps)')
ax4.set_ylabel('Temp. (K)')
ax4.plot(df["Time[ps]"], df["T[K]"], color="red",alpha=0.5)
ax4.set_ylim([500., 1700])
fig.suptitle("Fig.6-1b. Time evolution of total, potential, and kinetic energies, and temperature.", y=0)
#plt.savefig("6-1_liquid-Al_NVE_1.0fs_test_E_vs_t.png") # <- Use if saving to an image file is desired
plt.show()
1-5節 でも出てきましたが、系の全エネルギー \(E\) (Total energy) は、運動エネルギー \(K\) (Kinetic energy)と、ポテンシャルエネルギー \(V\) (Potential energy)で表されます。
このうち、原子の運動エネルギー\(K\)は古典力学の表式を用い以下のように計算できます。
ここで \(m_i, \mathbf{v}_i, \mathbf{p}_i\) はそれぞれ各原子の質量 (mass) 、速度 (velocity)および運動量 (momenta, \(\mathbf{p}=m\mathbf{v}\))です。ここで系の温度と運動エネルギーはほぼ同義で、以下のような関係で定義されます。
今回のケースでは初期温度を1600 Kと設定し、あとは古典力学の運動方程式に則って自然に系を時間発展させています。従って特に外力もないのでエネルギー保存則が守られ全エネルギーが数値誤差の範囲で一定に保たれます。そしてかなり早い段階で温度が低下し、運動エネルギーとポテンシャルエネルギーのそれぞれにエネルギー分配則が働き、最終的に温度は初期温度の半分程度に落ち着きます。自然界では実際に物質はそのようにふるまい、今回のような最もシンプルなMDシミュレーションではそういった状況を再現します。このような計算対象である系の状態分布をNVEアンサンブルと言います。（または小正準集団アンサンブル。NVEのNは原子数、Vは体積、Eは全エネルギーが一定であることを表します。）ちなみに、ここでいうアンサンブルとは統計力学の考え方で物質の状態分布のことを指します。
NVEアンサンブルにおけるMDシミュレーションの設定¶
さて、上記のNVEアンサンブルのMDシミュレーションでは初期温度や計算時間のような計算条件を設定する必要があります。それ以外にも比較的自明ではないパラメーターとしてMDのタイムステップを設定する必要があります。このMDのタイムステップはどのように設定したらよいでしょうか？
答えは求めている計算精度に依存します。NVEアンサンブルでは2階の常微分方程式である古典の運動方程式解くわけですが、この積分プロセスの精度がタイムステップの大きさで決まります。（積分手法の詳細については本節の最後に記載します。）タイムステップサイズを変えた時の全エネルギーがどのように変化するか見てみましょう。計算時間は仮に1 nsecとします。（古典MDなどで良く用いられる時間スケールです。）
Fig.6-1c. Time evolution of the total energy with respect to time step size.
上記はタイムステップを0.5から5.0 fsecまで変化させたときの計算結果です。本来一定であるはずの全エネルギーが、タイムステップが大きいときに時間の経過とともに大きくぶれてしまっていることが確認できます。
全エネルギーの時間依存性確認のため、上記のMDシミュレーションスクリプトのtime_step 引数を変えて実行し、そのエネルギープロファイルを比較してみてください。 可視化にあたっては、上記Fig.6-1b の可視化コードを参考にしてください。
上記の計算結果からMDのタイムステップサイズを変化し、横軸をステップサイズ、縦軸を全エネルギーの揺らぎから得られたRMSEをプロットすると以下のようになります。
Fig.6-1d. Total energy ERMSE as a function of MD step size.
縦軸は原子ひとつ辺りのエネルギーです。Log-logプロットで概ね直線上にのることが見て取れます。単位原子当たりの値にして、0.25 fsで1.2e-6 eVから5.0 fsで6e-5 fsまで変化することが分かります。この場合原子数は108原子なので、系全体で1e-4 eVから1e-3 eV変化しうることが確認できます。
エネルギーの観点では時間ステップを1 fsecとすると、1 nsecのシミュレーションを100原子系で0.5 meV程度、1000原子系だと約5 meV程度の誤差が予想されます。もちろん計算したい物性などによりますが、たいていの計算では十分な精度ではないでしょうか？これが例えば計算時間を稼ぎたいからとタイムステップを5 fsecとするとほぼ一桁誤差が大きくなります。今後MD計算を活用するうえで、この計算精度がどの程度かというのは頭の片隅に止めておく必要があります。
その他のアンサンブルとMDシミュレーション¶
本節の最後にNVEアンサンブル以外のMDシミュレーション法について触れておきます。
NVEアンサンブルのMDシミュレーションはシンプルかつパワフルな計算手法ですが、我々が一般に観察したい現象を再現するには不都合な点が多く存在します。 NVEのシミュレーションは計算対象のセルそれそのものが全てであり、外界というものが存在しません。しかし我々が実際に観察したいミクロな世界にはそれを取り巻く環境が存在します。外部環境の温度や圧力の影響により、計算対象の状態が変化しますが、これはNVEアンサンブルの範疇では再現が出来ません。
このため統計力学の世界ではcanonicalアンサンブル（NVT）、isothermal-isobaricアンサンブル（NPT）、等をはじめ様々な状態分布があり、それらを活用することで特定の条件下での計算が可能になります。次節以降ではこれらの統計力学アンサンブル下での計算方法について学びます。
参考文献¶
[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] 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
[3] 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#