リアルタイム三次元可視化分子動力学シミュレータMolecula Numericaの開発

阿部 博史, 長嶋 雲兵


Return

1 はじめに

科学技術計算分野において,その主たる開発および実行環境は,ハイエンドのスーパーコンピュータで,開発言語の多くはFORTRANであった.習得が楽で,最も容易にスーパーコンピュータの性能を引き出せるようにコンパイラやライブラリが用意されていた事がこの状況を作っていた.
一方その他のC,C++などの言語の国産スパコン上での性能が芳しくない事から,これらの科学技術計算分野での普及は停滞していた.しかしながら汎用CPUを用いたワークステーションやPC上ではその実行性能は遜色の無いものだった.
この科学技術計算分野におけるFORTRAN寡占の状況は近年の汎用CPUを用いた分散コンピュータが脚光を浴びてから徐々に変化し始めている.実行性能に問題がなくなったため,そのオブジェクト指向言語,そしてテンプレート対応といった面が重要視され始めている.既に,欧米では多くの科学技術系商用ソフトウェアがC++で開発されているし,日本でもC++を用いた開発プロジェクトが立ち上がっている[1, 2].また,そのほとんどがC,C++で提供されるGUIライブラリとの整合性の良さもFORTRAN系言語には無い特徴である.
さて,これらのソフトウェアを動作させるプラットフォームは,ハイエンドでは依然専用CPUマシンも存在するが,多くの研究者が用いるのは汎用CPUベースのクラスターコンピュータであり,その個々のコンピュータはいわゆるPCである.最近では64 bitメモリ空間を使えるOSもPC向けに存在し,メモリも4 GByteを超えて搭載できるPCも存在する様になった.
また,三次元グラフィックス能力もCPUに劣らず進化を遂げ,リアルタイムで三次元表示する能力を充分に有するようになっている.画像表示に用いる座標変換などの演算をGraphics Processing Unit ( GPU ) 上で並列に行えるため,CPUに負荷無く処理が行える.
このような状況を踏まえ,筆者らは高性能PCをハードウェアプラットフォームとして,リアルタイム/対話型処理を可能とし,また,C++のオブジェクト指向型テンプレート言語としての特徴を活かした技術をMDに応用したソフトウェアMolecula Numericaを開発した.
ソフトウェアを設計する際に,「ポータビリティ」,「高速性」,「拡張性」に留意した.この分野のプログラム利用者はUNIX利用者が最も多いと思うが,異なるOSでも同じソフトウェアが動作するのは利用者にとって利便性が高い.計算処理をでき得る限り速くするという事は時間の有効利用につながる.また,新たな力場を導入したいときや,機能を追加したい時に,低コストでプログラムの変更ができる事も重要である.
本稿では,現在の高性能PCの科学技術計算用プラットフォームとしての分析を行い,その有用性を確認した後,オブジェクト指向,テンプレート言語であるC++を用いて分子動力学 ( MD ) シミュレータを実装した事例について紹介し,最期に計算例について紹介する.

2 科学技術計算プラットフォームとしてのPC

近年までPCに用いられるCPUは,その動作クロックの向上によってその性能の向上を図ってきた.昨今では,動作クロックの高速化ではなく,1つのCPUに複数のプロセッサを実装する,マルチコア化により高速化を進めている.
2008年9月現在,最も多くのコアを持つ汎用CPUは東芝などが開発するCELLで,9個のコアを持つ.また,Intel社のXeonが最大4つのコアを持つ.Xeonの理論的な計算性能は単精度浮動小数点演算で最大51.2 GFLOPS ( 4コア× 4演算× 3.2Ghz ) という値となる.この値は,NEC社製SX-9の118.4 GFLOPSとそれほど遜色無い.これには,拡張計算ユニットの存在が大きい.
一般にCPUの実行性能でそのコンピュータの性能を評価しがちであるが,メインメモリとCPUの間でのデータ転送速度も重要である.CPUの計算性能がいくら高くてもデータをメモリからロードするのが遅ければそこがボトルネックとなるからである.スーパーコンピュータとPCとの差別化は,それらのCPUの演算速度が近づいている今,主にCPUと主記憶装置間のデータ転送速度で行われている.既出のSX-9では256 GB/sのデータ転送速度を持ち,一方のXeonでは,12.8 GB/sでその差は20倍ある.これがスーパーコンピュータのスーパーたる所以であるが,値段の差から考えれば納得せざるを得ない.CPUの高速性を活かすためには,メモリからロードしたデータをできるだけキャッシュ上で使い回す事が必要だが,今回は特に行わなかった.
三次元表示ゲームソフトからの要求に応じる形で,三次元グラフィックスの表示能力の向上は凄まじいものがあり,また,その汎用化により安価になった.またそのバスもPCI Expressインターフェースが用いられ始め,着々と高速化が進められている.ハイビジョンクラスの画像データを取得しようとすると,一枚の画像あたり,1920×1080×24 = 6.2 Mbyteのデータとなる.現在のPCI Express 2.0規格での転送速度は,1レーンあたり500 MB/s,一般に複数レーンが用いられるので16レーンだと総転送速度は8 GB/sに達する.そのため,たくさんの描画命令を送信したり,画像のピクセルデータを取得したりするためのコストは非常に低いことになる.
以上の分析にあげられた鍵となる要素について個別に説明する.

2. 1 拡張計算ユニット

PentiumなどのIntel社製プロセッサにおけるSSE, PowePCにおけるAltivec, CellにおけるSPEなど,現在の多くの汎用プロセッサには拡張計算ユニットが内蔵されている.ここにあげたCPUでは全てSingle Instruction Multiple Data ( SIMD ) 型の並列ユニットが拡張されている.SIMDでは,複数のデータを一度の命令で処理をする.そのため単純にデータの並列数倍の高速化が行える.先に例としてあげた全てのSIMDユニットは,128 ビットのデータを同時に扱え,単精度浮動小数点演算では同時に4つの演算が行える.例えばIntel社製のXeon X5470では,3.3 GHzの動作クロックなので, 理論上1 CPUコア当たり, 4×3.3 = 13.2 GFLOPSの処理速度がある事になる.
SIMD命令を使うと速くなるが,機種依存性も高くなる.ポータビリティを重要視すればある程度速度を犠牲にする事も必要となる.Moleculaでは物理量を格納する配列クラスHvectorの中にこのCPUごとに異なるSIMD命令を隠蔽し,クラスの外からは一つのインターフェースで扱えるようにした.

2. 2 共有メモリマルチコア

複数のプロセッサコアでメインメモリを共有するこのアーキテクチャは,科学計算の分野ではNEC社製SX-4をはじめ大きな成功がある.
PC市場でも,既出のXeonでは,4コアの物が流通しており,PCにはこれが2つまで搭載されている.つまり1台のPCに合計8つのCPUコアが搭載されメインメモリを共有しているという共有メモリマルチコアシステムが市販されている.
こういった共有メモリマルチコアのシステムでの高速化は複数の処理単位を同時に複数のコアで処理する事によって行われる.これはマルチスレッド処理と呼ばれている.

2. 3 高性能グラフィックス

よりリアルで,そしてよどみない動きを要求される三次元ゲームソフトウェアからの要求から,GPUも,データ転送速度も高速化の一途をたどり,今や非常に複雑な形状のレンダリングも瞬時に行えるようになった.
バス速度も非常に速いため,分子動力学に特有な多数の原子の描画命令をGPUに送信するのに適している.
分子動力学分野での主な三次元表示は,原子の位置表示である.この処理はCPU側での処理が殆どなく,描画に関わるほとんどの処理をGPU側に任せる事ができる.そのため,この種の三次元処理はシステム全体では非常に効率よく行える.

3 プログラム構成

Figure 1にMoleculaを構成する主要なクラスのクラス図を示す.


Figure 1. Class Diagram of the main classes of Molecula Numerica.

左の四角で囲んだのが,プログラムのデータ構造に関する重要なクラスである原子/分子のモデルデータを定義したクラスである.抽象クラスであるElementクラスから単原子の定義であるAtomクラスの系列,また複数の原子 ( site ) からなるMoleculeクラスの系列がある.白抜き矢印の元のクラスは,矢印先のクラスの機能を継承している.各々のクラスはそれぞれの2体間力のパラメータにアクセスするための関数 ( get_ ) が用意されている.ここでLocusは位置座標クラスである.
右はプログラムの実行に関する重要なクラスを示した.GUIとのインターフェースと表示に関する処理をつかさどるBaseControllerと,原子や分子の重心座標など物理量を保持するWorldクラス,そしてThreadや計算処理を制御するWorldControllerである.Hquatは四元数クラスでデータと演算が定義されている.これは分子の回転運動だけでなく,三次元画像の回転にも用いている.Hvectorはベクトル配列クラスで高速なベクトル演算を提供する.図中で黒菱のついたクラスは矢印先のクラスの中に内包されているという関係,単なる線は関連の存在を示している.

3. 1 Graphical User Interface/Framework

GUIのためのライブラリは多く存在する.最も古い歴史を持つX11ライブラリを始め,Qt,KDEなどがよく知られている。そのよく知られているライブラリの中で,Linux, MacOS X, MS Windowsに対応しているwxWidgets[3]を採用した.
ライセンス条項が非常に緩く,また開発が非常に活発である事や,対応する全てのプラットフォームの開発の足並みがよく揃っている事が採用の大きな理由である.
ポータビリティという観点から,プログラムの構成のうち,GUIに関わる部分をいくつかの限定されたクラスに留め,ライブラリの変更にも最小限で対応できる構成とした.

3. 2 三次元可視化

三次元可視化を行うには三次元グラフィックライブラリを用いるが,現在ではOpenGLに統一されていて,ほとんどのプラットフォームに移植されているためポータビリティが高い.また描画性能も各GPU開発会社がライブラリを最適化しているため非常に高性能である.
MoleculaではGUIライブラリであるwxWidgetsのOpenGLラッパークラスであるwxGLCanvasを親クラスとするView3Dクラスを定義し,その中で三次元処理を行った.GUIライブラリのクラスを三次元処理に用いる事により独立性が損なわれる事になるが,ラッパークラスであるのでそれほど損なわれないと判断した.
また計算の時間発展にあわせて表示を更新するメカニズムは,wxWidgetsのwxFrame::OnTimer関数を用いた.これはある一定の時間間隔で周期的に実行される関数で,表示の更新によく用いられる.ただし,実際の更新はシステムに依存する.
動画像ファイル作成のための,ある時間間隔での表示画像をデータファイルに出力する機能の実装には,OpenGLの関数,glReadPixelsを用いた.以前のPCではこの関数のコストが大きかったために使いづらかったが,現在のPCでは全く気にならない.

3. 3 マルチスレッド

共有メモリ型マルチコアプロセッサで並列化する場合,最もよく用いられるのがマルチスレッド化である.
スレッドとは,1つのプログラムの実行処理 ( プロセス ) の中での一つの処理単位である. プロセスには最低一つのスレッドが含まれる.ここで必要となるマルチスレッドとは,プロセスの中で,同時かつ複数のスレッドが処理されることを言う.
ウェブブラウザのような汎用ソフトでは,データのダウンロード処理などを一つのスレッドとし,複数のダウンロードを同時に行ったりするためにマルチスレッドを用いる.
科学技術計算では高速化のためにこのマルチスレッドを用いる.そのためには,計算負荷の大きい部分を並列化し,それをマルチスレッド化する事になる.
マルチスレッド化するには,ハード専用のコンパイラや,OpenMPやPOSIX Thread ( Pthread ) などのライブラリを用いる.Moleculaでは,ポータビリティのためPthreadを採用した.Microsoft社のWindowsではPthreadが標準では用意されていないが,オープンソースで用意されている[4].
Moleculaでは力の計算部分 ( WorldController::CalcForce ) と時間発展部分 ( WorldController::Advance ) がマルチスレッド化されている.各並列スレッドはCalcForce,Advanceのそれぞれの前後で同期処理される.これによって各スレッドで処理されるデータの整合性が保たれる.
スレッドの実際の処理はOSのスケジューリングにまかされる事になり,ユーザのコントロールは効かないが,スレッドの個数はプロセッサコアの数程度が最も効率が良いようである.

3. 4 C++実装

オブジェクト指向言語として知られるC++であるが,もう一つの大きな特徴としてテンプレートがある.これは,変数の型に対して抽象化を行い,プログラムの形 ( アルゴリズム ) だけを定義するための機能である.この機能は利用法が発達しており,自動的にソースコードを生成する事ができ,プログラミングをコンパイラに肩代わりしてもらう事ができる.
以下の節では,Moleculaで用いたC++実装の特徴的な物について紹介する.

3. 4. 1 分子クラス

文献[2]にあるように,分子クラスを定義するのがオブジェクト指向デザイン的であると言える.Moleculaではリジッド分子モデルを採用しているため,力の作用点であるサイトが一つの単原子とサイトが複数ある複数原子分子という分類と,Lenard-Jones ( LJ ) 力やCoulomb力などの組み合わせによる分類を組み合わせたクラスを定義した.
例えば,LJ力と単原子の組み合わせであるAtomクラス,それにCoulomb力が加わったC_Atomクラス,複数原子分子のLJ力のみの場合はMoleculeクラス,などである.例えば水はSPC/Eモデルの時はC_Moleculeクラスである.また個別の特性を持つもの,例えばNa+モデルを考える時,イオン同士ではBMH-Tosi-Fumi ( BMHTF ) ポテンシャル,その他ではLJ,Coulomb力の様にしたいときには,C_Atomクラスを親としてNa+特有のクラスを定義することで対応できる.
こうする事によって,共通の力で定義できる原子・分子は,共通のクラスで定義される異なるパラメータを持つ別個のデータとして,また特殊な相互作用をする原子/分子の場合は新たなクラスとして定義する.新たな原子や分子の計算ができるようにするには、単にLJ,Coulomb力のパラメータが異なる原子や分子であればプログラミングレベルでの変更は不要である.また新たなポテンシャルモデルの原子分子を追加する場合は、その記述および後述するTYPELISTへの登録が必要となる。
ここで,原子/分子クラスの実装の中には,ポテンシャルのパラメータ,例えばLJポテンシャルのσやε,は保持されているが,力の計算そのものは含まれていない.力の計算部分は,次節で説明するダブルディスパッチというメカニズムにより,二つの原子/分子クラスの組み合わせ対して定義された相互作用が,その原子/分子クラスの型情報からコンパイラによって自動的に割り振られる.そのため,手続き型言語では必要になる原子/分子の組み合わせによる条件分岐をプログラミングする必要が無くなるのである.

3. 4. 2 ダブルディスパッチ:2体間力の実装

分子動力学において基本となる力は2体間力であろう.しかしながらこの力は,原子/分子の組み合わせで異なるのが一般的である.
このような組み合わせで異なる力は,科学技術計算分野でよく用いられるプログラム言語による実装では,条件分岐によっているのがほとんどだと思われるが,このような場合C++では仮想関数を用いたダブルディスパッチ法がよく用いられる[5].この方法では,2体間力は原子/分子クラスのメンバ関数として定義される.しかし,この方法だと,新たに原子/分子クラスを追加すると,全ての既存の原子/分子クラスに新たに関数を追加する必要が出てくるため,変更コストが高い.
例えば,Na+,Cl-,水の3種類のクラスが登録されている状態で新たな原子Xを登録するとする.相互作用は各々のクラスのメンバ関数Interactで行う.ユーザが行うのはXとの相互作用を計算するメンバ関数Interactの各クラスへの追加である.新規クラスでは既存クラスとの相互作用の計算のためのメンバ関数の定義が必要となる.
    Class X {
                    ...
                    void Interact( X& a1 ) {...}
                    void Interact( Water& a1 ) {...}
                    void Interact( Na& a1 ) {...}
                    void Interact( Cl& a1 ) {...}
    };
また,既存クラスでは,
    Class Na {
                    ...
                    void Interact( X& a1 ) {
                                            // Xとの相互作用の計算
                    }
    };
このように,全ての既存のクラスにメンバ関数Interact ( X& ) を追加しなければならない.
また手続き型言語では,相互作用を計算する関数内部で,その2つの原子/分子の組み合わせに対して条件分岐を行うことになり,多くの種類の原子/分子を扱う時,プログラムが長く見通しが悪くなり変更コストが高くなる.
そこで,Loki[6]と呼ばれるC++ライブラリで実装されているダブルディスパッチ処理用テンプレートクラスStaticDispatcherを用いて,この分子間力の部分の実装を,原子/分子クラスから分離する事で変更コストを下げる事にした.
実行コードにおいては条件分岐を行っているのだが,その条件分岐の部分をコンパイラが自動的に作成してくれるため,手作業の部分を大幅に減らせる事になる.
このメカニズムは,Lokiの二つのTemplateクラス,LOKI_TYPELIST_nとLoki::StaticDispatcherを用いて実装されている.宣言の部分は,
    typedef Loki::StaticDispatcher<
    Interactor,
    Element,
    LOKI_TYPELIST_3(Na,Cl,C_Molecule),
    false,
    Element,
    LOKI_TYPELIST_3(Na,Cl,C_Molecule),
    void
    Dispatcher;
となり,TYPELISTに登録されているクラス間で相互作用が計算される.
相互作用部分は,Interactorクラスのメンバ関数Fireで以下のように定義される.
    inline void Fire( Atom1& a1, Atom2& a2 ) {
            ...
            LJ::calc_force( ... );
            ...
            BMHTF::calc_force( ... );
            ...
            Coulomb::calc_force( ... );
            ...
    }
このFireは全ての原子/分子クラスの組み合わせに対して実装されなければならない.いくつかのクラスに対しては,Templateを用いる事で実装を省略する事が出来る.例えばAtomクラスにおいてはLJ力だけが相互作用なので,以下のように出来る.
    template<class ATOM>
    inline void Fire( Atom& a1, ATOM& a2 ) {
            ...
            LJ::calc_force( ... );
            ...
    }
ATOMにはコンパイラが自動的に必要とされるクラスを代入してソースを補完してくれる.
 このように,プログラムソースから,従来のプログラムシステムでは必要だった原子/分子の種類による条件分岐を省略でき,いろいろなポテンシャルが簡単に組み込めるところに本プログラムシステムの優位性がある.

4 実装されたアルゴリズムとポテンシャル

本章では,Moleculaに実装された数値計算のためのアルゴリズムと,塩の結晶と水の相互作用のために実装されたポテンシャルについて説明する.

4. 1 並進運動

原子/分子の並進運動は,その重心についての運動方程式を,Leap-Frog法により離散化した.
また,適当なステップ数で系全体の並進運動を除去し,系全体の平行移動による,実質的な温度の低下を防いでいる.

4. 2 回転運動

Moleculaでは分子構造が固定されたリジッドモデルを採用しているため,その剛体回転運動を解く必要がある.ここでは回転運動を示す支配方程式に四元数表記を用いた[7, 8].アルゴリズムはFinchamのExplicit predictor-correctorスキームを用いた[8].

4. 3 ポテンシャル

先にあげたように2体間力の実装は容易に行える.ここでは,基本的なLJポテンシャルと,CoulombポテンシャルとBMHTFポテンシャルを実装した.
水モデルはSPC/Eモデルとし,イオンとの相互作用はLJ+Coulomb力とした.ナトリウムイオンと塩素イオンの相互作用はBMHTFモデルとした[9].各力モデルのパラメータは,文献[9]のTable I, IIのデータをそのまま用いた.
2体間力の実装はInteractorクラスで行うが,今回の場合,水,Na+,Cl-の3種類のクラスのうち2つの組み合わせとして,全部で9種類のFire関数を実装した.この組み合わせは可換なので6種類まで減らせる.

4. 4 温度制御

温度制御は,NVEアンサンブルのための温度制御を行わない物と,NVTアンサンブルのためBerendsen制御といくつかの単純な温度スケーリングによる制御を実装した.

5 塩の水への溶解

Moleculaを用いた,塩の水への溶解の様子のシミュレーションを紹介する.
初期条件としては,塩の結晶を中心に置き,その回りに水分子を格子状に配置して囲んだ ( Figure 2 ) .塩の結晶は,Na,Clそれぞれ32個ずつ,水分子は32個の固まりが4つで計128個である.水の格子間隔は0.3 nm,塩の格子間隔は0.291 nmである.温度はBerendsen制御で,300 Kとした.


Figure 2. Initial distribution of the crystal of NaCl ( Na:big, Cl:small sphere ) and water ( surrounding NaCl ) .

水の中にNaイオン1つが溶け出したところをFigure 3に示す.水は計算上では存在するが,塩の様子を見るために,表示を不可視にしている.


Figure 3. A Na ion is dissolving into water. Water molecules are invisible but exist surrounding the crystal.

なお,このシミュレーションの動画は動画サイトYoutube[10]で見る事が出来る.

6 まとめ

科学技術計算用プラットフォームとして,現在の高性能PCは,ミニチュア版スーパーコンピュータとも言うべき性能を持っている.また,三次元可視化のためのグラフィックスもハード/ソフトともに充分に成熟している.
これらの機能を踏まえ,三次元リアルタイム可視化機能を備え,C++テンプレート機能を利用して,分子間相互作用の実装コストを低減させた,マルチプラットフォームで動作する,Molecula Numericaを開発した.
本プログラムは、MDの初学者ならびに,ポテンシャルを新たに作成して新たな系のMDを実行されたいとお考えの研究者の利用を念頭に置いている.
このソフトウェアのバイナリはMacOS X 10.4用と,Windows XP用がVector Inc.[11]などのダウンロードサイトから無償でダウンロードできる. またソースコードも非営利利用に対し,条件付きで無償提供される.興味のある方は筆者まで連絡願いたい.

参考文献

[ 1] Hiroshi Abe, Comput. Accel. Phys. Conf., 1996, 41.
[ 2] 佐々木徹,長嶋雲兵,塚田捷, J. Chem. Software, 6, 165 (2000).
[ 3] Cross-Platform GUI library :
http://www.wxwidgets.org/
[ 4] Pthread for Win32:
http://sourceware.org/pthreads-win32/
[ 5] Scot Meyer, More Effective C++ (1998), 218.
[ 6] Andrei Alexandrescu, Modern C++ Design (2001), 277.
http://loki-lib.sourceforge.net/
[ 7] 上田顕, 分子シミュレーション (2003).
[ 8] David Fincham, CCP5 Newsletter, 2, 6 (1981).
[ 9] E. Sanz and C. Vega, J. Chem. Phys., 126, 014507 (2007).
[10] http://www.youtube.com/watch?v=NeDmsU3pWFU
[11] http://rd.vector.co.jp/vpack/browse/person/an015761.html


Return