分子動力学シミュレーション

TL;DR

• MDシミュレーションでは古典力学の運動方程式に則った計算方法で、全ての原子の位置と速度の時間発展を露わに観察できる。

• 時間スケールの制約があり一般的に数nsから数十nsぐらいが上限。その時間内に発現しない現象は別の計算手法を検討する必要がある。

• MDシミュレーションには再現したい状態に応じて様々な状態（アンサンブル）があり、最も単純なアンサンブルはNVEアンサンブルという。

• NVEアンサンブルでは全エネルギー保存則が成立し、計算精度および計算時間の観点から、時間ステップサイズを適切に設定する必要がある。

この章では系の時間発展をシミュレーションする、 分子動力学シミュレーション（Molecular dynamics, MD） について学びます。

MDシミュレーションでは個々の原子の軌道の時間発展を露わに取り扱うシミュレーションで、古典力学の運動方程式を積分することによって計算対象である原子の座標と速度を逐次計算していく手法です。この計算手法自体は原子間に働く力やエネルギーのモデルとは独立した理論で、古くから分子シミュレーションの分野で用いられてきました。従って、理論的な背景や事例等は多くの書籍や文献が多く存在するのでそれらを確認していただければと思います。（参考文献：[1-3]） 本チュートリアルではあくまでMatlantisを使って実践的にこれらの計算を実行するために必要な知識を習得するところが目的となっています。

早速ですが具体的な事例を通して、どんなことがMDシミュレーションで観察出来るかみてみましょう。

• https://wiki.fysik.dtu.dk/ase/tutorials/md/md.html

事前準備 - 必要なライブラリのインストール

本章ではところどころ、ASE上で簡易に利用できる古典力場であるASAP3のEMT力場を利用しています。古典力場であるため精度や用途はかなり限定的ですがシンプルかつ高速で実行可能なのでチュートリアルとして事例を示すには十分な機能を持っています。ちなみにASAP3-EMTで利用可能な元素はNi、Cu、Pd、Ag、Pt、Auとこれらの合金に限定されます。インストールの仕方は以下のとおりです。

!pip install --upgrade asap3

Looking in indexes: https://pypi.org/simple, http://pypi.artifact.svc:8080/simple
Requirement already satisfied: asap3 in /home/jovyan/.py39/lib/python3.9/site-packages (3.13.4)
Collecting asap3
  Downloading asap3-3.13.5.tar.gz (849 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 849.3/849.3 kB 46.4 MB/s eta 0:00:00
  Installing build dependencies ... done
  Getting requirements to build wheel ... done
  Installing backend dependencies ... done
  Preparing metadata (pyproject.toml) ... done
Requirement already satisfied: ase>=3.23.0 in /home/jovyan/.py39/lib/python3.9/site-packages (from asap3) (3.23.0)
Requirement already satisfied: numpy>=1.18.5 in /home/jovyan/.py39/lib/python3.9/site-packages (from ase>=3.23.0->asap3) (1.26.4)
Requirement already satisfied: scipy>=1.6.0 in /home/jovyan/.py39/lib/python3.9/site-packages (from ase>=3.23.0->asap3) (1.13.1)
Requirement already satisfied: matplotlib>=3.3.4 in /home/jovyan/.py39/lib/python3.9/site-packages (from ase>=3.23.0->asap3) (3.9.0)
Requirement already satisfied: contourpy>=1.0.1 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (1.2.1)
Requirement already satisfied: cycler>=0.10 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (4.53.0)
Requirement already satisfied: kiwisolver>=1.3.1 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (1.4.5)
Requirement already satisfied: packaging>=20.0 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (24.1)
Requirement already satisfied: pillow>=8 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (10.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (3.1.2)
Requirement already satisfied: python-dateutil>=2.7 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (2.9.0.post0)
Requirement already satisfied: importlib-resources>=3.2.0 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (6.4.0)
Requirement already satisfied: zipp>=3.1.0 in /home/jovyan/.py39/lib/python3.9/site-packages (from importlib-resources>=3.2.0->matplotlib>=3.3.4->ase>=3.23.0->asap3) (3.19.2)
Requirement already satisfied: six>=1.5 in /home/jovyan/.py39/lib/python3.9/site-packages (from python-dateutil>=2.7->matplotlib>=3.3.4->ase>=3.23.0->asap3) (1.16.0)
Building wheels for collected packages: asap3
  Building wheel for asap3 (pyproject.toml) ... done
  Created wheel for asap3: filename=asap3-3.13.5-cp39-cp39-linux_x86_64.whl size=4432248 sha256=31f8e8a60ef879b55c6cd93d7518f8fd0183a6632500d1a2f34ac3976ea5cfe3
  Stored in directory: /home/jovyan/.cache/pip/wheels/86/c6/6c/ad1816075e4617db75293ef31ab3e288466842bcc0727e2fdc
Successfully built asap3
Installing collected packages: asap3
  Attempting uninstall: asap3
    Found existing installation: asap3 3.13.4
    Uninstalling asap3-3.13.4:
      Successfully uninstalled asap3-3.13.4
Successfully installed asap3-3.13.5

[notice] A new release of pip is available: 24.2 -> 24.3.1
[notice] To update, run: pip install --upgrade pip

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実行しますが、数十秒ほどで計算が完了します。

%%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.144732392     9.717662136    22.427070256   696.10
   2000     32.144595676     10.384753227    21.759842449   743.89
   3000     32.144485776     10.292545081    21.851940695   737.28
   4000     32.144656501     9.569029340    22.575627161   685.46
   5000     32.144223666     11.260312865    20.883910801   806.61
   6000     32.144545079     10.854854065    21.289691013   777.56
   7000     32.144674572     9.749442021    22.395232551   698.38
   8000     32.144463479     10.658602272    21.485861207   763.51
   9000     32.144319804     10.387294028    21.757025776   744.07
   10000     32.144629878     10.025983768    22.118646110   718.19
   11000     32.144361501     11.236237078    20.908124423   804.88
   12000     32.144282709     11.138136472    21.006146236   797.86
   13000     32.144353906     10.941409697    21.202944209   783.76
   14000     32.144429188     10.722515666    21.421913523   768.08
   15000     32.144375555     10.939219269    21.205156286   783.61
   16000     32.144592848     9.409488230    22.735104618   674.03
   17000     32.144024018     11.224044353    20.919979664   804.01
   18000     32.144368456     10.108417246    22.035951209   724.09
   19000     32.144305600     10.668868800    21.475436800   764.24
   20000     32.144225902     10.733945978    21.410279924   768.90
   21000     32.144479804     9.741119218    22.403360586   697.78
   22000     32.144046753     10.898466211    21.245580542   780.69
   23000     32.144622226     9.356067360    22.788554866   670.20
   24000     32.144735318     8.279816397    23.864918921   593.11
   25000     32.144262614     10.770595178    21.373667437   771.53
   26000     32.144587780     10.012490263    22.132097518   717.22
   27000     32.144457659     10.548355368    21.596102291   755.61
   28000     32.144223772     11.000952065    21.143271706   788.03
   29000     32.144108672     10.349540962    21.794567710   741.37
   30000     32.144032797     10.306106856    21.837925940   738.26
   31000     32.144272984     10.845360184    21.298912801   776.88
   32000     32.144080434     11.683899064    20.460181371   836.95
   33000     32.144326494     10.522818948    21.621507546   753.78
   34000     32.144349760     10.829084255    21.315265504   775.72
   35000     32.144530782     9.748968613    22.395562168   698.35
   36000     32.144264785     11.057956331    21.086308454   792.11
   37000     32.144231508     10.737533644    21.406697864   769.16
   38000     32.143824362     11.202989634    20.940834728   802.50
   39000     32.144470371     9.805347796    22.339122575   702.39
   40000     32.144075224     11.038327085    21.105748139   790.71
   41000     32.144660333     10.122700133    22.021960200   725.12
   42000     32.144399470     10.471717961    21.672681509   750.12
   43000     32.144301319     10.889584913    21.254716406   780.05
   44000     32.144712144     10.040503091    22.104209052   719.23
   45000     32.144470171     10.335349862    21.809120309   740.35
   46000     32.144498444     10.831459795    21.313038648   775.89
   47000     32.144480014     10.745521376    21.398958638   769.73
   48000     32.144208898     10.811019456    21.333189442   774.42
   49000     32.144394973     10.520276617    21.624118357   753.60
   50000     32.143711391     11.693195391    20.450516000   837.62
   51000     32.144216072     10.794231584    21.349984488   773.22
   52000     32.143954477     11.167702465    20.976252011   799.97
   53000     32.144197026     10.761047389    21.383149636   770.84
   54000     32.144692097     10.650248738    21.494443358   762.91
   55000     32.144292435     11.069017237    21.075275198   792.91
   56000     32.144451534     10.818961116    21.325490419   774.99
   57000     32.144159985     10.436902720    21.707257265   747.63
   58000     32.144409863     10.518871169    21.625538694   753.50
   59000     32.144290785     10.147827359    21.996463427   726.92
   60000     32.144252442     10.561857643    21.582394799   756.58
   61000     32.144211404     11.237046923    20.907164481   804.94
   62000     32.144335612     10.234587556    21.909748055   733.13
   63000     32.144553693     9.734651928    22.409901765   697.32
   64000     32.144056336     11.331668716    20.812387620   811.72
   65000     32.144544495     10.073495649    22.071048846   721.59
   66000     32.144261257     11.553758409    20.590502848   827.63
   67000     32.144400715     10.388117319    21.756283396   744.13
   68000     32.144180052     11.049408118    21.094771934   791.50
   69000     32.144331305     10.787515820    21.356815485   772.74
   70000     32.144363121     9.757896544    22.386466577   698.99
   71000     32.143995748     12.113379277    20.030616471   867.72
   72000     32.144360474     9.849801292    22.294559182   705.57
   73000     32.144573361     10.621878850    21.522694511   760.88
   74000     32.144018243     11.163225758    20.980792485   799.65
   75000     32.144376982     10.221662849    21.922714133   732.21
   76000     32.144588102     9.451493371    22.693094731   677.04
   77000     32.144417191     10.553128007    21.591289185   755.95
   78000     32.144597951     10.231912256    21.912685695   732.94
   79000     32.144830327     9.925136042    22.219694285   710.97
   80000     32.143996211     10.923923218    21.220072993   782.51
   81000     32.144574158     9.802931001    22.341643157   702.21
   82000     32.144521087     10.157138062    21.987383025   727.58
   83000     32.144238473     10.784705296    21.359533177   772.54
   84000     32.144152374     11.426586870    20.717565505   818.52
   85000     32.144128400     10.351754495    21.792373904   741.53
   86000     32.144538231     10.954848440    21.189689791   784.73
   87000     32.144430968     10.949426732    21.195004236   784.34
   88000     32.144274375     10.209803957    21.934470418   731.36
   89000     32.144706378     9.682925565    22.461780814   693.62
   90000     32.144748840     10.276817862    21.867930978   736.16
   91000     32.144569186     10.412884902    21.731684284   745.90
   92000     32.144353339     10.693038177    21.451315163   765.97
   93000     32.145033728     9.767317058    22.377716670   699.66
   94000     32.144224948     11.741172068    20.403052880   841.05
   95000     32.144571344     10.235993503    21.908577841   733.23
   96000     32.144506261     11.768195101    20.376311159   842.99
   97000     32.144586363     10.502267187    21.642319176   752.31
   98000     32.144658691     10.422999129    21.721659562   746.63
   99000     32.144269295     10.908635273    21.235634022   781.42
   100000     32.144367925     10.044230000    22.100137925   719.50

Normal termination of the MD run!
CPU times: user 44.9 s, sys: 355 ms, total: 45.3 s
Wall time: 1min 13s

プログラムの流れはスクリプト内のコメントを見ながら流れを追っていけばだいたい理解できるとは思いますが重要なポイントのみについて解説します。

（１）初速度分布

構造作成、計算パラメーターの設定が終わったのち、各原子の初速度は指定の温度に対応した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実行後、軌跡の可視化は以下のように行なえます。

from ase.io import Trajectory
from pfcc_extras.visualize.view import view_ngl


traj = Trajectory(traj_filename)
view_ngl(traj)

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:   27.4s
[Parallel(n_jobs=16)]: Done 101 out of 101 | elapsed:  1.6min finished

NVEアンサンブルにおける物性値のヒストリー

各原子の速度（すなわち運動量）も分かるため各種エネルギーの履歴を追うことも出来ます。例えば、以下は上のfcc-Alの溶融の際の全エネルギー（Tot.E.）、ポテンシャルエネルギー（P.E.）、運動エネルギー（K.E.）、そして温度（Temp.）の時間プロファイルです。

import pandas as pd

df = pd.read_csv(log_filename, delim_whitespace=True)
df

/tmp/ipykernel_49555/510401485.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(log_filename, delim_whitespace=True)

	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.427	9.718	696.1
2	2.0	32.145	21.760	10.385	743.9
3	3.0	32.144	21.852	10.293	737.3
4	4.0	32.145	22.576	9.569	685.5
...	...	...	...	...	...
96	96.0	32.145	20.376	11.768	843.0
97	97.0	32.145	21.642	10.502	752.3
98	98.0	32.145	21.722	10.423	746.6
99	99.0	32.144	21.236	10.909	781.4
100	100.0	32.144	22.100	10.044	719.5

101 rows × 5 columns

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)で表されます。

\[E = K + V\]

このうち、原子の運動エネルギー\(K\)は古典力学の表式を用い以下のように計算できます。

\[K = \sum_{i=1}^{N} \frac{1}{2} m_i {\mathbf{v}}_i^2 = \sum_{i=1}^{N} \frac{{\mathbf{p}}_i^2}{2 m_i}\]

ここで \(m_i, \mathbf{v}_i, \mathbf{p}_i\) はそれぞれ各原子の質量 (mass) 、速度 (velocity)および運動量 (momenta, \(\mathbf{p}=m\mathbf{v}\))です。ここで系の温度と運動エネルギーはほぼ同義で、以下のような関係で定義されます。

\[K = \frac{3}{2} k_B T\]

今回のケースでは初期温度を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#