PLAINセンターニュース第97号
page 1

数値シミュレーションエンジン

(2) 専用エンジンによる数値シミュレーションの方法

薄田竜太郎
理化学研究所 計算科学技術推進室

図1 専用エンジンを用いたシステムの構成

 図1のように専用エンジンはPCやワークステーションであるホストと組み合わせて使用します。専用エンジンのプロセッサーを搭載したボードをホストの外部バスに接続します。この接続方法は例えばPCやワークステーションの外部バスにネットワークカードを接続するのと同様です。重力多体シミュレーション専用エンジンを使うときはホストは各質点の位置と質量を専用エンジンのメモリに書き込みます。これをもとに専用エンジンは各質点に働く重力を計算し、ホストがその計算結果を読みとります。分子シミュレーション専用エンジンを使うときはホストは各粒子の位置、電荷、種類を専用エンジンのメモリに書き込みます。これをもとに専用エンジンは各粒子に働くクーロン力や分子間力を計算し、ホストがその計算結果を読みとります。これらの力の計算以外のシミュレーションの計算はホストが行います。シミュレーションプログラムの開発と実行は汎用コンピューターの場合と同じようにホスト上で行います。CやFortranのプログラムから専用エンジンにアクセスするためのライブラリがあり、シミュレーションプログラムはこれを使って専用エンジンにアクセスします。専用エンジンが接続されていないコンピューターでもプログラム開発ができるように専用エンジンのエミュレーションライブラリがある場合もあります。

 具体例として、分子シミュレーション専用エンジンMDGRAPE-2を使って粒子間のクーロン力を計算する場合を説明します。クーロン力を重力に、電荷を質量に読みかえると質点間の重力の計算になります。MDGRAPE-2ボードはPCIカードですので、PCやワークステーションのPCI拡張スロットに接続します。N粒子に働く、自分以外のN-1粒子からのクーロン力をMDGRAPE-2で計算するCプログラムは次のようになります。

#include "m2_unit.h"

double r[N][3]; /* 各粒子の位置のx、y、z成分 */
double q[N]; /* 各粒子の電荷 */
double f[N][3];/* 各粒子に働くクーロン力のx、y、z成分 */
M2_UNIT *mu; /* MDGRAPE-2へのポインタ */
int i;
... /* rとqに各粒子の位置と電荷を代入する */
mu = m2_allocate_unit("grav.table", M2_FORCE, 0, 100, NULL_INT);
/* 1辺100の立方体の中で粒子に働くクーロン力を計算するようにMDGRAPE-2を初期化する */
m2_set_positions(mu, r, N); /* MDGRAPE-2に位置を書き込む */
m2_set_charges(mu, q, N); /* MDGRAPE-2に電荷を書き込む */
m2_calculate_forces(mu, r, N, f); /* MDGRAPE-2で電場を計算し、fに返す */
for(i = 0; i < N; i ++){
int k;
for(k = 0; k < 3; k ++){
f[i][k] *= q[i];
 }
}

 m2_から始まる関数はMDGRAPE-2にアクセスするライブラリの関数です。MDGRAPE-2は正確には各粒子の位置における電場を計算しますので、最後に各粒子の電荷を掛けてクーロン力を求めています。実際の分子シミュレーションではこれ以外にいろいろな計算が必要ですが、これらはホストで行いますので、その部分は汎用コンピューター上のプログラムと同様になります。このプログラムをホスト上でコンパイルしてMDGRAPE-2にアクセスするライブラリをリンクし、実行します。

 ホストにMDGRAPE-2ボードを2枚以上接続して使う場合もプログラムは同じです。汎用CPUでは2つ以上のCPUで並列に計算する場合は一般にプログラムを変更する必要があります。しかしMDGRAPE-2ではボードを2枚以上使う場合はライブラリの内部でボードを並列に動作させますので、プログラムを変更する必要はありません。

 MDGRAPE-2を使った分子シミュレーションの例を図2に示します。これはNaClの液体が冷えて結晶化する様子をシミュレートしたものです。最初、液体のときはNaイオンとClイオンが不規則に混ざり合っていますが、融点以下に冷やすと結晶化して規則正しく並ぶ様子がわかります。イオン間にはクーロン力以外に分子間力も働きますので、MDGRAPE-2で両方の力を計算しています。

 クーロン力を直接計算せず、ほかのアルゴリズムを使う場合には、アルゴリズムの中で粒子間の中心力として表される部分を専用エンジンで行い、それ以外の部分はホストで行います。先の例ではクーロン力を直接計算していますが、周期境界条件の場合や粒子数が多い場合には直接計算は計算量が多いので、計算量を減らすためにいろいろなアルゴリズムがあります。周期境界条件のときにクーロン力の計算量を減らすためにEwald法というアルゴリズムがあります。Ewald法はクーロン力を実空間と波数空間の2つの部分に分けて計算します。実空間部分はある中心力の形に表されます。分子シミュレーション専用エンジンは任意の中心力を計算できますので、実空間部分は専用エンジンで計算します。プログラムはクーロン力の場合と似たものになります。波数空間部分はホストで計算します。波数空間部分を高速に計算する機能をもった専用エンジンもあり、その場合は波数空間部分も専用エンジンで計算します。粒子数が多い場合にクーロン力の計算量を減らすために高速多重極子法などのアルゴリズムがあります。専用エンジンを使う場合にはこれらのアルゴリズムを専用エンジン向きに一部変更したものを用います。これらのアルゴリズムの中には、直接計算するより計算量は少ないですが粒子や疑似粒子間のクーロン力の計算をする部分があります。この部分を専用エンジンで計算します。アルゴリズムの中の他の部分はホストで計算します。現在の専用エンジンの性能では粒子数がおおむね数万よりも多い場合は直接計算するよりも高速アルゴリズムを用いた方が速くなります。汎用CPUの場合と同様に高速アルゴリズムを使うと直接計算よりもプログラムは複雑になります。

 よく使われるシミュレーション用のソフトウェアパッケージを専用エンジンを使ったシステムの上で動くように移植したものもあります。この場合専用エンジンはパッケージのバックエンドとして働きます。利用者は、汎用コンピューター上のパッケージと同じインターフェースで使うことができます。

(左)初期状態の液体状態

(右)融点以下に冷却した後の結晶 イオン数13824

図2 : 分子シミュレーション専用エンジンMDGRAPE-2を使ったNaClの結晶化の
シミュレーション(古石貴裕氏等による) 赤 Naイオン 緑 Clイオン



この号の目次へ

大計算機運営委員会・議事録(抄)へ
0
Mail to CPISPLAINnews HOME
この号の目次へ


(145kb/4ページ)

第98号へ
第96号へ
バックナンバー一覧へ
著者リスト一覧へ