数学・数理科学からのアプローチによる分子動力学計算

Abstract

In the present report, we introduce our new three approaches based on the mathematics / mathematical sciences to the classical molecular dynamics calculations; 1) the approach by the analytical mechanics instead of the classical mechanics; 2) the approach to the periodic boundary condition by the torus model; 3) the approach by the mathematics-based programming substituting for the procedural (imperative) programming. Then, we found that these approaches are very effective to make calculations simpler, more compact, steadier, and firmer. Thus, it was concluded that the mathematical / mathematical scientific approaches are the promising ones to overcome various problems that we confront in the computational chemistry.

Translated Abstract

In the present report, we introduce our new three approaches based on the mathematics / mathematical sciences to the classical molecular dynamics calculations; 1) the approach by the analytical mechanics instead of the classical mechanics; 2) the approach to the periodic boundary condition by the torus model; 3) the approach by the mathematics-based programming substituting for the procedural (imperative) programming. Then, we found that these approaches are very effective to make calculations simpler, more compact, steadier, and firmer. Thus, it was concluded that the mathematical / mathematical scientific approaches are the promising ones to overcome various problems that we confront in the computational chemistry.

1 序論

現在，コンピュータ化学は，材料産業分野において広く実活用化されつつあるものの，その高精度ゆえの実活用限界の低さ，計算プログラムの膨大化・複雑化，実装技術の古さ，手続き(命令)型プログラミングの限界，他，様々な課題・問題も顕在化している [1,2,3]．我々は，現状からのブレークスルーを探るべく，コンピュータ化学への数学・数理科学的なアプローチを試みている．本報告では，分子動力学(molecular dynamics: MD)計算 [4]を題材とした，新たな3つのアプローチについて紹介する．まず第1は，古典力学(classical mechanics: CM)に代わる解析力学(analytical mechanics: AM)によるアプローチ，第2は，周期的境界条件(periodic boundary condition: PBC)へのトーラスモデルによるアプローチ，そして第3は，手続き型プログラミングに代わり，数学・数理科学の知識と技術を積極的に関数型プログラミングに取り入れた数理ベースプログラミング [5]によるアプローチである．以下，それぞれについて説明する．

1.1 解析力学

CMにおいては，加速度と力の関係を表す運動方程式を微分方程式で表現し，その微分方程式の解を求めることにより，運動をシミュレートする．また，運動を記述する変数は，空間での静止直交座標系による位置座標を表し，1次元，2次元，または3次元空間の方程式で運動をシミュレートする．一方，AMにおいては，使われる変数は静止直交座標系で直接表す必要はない．運動を記述する方程式として，例えば，運動エネルギーと位置エネルギーの差を表すラグランジュ(Lagrange: L)関数を構成し，運動方程式は，L関数の積分から定まる作用汎関数に対する最小作用の原理から導かれる．本報告では，AM的アプローチによりMD計算を行う新たな手法を提案する．

1.2 トーラスモデル

コンピュータ化学において，PBCでは，複数の粒子間の運動の計算を格子状に分割した空間を想定して行う．そのとき，各格子内の粒子の状態は相対的に同じであると仮定して，格子内，および隣接する格子にある粒子間の運動をCMによる運動方程式で記述し，時間発展させることによりシミュレーションを行う．固体材料において，結晶は繰り返し構造になるので，格子状に分割された隣接する格子内の状態は相対的に同じであるという仮定は自然である．一方，シミュレーション計算は，操作的であり，一つの格子内の粒子を隣接する粒子と捉えてエネルギー計算を行うなど，計算上の便法を使っており，系全体の数学的な定式化は存在しない．シミュレーションアルゴリズムとして，様々な解は求められているが，その解法は，操作的，かつ恣意的である．

今回，我々は，一つの格子が繰り返す構造を表現する座標系で状態を表現し，その座標系内の粒子の運動を表すL関数を構成し，AMの考え方に従い運動方程式を求め，その方程式の解を微分方程式の数値解法により計算する手法を新たに提案する．繰り返し構造を表現する座標系としては，n次元トーラスを用いる．本手法の特徴は，繰り返し構造が完全に数学的に微分方程式で定式化されている点と，シミュレーション計算について，純粋に微分方程式の数値解法アルゴリズムを用いている点である．この手法，すなわち，問題の定式化とシミュレーション計算を分離することにより，プログラム実装時のミスを減らすことも検討している．

1.3 数理ベースプログラミング

現在のコンピュータ化学分野における主流である手続き型プログラミングについては，それ自体が計算プログラムの膨大化・複雑化を招いている主要因でもあり，増大する一途のバグ・エラーの解決は困難を極める．我々は，上記の問題を克服し得るものとして，数学・数理科学を積極的に関数型プログラミングに取り入れた数理ベースプログラミングに着目し，これをMD計算のプログラム作成に適用した．そして，同プログラミング，およびAM的手法，トーラスモデルを用いることにより，計算のシンプル化，コンパクト化，安定化，強固化，そしてバグ・エラーの大幅な削減が図れることを見出した．

2 計算方法

Figure 1はPBC付きの1次元空間$ℝ$を表している．ここで，Figure 1に現れる正の実数$r$が$r=1$の場合を考える．この場合, PBC 付きの 1 次元空間$ℝ$は, $n$を整数とするとき, $2n\pi$ から$2\left(n+1\right)\pi$までの長さ$2\pi$の区間に分割される. 各区間に現れる粒子$m_i$の座標の集合を考えると$\text{{}{x}_i+2n\pi \text{|} n\in \mathbb{Z}\}$ ($0\le x_i<2\pi$) となる. 粒子間の距離は区間内の距離とは限らない. それらは最も近い粒子間の距離で決まる ので, PBC下では次のような条件分岐による計算が必要となる   

$m_i と m_j の距離=\{\begin{array}{c}\begin{array}{c}\left|x_i-x_j+\pi \right| (x_i-x_j<\pi )\\ \left|x_i-x_j\right| \left(\left|x_i-x_j\right|\le \pi \right)\end{array}\\ \left|x_i-x_j-\pi \right| (\pi <x_i-x_j)\end{array}$(1)

Figure 1.

 The Schematic figure of the three particles (m₁, m₂, m₃) on the one-dimensional periodic lattice model.

Figure 2は「PBC付きの1次元空間$ℝ$に配置された粒子」と「円周上に配置された粒子」が対応することを示している．二次元空間上の空間$S=\left\{\left(\cos \theta ,\sin \theta \right)\in {ℝ}^2:\theta \in ℝ\right\}$は一次元トーラスと呼ばれるが，ここでは1次元周期空間と呼ぶ．今，Figure 1に現れる正の実数$r$が$r=1$の場合を考えているので，1次元周期空間$S$の半径も$r=1$であることに注意する．そして「位置$m_i$に配置された$ℝ$の粒子」は「$S$に配置された円上の位置$\left(\cos {\theta }_i,\sin {\theta }_i\right)$」に配置される．ここでは$r=1$から$m_i={\theta }_i$であり，実際に「PBC付きの1次元空間$ℝ$に配置された粒子の位置」と「円周上，つまり，一次元周期空間$S$に配置された粒子の位置」が対応していることがわかる．また，「$ℝ$に配置された粒子の運動」と「$S$に配置された粒子の運動」も対応する．

Figure 2.

 The Schematic figure of the three particles (m₁, m₂, m₃) on the one-dimensional torus model.

時刻$t\in ℝ$における，ある空間上での点を$p\left(t\right)$で表す．そして，$p\left(t\right)\in ℝ$の場合は「$p\left(t\right)$の$ℝ$での運動」を表し，$p\left(t\right)\in S$の場合は「$p\left(t\right)$の$S$での運動」を表すことにする．$ℝ$上に周期的に配置された粒子たちの運動が$ℝ$上のL運動方程式の解であるとき, 対応する1次元周期空間$S$上に配置された粒子が満たす$S$上のL運動方程式を考察する．そして, その1次元周期空間$S$上の運動方程式を解き，その解を$ℝ$上に周期的に配置することによって，求める粒子たちの運動を観察する様な方法を提案する．L関数$L$は，ポテンシャル関数$T$と運動エネルギー関数$V$の差 $L=T-V$で表現される．そして，L運動方程式は，次の形によって与えられる．   

$\frac{d}{dt}\left(\frac{\partial L}{\partial {\theta }_i}\right)-\left(\frac{\partial L}{\partial \dot{{\theta }_i}}\right)=0$(2)

ここで，時刻$t\in ℝ$で位置${\theta }_i\left(t\right)$と位置${\theta }_j\left(t\right)$にある粒子たちの距離$r_{ij}$ は，PBC付きの1次元空間$ℝ$(つまり，Figure 1)で与えた様な分岐を伴う距離とは異なり，次の様な形によって一意的に与えることが出来る．   

$r_{ij} =\left|2\cdot {\cos }^{-1}\sqrt{1-{\left(\frac{\sqrt{{\left(\sin \left({\theta }_i\right)-\sin \left({\theta }_j\right)\right)}^2+{\left(\cos \left({\theta }_i\right)-\cos \left({\theta }_j\right)\right)}^2}}{2}\right)}^2}\right|$(3)

ここではポテンシャル関数をLennard-Jones ポテンシャル [4]から与えることにする．つまり，L関数を構成するポテンシャル関数と運動エネルギー関数が，各々，次の$T$と$V$によって与えられる．ここで粒子の質量は$1$である．   

$T=\sum _{i<j}\frac{1}{4}\left(\frac{1}{r_{ij}^{12}}-\frac{1}{r_{ij}^6}\right), V=\sum _i{\dot{{\theta }_i}}^2$(4)

3 計算結果・考察

関数$\varphi :ℝ\to S$を$\varphi \left(x\right)=\left(\cos x,\sin x\right)$とする．また，関数${\varphi }^{-1} :S\to 2^{ℝ}$を${\varphi }^{-1}\left(x\right)=\left\{x+2n\pi :n\in \mathbb{Z}\right\}$(ただし，$\left(\cos x,\sin x\right)\in S$である)とし，関数$\overline{{\varphi }^{-1}} :2^S\to 2^{ℝ}$を$\overline{{\varphi }^{-1}}\left(X\right)=\cup \left\{{\varphi }^{-1}\left(s\right):s\in S\right\}$とする．このとき，$S$上のL方程式の解$P=\left\{p_i\left(t\right) :t\in ℝ\right\}$に対して$\overline{{\varphi }^{-1}}\left(P\right)$がPBC付きの1次元空間$ℝ$上でのL方程式の解となるように，$S$上のL方程式を定めることが重要である．特に，Figure 3は次の初期条件から計算された粒子の位置の変化を表している．また，正の実数は$r=1$としており，各粒子の質量は$1$としている．以下に，m_iに関する初期値を示す．   

$m_1=0.23, m_2=2.09333, m_3=4.18666$

  

$\dot{m_1}=-0.0000260, \dot{m_2}=0.0000330, \dot{m_3}=-0.0000160$

Figure 3.

 The results of MD calculations for the three particles (m₁, m₂, m₃) on the one-dimensional torus model. The horizontal and vertical axes give the time steps and positions of three particles, respectively.

Figure 3が示す様に，今回与えた新しいモデルによる計算結果が妥当な結果を与えていることを確認出来た．

4 結論

本報告では，我々の，コンピュータ化学への数学・数理科学的なアプローチによる試みを紹介した．MD計算について，CMに代わるAMによるアプローチ，PBCへのトーラスモデルによるアプローチ，そして手続き型プログラミングに代わる数理ベースプログラミングによるアプローチの結果，同計算のシンプル化，コンパクト化，安定化，強固化，そしてバグ・エラーの大幅な削減を図ることに成功した．参考までに，Pythonで作成されたMD計算のプログラム [4] の書き換えでは，配列変数の統合と高階関数の使用によるごく軽いレベルの書き換えでも，バグ・エラーの原因となりやすい「forループ」の数の約89% の削減による安定化，強固化，プログラム行数で約9%のシンプル化，コンパクト化，そして，一般的なパーソナルコンピュータでの324分子のMD計算において，平均約12%の高速化を実現した．

本手法の特徴は，ユークリッド空間において周期解を求めるのではなく，トーラスモデルで解を求めるという発想の転換を行った点である．さらに，トーラスモデルでの運動方程式をAMの観点から具体的な方程式によって定めた点である．すなわち，シミュレーションプログラムの繰り返し計算が完全に数学的に微分方程式で定式化され，純粋に周期空間上の微分方程式の数値解法アルゴリズムに基づく計算のみである点である．その利点として，周期性判定のための条件判断等の例外計算が不要となり，プログラムが簡素化され, 理論的な計算量は落ちていることが期待される．トーラスモデルは厳密に数値解を計算しているが，運動方程式が複雑で従来法の様な具体的な高速計算技法が未成熟であるという点が課題である．しかし，モデル依存の繰り返し計算ではなく，微分方程式の数値解法に還元されているので，従来法に依存しない改善も期待される．また，数の計算のための関数だけでなく，数式処理のためのアルゴリズムを実現した関数(例えば，Mathematica言語の求解部分のNDSolve関数)も多々開発されている．さらには，単純な数式変換ではない，新たな等式や式変形の理論が，別の視点から，数学者らにより研究され証明されている．これらを活用した数だけの計算だけではなく，数式の計算も意識したプログラミング作法が，数理ベースプログラミングである．

今後の課題として，具体的なコンピュータ化学のライブラリを関数型プログラミングによるライブラリとして再構築し，数理ベースのプログラミングの事例を拡張することも検討している．

本研究成果は，AIMaP事業 [6]の一環によりコーディネートされた共同研究から得られたものである．

参考文献

• [1]   S. Higai, Ouyou-suuri, Bulletin of Jpn. Soc. for Indus. and Appl., Math., 29, 170 (2019). (in Japanese)
• [2]S. Higai, Frontiers of Computational Sciences, S. Nakamura, A. Makinouchi, T. Adachi, M. Sugimoto supervised, Nakamura Labo. ed., Kindaikagakusha, 193 (2019) (in Japanese).
• [3]S. Higai, H. Kamisaka, CREST-SAKIGAKE-AIMaP Joint Symp. 2019, Tokyo, Japan, 10, 11 Mar. (2019) (in Japanese).
• [4]   R. H. Landau, M. J. Paez, C. C. Bordeianu, Computational Physics: Problem Solving with Python, John Wiley & Sons, Chap. 16 (2015).
• [5]   S. Higai, Y. Mizoguchi, R. Fukasaku, Ensemble, Bull. Mol. Sim. Soc. Jpn, (2021). (to be published) (in Japanese)
• [6]   Advanced Innovation powered by Mathematics Platform (AIMaP), MEXT, Jpn., Official website;https://aimap.imi.kyushu-u.ac.jp/wp/