分子動力学法入門
能勢修一
慶應義塾大学理工学部・大学院理工学研究科

i
はじめに
この原稿は、分子動力学法入門のテキストとして準備していたものの中で、プログラミ
ング実習の部です。このテキストは未完成で、図、表、はついていませんし、未完のコラ
ム等への引用を含んでいます。使用の場合、この点に注意しながら、利用してください。
プログラミング実習の部は、４年卒業研究で分子動力学法を用いる場合に利用していました。この部分で焼く２ヶ月でやりとげ、より高度のシミュレーションへ進むことができます。

目次
第 1 章 粒子配置をつくる

1

第 2 章 粒子間距離の分布を調べる。

5

第 3 章 Verlet 法を用いて調和振動子の運動方程式を積分する。

9

第 4 章 予測子・修正子法を用いて調和振動子の運動方程式を積分する。

18

第 5 章 レナード・ジョーンズ相互作用で相互作用する２つの粒子の系の運動

23

第 6 章 多粒子系のシミュレーションを行なう。

26

1

1

第 1 章 粒子配置をつくる
MD シミュレーションでは、多数の粒子の運動をニュートンの運動方程式を数値的に積分することによりコンピューターの中で再現していく。よく知られているように、古典力学
では、ある時刻 t での粒子配置（N 粒子すべての座標 r 1 (t), r 2 (t), · · · , r N (t) ）および粒子の速度（ṙ 1 (t), ṙ 2 (t), · · · , ṙ N (t)）を定めるとその後の各粒子の運動は一通りに決まってしまう（この性質より力学法則は決定論的であるという）。
シミュレーションを始めるためには、まず出発点となる粒子の配置をつくることが必要となる。日常生活で接する物質は、アボガドロ数程度の非常に多数の粒子の集合体であるが、我々がシミュレーションで調べることができるのは、コンピューターの速度、記憶容量などの制限のため、通常数百ないし数万個程度である（最近、百万個を越える例もい
くつか報告されている）。統計力学的なシミュレーションの目的は、非常に多くの粒子の集まった巨視的な体系（統計力学では密度を一定に保ったまま粒子数を無限大にする極限（熱力学的極限と呼ぶ）を考える）の示す性質を、限られた数の粒子を用いて微視的な観点から精度よく推測することである。
ある限られた大きさの領域（これを MD セルと呼ぶ。平行６面体の形に選ぶことが多い。）の中での粒子の運動を考える。粒子配置はこれから調べようとする状態にふさわし
いものにしておくことが望ましいが、液体状態のように乱れた配置をいきなり作り出すことは難しい。シミュレーションを用いて正しい巨視的性質を得るためには、統計平均を行
いいろいろの物理量を求めるシミュレーションの本体部分を行う前に、調べようとする条件に合致する状態を作り出すこと、そしてそれが平衡状態となっていることを確認するという準備段階が必須である。液体状態の粒子配置もこの準備段階でつくるとよい。
粒子の重なり合いを避け、ある一定密度の粒子配置をつくるには、規則的に粒子が並んだ結晶配置が便利である。結晶では、ある基本的な単位（ユニットセルまたは単位胞と呼ぶ）の中での粒子配置を（３次元の場合）３個の方向へ周期的に並べていくことにより結晶全体の粒子配置がつくられる。
ユニットセルの形を a1 , a2 , a3 の３個のベクトルを辺とする平行六面体とする。このユニットセル内の粒子の座標を ri (i = 1, · · · , n) とすると完全結晶（粒子の配置に乱れがな

第1章

2

粒子配置をつくる

く、粒子が運動を止める温度 T = 0K の状態に対応する）の場合、すべての粒子の座標は

r i + m1 a1 + m2 a2 + m3 a3(m1 , m2 , m3 = 0, ±1, ±2, · · · , ±∞),(i = 1, · · · , n)(1.1)

の形で表される。
課題 1.1
ユニットセルを定めるベクトル a1 = (a1x , a1y , a1z ), a2 , a3 およびユニットセルの中での粒子の座標 ri = (xi , yi , zi ) (i = 1, · · · , n) が与えられたとき a1 , a2 , a3 方向へそれぞれ

m1 , m2 , m3 個のユニットセルが並んだ粒子配置をつくり大きさ (3, N ) の配列 R(i, j) に書き込むプログラムをつくりなさい。ただし、i = 1 は X 座標、i = 2 は Y 座標、i = 3 は

Z 座標、j は粒子の番号を示す。また、全粒子数 N は N = n × m1 × m2 × m3 である。
注意点： このプログラムの作成にはユニットセル内に含まれる粒子数 n、ユニットセルの数 m1 , m2 , m3 の値がパラメターとして必要だが、プログラムの主要部分はこのパラメターの値に依存しない形で書くことが望ましい。パラメターの値は parameter 文で定義するか、プログラムの最初の部分で入力データとして読み込むようにする。ここでは、プログラムの変更なしに種々のテストができるようにパラメターを入力データとして読み込む形でつくってみなさい。この場合、配列の大きさを決めるためパラメターの値の取り得る範囲についてあらかじめ定めておく必要がある。最初のテストとしてはそれぞれの最大値を４としておけば十分である。このとき粒子数 N の最大値は 256 となる。
座標の値をどの単位系を用いて表すかは重要な問題である。コンピューターで取り扱うことのできるのは数値のみであり、次元を持った物理量を直接扱うことはできない。コン
ピューターで扱うデータはすべてある単位系を用いて規格化し無次元量とした数値と考えねばならない。この時、MKS 単位系などの実スケールの表示を用いると原子レベルでは数値が大変小さな量となる。例えば、1A=10−10 m である。このような量を用いると、例
えば高次の巾の計算の際、計算機による実数表示の限界を越え困難を生ずることがある（コラム

LJ 系がブラックホールに？ 参照）。これを避けるためには、原子レベルの量を単位にして表示するとよい。長さについては、A を単位に用いるかその体系の代表的な長さ（レナード・ジョーンズ相互作用の粒子半径 σ など）で規格化して表すとよい。

課題 1.2単純立方格子、体心立方格子、面心立方格子のいずれの構造もつくることができるプ
ログラムをつくりなさい。MD セルの形は３辺が互いに直行している直方体の形に選ぶ。

第1章

3

粒子配置をつくる

単純立方格子は、ユニットセルを辺の長さ１の立方体としたとき、この中に１つだけ粒子が含まれるものである。この座標を (1/2, 1/2, 1/2) とする。体心立方格子は立方体内に２つの粒子が含まれ立方体の頂点および重心に粒子が位置する。体心立方格子の単位胞内の粒子位置としては、(1/4, 1/4, 1/4) と (3/4, 3/4, 3/4) を選ぶとよい。面心立方格子では立方体に４つの粒子が含まれ、立方体の頂点および面の中心に粒子が位置する。粒子の位置
は (1/4, 1/4, 1/4)、(1/4, 3/4, 3/4)、(3/4, 1/4, 3/4)、(3/4, 3/4, 1/4) と選ぶとよい。
注意点：

格子の種類の選択はパラメターの入力で決めるとよい。上の結晶構造のユニットセル内の粒子配置は、MD セルの境界面上に粒子が並ばないように、全体を平行移
動して粒子が MD セルの内部に配置されように選んでいる。境界上に粒子が位置すると粒子を運動させた時、すぐに MD セルから出てしまう。これを避けておくほうが、後程述べるプログラムのチェックのときに有利である。

課題 1.3粒子が完全結晶の位置より少し乱れた配置をつくる。課題 1-2 でつくった粒子配置を元に結晶格子点の近傍で粒子を少し移動させ乱れた配置をつくる。これは温度が有限で熱揺らぎがある時の結晶の構造に対応している。
注意点： 通常コンピューターのシステムに組み込まれている一様乱数のルーチンを用
い、変域が −1 ≤ ξ ≤ 1 の乱数 ξ を逐次生成していく。各座標ごとに、異なる乱数を用いて

xi = xi + ξδ

(1.2)

の様に、完全結晶の格子点の座標 xi を乱れた値 xi に変えていく。ここでは、変位の最大幅を δ としている。この値は最近接粒子間距離の 10 % 以下の値にするとよい。乱数の変域が 0 ≤ ξ  ≤ 1 の場合には、ξ = 2ξ  − 1 とするとよい。コンピュータで使用する乱数は通常擬似乱数とよばれるものである。これが本当に乱数としての条件を満たしているかど
うか十分な吟味が必要なこともあるが（乱数を多く用いるモンテカルロシミュレーション
では特に重要）、MD 法の初期配置をつくる場合には、乱数の良否はそれほど問題にならない。

擬似乱数の特徴は初期値となる入力データに同じものを用いると全く同じ乱数列をもう一度再現することができることである。念のため、乱数の初期値の値を記録しておくとよい。計算に問題が生じた時、再実行により誤りを探すことができる。

第1章

粒子配置をつくる

4

整数型演算のオーバーフローをチェックしないコンピューターでは、32bit の整数型変数を用いて、次のような文定義関数を用いてもよい。この時乱数列の周期は 232 となる。
これは、混合合同法と呼ばれる擬似乱数の生成法である。
プログラムの先頭でIR(IY)=513*IY+495485483
R(IY)=DFLOAT(IY)/2147483648.D0の２つの文関数を定義しておく。関数 IR(IY) では、積および和の演算により -2147483648
から +2147483647 までの整数値が得られる。IR(IY) の定義のうち、和の数 495485483 は
奇数であれば何に置き換えてもほとんど影響ないが、積の数 513 は重要であり、この値により乱数の良否は大きく変わる。ここで用いた 513 は著者が頻度テストなどの簡単なチェックを行ったものである。これを 129 や 257 に置き換えると乱数としては不満足なものとなる。関数 R(IY) により、整数 IY を -1.0 から 1.0 の間の値の実数に変換している。
この文定義関数を用いる場合は、まず乱数の初期値として整変数 IY の値を入力する。
乱数を生成するときに、毎回 IY の値を IY=IR(IY) と更新してから関数 R(IY) を用いる。

IY=IR(IY)
X=R(IY)
これで、X は −1 ≤ x ≤ 1 の一様乱数となる。
乱数には、正規乱数を用いる方が望ましいと考える人があるが（特に速度の初期値をマックスウェル分布にするため）、MD シミュレーションでは望ましい条件の状態を作り出すために必ず準備段階が必要なことを考えると、一様乱数を用いて特に困ることはない。
逆に正規乱数を用いた場合には不利な点が２つある。それは、正規乱数の計算には一様乱数より少し余分に計算時間がかかることと、まれではあるが 4σ （σ は標準偏差）を越えるような非常に大きな値が現れる可能性（約１万分の１）があり、変位が予想される限界を越えて粒子が重なり合ってしまう等の困難を起こしてしまう場合があることである。

5

第2章

粒子間距離の分布を調べる。

いろいろな状態で、粒子配置がどのようになっているかを特徴づける代表的な量が動径分
布関数 g(r) である。これは、ある粒子から距離 r 離れた点に他の粒子がどのくらいいるのかを表すものである。より正確には、r と r + dr で挟まれた区間にいる他の粒子の数が ρ4πr2 g(r)dr に等しくなるように定義されている。球殻の体積が 4πr2 dr であり、系全体の密度が ρ であるから、粒子が全く独立に一様に分布するとき g(r) = 1 となるように規格化されている。g(r) の特徴を図??

に示している。結晶では、特定の距離に鋭い山を持つ。この山の形は、絶対 0 度では δ 関数的だが、有限温度では熱振動のために幅を持つようになる。液体では、中心の粒子と重なるため他の粒子は r = 0 付近にくることができないために r が小さな場合 g(r) = 0 となる。最近接距離付近に鋭い山を持つが、それより遠くになるにつれ構造の特徴がだんだん消え一様分布に近づいていく。気体では、近距離の部分を除き、すぐ一様分布となってしまう。
 rn(r) =
ρ 4πr2 g(r)dr

(2.1)

0

は、距離 r までにいる他の粒子の数を表しており、これもぜひ計算しておくべき量である。

g(r) は、X 線回折、中性子回折などの実験により求められ、実験とシミュレーションの直接の比較ができる重要な量となっている。
課題 2.1
課題 1-2 または 1-3 で作成した粒子配置より動径分布関数および n(r) を求めるプログラ
ムをつくれ。また、プリンターを用いて、g(r) を描いてみよ。

注意点： g(r) を求めるには距離が r と r + dr の区間に含まれる粒子対の数 ρ4πr2 g(r)drを計算するとよい。g(r) を計算する最大の距離 R （通常は相互作用のカットオフの距離

rc に選ぶ）までの区間を刻み幅 dr で区分けする。そして、この区間に入る粒子間距離が計算中に何回現れたかを数え上げる。区間の数より大きな整数型の一次元配列（例えば、

NG(0:1000) とする）を用意し、初期化でその値をすべて 0 にしておく。粒子対の距離 rijを求めた時点で NN=INT(rij /dr) は (NN-1)*dr と NN*dr の間の区間に NN というラ

第2章

粒子間距離の分布を調べる。

6

ベルを付けたことになる。そこで、NG(NN)=NG(NN)+1 と NG の NN 要素に 1 を加えておく。区間のラベルを付けるとき NN=INT(rij /dr + 0.5) とすると g(r) 等の計算を区間の中心の値 NN*dr で代表することになる。すべての計算終了後、NG(NN) の値を

N/2*NT で割ると ANG(NN)=4πr2 g(r)dr が求まる（N は粒子数、NT は計算に用いた
粒子配置の数である）。ここでは粒子対の距離の計算で２重に計算していないとする。NT
でも割るのは、異なる粒子配置についての平均より構造を求めるときに適用できるように記述しているからである。一つの粒子配置に限っている場合には NT=1 である。累積粒
子数 n(r) は ANG(NN) を距離 0 に対応する区間から順次加えていくと求まる。g(r) はANG(NN)/4πr2 drρ で求められる。
g(r) を計算するのは、物理的に意味のある範囲に限らなければならない。相互作用のカットオフを用いる場合は、通常、これ以上の距離の粒子間距離を求めないので、この距離が最大値となる。周期境界条件を用い、かつ minimum image convention を用いる場合（この後で説明する）には、MD セルの辺の長さの最小値の半分（rh とする）を越えた距離については、粒子対が正しく計算されない。この場合にも、rh を越えた距離についての g(r) は正しい値とはならないので、計算してはいけない。
完全結晶の場合には、ある１つの粒子を選び出し、そのまわりの他の粒子の距離を調べてもよいが、実際のシミュレーションでは、粒子の位置は乱れており、また時々刻々変って行く。このため、すべての粒子について、そのまわりの様子を調べ、g(r) を計算する必要がある。
すべての粒子対について考慮するには、周期境界条件を用いなければならない。これを用いず、壁に囲まれた領域を考えると、壁に近い粒子では、壁の方向には他の粒子がいないため粒子に取り囲まれた巨視的状況からのずれが大きくなる。周期境界条件を用いると、同じ粒子配置を持つ MD セルが周期的に配列しているのと同じことになる。ある粒子 i と他の粒子 j の対で考える場合、同じ MD セル内にいる粒子 j だけでなく、粒子j の他の MD セルに含まれる鏡像 j  も考慮しなければならない。場合によっては、i − jの対の距離より i − j  の対の距離の方が短いこともある。g(r) を計算する最大距離 R を

MD セルの最短の辺の長さの半分より小さくする場合には、j および j の鏡像のうち、 i
から見て距離 R の範囲に入るものは高々１つとなる（これを minimum image conventionという）。この時、鏡像の中から最短距離にあるものを選び出すには次のようにするとよい。x 方向の MD セルの辺の長さを Lx とする。
（同様に、Ly 、Lz を定める）。ここでは、

x 方向のみについて示す。辺の長さ Lx を変数 CLX、その半分を CLXH=CLX*0.5d0 とする。プログラムでは次のようにする。

第2章

粒子間距離の分布を調べる。

7

DX=R(1,I)-R(1,J)
IF(DX.GT.CLXH) DX=DX-CLXIF(DX.LT.-CLXH) DX=DX+CLX
同一の MD セルに含まれる粒子の間の座標の差 DX が、Lx /2 より大きい（−LX /2 より小さい）ときには x の負（正）方向に隣あう MD セルに含まれる鏡像の方が近くなる。
それを判定しているのが２つの IF 文である。この操作は、FORTRAN の組み込み関数

INT(X) や MOD(X,1) を用いても行うことができる。
DX=DX-CLX*INT(DX/CLXH)
この場合には、判断が含まれないためベクトル演算が導入された当初はベクトル化が可能となり IF 文を含む場合より優れていた。現在では、この程度の判断では同様のベクトル化が可能である。実際に演算を行うと IF 文を含む方が速いようである。これは、IF の判定の後実際にその後の操作を行う確率が 1/4 であることによる。
簡単に g(r) を図示するには通常の文字入力を利用するとよい。１行に出力可能な文字数により図の精度が決まる。例えば 40 文字を用い g(r) の値を 0.1 刻みで表す場合は次のようになる。

CHARACTER*1 A(40)/40*’ ’/
CHARACTER*1 ASTER/’*’/,ONE/’1’/,ZERO/’0’/,BLANK/’ ’/g(r) の値を GR という変数で表すことにすると
DO 200 J=1,NR
GR=G(J)
IK=INT( GR*10.0D0 + 1.5D )IF(IK.GT.40) IK=40
A(1)=ZERO
A(11)=ONE
A(IK)=ASTER
WRITE(*,110) X(J),G(J),N(J),A

第2章

粒子間距離の分布を調べる。

8110 FORMAT(1H ,3F6.2,3X,40A1)
A(IK)=BLANK
200 CONTINUE
注意点：

WRITE 文で距離（X(I)）等の値を書くとき書式（FORMAT）文を用いて、固定小数点表示で出力しないと、その後に出力する g(r) のグラフが行ごとにずれてしまう。上のプログラム例で 10.0D0 をかけているのは g(r) の値を 0.1 の精度で出力（0.1 は10 の逆数）することに対応している。1.5d0 を加えているのは g(r) の値が 0.95 から 1.05の時 1 になるようにし、また１カラム目を 0 に、11 カラム目が 1 に対応するようにしているからである。g(r) = 1 の線を引くために 11 カラム目に文字の 1 を出力している。この程度のグラフでも概形を知るためには十分である。もっと精度のよいグラフを描くにはグラフィクスを用いる必要がある。

9第 3 章 Verlet 法を用いて調和振動子の運動方程式を積分する。
MD 法では古典力学の運動方程式を数値的に積分して粒子の運動を調べているが、全エネルギーが保存される保存系を扱う場合が多い。この時、エネルギー保存則が満たされているかどうかを見ることにより計算が正しく行われているかどうかの判定をすることができる。これは、MD 法の利点の一つであり、これを充分に活用することにより、プログラムの誤りを見つけ、比較的簡単に正しいプログラムを作り上げることができる。ここでは、
まず、解析的に解を得ることのできる調和振動子について考え、計算が正しく行われていることをどのように確認するかという手続きを述べる。
一次元調和振動子の運動方程式は座標（バネの変位）を q 、運動量を p とすると

dq
p
=
,
dt
m

dp
= −mω 2 q ,
dt

または

m

d2 q
= −mω 2 q
dt2

(3.1)

1
p2
+ mω 2 q 2 は一定値2m 2に保たれる。全エネルギーのように、運動の間一定に保たれる量を保存量という。長さ、である。この運動方程式に従う運動では、常に全エネルギー E =時間、エネルギーについて適当な規格化を行うと運動方程式は

dQ
=P,
dt

dP
= −Q ,
dt

または

d2 Q
= −Q
dt2

(3.2)

1 2 1 2
P + Q と表される。厳密に
2
2は、P は運動量だが、ここでは、簡単のため P を速度と呼ぶことにする。t = 0 での、の簡単な形に変形できる。この時、全エネルギーは E =P 、Q の値を P0 、Q0 とした時の微分方程式の解（解析解と呼ばれる）は
P = P0 cos t − Q0 sin t ,Q = P0 sin t + Q0 cos t

(3.3)

である。
数値計算では、時間に対し、連続的に解を求めていくことはできない。刻み幅を ∆t として、この整数倍の時刻 tn = n∆t での座標 Q、速度 P のみを考える。これを Qn 、Pnと表す。そして、微分方程式を差分方程式に置き換え、すでに得られている値（例えばQn 、Pn ）より、もう一つ先のステップでの値（Qn+1 、Pn+1 ）を順次求めていく。このよ第3章

Verlet 法を用いて調和振動子の運動方程式を積分する。

10

うにして得られた数値解を解析解と比較することにより数値積分の精度（計算の誤差）を知ることができる。計算の誤差には、大きく分けて、微分を差分に置き換えて数値積分を行うという計算方法に由来するもの（これを打切り誤差という）と、計算機の内部での演算がある一定の有効桁数の範囲で行われていることに起因するもの（これを丸め誤差というがある。一般に、時間刻みの大きさ ∆t を小さくすると、打切り誤差は、小さくなる。
丸め誤差は、逆に大きくなっていく。これは、一定時間 t = n∆t まで計算する場合、∆tを小さくすると、計算回数が増えるためである。計算によって生ずる誤差に注目することにより、期待通りに計算が行われているかどうかを判断することができる。
調和振動子の場合には、式 (3.3) のように、解析的に解が求まり、これと数値解を比較することにより誤差を知ることができるが、このように解析的に解が得られるのは非常に限られた場合だけである。我々が通常興味を持ちシミュレーションを行うのは解析解が求められない場合である（だからこそシミュレーションを行う意義があるのである）。この場合にも適用できる誤差の評価法を、調和振動子という簡単な例を通して考えていく。

MD 方で用いられている代表的な数値積分法である Verlet 法（特に leap frog の形）では、運動方程式の積分を次のような手順で行う。計算を１ステップ進めるために必要な既知のデータは tn = n∆t での座標の値 Qn とそれより半刻み前の時刻 tn−1/2 での速度

Pn−1/2 である。
(1) n ステップでの座標の値 Qn より、このステップで粒子に働く力の値 Fn を求める。
調和振動子では Fn = −Qn である。

(2) 微分方程式 dP/dt = F を差分式
Pn+1/2 = Pn−1/2 + Fn ∆t(3.4)に置き換え、tn+1/2 での速度 Pn+1/2 を求める。調和振動子では

Pn+1/2 = Pn−1/2 − Qn ∆t(3.5)となる。ここで、現れる速度の時間が座標に比べて半刻みずれているのが leap frog
（日本では馬とびと呼んでいる屈んだ人の背中に手をついて飛び越す遊び 飛ぶ人と馬になる人が交代しながら進んでいく）という名前の由来である。

(3) n + 1 ステップでの座標 Qn+1 は微分方程式

dQ= P を差分化した
dt

Qn+1 = Qn + Pn+1/2 ∆t

(3.6)

第3章

Verlet 法を用いて調和振動子の運動方程式を積分する。

11で求められる。

(4) Qn と同時刻の速度 Pn が必要な場合（全エネルギーの計算など）には Pn−1/2 と Pn+1/2の平均値が Pn であるとする。すなわち

1
Pn = (Pn+1/2 + Pn−1/2 )
2

(3.7)

を用いる。

計算を開始するためには、初期条件として Q0 と P−1/2 を与えるのだが、この２つの量はおなじ時刻の量ではないので、解析解との比較が難しくなる。式 (3.4) および (3.7) より得られる Pn−1/2 と Pn の間に成り立つ関係式

Pn = Pn−1/2 + Fn

∆t
2

(3.8)

を用い、n = 0 と置くと、同時刻の初期値 Q0 、P0 より

P−1/2 = P0 − F0

∆t
2

(3.9)

と定まる。
課題 3.1

Verlet 法により調和振動子の運動方程式を積分せよ。
∆t、Q0 、P0 を与えて、Qn 、Pn を順次求めるプログラムをつくる。この時、保存量と1なるべきエネルギー En = (Pn2 + Q2n ) と Verlet 法を用いて調和振動子の問題を解いたと
2
きの真の保存量（これについてはシンプレクティックな数値積分法の項を参照のこと）



1
1
2
2 2
Sn =
(3.10)
P + 1 − (∆t) Qn2 n4の値も計算しておく。
初期値を Q0 = 1.0, P0 = 1.0 、時間刻みを ∆t = 0.01 とし、20 ステップまで１ステッ
プごとに Qn , Pn , En , Sn を出力して、解析解の値 P̂n , Q̂n , Ên と比較してみよ。En と Snの値を表?? に示す。保存量 En の値がわずかずつ変化するのに対し、Sn の値は厳密に保存されるはずである。

第3章

Verlet 法を用いて調和振動子の運動方程式を積分する。

12

注意点： 計算は倍精度演算を用いて行うこと。計算が正しく行われているかどうかは打切り誤差の振舞いにより判定する。丸め誤差の影響を小さくするために倍精度演算が必要である。一旦プログラムの正しいことが確認できたなら、計算精度をあまり要求しない場
合には単精度演算を用いてもよい。
はじめてシミュレーションのプログラムを書き上げるときの関門はなかなか正しいプロ
グラムにならないことである。複雑な計算を一気にまとめて仕上げようとするとどこが誤っているのか判断が難しい。この本ではできるだけ要素に分け、部分ごとに確認しながらプログラムを作り上げるように構成している。
上の簡単な積分で保存すべき量 Sn （少なくとも 10 桁程度は保存すべき）が変動した
り、En の値が表に示した値と異なっている場合には、作成したプログラムにエラーがあることを示している。プログラムのエラーを見つける一般的な注意は？？を参照のこと。
もしこの時点で、正しい結果が得られない場合はプログラムチェックの貴重な機会と考えて、よく調べることを勧める。

課題 3.2
課題 3.1 と初期条件を同じにし、時間の刻み幅 ∆t を 0.1, 0.05, 0.02, 0.01, 0.005, 0.002,

0.001 と変え、t = 0.2 の時の Q, P, E の誤差を求めてみよ。計算に必要なステップ数は刻み幅 ∆t により異なり、t/∆t である。この時、途中のステップの値を、出力したり、デー
タとして取っておく必要はない。
誤差 δx は数値計算で求めた値 x と解析解の値 x̂ の差 δx = x − x̂ と定義される。ここで x は座標 Q、速度 P 、またはエネルギー E である。それぞれの誤差を刻み幅 ∆t に対して両対数プロットをしてみよ。数値計算法に由来する打ち切り誤差は時間の刻み幅 ∆tに対し、δx = A(∆t)k のように冪的に依存する場合が多い。この場合、両対数プロットを行うと ln δx = ln A + k ln ∆t と直線になる。この直線の傾きより、誤差の次数 k を決めることができる。この誤差の次数が、用いた数値計算法で予想される値に一致した時、数値積分の（または MD シミュレーションの）プログラムの主要部分が正しく作成されていると確認できる。もちろん、エネルギー保存に影響しないエラーはあるのでこれですべてのエラーに対してチェックできているわけではないが。単に、エネルギーの変動が小さいというのは正しい判定法ではない。プログラムにエラーが含まれている場合でも、時間の刻み幅 ∆t を小さく選ぶとエネルギーの変化も小さくなるが、この場合には、１ステッ
プ当たりの誤差は ∆t に比例し、一定時刻で誤差の比較をした場合には ∆t にはほとんど

第3章

Verlet 法を用いて調和振動子の運動方程式を積分する。

13

依存しなくなる。通常我々は、１ステップ当たりの誤差が ∆t よりも高冪になる数値積分法を用いるので、計算の良否の判定にはこの誤差の次数が再現されることを利用するとよい。

Verlet 法で一定の時間間隔について刻み幅 ∆t を変えて比較したときの誤差の次数は座標、速度、エネルギーともに 2 である。

課題 3.2 では、一定の時間間隔をとり本質的に同じ運動を考えているので、非常に厳密な誤差の比較ができるが、∆t を小さな値にすると多数の計算ステップが必要となり、計算に要する時間が膨大になる。誤差評価には、次のような簡便法も利用できる。

課題 3.3
課題 3.2 と同じ初期条件に対し、∆t = 0.001, 0.0005, 0.0002, 0.0001, 0.00005, 0.00002,

0.00001, 0.000005, 0.000002, 0.000001 のそれぞれの場合について、20 ステップ目でのQ, P, E の誤差を求め、両対数プロットを行い、誤差の次数を求めよ。極端に ∆t を小さくすると直線から外れるようになるが何が原因か考えてみよ。

一定の時間で比較する場合、t/∆t の計算回数が必要となる。したがって、一定回数で誤差評価を行った場合、一定時間で行ったものより、1 だけ大きい誤差の次数が得られる
ことが多い。（この場合、Verlet 法では、誤差の次数は 3 となる。後で述べる予測子・修
正子法ではどちらの評価でもおなじ次数が得られる）。同一ステップ数で比較する場合、比べている運動の範囲が異なっているため、あまり大きな ∆t をとると結果が大きく食い違ってしまう。
刻み幅を小さく取ったときには、数値積分法に起因する打切り誤差ではなく、計算機の中で限られた有効桁数を用いて計算を行っていることに由来する丸め誤差の方が大きくなる。浮動小数点を 8 バイトで表現する倍精度計算の場合、有効桁数は約 1014 桁である
（値の内部表現の違いによりこの桁数は計算機により異なることがある）から、10−14 より
小さな誤差については評価することはできない。このため、∆t が非常に小さいときの誤差は 10−14 程度の値になってしまう。一定の時間間隔で誤差の比較をする場合、∆t をあまり小さくすると計算ステップ数の増加のため丸め誤差が増加し、返って大きな誤差が現れるようになる。

第3章

Verlet 法を用いて調和振動子の運動方程式を積分する。

14∆t に対する誤差の変化の概念図を図？？に示す。シミュレーションでは、許される限り大きな ∆t を取って、計算を効率的に行うことが必要になる。したがって、どの程度の誤差まで許すかによって、∆t の大きさが決まってくる。
解析解が知られていない場合の誤差評価
調和振動子のように、解析的に解が求まる場合、数値的な結果を解析解と比較することにより、計算の精度を直接確認することができる。しかし、シミュレーションは通常、解析解が知られていない場合に行う。このようなときにも計算がどの程度正しく行われているか判定できないと困る。保存系においては、全エネルギーの保存の様子を調べるとよい。長時間での振舞いも重要だが、１ステップで現れる違いに着目して誤差を評価する方が計算時間が短くてすみ効率的である。
誤差の評価には、１ステップ進めたときのエネルギーの変化の大きさの平均値を用いる。
1 ∆=|Ei+1 − Ei |(3.11)N i
ここで、エネルギーの変化量の絶対値を加えて行くのは、計算法によっては（Verlet 法はその代表例）長時間にわたるエネルギーの変化が正負両方に交互に起こり、打ち消し合いがあるためである。時間間隔がたまたま打ち消し合いにより誤差が 0 に近くなる値であった場合、誤差の ∆t 依存性をうまく見積もれなくなる。絶対値を用いた評価では、打ち消し合いがないのでこの困難は生じない。
この場合でも、おなじ時間間隔まで積分した誤差を比較した方がよいが、プログラムのチェックのためには、おなじステップ数だけ行うという簡便法でも十分である。この場
合、誤差の次数はそれほど正確には求まらないが、プログラムにエラーがあるかないか手軽に判定できるためデバッギングの段階では非常に便利である。
計算法によっては（後で述べる予測子・修正子法はその代表例）、初期条件が正しく設定されていないと、初めのいくつかのステップで非常に大きな誤差が現れるものがある。
この場合、初めの 10 ∼ 20 ステップの結果を捨て、その後の 10 ∼ 20 ステップについてのエネルギーの変化量から誤差を見積もればよい。

課題 3.4

1 ステップ当たりのエネルギーの変化量 ∆ を誤差の目安として、課題 3.3 と同じ条件で計算を行いエネルギーの誤差の次数を求めよ。計算は、各 ∆t について 20 ステップづつ
行う。

第3章

Verlet 法を用いて調和振動子の運動方程式を積分する。

15

課題 3.5
初期条件を変え計算の誤差評価を前と同様に行なえ。特に、Q0 = 1.0, P0 = 0.0 とした時、誤差の次数が違って見えるのはなぜか考えて見なさい。

速度が 0 の静止状態や、変位が 0 などの特殊な初期条件の場合、誤差の次数が変ってくることがある。プログラムのチェックのためにはできるだけこのような状態を避けた方がよい。

課題 3.6調和振動子の運動を Verlet 法を用いて解く問題は、代数的に解を求めることができる。
次の手順に従い、解析解を求めてみよ。
座標 Qn 、運動量（速度）Pn を用いると、Verlet 法は２つの漸化式

Pn+1/2 = Pn−1/2 − Qn ∆t(3.12)

Qn+1 = Qn + Pn+1/2 ∆t(3.13)

で表される。式（3.13）で n を n − 1 に置き換えた式

Qn = Qn−1 + Pn−1/2 ∆t(3.14)

を用い、この３つの式より、Pn+1/2 , Pn−1/2 を消去すると、Qn+1 , Qn , Qn−1 ３項の間の漸化式


Qn+1 − 2 − (∆t)2 Qn + Qn−1 = 0(3.15)

が得られる。この漸化式の解は、Qn = Aλn の形に解を仮定して式に代入することにより得られる λ についての方程式
λ2 − 2 − (∆t)2 λ + 1 = 0(3.16)

Qn = Aλn1 + Bλn2

(3.17)の解 λ1 , λ2 を用いてと表される。

第3章

Verlet 法を用いて調和振動子の運動方程式を積分する。

161
cos α = 1 − (∆t)2 と置くと、λ1,2 = e±iα と表すことができることを示せ。な2
お、∆t が小さいとき、

問 (1)α = ∆t +1
3(∆t)3 +(∆t)5 + · · ·
24
640

である。
時刻 tn での運動量 Pn はPn =

 Qn+1 − Qn−11
Pn+1/2 + Pn−1/2 =2
2∆t(3.18)

で求められる。
問 (2)

初期値 Q0 , P0 を用いて Qn , Pn を表すと
∆tQn = Q0 cos nα +P0 sin nα
sin α


sin α
Pn = P0 cos nα −Q0 sin nα∆t(3.19)(3.20)となることを示せ。
問 (3)

1
P2 +
Sn =
2 n



sin α
∆t

2

Q2n

(3.21)

が保存量となっていることを示せ。この式に現れる

2
1
sin α
= 1 − (∆t)2∆t4であるから、Sn と調和振動子の全エネルギー En との差は − 18 (∆t)2 Q2n と (∆t)2 の程度である。

課題 3.7
シミュレーション解法になっている Verlet 法の特徴をよく理解するには、長時間にわた2π2πる性質を調べてみるとよい。時間刻みを ∆t =、および ∆t =に選び、振動の100
1000
一周期にわたる Q, P, E の誤差 δQ, δP, δE, δX = (δQ)2 +)δP )2 を求め、時間とともにどのように変化するか図に描いてみよ。この結果からどのようなことが読み取れるか。

第3章

Verlet 法を用いて調和振動子の運動方程式を積分する。

17

課題 3.82πに選び１周期ごと（ここでは 100 ステップ）に誤差を求め、長100時間での δQ, δP, δE, δX の変化を調べてみよ。（20 周期程度でよい）。

時間刻みを ∆t =課題 3.9
課題 3.8 をいくつかの異なる ∆t で行い、δE の ∆t 依存性を調べてみよ。シンプレクティッ
ク解法の特徴として、計算を長時間行なってもエネルギーの誤差は (∆t)k のオーダーにとどまる。ただし、δQ, δP はどんどん大きくなっていく。エネルギーの誤差は１周期ではほとんど打ち消しあっており (∆t)4 程度でしか現れなくなる。

課題 3.10
時間刻みを ∆t =2πと大きく選び、積分してみよ。
4

これは、調和振動子の１周期の振動をたった４点で解いていることになるが、厳密解法と2πしての Verlet 法の安定限界は調和振動子の場合 ∆t = 2 であり、これより小さな ∆t =4では解は安定であり、単振動的な振舞いが得られる。ただし、この場合、点の間隔があまりにも開いているので一目では単振動であることを確認することは難しい。(Qn , Pn ) の点を次々に平面上に書いていくと、それが１つの楕円の上にのっていることが確かめられる。この楕円を Q 軸に投影した解 Qn は曲線

Q(t) = Q0 cos(αt) +P02

1 − π16sin(αt)の上にのる。ただし、α は



2
π2
−1
α = cos1−= 1.150167 · · ·π8である。正しい解析解は、Q̂(t) = Q0 cos t + P0 sin tである。図？？では、この 2 つの曲線と、Verlet 法で解いた点が示されている。単振動的な振舞いをしている点では正しい解に似ているが、振動の振幅および周期にかなりの違いが見られる。この時、エネルギーは 0.69157 ≤ E ≤ 1.80497 の範囲で振動する。時間刻
み ∆t を大きく取った場合には、長期的にエネルギーの値が安定していても解析解との誤差はどんどん大きくなっていることに注意しておく必要がある。

18

第 4 章 予測子・修正子法を用いて調和振動子の運動方程式を積分する。
Verlet 法の次に、MD 法でよく用いられる数値積分法は予測子・修正子法である。方法の詳細については数値積分の章を参考にしていただきたい。
ここでは、２階微分方程式に対し、５階微分の項まで考慮する方法（５次の方法と呼ぶ）により説明する。現実の応用では、この次数での計算が最もよく用いられるが、これは、予測子・修正子法を用いて行った最初の重要な MD シミュレーション (Rahman とStillinger による水の研究）でこの次数が用いられたためである。
この方法で、運動方程式md2 r= F (r)dt2(4.1)

を積分するためには、時刻 t での座標 r(t) およびその導関数 r(1) , r(2) , r(3) , r(4) , r(5) を必要とする。通常は計算機での処理に便利なように

Xk =1(∆t)k r(k) ,k!

(k = 0, 1, 2, 3, 4, 5)(4.2)の形でデータを保存しておく。
予測子・修正子法では予測および修正の２つの過程を繰り返して実行していく。予測の段階では時刻 t の値 Xk (t) より、次の時刻 t + ∆t での値 Xk (t + ∆t) をテイラー展開の形で近似する。

X0 (t + ∆t) = X0 (t) + X1 (t) + X2 (t) + X3 (t) + X4 (t) + X5 (t)(4.3)

X1 (t + ∆t) = X1 (t) + 2X2 (t) + 3X3 (t) + 4X4 (t) + 5X5 (t)(4.4)

X2 (t + ∆t) = X2 (t) + 3X3 (t) + 6X4 (t) + 10X5 (t)(4.5)

X3 (t + ∆t) = X3 (t) + 4X4 (t) + 10X5 (t)(4.6)

X4 (t + ∆t) = X4 (t) + 5X5 (t)(4.7)

X5 (t + ∆t) = X5 (t)

(4.8)

第4章

予測子・修正子法を用いて調和振動子の運動方程式を積分する。

19

展開に用いた項が有限なので、ここで求めた新しい座標 X0 (t + ∆) を用いて計算した力

F (X0 (t + ∆)) と加速度 X2 (t + ∆) より、
δ=

1 F (X0 (t + ∆))(∆t)2 − X2 (t + ∆)
2
m

(4.9)

を計算しても 0 にならない (運動方程式を満たさない）。
そこで、修正の段階では、予測で求めた値 Xk (t + ∆) を運動方程式に現れた差 δ に比例する値だけ変更する。

Xk (t + ∆) = Xk (t + ∆) + Ck δ

(4.10)

この修正の式に現れる係数 Ck は、運動方程式を満たし、かつ、計算誤差ができるだけ小さくなるように定められる。Ck を定める手順は一階微分の場合に Nordsieck が与えている。２階以上の微分方程式については、Gear が拡張した。係数 Ck の値は、Gear の本に与えてある。２階微分方程式で５次の場合には

C0 =

3
251
11
1
1
, C1 =, C2 = 1, C3 = , C4 = , C5 =
20
360
18
6
60

(4.11)となる。

課題 4.1５次の予測子・修正子法を用いて、調和振動子の運動方程式

d2 Q
= −Q
dt2

(4.12)

を積分せよ。

∆t, Q0 , P0 = Q̇0 を与えて、Xk (n∆t) を順次求めるプログラムをつくる。まず、t = 0
で、高次の微係数の値を正しく与えて計算する場合を考える。調和振動子の場合はその値は簡単に求まり、X0 (0) = Q0 ,X4 (0) =1
1
X2 (0) = − Q0 (∆t)2 , X3 (0) = − P0 (∆t)3
2
6
1
X5 (0) =(4.13)
P0 (∆t)5120

X1 (0) = P0 ∆t,

1
Q0 (∆t)4 ,24である。
初期値を Q0 = 1.0, P0 = 1.0, ∆t = 0.05 とし、20 ステップまで１ステップごとに座標、速度、エネルギーの予測値、X  (n∆t) , En
X0 (n∆t), 1∆tX1 (n∆t)X0 (n∆t),∆t

第4章

予測子・修正子法を用いて調和振動子の運動方程式を積分する。

20

および、修正値
X1 (n∆t)
X1 (n∆t), En X0 (n∆t),
X0 (n∆t),∆t∆tを出力して、解析解の値 Q̂n , P̂n , Ên と比較してみよ。

注意点： ここでは、予測の段階で求めた座標、速度を用いて計算したエネルギー、力などの値と、修正を終えた段階で求めた値を比べている。厳密に言えば、いろいろの量の計
算には修正後の値を用いるべきだが、予測の段階で求めた結果も修正の段階を経たより正
確な値のまわりを揺らぐだけであり、統計的な誤差は現れないことを認識しておくことは重要である。多粒子系の場合、力、エネルギーの計算の部分がシミュレーションの大部分の時間を占める。この時間を要する部分の計算を予測、修正の２つの段階で行うのではな
く、予測の１段階だけにすることにより、計算速度を約２倍にすることができる。

課題 4.2
課題 4.1 と初期条件を同じにし、時間刻みを、0.2,

0.1, 0.05, 0.04, 0.0625, 0.02,

0.01, 0.00625, 0.005, 0.004 と変え、t=1.0 での Q, P, E の誤差を求めてみよ。誤差の次数は、それぞれいくらになるか、予測値で評価した場合と、修正値で評価した場合で違いが見られるか。

課題 4.3
予測子・修正子法では、座標、速度などが連続的に変化すると仮定している。もし、この仮定が破れた場合、非現実的な振動が見られる。しかし、この振動は数ステップのうちにどんどん小さくなり、運動方程式を満たす連続的な解へ近付いて行く。ただし、それは、初期条件で定めた運動とは、異なるものとなっている。
初期条件として、座標の高階の微係数に正しい値を与えるのは複雑な相互作用では困難となる。初期値として、２階微分以上には 0 を与え、課題 4.1 と同様の条件で計算を行い、１ステップごとにどのようにエネルギーの値が変化するか調べよ（20 ステップ程度）。
計算ステップが進むとともにどのようにえねるぎーの変動が小さくなるか評価してみよ。
エネルギーの値を書き出して、数値が変わらない桁数がどのように増していくかを調べる
か、エネルギーの変動 ln |∆E| = ln |En+1 − En | をステップ数に対してプロットしてみるとよい。また、初期条件で設定したエネルギーの値と、最終的に落ちついたエネルギーの第4章

予測子・修正子法を用いて調和振動子の運動方程式を積分する。

21

値の差はどの程度か、∆t を変えて調べてみよ。

課題 4.4
初期条件として、高次の微係数に正しい値を設定していない場合について、計算の誤差を評価せよ。この場合、初期条件で定めた解からはずれてしまうため、解析解との比較は、あまり意味がない。各ステップ毎に、エネルギーがどの程度変動するかを目安として誤差の評価を行え。予測子・修正子法では、正しく初期値の設定がなされていないと初め大きな変動が現れる。初めの 20 ステップの結果を捨て、次の 20 ステップの結果より、誤差を見積もれ。

課題 4.5
予測子・修正子法では、力の部分に速度が含まれる場合でも、本質的な変更なしに数値積分を実行することができる（Verlet 法では、差分方程式を変えなければならない）。速度依存の力が働く場合、非保存力であることが多く、誤差の評価が難しい。ここでは、解析的に解が求まる減衰振動子

d2 Q
dQ= −Q − 2a2dtdt(4.14)

を用いて、力の項に速度が含まれた場合、誤差の次数に差が出るかどうかを調べてみよ。
この方程式を、予測子・修正子法を用いて計算し、座標 Q、速度 Q̇ = P 、エネルギー
1
E = (P 2 + Q2 ) に現れる誤差を評価せよ。厳密な比較を行うためには、高階の微係数に2ついても初期値を与えておく必要がある。t = 0 で Q(0) = Q0 , Q(1) (0) = P0 の時

Q(2) (0) = −Q0 − 2aP0
Q(3) (0) = 2aQ0 + (4a2 − 1)P0
Q(4) (0) = −(4a2 − 1)Q0 + 4a(1 − 2a2 )P0
Q(5) (0) = 4a(2a2 − 1)Q0 + (1 − 12a2 + 16a4 )P0(4.15)である。
また、解析解は



√
√
+
aQ
P
0
0
Q(t) = esin( 1 − a2 t)
Q0 cos( 1 − a2 t) + √1 − a2


√
√
+
aP
Q
0
0
−at
P (t) = Q̇(t) = esin( 1 − a2 t)
P0 cos( 1 − a2 t) − √1 − a2−at(4.16)(4.17)

第4章

予測子・修正子法を用いて調和振動子の運動方程式を積分する。

22となる。計算では a = 0.1 を用いてみよ。計算の始めのステップでは通常の解析より誤差が大きくでるようである。計算を少し進めた時点（１０ステップ程度）の座標・速度の値
より解析解の形を決め直すとより正確に誤差評価を行うことができる。

課題 4.6
エネルギー保存則が成り立たない場合でも、エネルギー変化の式より保存される量を求1
め、計算の精度チェックに用いることができる。減衰振動子の場合、E = (P 2 + Q2 ) の
2
時間変化は

dP
dQ
dE
=P
+Q
= P (−Q − 2aP ) + QP = −2aP 2
dt
dt
dt

(4.18)である。これを積分して得られた式

E(t) − E(0) = −2a t

P (t )2 dt

(4.19)P (t )2 dt

(4.20)

0より、初期条件でのエネルギーの値 E(0) を t
E(0) = E(t) + 2a0と書くと、計算精度を評価するために利用できる。ただし、この式では、数値積分の部分をどのように求めるかにより、全体の精度が決まってしまうため、微分方程式を数値的に解いた時の誤差の次数より小さな次数しか求まらないことが多い。しかし、プログラムエ
ラーがあるかどうかの判定には、この程度のもので十分有用である。
 t
数値積分A(t )dt の計算に最も簡単な台形則01
1
A(t )dt = ∆t A(0) + A(∆t) + A(2∆t) + · · · + A((n − 1)∆t) + A(n∆t)
2
2
0

 t







を用いると、誤差は (∆t)2 程度である。

(4.21)

23

第 5 章 レナード・ジョーンズ相互作用で相互作用する２つの粒子の系の運動
粒子間の相互作用がレナード・ジョーンズ型

φLJ (rij ) = 4

σ
rij



12
−

σ
rij

6
(5.1)

の場合の粒子運動を求める。ただし、rij = |r i − rj | は、粒子 i と j の距離である。２粒子の場合、全エネルギーはH=m 2 m 2ṙ + ṙ 2 + φ(r12 )
2 1
2

(5.2)

であり、運動方程式は∂φ(r12 )∂φ(r12 ) ∂r12
∂φ(r12 ) r 1 − r2
=−
=−
∂r1
∂r12 ∂r1∂r12
r12∂φ(r12 )∂φ(r12 ) ∂r12
∂φ(r12 ) r 2 − r1
mr̈2 = −=−=−∂r2∂r12 ∂r2∂r12
r12
mr̈1 = −

(5.3)
(5.4)

となる。もし必要ならば、アルゴンの値

m = 39.9mu = 6.63 × 10−26 kg= 120kK = 1.65 × 10−21 J
σ = 3.4A = 3.4 × 10−10 mを用いるとよい（mu は原子量単位）。
座標、速度などの値を実単位（例えば上記の m, , σ についての MKSA 単位系での値）で表すことは勧められない。計算機の中での数値の表現として非常に大きい、または非常に小さな数が現れ困難を生ずることがある（コラムを参照のこと）。原子レベルの大きさを持つ量で規格化し無次元の量として表すことが望ましい。レナード・ジョーンズポテン
シャルの場合には、長さを σ 、質量を m、エネルギーをを単位として表すと、すべて

第5章

レナード・ジョーンズ相互作用で相互作用する２つの粒子の系の運動24の量が物質に依存するパラメターを一切含まない形で調べることができる。これは、多くの物質に共通する普遍的な性質を調べる場合にふさわしいやり方である。実験との比較を重視して、ある特定の物質の示す性質を調べる場合には、 , σ などで規格化すると、数値的な比較が難しくなる。この場合、実単位を基に単位の大きさを適切に選び数値が 1 程度の値となるようにするとよい。例えば、長さについては A を、質量については原子量単位または g/mol を、圧力については MPa, GPa または気圧を単位として用いると実験値との比較が容易である。

課題 5.1
レナード・ジョーンズ相互作用を持つ２つの粒子について運動方程式を積分するプログラ
ムをつくれ。Verlet 法および予測子・修正子法の両方について行え。
まず、プログラムに誤りがないことを確認するために誤差の評価を行う。初期配置として、２つの粒子が 2σ 程度離れて静止している状態を用いるとよい。解析解は得られないのでエネルギーの変化に着目して誤差を調べよ。特に、予測子・修正子法については、高次の導関数について、正しい初期値を設定せずに誤差の評価を行ってみよ。

課題 5.2
時間刻みを、系の代表的な時間 τ =mσ 2（アルゴンの場合 2.1 × 10−12 s）の 100 分の 1
程度（規格化した座標系では 0.01τ 、実時間を用いる場合は 0.02 ps）に選び 1000 ステッ
プ程度計算し、エネルギーの値がどのように変化するか調べてみよ。
エネルギーの変化が大きいのはどのようなときか。初期条件を変えて調べてみよ。全エ
ネルギーが負となるような状態を選ぶ。また、全体としての粒子の動きがないように、全運動量 m(v 1 + v 2 ) が 0 となるように選んでおく。

注意点： 粒子の速度により、計算の精度（全エネルギーの変化量）が大きく変わることに注意しておくこと。特に、粒子が接近し、反発してまた離れていく過程で誤差が大きく現れる。誤差の制御は、この誤差が一番大きくなる部分に注目して ∆t の大きさを選ぶことによって実行できる。場合によっては、速度が大きくなると、∆t の大きさを小さくすることも必要となる（万有引力で相互作用する星の運動の場合には必須である）。

課題 5.3
課題 5.2 を時間刻みをもっと小さな値に取って繰り返してみよ。全エネルギーの変化量は

第5章

レナード・ジョーンズ相互作用で相互作用する２つの粒子の系の運動

25

どの程度になるか。

課題 5.4計算の最後のデータをファイルに保存し、さらにこのファイルを読み込むことにより、計算を継続できるようなプログラムをつくれ。計算を再開する際、必要であれば、時間刻みの大きさ、粒子の速度の変更ができるようにする。プログラムの先頭で、パラメターを読み込むことにより、初期配置から始める場合と、計算を継続する場合が選択できるように
する。

26

第 6 章 多粒子系のシミュレーションを行
なう。
周期境界条件を考慮し、レナード・ジョーンズ相互作用を持つ多粒子系のシミュレーショ
ンを行なうプログラムを作る。ステップ５で行なった２粒子の場合と比較すると、多粒子系では、周期境界条件と相互作用のカットオフを付け加える必要がある。このうち、相互作用のカットオフは誤差の評価を難しくする。なぜなら、粒子間距離がカットオフの距離

rc となったとき、位置エネルギー、力に不連続な変化が生じ、エネルギー保存が満たされなくなる。カットオフによる誤差は時間刻みの大きさ ∆t にあまり依存せず、またその大きさも通常の誤差評価で調べる値よりも大きいため、カットオフがあると今まで述べた
誤差評価法はそのままでは用いることができない。プログラムにエラーがあるかどうかの判定は、カットオフの影響が現れないようにして行なう必要がある。このためには、次のような手法が有効である。

1. ステップ５で調べたような数個の粒子からなる孤立系を考える。この時、粒子対の相互作用をすべて求めるためカットオフは必要ない。周期境界条件と両立させて、この状態をつくるためには、一辺の長さを粒子間距離の 10 倍程度と大きく取ったシミュ
レーションセルの中心付近に数個の粒子を置くとよい。カットオフの距離は粒子間距
離より十分大きく取る。この配置により、数値積分法が正しくプログラムされている
か、ポテンシャルと力の計算が整合しているか等を確かめることができる。

2. 周期境界条件をあらわに考慮した場合のチェック法。fcc 等の完全結晶の粒子配置をつくる。相互作用のカットオフの距離を rc を g(r) が 0 となる領域に選ぶ（例えば、最近接粒子間距離 r1 と第２近接粒子間距離 r2 の平均値（ 12 (r1 + r2 ) にとる）。小さな初期速度を与え、時間刻み ∆t を小さく選んでおくと、計算を開始してしばらくの間は粒子間距離がカットオフ距離 rc に達することがないようにすることができる。
この短い時間の間に誤差の次数を計算する。この配置により、周期境界条件に関連した部分のプログラムのチェックができる。クーロン力をエバルト法を用いて計算する
場合にも、このチェック法は利用できる。

第6章

27

多粒子系のシミュレーションを行なう。

本格的にシミュレーションを行なう場合にはカットオフに伴うエネルギーの揺らぎをそれほど気にする必要はない。原因のわかっている誤差はあってもよい。これは、計算を実行可能にするために行なったやむを得ない近似である。ただし、その誤差が予想される程度であることを常に確認しておくとよい。また、結果の信頼度を考えるときこの誤差があることを考慮する必要がある。
カットオフに伴うエネルギーの跳びを見かけ上少なくするには、カットオフの距離 rc
でポテンシャルが 0 となるようにエネルギーの原点を変えたシフトポテンシャル

φsh (r) = φ(r) − φ(rc )(6.1)

を用いるとよい。カットオフおよびシフトのない元の体系と比較する必要がある場合には、rc 以内に含まれる平均粒子数 nc （これは、動径分布関数を求めるためのヒストグラ
ムから求めることができる）を用い、

1
E = Ec + φ(rc )N nc + 2πN ρ2 ∞φ(r)r2 dr

(6.2)

rc

とすればよい。右辺第２項は、ポテンシャルのシフトに対する補正、第３項は、カットオ
フに対する補正である。
先に述べた２つのチェック法が利用できるように、粒子の初期配置として、

1. 大きなセルの中心付近に数個の粒子を置くもの、
2. fcc などの結晶配置、3. これからシミュレーションを行なう予定の結晶構造を入力により選択できるプログラムとする。液体、気体など乱れた配置を持つ状態でもまず結晶状態に粒子を並べてから行なう。これは、液体、気体に現れる乱れた粒子配置を直接作り出すのは難しいからである。シミュレーションを行なう状態とおなじ密度を持つ結晶構造をつくり、それを融点の２倍程度の高温にして、融解させ、その後に望みの温度に戻すことにより、液体、気体状態に対してもそれほど困難なしにシミュレーションを行なう初期状態をつくることができる。
温度を変更するためには速度のスケーリングを用いる。温度 Teq は運動エネルギーの平均値と関係している。

1
i

2


mi v 2i

=

3N
kTeq
2

(6.3)

第6章

多粒子系のシミュレーションを行なう。

28

運動のエネルギーの値と

1
i

2

mi v 2i =

3N
kT
2

(6.4)

の関係で定義される量 T を瞬間の温度と考えると、望みの瞬間の温度 T  に変更するには、すべての粒子の速度をv i = sv i

(6.5)

1
3N   1
3N
kT =mi v i 2 = s2
mi v 2i = s2
kT
2
2
2
2
i
i

(6.6)

と s 倍にするとよい。この係数 s は

より


s=

3N
kT 
 21
2
i 2 mi v i(6.7)

である。ただし、このスケーリングで変更できるのは瞬間の温度（運動エネルギー）の値
であり、本当の温度が T  になったわけではない。本当の温度を T  に近くするためには、すくなくとも速度スケーリングを 10 ステップ程度の間隔で 5∼ 10 回程度行なわなければならない。
シミュレーションは次の３つの段階に分けて行なう必要がある。本当のシミュレーショ
ンは第３段階である。
第１段階はシミュレーションを行なう状態を準備することである。特に、液体、気体状
態では、先に述べた手順により、調べる状態にふさわしい配置をつくっておく。結晶状態では調べようとする粒子配置を初期配置としてつくる。調べる温度にふさわしい速度の分布を与え、その後、速度スケーリングを用いながらシミュレーションを行い、指定した温度の近くになるようにする。（残念ながら、通常の MD 法では厳密に指定した温度にすることはできない。温度は、運動エネルギーの平均よりシミュレーションの後で正確な値がわかるものである。）
第２段階では、シミュレーションを続けながら、圧力、位置エネルギーなどいくつかの量の時間変化を調べ、平均の値が時間とともに変化していく傾向がなくなり、系が平衡状態に達するまで待つ。次の第３段階のシミュレーションで求める統計平均の結果の信頼度を高くするためには、この第２段階で系が平衡に達していること確認しておくことが重要

第6章

29

多粒子系のシミュレーションを行なう。

である。
第３段階がシミュレーションの本体部分であり、統計平均はこの部分のデータのみを用いて行なう。各ステップで求められた運動エネルギー、位置エネルギー、圧力などの値を計算ステップ数だけ加えていき、平均量を求める。平均を計算する場合には、結果の信頼度が評価できるように、できるだけその揺らぎまでもとめておく。これらの量の計算では、統計平均を開始する直前（または最初）の値を控えておき（これを A0 とする）、この値との差を求めてから和を求めると計算に伴う丸め誤差を小さくすることができる（仮平均の方法）。
平均では Ai − A0 、揺らぎでは (Ai − A0 )2 の和を求める。真の平均値は1 
Ai  = A0 + Ai − A0  = A0 +(Ai − A0 )
M i=1
M

(6.8)

真の揺らぎの分散は2 (δAi )2  = (Ai − Ai )2  = Ai − A0 − (Ai  − A0 )
2
= (Ai − A0 )2  − Ai − A0 
M

2
1 
=
(Ai − A0 )2 − Ai − A0 
M i=1

(6.9)

で求めることができる。M は計算のステップ数である。たとえ他の部分での計算が単精度でおこなわれていたとしても、統計平均を計算する部分は倍精度演算で行なわれるべきである。
実際のシミュレーションでは、第１、第２、第３のそれぞれの段階でどのように実行して行くのか定め、それをプログラムに入力データとして与えなければならない。これをジョブのスケジューリングという。いろいろの状況すべてに対応できるようにと考えるとプログラムは複雑なものとなっていく。場合によっては、シミュレーションの実行のたびにソースプログラムを変更していくことになる。しかし、ソースプログラムの変更はできるだけ行なわない方がよい。一部を変更する際、他の部分まで変えてしまい新たにエラーが生ずる危険がある。完成したプログラムはできるだけそのままの形で用いるべきである。
筆者の場合は複雑なジョブのスケジューリングはプログラムで行なわず、ジョブを短く切り、それを続けて行なう形で処理している。計算の継続はシミュレーションの最後で出力するデータファイルを介して行なっている。この場合、プログラムに入力すべきスケ第6章

多粒子系のシミュレーションを行なう。

30

ジューリングに関するデータは計算ステップ数、温度、時間刻みの大きさ、速度スケーリ
ングの実行/非実行、速度スケーリングの間隔、データの出力の間隔、程度になる。
例えば、初めて、液体状態でのシミュレーションを行なう場合、次の４つのジョブを続けることになる。

• ジョブ１ 粒子の初期配置として結晶構造をつくる。温度を融点の２倍程度の高温に設定し、速度スケーリングを行ないながら、結晶を融解させる。

• ジョブ２ 温度を最終的なシミュレーションを行なう温度に変更し、速度スケーリ
ングを行なって、この温度になるようにする。ここまでが第１段階に相当する。

• ジョブ３ 速度スケーリングは行なわずにジョブ２の計算を継続する。ある時間間隔でモニターする量を出力し、平衡状態に近づいているかどうか判定する。（第２段階）• ジョブ４ 長いシミュレーションを行い、統計平均を求める。場合によっては、途中のステップでの粒子の座標、速度をファイルに出力し、後で、もっといろいろの解析ができるようにする。（第３段階）

特に、その体系について初めてシミュレーションを行なう場合や、時間的に余裕のある場合には、それぞれの段階の結果を見て、次のステップを進めるほうがよい。計算になれてくると、見込みでスケジューリングを行い、第３段階まで一気に行なってしまうが、第２段階で十分平衡に達していないことがわかり、結局、全体が無駄な計算になってしまう
場合が出てくる。

課題 6.1
レナード・ジョーンズ相互作用を持つ多粒子系について、次に記す条件を満たすプログラ
ムをつくる。数値積分には、Verlet 法か、予測子・修正子法を用いる。誤差の評価を行なって、プログラムのエラーを見つけ、正しいプログラムにする。

1) 課題 5.4 と同様に、初期配置から始める場合と、以前の計算の継続がえらべるようにする。
計算の継続ができるように最終配置のデータをファイルに保存する。データの種類によって別のファイルに保存する必要がある。FORTRAN を用いる場合には入出力番
号をいつでも決めておくと間違いが少なくなる。筆者を用いている例では、計算の継

第6章

31

多粒子系のシミュレーションを行なう。

続可能なデータが１０番、計算の出力する座標・速度が９番、エネルギーなどの変化をグラフィックに出すためのデータが２番となる。
初期配置は、プログラムチェックのために数個の粒子だけを考える場合と N = 256 (4×4 × 4 × 4) までの fcc 結晶が選べるようにする。シミュレーションのセルの形は辺の長さがそれぞれ Lx , Ly , Lz の直方体とする。

2) 速度スケーリングを用いて、指定した温度に変更できるようにする。
3) 計算の再現が可能となるように計算の開始に用いたデータを記憶しておく。計算に問題点があった場合、まったく同じ計算をもう一度行なうことが必要となる。計算機で用いる乱数は通常擬似乱数であるから、乱数の種（初期値）を記録しておくと乱数についても計算の再現が可能となる。

4) 平衡状態量としては少なくとも、運動エネルギー
K(v) =

1

mi v 2i(6.10)φ(rij )

(6.11)

2i

位置エネルギー

U (r) =

i<j

全エネルギーE = K(v) + U (r)(6.12)

圧力（厳密にはこの量の平均が圧力）

P =

1
3V



mi v 2i −

i


i<j

rij

dφ(rij )drij(6.13)の平均および揺らぎの分散を求める。第１、第２段階では、これらの量の変化により平衡に達しているかどうかの判断をする。
シミュレーション終了後、温度

T =

2
K
3N k

(6.14)

第6章

多粒子系のシミュレーションを行なう。

32

定積熱容量

CV =

k
(δK)2 
2
−
3N
K2

(6.15)を求めよ。

5)構造を表す量として、動径分布関数 g(r) を求めよ。ステップ２での結果を利用すること。g(r) を求めるためのヒストグラムは、力の計算を行なう部分で粒子対の距離を求めたときに計算しておく。この力の計算部分は、MD シミュレーションで最も計算労力を必要とする部分であり、この部分の計算の効率化が全体の計算速度を決める重要な部分であるが、まず一番初めは、すべての粒子対についてエネルギー、力を計算する形で行なってみるとよい。例としては次のようになる。x 方向のセルの長さ

Lx を変数 alx、その半分を alxh、粒子の座標は (3, N ) 次元の配列 x、カットオフの
長さ rc の２乗を rcut2 とする。

do 2000 i=1,N-1
do 1000 j=i+1,Nxx=x(1,i)-x(1,j)
if(xx.gt.alxh) xx=xx-alx
if(xx.lt.-alxh) xx=xx+alh
Y、Z 成分も同様
rr=xx*xx+yy*yy+zz*zz
if(rr.gt.rcut2) go to 1000
r=dsqrt(rr)g(r) のためのヒストグラムの計算
エネルギー、力の計算

1000 continue
2000 continue
スーパーコンピューターを用い、高速計算を考える場合 g(r) の計算が上の例のような位置にあると、この最も計算時間を要する部分のベクトル化ができなくなる。この時は、DO ループの中ではヒストグラムの計算はせず、粒子対の距離を配列に保存しておく。そして、DO ループ終了後、この配列のデータを読んでベクトル化不可能なヒストグラムの計算だけを別に行なうようにする。こうすると、力の計算の主要部分

第6章

多粒子系のシミュレーションを行なう。

33

はベクトル化できるようになる。

6) 動的な振舞いを調べるために、計算の直前（または最初）のステップの座標 ri (0) を配列として取っておいて、平均２乗変位

δ(t) =

1 |r i (t) − r i (0)|2N i(6.16)を求める。長時間ではこの値は 6Dt に近づいていく。できれば、x, y, z 方向ごとに計算しておくと３つの値の違いより、信頼性を調べることもできる。
数値積分で粒子を動かしシミュレーションセルから出た時、周期境界条件により、元のシミュレーションセルに戻すが、この時、r i (0) の値も同様の変換をしておくと変位 r i (t) − r i (0) は時間について連続的になる。
プログラム例

if(x(1,i).gt.alx) thenx(1,i)=x(1,i)-alxx0(1,i)=x0(1,i)-alx
endif
if(x(1,i).lt.0.0d0) thenx(1,i)=x(1,i)+alxx0(1,i)=x0(1,i)+alx
endif
課題 6.2
課題 6.1 でつくったプログラムを用いて、N = 108 の系の液体状態をつくる。

1) fcc 結晶から始め、温度を融点 (アルゴンの場合、実単位で 84K、 /k 単位では 0.7）の２倍程度の高温に設定する。速度スケーリングを行って、結晶構造が壊れ、乱れた配置となることを確かめる。時間刻み ∆t は 0.02ps （実単位）または 0.01（規格化したとき) 程度にする。100 ステップ程度でよい。もし、プログラムの暴走（運動エネルギーが急激に増加して数値計算不可能になること）が起これば、もっと小さな∆t にして再度試みてみる。ただし、暴走が起こる一番の原因はプログラムにエラーがあることである。

第6章

多粒子系のシミュレーションを行なう。

34

2) 乱れた配置となったことは、動径分布関数の形、平均２乗変位が時間に比例して増加すること、または、すべての粒子の座標をある平面に投影した時、規則性が消えることよりわかる。

3) 温度を、融点に変え、速度スケーリングを行いながら、100 ステップ程度計算を行う。
運動エネルギー、位置エネルギーがどのように変化するか調べてみよ。特に、速度ス
ケーリングを行った後、どのように変化するか。
もし、速度スケーリングを行う前の状態が平衡状態であれば、ルシャトリエの原理（平衡からはずれるような変化が系に加えられると、それを妨げる方向に系は変化する）により、変化の方向がわかる。温度を下げたときには、温度を上げる方向に変化が起こる。ただし、その変化量は、始めに加えられた変化を打ち消すほど大きくはない。

課題 6.3

N = 108 の系について、融点近傍の温度で平衡になっていることを確かめた上で、2000
ステップ程度のシミュレーションを行ない、平衡状態量、構造、拡散係数を調べてみよ。
平均については、全体での平均だけでなく、500 ステップ程度毎の部分平均および揺らぎも計算できるようにしておくと、平均値の信頼度を見積もることができる。
全体が M ステップとする。部分平均を m ステップ毎に M/m 回行ったとする。M ス
テップの平均量を AM 、m ステップでの平均量を Am とする。これらの平均を多数回行ったとき、異なる部分平均に含まれるデータが互いに独立と考えてよい程各部分平均のス
テップ数が長いと仮定すると、平均値の平均と平均値の揺らぎはAM  = Am m(δAM )2  =(δAm )2 Mとなる。これから、M ステップで平均した値の信頼度はm
(δAm )2 
AM ±Mである。もし、平均に用いるデータがすべて独立であれば、
1AM ±(δA)2 M(6.17)(6.18)(6.19)(6.20)のように揺らぎの分散によって誤差が見積もれる。しかし、実際のデータは、近い時間については独立でないので、上の関係は厳密ではないが、信頼度について目安を与える。計第6章

多粒子系のシミュレーションを行なう。

35

算結果を

AM ±

(δA)2 (6.21)のように、平均と、標準偏差を与える形で示しておくとよい。これにより、結果を見た人は実際の誤差は式 (6.20) 程度だと、理解することができる。
拡散係数 D は、平均２乗変位 δ(t) を時間に対しプロットし、δ(t) が直線的になった部分の傾きが 6D となることから求めるのが最もよい結果を与える。本来の定義である、
δ(t)limの収束は遅い。
t→∞ 6t

課題 6.4
課題 6.1 でつくったプログラムでは、力の計算で、粒子数が N のとき 12 N (N − 1) の粒子対についての和が現れる。これは、粒子数 N の２乗に比例して増加して行くため、粒子数を大きく取ると、計算時間が急激に増大して行く。全体の粒子数を増やしても、カット
オフの距離を一定値 rc にしておくと、１つの粒子から見て rc 以内にある粒子の数はほぼ一定であるため、実際に相互作用している粒子対の数は N に比例する数程度でしか増加しない。このことを、うまく計算に取り入れると、粒子数が多い場合（100 個以上になると有意の差がでる）計算時間を短くすることができる。
各粒子について、近接粒子の表をつくり、その表に記録されている粒子についてのみ、相互作用の計算を行うようにする。こうすると、計算労力が N 2 に比例するのは、表をつくる時のみで、力の計算部分は N に比例するようにできる。表を新しいものに置き換える頻度を 10 ∼ 20 ステップに１回程度にすると、表をつくるための労力は全体に対し小さくなる。
表をつくる際、注意すべき点は、表に含まれていない粒子対の距離が近くなって相互作用の範囲内に入ってくるということが起こらないようにすることである。このために、近接粒子の表をつくるときには、rc より少し大きな距離 rt を用いる。表をつくる間隔を大きく取ると、表をつくるための労力は減るが、rt をより大きく選ばなければならなくん
り、力の計算部分での手数が増えてしまう。
記憶容量の点で不利だが理解がしやすいので、まず２次元の配列 ntable(j,i) を用いる
場合を説明する。この時には、粒子 i 毎に、近接粒子が記録されている（その j 番目がntable(j,i) ）。粒子 i について記録されている近接粒子の数を nnum(i) とする。また、rtの２乗を rtab2 とする。

第6章

多粒子系のシミュレーションを行なう。

表をつくるためのプログラムは次のようになる。

do 2000 i=1,n-1
nn=0
do 1000 j=i+1,nxx=x(1,i)-x(1,j)
if(xx.gt.alxh) xx=xx-alx
if(xx.lt.-alxh) xx=xx+alx
y、z 成分も同様に
rr=xx*xx+yy*yy+zz*zz
if(rr.gt.rtab2) go to 1000
nn=nn+1
ntable(nn,i)=j1000 continuennum(i)=nn
2000 continue
ここでつくった表を用いて力の計算を行う部分のプログラムは次のようになる。

do 4000 i=1,n-1
nmax=nnum(i)if(nmax.eq.0) go to 4000
do 3000 nn=1,nmax
j=ntable(nn,i)xx=x(1,i)-x(1,j)
if(xx.gt.alxh) xx=xx-alx
if(xx.lt.-alxh) xx=xx+alx
y、z 成分も同様に
rr=xx*xx+yy*yy+zz*zz
if(rr.gt.rcut2) go to 3000
力、エネルギーなどの計算

3000 continue
4000 continue

36

第6章

多粒子系のシミュレーションを行なう。

37

配列 ntable(j,i) の j の部分の大きさは、rt 以内の粒子数の平均にある程度の余裕をみて定めておく。 j > i の場合についてのみ表をつくったので、i が大きくなると表に含まれる粒子が少なくなる。２次元配列を用いた場合、その少なくとも半数は空となっており、記憶領域を無駄にしている。
表を１次元配列にすると記憶容量を約半分にすることができる。この時には、粒子の表

ntable(i) と i-1 番目の粒子の表の終りの位置を示す npoint(i) の２つの配列を用いる。
この場合のプログラム例は、

nn=0
do 2000 i=1,n-1npoint(i)=nn
do 1000 j=i+1,n···rr=xx*xx+yy*yy+zz*zz
if(rr.gt.rtab2) go to 1000
nn=nn+1
ntab1(nn)=j1000 continue
2000 continue
npoint(n)=nn
do 4000 i=1,n-1
nstart=npoint(i)+1
nend=npoint(i+1)
if(nened-nstart.lt.0) go to 4000
do 3000 nn=nstart,nend
j=ntab1(nn)
···
3000
4000

continue
continue

となる。

第6章

多粒子系のシミュレーションを行なう。

38

npoint(i) は i-1 番目の粒子の表の最後の位置を示している。従って、i 番目の粒子に接近している粒子の数は ntab1(i) の i=npoint(i)+1 から npoint(i+1) までに収められている。
もし、計算に用いる粒子数が極端に多く（１万個以上）なると、１つの粒子と相互作用する粒子の数がほぼ一定の場合でも、相互作用する粒子対の表の配列が非常に大きくなり計算機の記憶容量の制限に引っかかってしまう場合がある。この場合は、シミュレーショ
ンセルをブロックに分けて表をつくることにより計算を効率的にできる。シミュレーショ
ンセルを x, y, z 方向に、それぞれ m 分割し、全体で m3 のセルに分ける。この時、セルの辺の長さを、カットオフの距離より少し大きくなる程度の値にとると、ある１つのセルに含まれる粒子と相互作用する粒子は、必ずそのセルに隣接する 27 個のセルの中にある。この場合には、各ブロックに含まれる粒子の表を作っておくだけで、力の計算ができるようになる。

課題 6.5
課題 6.3 N = 256 まで計算できるプログラムについて、力の計算を粒子対の表を用いて行う形に変更しなさい。相互作用のカットオフの距離を N = 108 の時のシミュレーショ
ンセルの辺の長さの半分とした時、表を用いない場合に比べどの程度の計算時間の短縮が可能か、N = 108, N = 256 の場合について調べなさい。

CPU 時間を測るためにはシステムに内蔵されているルーチンを用いるとよい。初期設
定の部分を省き、新しい粒子配置を作り出す部分の CPU 時間を比較する。時間の比較が目的だから、多くのステップ数の計算を行う必要はないがステップ数が短すぎると計算が速すぎて、精度よく CPU 時間が計測できないことがある。(特に、秒単位でしか測定できないときには注意が必要である）。