分子シミュレーション

(1)

ー

分子シミュレーション

(25ページ参照）

(2)

2 第1

章分子シミュレーション

1.1

本書の目的

分子シミュレーションの方法は古典系については古典力学，統計力学を基礎としている．量子系にはさらに量子力学が加えられる．

多数の粒子を体精一定，エネルギー一定のもとで運動させるとどんなマクロの性質があらわれるのか，計算機を用いて運動方程式を数値積分し，力学鼠の時間平均をとおしてマクロな性質を調べる分子動力学法

(MD

法）と，統計力学の集団理論に基づき多数の粒子の空間配置を確率法則によって生成させるモンテカルロ法

(MC

法）とによる分子シミュレーションは前世紀の

50

年代に始められた．力学法則による結果と統計法則による結果か一致するこ

と，および剛体球系の固相一流動相転移の存在が明らかにされ，レナード・

ジョンズポテンシャルの系にも応用されてその有効性が明らかになり，計算機の発達に伴って，対象も球対称の単純な粒子系から，イオン系，水，高分子，…と構造をもつより複雑な体系に広がっている．そして物理，化学，

さらには生体高分子など生物の分野でも実験，理論に加えて物質科学の有力な研究手段として，その地位を確立してきた．

それにいたるまでには方法論の開発・発展が大きく貢献してきた．本書の目的の一つは，古典系のシミュレーションについて種々の方法の物理的数学的根拠を明確に示し，そのうえでアルゴリズムを説明し，具体的にプログラムができるようにすることにある．

原子，分子の電子構造を非経験的に求める道は分子軌道法として化学の分野で開拓されてぎた．たとえば，水

2

分子の相互作用は，電荷，水素イオン，酸素原子が剛体的に拘束された構造にモデル化され，パラメータを含むポテンシャルエネルギーの形に表現されている．パラメータは実験との比較によって決められる．この経験的に相互作用をモデル化する方法に加えて，

非経験的分子動力学法（カー—パリネロの方法， CP 法）が 80 年代半ばに開発

され，原子，イオン間の力の媒介をする電子雲の動きを直接調べることが可

(3)

2

分子間力，分子モデルと計算機実験

原子・分子間に働く力をどのようにモデル化するかは分子シミュレーションにおいてもっとも重要な課題である．対象とする物質をどうモデル化するかは，その何を調べようとするかによって異なる．いろいろな見方で，

単純な物質からより複雑な物質へ，いくつかの例をあげ

る．

(4)

8 第2章分子間力，分子モデルと計算機実験

2. 1

はじめに

シミュレーションの対象とする系のマクロな性質（物性，構造）は構成要素の原子，分子の示すミクロの運動を平均したものである．その運動の特徴は原子，分子のもつ属性と分子間力をとおしてあらわれる．分子間力は原理的には電子の多体問題を量子力学的に取り扱うことによって近似的に求められる．これは本書の後半の問題である．本章では分子間力を経験的情報からモデル化してきた経緯をいくつかの例について述べる．

希ガス（不活性ガス）原子間に働く力は

2

体相互作用（対ポテンシャル，

pair potential)

によってよく記述される．対ポテンシャルは希ガスの他に単純な金属やイオン系によく用いられている．しかし構造をもつ分子では

3

体力，

4

体力が無視でぎなくなる．このような分子の分子間力をどう表現するか，

対象をどのようにモデル化するか，系が複雑になってくるはどモデルの特徴のどこがよくて，何がわるいのかその見分けが困難になってくる．ファン・

デル・ワールスの状態方程式から話を始め，単純な原子モデルからより複雑な分子モデルヘ，計算機実験の方法にも触れながら後章への準備とする．

2.2

ファン・デル・ワールスの状態方程式

原子，分子間には数オングストロームの近距離で分子間力が働く．イオンなどではさらに長距離のクーロンカが働く．分子の形を最も単純な球状分子の場合に限ると，分子間力の斥力は分子の大ぎさと対応させることができるであろう．

ニュートンは気体が圧力を示すのはアトム間に斥力が存在するからだと仮

定し，その斥力を

F(r)=

1 /

r, r S:: r1 ; F(r) =

0 ,

r > r1

いは最近接

アトムとの距離）とおいた．マクスウェルは気体運動論において罪離の逆ベ

きの斥力を論じ，とくに

r‑5

の斥力（ポテンシャルではない）の場合，気体の

(5)

2. 2 ファン・デル・ワールスの状態方程式

︐

分布関数を知らなくても粘性係数など輸送係数を計算できることを示した．

この分子モデルは彼の周囲の人々は別として，当時のヨーロッパ大陸の科学者には受け入れられなかった．その理由の

1

つに引力が含まれていないことがあったそうである（文献

[2.l]).

分子間力が斥力と引力からなることは，ファン・デル・ワールス

(vander Waals)

の状態方程式によって示された

(1873

年） .

1

モル当たりの式は

(P+~2)(V-b)=RT

である.

a, b

は定数．この状態方程式を

P = RT V‑b―v

^a

2

(2.1)

( 2 . 2 ) と変形すると，第 1項の (V‑b )は，分子が大きさをもっため分子の運動でできる空間の容積が

b

だけ減少したためとみることができる．つまり分子間に斥力が働くため，ある距離以上は互いに近付けないことを意味する．圧力は第

2

項の分だけ弱められる．これは分子間に引力が働いて引ぎ合うため，

p o

T>T, T=T,

}r<T,

＞

図2.1 ファン・デル・ワールスの状態方程式

(6)

4 データ解析

エネルギー一定の系の長時間乎均，温度一定の系のカ

ノニカル平均について，例として状態方程式，動径分布

関数，静的構造因子，および自己拡散係数，速度相関関

数の求め方を述べる．

(7)

88 第4章データ解析

4. 1

はじめに—長時間平均とカノニカル平均

粒子数

N,

体積

V

で相互作用ポテンシャル

U(ri,r2, …

，

r

りの質点系を考えよう．この系はエネルギー

E

が一定の保存系であるから，前章で述べた

M D

法は近似的にエネルギー

E

を一定に保つとみてよい．そこで粒子数

N,

体積

V

およびエネルギー

E

を一定に保つ

M D

法を

NVE‑MD

法と略記することにする．今後，同様の略記法を圧カ一定，温度一定の場合にも用いる．

初期条件を与え，

NVE‑MD

法によって力学系の時間発展を追跡していくと，各粒子の位置と運動量が時間の関数として与えられる．運動エネルギーやポテンシャルエネルギーなどの力学鼠を時間の関数としてグラフに描く

と，力学量はゆらぎながら変化していき，ある時間（緩和時間という）経過すると一定値のまわりで変動するようになる．いろいろの力学量が一定値のまわりでゆらぎ始めたとき，力学系は熱力学でいう熱平衡状態に入ったとみなすことにする．系の熱力学量，温度，内部エネルギー，圧力などは，それぞれ対応する力学量を，この熱力学状態に含まれる各粒子の位置と運動量で決まるミクロ状態について平均したもので与えられると考える．

この平均は統計力学でいう長時間平均である．粒子の座標と運動量（あるいはその一部分）の関数である力学量を時間

t

の関数として

A(r(t),p(t))

=

A(t)

と記すと，

A(t)

の長時間平均

<A(t)〉t

は

〈A(t)〉1= lim+ _T_{→ ②}

1

^A(t)dt ⁽⁴^.¹⁾

で定義される．この長時間平均はミクロカノニカル集団に関する平均（統計

平均）と等しいと考える．一方，統計力学の集団理論によれば，ミクロカノ

ニカル集団についての平均は，熱力学的極限（条件

N/V=

一定のもとで

粒子数

N

→

00

とする極限）においてカノニカル集団，グランドカノニカル

集団についての統計平均と一致する．

(8)

4. I はじめにー長時間平均とカノニカル平均 89

モンテカルロ

(MC)

法において最も璽要な統計集団はカノニカル集団である

. N

個の同種粒子からなる古典系について，カノニカル平均は次のように定義される．カノニカル集団に属する体系は体積 (V)一定，温度 (T)一定のミクロの状態をとる．この体系がエネルギー

H(p

尺

rN)= K(pN)

+

U(r

り＝芦ご+

U(n, n,

… ，

r

砂 ( 4 . 2 ) をもつミクロ状態を体積要素内

dridr2

…

drNdp1dP2・・・dpN

=

drN dpN

に見出す確率は

1

ZNN!h sN exp(‑/3H) drN dp

尺

_/₃= _ksT¹

⁽ ⁴ ^. ³ ⁾ で与えられる．ここに

exp(‑SH)

はボルツマン因子である．力学量

A

を粒子の位置座標の関数とすれば，そのカノニカル平均は

〈

A〉NVT= 1

N! h3N zN

I I

^Aexp(‑SH) dpN drN = JN

I

^Aexp(‑SU) drN

( 4 . 4 ) で与えられる．ここで

ZN(V, T) = N! ¹h3N

I I

exp(‑SH) dpN drN

=(三）

³^N

J ! ^仏

⁽^V^,T)

1 1

＝心 N!QN(V, T)

を分配関数

(partitionfunction)

と呼ぶ

. A

は

A = h

✓

21rmkBT

(4.5

a )

(4. 5 b)

で定義され，熱ド・ブロイ波長

(thermalde Broglie wave length)

という長

さの次元をもつ量である．因子

(1/N!)

は同種粒子の識別不可能性から由

来する．また

(9)

7 解析力学の復習

一般化座標の時間変化を記述するラグランジュ

(1736‑ 1813)

の運動方程式が導かれる．この方程式は作用が停留値をとる，つまり（配位空間における）作用の変分

=0

というハミルトンの変分原理からも導かれる．これはオイラーの変分法において変分関数にラグランジアンをあてはめたものである．さらに一般化運動鼠を一般化座標とは独立な変数とし，ラグランジアンからハミルトニアンにルジャンドル変換すると，ハミルトン

(1805‑1865)

の正準方程式が導かれる．この方程式は位相空間におけ

る作用の変分

= O

という条件からも導かれる．この条

件をみたす正準変数の変換が正準変換である．

(10)

144 第7章解析力学の復習

7.1

はじめに

ニュートンの運動方程式は

3

次元の直交座標系で表されている．しかし質点系の位置を直交座標系の座標を用いて表すのが常に便利であるとは限らない．要は質点系の位置を一義的に指定でぎる変数であればよいわけである．

一般化座標を導入しニュートンの方程式をかき換えたのがラグランジュの運動方程式である．ラグランジュ形式で表した運動方程式はとくにホロノーム拘束系の場合に有効に応用される．また粒子の位罹座標以外の変数を新しい自由度として導入することができる．圧カ一定の

M D

法では系の体積が座標として導入され，温度一定の

M D

法では系と熱浴との熱エネルギーの授受を担当する変数が力学座標として導人される．このように本来の力学系の

自由度の他に新しい自由度が導入された系を本書では拡張系と呼ぶ．

ラグランジュの運動方程式はニュートンの方程式と同様に，座標を従属変数とする

2

階の微分方程式である．一般化速度は一般化座標とその時間微分の関数で与えられる．ハミルトン形式では，新たに一般化運動量を定義し，

一般化運動量と一般化座標を互いに独立な従属変数とする

1

階の連立微分方程式にかき換えられる．こうして導かれたハミルトンの正準方程式は形式的に美しい形をとり，正準変換に対して不変であって，系の時間発展は時間発展演算子を用いた見通しよい形式で表される．この時間発展演算子を差分化して分解した差分法はシンプレクティック差分法として最近注目されている．

このように分子動力学の方法は解析力学を基礎として著しく発展してぎ

た．本章ではその基礎の部分を復習する．

(11)

7. 2 ラグランジュの運動方程式 145

7.2

ラグランジュの運動方程式

7.2.1 一般化座標

拘束条件

(5.2)

のもとで，

N

個の質点から成る質量叩の質点系の運動方程式を

m~!°/ ⁼

^FJ ⁺

^R^J^, ^j⁼

¹ ^, ² ^, ^… ^， ^N, ⁽ ⁷ ^. ¹ ⁾

Rj=一区入a l

0 ( ] )

a

a~I Orj' a=l,2,・・・,l (7.2)

とする．この方程式を一般化座標で表した運動方程式，すなわちラグランジュの運動方程式

(Lagrange'sequation of motion)

にかき換えることを考えよう．

この質点系の位置が，

r

個の変数小，

q2,

…，のの

1

価連続関数によって次式で指定されるとしよう．

r,=r,(q1,q2,・・・,qr,t), i=l,2,・・・,N. (7.3)

この関数の具体的な形は系の構造に関係するばかりでなく，変数の選び方にも関係する．この変数

qr,qz,

… ，

qr

を一般化座標

(generalizedcoordi‑ nates)

と呼ぶ．その個数

r

も一般化座標の選び方に関係し，一般に

qr,q2,

…，小は独立とは限らないので，系には

P

個の拘束条件

如

(q1,q2,・・・,qr,t)=O, a=l,2,・・・,P (7.4)

が課されているとする．直交座標での［個の条件式 ( 5 . 2 ) を一般化座標でか

き表したのが

P

個の条件式

(7.4)

であるとする.

p

<

l

と考えてよい．このように座標の関数で表される拘束をホロノーム拘束と呼ぶことは第

5

章で述べた．さらに時間

t

を陽に含むときレオノーム拘束

(rheonomouscon‑ straint),

含まないときスクレロノーム拘束

(scleronomousconstraint)

とい

う．また式

(7.4)のように一般化座標の関数として表せない拘束を非ホロノ

ーム拘束と呼ぶ．拘束条件式に一般化座椋の時間微分（一般化速度という）

(12)

10 長距離力の計算

エワルドの方法ではクーロンポテンシャルを近距離と

遠距離に分け，近距離部分は実窄間で，遠距離部分はフ

ーリエ変換した波数空間で計算することによって計算の

負荷を少なくする．こうしてポテンシャルエネルギーと

カの公式が導かれる．この方法はポテンシャルが距離の

逆べきの場合にも拡張できる．ついで，粒子系の計算負

荷を減じる方法として多重極展開と空間分割を述べる．

(13)

218 第 10章長距離力の計算

10.1

はじめに

電解質水溶液，固体電解質など荷電粒子間にはクーロンカが働く．クーロンカは

LJ

ポテンシャルに比べてエネルギーの値が大きいばかりでなく，カの作用範囲がはるかに遠距離まで及ぶ．このため力の作用距離を切断する近似はよくない．クーロンカを精度よく計算する方法にエワルドの方法（略してエワルド法）

[10. l]

がある．これはイオン結品の凝集エネルギーの計算に開発された方法であって，結晶の単位セルをシミュレーションの基本セルに対応させ周期境界条件を課する．クーロンポテンシャルを近距離部分で効くポテンシャルとその残りの部分に分け，残り部分のポテンシャルは逆格子空間における級数で表し，力の作用距離を切断することなく遠い距離からのカを取り入れるのである．この方法では力の計算回数が粒子数

N

に対して

O(N

りで増加する．エワルドの方法は距離の逆べきポテンシャルの場合に拡張できる

[10.2].

これに対して粒子に働く長距離力として重力（星を粒子とみる）の場合も含め，遠距離にいる粒子からの力を粒子群からの力としてまとめて取り扱い，

計算時間を

O(NiogN)

あるいは

O(N)

に短縮しようとする方法が最近開発されてきた．遠距離粒子群のまとめ方に，重力の場合の

8

分木構造による空間分割

[10.3,10.6],

クーロンカの場合には空間のセル分割と多重極展開を併用したものがある

[10.4,10.5].

計算機の能力上昇に伴って百万個の粒子系を取り扱うことも可能になってきた．この方法は水溶液中の

1

個のタンパク質など周期境界条件を課さない系にも有効である．なお重力多体系の数値計算については文献

[10.7]

にまとめられている．

この章ではエワルドの方法と，遠距離にある粒子群をまとめて取り扱う方

法の例を紹介しよう．

(14)

10.2 工ワルドの方法 279

10.2

工ワルドの方法

10.2.1 ポテンシャルエネルギーと力の公式

エワルド法では，荷電粒子の受ける力として基本セル内の他の粒子，

びすべてのイメージセル内の粒子からの力を取り人れる．

およ

体系は

N

個の正負イオンからなり，基本セルは体積

V =

いの立方体とし，

条件

イオン

j

の電荷を

Zje

とおく．系全体として電荷は中和しており，

N

~Z;e = 0

;~1 (10.1)

が成り立っているとする．エワルド法はこの条件がみたされなければ使えな

し‘• その理由は次節で述べる．格子欠陥など電荷が中和していない場合に応用するには，中和の条件が成立するように電荷

Z;e

を補正する必要がある．

体系のクーロンエネルギー

(/Jc

はイメージセルのイオンからの寄与を含めると，

1 N N

(})c =

I :

~,

I :

^Z^j^Z^k^e²

2 た 1k~1

^v47fcolL11+rk~ri (10.2)

で与えられる．ここでベクトル

II=(l.Jむ））Y, 叫；））x, l.Jy, l.Jz =

0 ,

土

1 ,

土

2 , … は基本セルおよびすべてのイメージセルの位置を示し， vについてはすべてのセルについて和をとる．

項は除かなければならない．

ただし基本セル

II= (0, 0, 0)

の場合，

このことを 2 いの＇で示す．

k=j

の

これから基本セル内のイオンの位置

r

を基本セルの一辺の長さ

L

を単位にして測り，

r =Ls= (Lsx, Lsy, Lふ），〇:'.':'.:Sx, Sy, Sz'.':'.: 1 (10. 3)

によって無次元化した位罹座標

S

を用いることにする．

れた公式は圧カ一定の

M D

法，ギー

(})c

は

M C

法にそのまま使える．

こうすると，得ら

クーロンエネル

(15)

13 多電子系の鼠子力学

ー非経験的方法への準備一

水索状原子のシュレーディンガー

(1887‑1961)

波動方程式の固有関数，固有値をまとめておく．多電子原子系の電了の運動を静止した原子核のクーロン場のもとで記述するボルンーオッペンハイマー

(BO)

近似を導入する．

波動方程式と等価なアプローチの変分原理は多電子原子系の近似解を得る見通しをたてやすい．ヘリウム原子について原＋軌道関数を例示し，多電子問題を

1

体問題に還元するハートリーの方法を述べる．これにはつじつまの合う場 (SCF)の方法が登場する．しかしこの方法では，波動関数を

1

電子波動関数の積で表したが，パウリの排他原理によって反対称な波動関数で表さなければな

らなくなる．

(16)

302 第13章多電子系の量子力学ー非経験的方法への準備

13.1

はじめに

原子，分子間の相互作用ポテンシャルを，

LJ

ポテンシャルのようにパラメータを含み，原子，分子の位罹の関数で表して物質の構造，物性を調べる方法を経験的方法という．第

2

章で述べた相互作用ポテンシャルは経験的方法によっている．他方この相互作用を原子，分子を構成する原子核と電子の集合として量子力学によって理論的に求め，電子，原子核の質量，電荷，スピン以外にパラメータを含まない相互作用で表現する方法を非経験的方法

(ab initio method)

と呼ぶ．とくにこの方法を取り入れた分子動力学を第

1

原理分子動力学と呼んでいる．そのなかでカー—パリネロの法 (Car-

Parrine‑ llo method)

はよく知られている．本章ではそれにいたる多電子系の星子力学的取り扱いの理解を助ける必要な事項を述べて，次章のハートリー—フォ

ック

(Hartree‑Fock,HF)

の方法への準備とする．

13.2

水素状原子

電荷

Ze

の原子核と電荷ー

e

の電子

1

個からなる水素状原子（ヘリウムイオン

He

，十リチュウムイオン

Li2+

なども含まれる）について，電子の座標を

re,

原子核の座標を

R

とかけば，シュレーディンガーの波動方程式は

［がふが— _2M — _2m

ふ十

V(re‑

R l ]

^砂⁽^r^e^,^R⁾=

E

^砂⁽^r^e^,^R⁾

( 1 3 . 1 ) である（△ はラプラス演算子）．ここで電子と原子核とのクーロンポテンシャルを

V(re ‑R) = ‑ ¹ ^Ze2

47fco Ire‑RI

⁽ ¹ ³ ^. ² ⁾

とおく.

Eo

は真空の誘電率である.

2

粒子系を解くには古典力学と同様に重

心

G

の運動と相対運動に分離する．重心の座標と相対座標を

(17)

13. 2 水素状原子 303

RG = 1

m + M (mre + MR), r = re ‑ R (13. 3)

とおき，系の質量

Mc=M+m,

換算質景

(1/μ=1/M + 1/m)

を用いると

1 ふ十 1 △

e= 1 1 M m Mc

^{ふ十}

μ

^△ であるから，式

(13.1)

は

[—贔ふー;~ ^△

^+ V(

パ ]

lfl'(r, Re)= E'

叩

(r,Re)

となる.

E'=Ec+E, lfl'(r,Rc)=cp(Rc)</!(r)

とおいて変数分離すれば，

重心および相対運動のシュレーディンガーの波動方程式

[~2位叫り(Re)=

^Ee

り

^(R

^心

⁽¹³^.⁴⁾

[ ̲lt̲

△十

V(r)]り(r)= E rp(r), r = I rl (13.5) 2μ

か導かれる．体積

V(=

じ）の立方体内で重心運動の固有関数は，周期境界条件

tf;(RG + Lv)

=炒

(R

砧

I.I= (J.Jx, レY,叫，

を課すると次式で与えられる．

lJx, lJy, lJz =

0 ,

土

1 ,

土

2 , …

伽

(RG)= r1 v exp(ik・RG), ^{k =}

²

^L

⁷ ^f

^I^I^.

(13.

6 )

( 1 3 . 7 ) 原子核の質量は電子に比べてはるかに重いので換算質量は電子の質量でおき換えても差し支えない．以下ではμ の代わりに電子の質量

m

とおく．

13.2.1 水素状原子の固有状態

この節では中心カポテンシャル

V(r)

のもとでの電子の波動方程式の固有

関数，固有値を公式の形でまとめておく．

(18)

15 第

1

原理分子動力学法

密度汎関数 (DF)理論の碁礎は，系の基底状態のエネルギーが電子密度の汎関数で与えられることにあり，変分原理によって

1

体問題のコーンーシャム

(KS)

方程式が導かれる．方程式のなかの交換相関ポテンシャルは局所密度(LD)近似のもとで一様電子ガスの交換相関ポテンシャルで近似される．全電子についての解がSCF法で得られる．原

f

を価電

f

と原子芯に分けると，価屯

f

に働く原子芯のポテンシャルは内殻電子で一部遮蔽され弱くなった正電荷の原子核のつくるポテンシャルである．

このポテンシャルとして擬ポテンシャルの考えが生れ，

価屯子の波動方程式の解は全電子の場合よりはるかに容易に求めることが可能になった．さらに，原子核の運動と電了雲の時間変化を同時に記述する方法かカーとパリネロによって開発された．

(19)

360 第 15章第 1原理分子動力学法

15.1

はじめに

HF法は多電子系の基底状態を調べる問題を 1電子問題に還元しようとする近似法であった．密度汎関数理論

(densityfunctional theory,

DFT)は多電子問題を

1

電子問題にかき換える基礎を与えるもので，コーン—シャム方程式

(Kohn‑Shamequation, KS)

はシュレーディンガーの波動方程式の役割をする

[15.8‑15.10].

とくに電子密度の局所密度近似

(local density approximation, LDA)

は固体や分子の電子バンド理論，凝集機構の説明に用いられ，成功を収めてきた．しかし

LD

近似では説明できない基底状態には一般化勾配近似

(generalizedgradient approximation,

GGA)が開発されるなど，発展過程にある．

DF理論のように電子，原子核の質量，電荷など粒子の属性に関する数値を与えるだけで，系のハミルトニアンにパラメータを含まない方法は非経験的方法，

abinitio

法，あるいは第

1

原理的方法などと呼ばれ，その対象は主に基底状態の静的構造の研究に絞られていた．この方法を原子核のダイナミックスが取り扱えるように分子動力学法に広げたのがカー—パリネロ

(Car‑Parrinello, CP)

の方法

[15.17‑15.19]

であって第

1

原理分子動力学法とも呼ばれる．

CP

法では原子核（あるいは原子芯）の間に働く力を，時間変化する電子雲

の分布から直接求めて原子運動のダイナミックスを取り扱う．本章では DF

理論の基礎的考えを述べた後，原子核の運動方程式を導<.

CP

法は結晶状

態よりも液体など秩序の乱れた状態に応用され，有効性を示してきた．系の

バルク状態については周期境界条件を課する．このとき系の基本セルは結晶

の単位細胞にあたるので，固体の電子構造の理論がお手本になる

[15.6, 15. 7].

また原子核のクーロン引力は内殻電子と原子核をまとめた原子心で

おき換える近似を行う．それには原子核のクーロンポテンシャルを原子芯の

擬ポテンシャル

(pseudpotential)

でお苔換える．

(20)

15.2 HF方程式から KS方程式へ 361

本章で立ち人っていない

GGA,DF

の統計理論などについては教科暑

[15. 2]および，

より専門的な教科書

[15.3]

を参照されたい．

15.2 HF

方程式から

KS

方程式へ

ハートリー—フォック (HF) の方法を整理しておこう．

N

電子系のハミルトニアンを次のように与える．

BO

近似のもとで

fl=~( 一；叫 +~vext(r;)

_z

+ミ上

z z<; ru

vexi(r)

=―こ

^Za

a lr‑R』・

ここで

vexl(r)

は原子核のクーロン引力場である．

ヽ‘~ヽーー／

1 2

．

5 5 1 1

︵

基底状態を記述する

HF

方程式は

N

電子波動関数に規格化されたスレーター行列式 1 J fを用い，

ルギーの期待値〈

IJfI fl I IJf

〉が停留値をとるという条件から導かれた．

で電子密度

エネそこ

n(r) =~

苓

I(f!k,!5(r)l2 =~n15(r),

および交換ポテンシャル ( 1 4 . 5 0 )

瓜

(r)rp(r) = ‑

f

^d^r^'ⁿ_Ir‑r^f⁽^r^,^r_'^'_I⁾^r^p⁽^r^）

(15. 3)

( 1 5 . 4 )

を導入して，

HF

方程式をかぎ換えると次のようになる．

[— t

^•

^{＋戸}

⁽^r⁾⁺

f

^d^r^'^I

;~';,1

^{+ 成(}^r⁾^]^!^{^!^k^び⁽^r^)= s^Z^!^f^!^k^,^<¹⁽^r⁾

( 1 5 . 5 )

N

電子系の全エネルギーの式 ( 1 4 . 2 9 ) は次のようにかぎ換えられる．

鱈＝図九+ t f f

^nl~>^n~~? ^dr dr'+ E

^翌

h;

〈 =

^g

― ‑ !

—• + vext(r) rp,

. 〉

︑

＼ノ

︶

6 7

．

5 5 1 1

分子シミュレーション

ー

分子シミュレーション

章分子シミュレーション

本書の目的

分子シミュレーションの方法は古典系については古典力学，統計力学を基 礎としている．量子系にはさらに量子力学が加えられる．

多数の粒子を体精一定，エネルギー一定のもとで運動させるとどんなマク ロの性質があらわれるのか，計算機を用いて運動方程式を数値積分し，力学 鼠の時間平均をとおしてマクロな性質を調べる分子動力学法

法）と，統 計力学の集団理論に基づき多数の粒子の空間配置を確率法則によって生成さ せるモンテカルロ法

法）とによる分子シミュレーションは前世紀の

年代に始められた．力学法則による結果と統計法則による結果か一致するこ

と，および剛体球系の固相一流動相転移の存在が明らかにされ，レナード・

さらには生体高分子など生物の分野でも実験，理論に加えて物質科学の有力 な研究手段として，その地位を確立してきた．

原子，分子の電子構造を非経験的に求める道は分子軌道法として化学の分 野で開拓されてぎた．たとえば，水

非経験的分子動力学法（カー—パリネロの方法， CP 法）が 80 年代半ばに開発

され，原子，イオン間の力の媒介をする電子雲の動きを直接調べることが可

2

分子間力，分子モデルと計算機実験

単純な物質からより複雑な物質へ，いくつかの例をあげ

る ．

はじめに

希ガス（不活性ガス）原子間に働く力は

体相互作用（対ポテンシャル，

によってよく記述される．対ポテンシャルは希ガスの他に単純な 金属やイオン系によく用いられている． しかし構造をもつ分子では

体力，

体力が無視でぎなくなる．このような分子の分子間力をどう表現するか，

対象をどのようにモデル化するか，系が複雑になってくるはどモデルの特徴 のどこがよくて，何がわるいのかその見分けが困難になってくる．ファン・

デル・ワールスの状態方程式から話を始め，単純な原子モデルからより複雑 な分子モデルヘ，計算機実験の方法にも触れながら後章への準備とする．

ファン・デル・ワールスの状態方程式

ニュートンは気体が圧力を示すのはアトム間に斥力が存在するからだと仮

定し，その斥力を

1 /

0 ,

い は 最 近 接

アトムとの距離）とおいた．マクスウェルは気体運動論において罪離の逆ベ

きの斥力を論じ， とくに

の斥力（ポテンシャルではない）の場合，気体の

︐

分布関数を知らなくても粘性係数など輸送係数を計算できることを示した．

この分子モデルは彼の周囲の人々は別として，当時のヨーロッパ大陸の科学 者には受け入れられなかった．その理由の

つに引力が含まれていないこと があったそうである（文献

分子間力が斥力と引力からなることは，ファン・デル・ワールス

の状態方程式によって示された

年 ） .

モル当たりの式は

である.

は定数．この状態方程式を

a

( 2 . 2 ) と変形すると，第 1項の (V‑b )は，分子が大きさをもっため分子の運動で できる空間の容積が

だけ減少したためとみることができる．つまり分子 間に斥力が働くため，ある距離以上は互いに近付けないことを意味する．圧 力は第

項の分だけ弱められる．これは分子間に引力が働いて引ぎ合うため，

p o

＞

4

デ ー タ 解 析

エネルギー一定の系の長時間乎均，温度一定の系のカ

ノニカル平均について，例として状態方程式，動径分布

関数，静的構造因子，および自己拡散係数，速度相関関

数の求め方を述べる．

はじめに—長時間平均とカノニカル平均

粒 子 数

体 積

で相互作用ポテンシャル

，

り の 質 点 系 を考えよう．この系はエネルギー

が一定の保存系であるから，前章で述 べた

法は近似的にエネルギー

を一定に保つとみてよい．そこで粒子 数

体 積

およびエネルギー

を一定に保つ

法を

法と 略記することにする．今後，同様の略記法を圧カ一定，温度一定の場合にも 用いる．

初期条件を与え，

法によって力学系の時間発展を追跡していく と，各粒子の位置と運動量が時間の関数として与えられる．運動エネルギー やポテンシャルエネルギーなどの力学鼠を時間の関数としてグラフに描く

この平均は統計力学でいう長時間平均である．粒子の座標と運動量（ある いはその一部分）の関数である力学量を時間

の関数として

=

と記すと，

の長時間平均

は

分子シミュレーションの方法は古典系については古典力学，統計力学を基礎としている．量子系にはさらに量子力学が加えられる．

多数の粒子を体精一定，エネルギー一定のもとで運動させるとどんなマクロの性質があらわれるのか，計算機を用いて運動方程式を数値積分し，力学鼠の時間平均をとおしてマクロな性質を調べる分子動力学法

法）と，統計力学の集団理論に基づき多数の粒子の空間配置を確率法則によって生成させるモンテカルロ法

さらには生体高分子など生物の分野でも実験，理論に加えて物質科学の有力な研究手段として，その地位を確立してきた．

原子，分子の電子構造を非経験的に求める道は分子軌道法として化学の分野で開拓されてぎた．たとえば，水

る．

によってよく記述される．対ポテンシャルは希ガスの他に単純な金属やイオン系によく用いられている．しかし構造をもつ分子では

対象をどのようにモデル化するか，系が複雑になってくるはどモデルの特徴のどこがよくて，何がわるいのかその見分けが困難になってくる．ファン・

デル・ワールスの状態方程式から話を始め，単純な原子モデルからより複雑な分子モデルヘ，計算機実験の方法にも触れながら後章への準備とする．

いは最近接

きの斥力を論じ，とくに

この分子モデルは彼の周囲の人々は別として，当時のヨーロッパ大陸の科学者には受け入れられなかった．その理由の

つに引力が含まれていないことがあったそうである（文献

年） .

^a

( 2 . 2 ) と変形すると，第 1項の (V‑b )は，分子が大きさをもっため分子の運動でできる空間の容積が

だけ減少したためとみることができる．つまり分子間に斥力が働くため，ある距離以上は互いに近付けないことを意味する．圧力は第

データ解析

粒子数

体積

りの質点系を考えよう．この系はエネルギー

が一定の保存系であるから，前章で述べた

を一定に保つとみてよい．そこで粒子数

体積

法と略記することにする．今後，同様の略記法を圧カ一定，温度一定の場合にも用いる．

法によって力学系の時間発展を追跡していくと，各粒子の位置と運動量が時間の関数として与えられる．運動エネルギーやポテンシャルエネルギーなどの力学鼠を時間の関数としてグラフに描く

この平均は統計力学でいう長時間平均である．粒子の座標と運動量（あるいはその一部分）の関数である力学量を時間

平均）と等しいと考える．一方，統計力学の集団理論によれば，ミクロカノ

一定のもとで

とする極限）においてカノニカル集団，グランドカノニカル

法において最も璽要な統計集団はカノニカル集団である

個の同種粒子からなる古典系について，カノニカル平均は次のように定義される．カノニカル集団に属する体系は体積 (V)一定，温度 (T)一定のミクロの状態をとる．この体系がエネルギー

に見出す確率は

⁽ ⁴ ^. ³ ⁾ で与えられる．ここに

を粒子の位置座標の関数とすれば，そのカノニカル平均は

J ! ^仏

の運動方程式が導かれる．この方程式は作用が停留値をとる，つまり（配位空間における）作用の変分

次元の直交座標系で表されている．しかし質点系の位置を直交座標系の座標を用いて表すのが常に便利であるとは限らない．要は質点系の位置を一義的に指定でぎる変数であればよいわけである．

法では系の体積が座標として導入され，温度一定の

法では系と熱浴との熱エネルギーの授受を担当する変数が力学座標として導人される．このように本来の力学系の

ラグランジュの運動方程式はニュートンの方程式と同様に，座標を従属変数とする

階の微分方程式である．一般化速度は一般化座標とその時間微分の関数で与えられる．ハミルトン形式では，新たに一般化運動量を定義し，