概要

本稿では，分子シミュレーションでの利用を想定し，計算効率が良い新しい疑似乱数生成器（Pseudo Random Number Generator: PRNG）を提案する．本手法による計算速度の向上を示すために，散逸粒子動力学（Dissipative Particle Dynamics: DPD）シミュレーション向けに既存の手法と比較を行う．DPD法は，分子の熱ゆらぎを乱数によって再現することで，メソスケールの物理現象を効率的に取り扱う事ができる．しかし，並列計算環境で運動量を保存しながら乱数を用いるには，暗号化などの追加計算が不可欠であった．本研究では，DPD法の計算工程そのものに高品質な乱数が内在していることを示し，それを活用することで暗号化を要しない効率的なPRNGを提案する．本手法は複雑な乱数生成アルゴリズムや暗号処理に依存しないため，最小限の計算コストで乱数を生成可能である．本手法は他の多くの物理・化学シミュレーションでも，高速な乱数生成を可能にする可能性がある．

1　 はじめに

様々な手法が存在するシミュレーションの中でも，近年特に，計算機の発達を受けてますます利用価値が高まっている手法の一つが，分子シミュレーションである．実験では空間的，時間的に観察することのできない分子スケールの現象解明によって，新しいバイオマテリアルのデザインなどを始めとし，分子生物化学，生態材料科学などに役立てられている．スーパーコンピュータを用いての計算は，例えば，一般的に難しいとされてきた，全原子シミュレーションでの細胞小器官全体の再現や，100 μsを超える長時間のタンパク質シミュレーションなど，大規模な系の計算を可能にしており，ますます多岐にわたる応用的な利用が議論されている［1,2］．分子シミュレーションの一つである散逸粒子動力学法（Dissipative Particle Dynamics: DPD）は，全原子分子動力学法と違い，数個の粒子を1つの粗視化粒子としてみなし，少ない計算量で比較的大きな現象を取り扱う手法である［3–5］．一般的な粗視化分子動力学法では，粗視化に伴って粒子の自由度が減少し，分子内摩擦を扱えていないために，正確な動的性質を再現できないという問題が指摘されてきた．一方で，DPD法では，粒子間力を算出する際に，保存力項に加え，散逸力項とランダム力項を加えることで，粗視化に伴い再現できなくなる分子の摩擦と揺動の効果を再現し，比較的小さな計算コストで数百nm，数十μsほどのスケールの動的性質の解析が可能となった．しかし，スーパーコンピュータ等を用いてDPD法による計算の大規模並列化を行う際，生じる問題がある．分子動力学法など，力が保存力のみで構成される場合は，粒子間距離から力を計算することができるため，効率的な並列計算が可能であったが，DPD法のように，力が保存力のみならず，乱数を用いて計算されるランダム項を含む場合は，各計算ノードごとにそれぞれ異なった疑似乱数生成を行うため，各粒子ペアに共通したランダム力を再現することが難しく，高い効率の並列計算を実行することは困難である［6–10］．そこで，各計算ノードが共有する情報を，暗号学的ハッシュ関数で変換することにより，余分なノード間通信を伴わず，粒子ペアごとに共通した乱数を生成する方法が考案されている［11］．しかし，暗号化を用いたPRNGはMersenne Twister（MT）などの一般的なPRNGに比べて高い計算コストを必要とする．本研究では，シミュレーションの計算そのものに内在する乱雑性を利用し，暗号化を用いない高速なPRNGを提案する．

2　 手法

2.1　 散逸粒子動力学法

DPD法は，複数の分子の集合を1つの粒子とみなして計算を行う．その粗視化粒子の挙動は次の運動方程式

  

$$\begin{array}{c}\hfill \frac{d{\mathit{r}}_i}{dt}={\mathit{v}}_i\end{array}$$(1)

  

$$\begin{array}{c}\hfill m_i\frac{d{\mathit{v}}_i}{dt}={\mathit{f}}_i\end{array}$$(2)

に従う．$m_i,{\mathit{r}}_i,{\mathit{v}}_i,{\mathit{f}}_i$はそれぞれ粒子$i$の質量，位置，速度，作用する力を示す．一般的に全ての粒子の質量は等しいとされ，$m_i=m\equiv 1$と定められる．ここで${\mathit{f}}_i$は次に示す3つの2体間の力の総和として

  

$$\begin{array}{c}\hfill m\frac{d{\mathit{v}}_i}{dt}={\mathit{f}}_i=\sum _{j\ne i}^{}({\mathit{F}}_{ij}^{\mathrm{C}}+{\mathit{F}}_{ij}^{\mathrm{D}}+{\mathit{F}}_{ij}^{\mathrm{R}})\end{array}$$(3)

と定められる．${\mathit{F}}_{ij}^{\mathrm{C}},{\mathit{F}}_{ij}^{\mathrm{D}},{\mathit{F}}_{ij}^{\mathrm{R}}$はそれぞれ粒子$j$から$i$に作用する保存力，散逸力，ランダム力を表す．保存力は，粗視化粒子の相互作用の力で用いられているものと同様，ソフトポテンシャルに起因する斥力として，

  

$$\begin{array}{c}{\mathit{F}}_{ij}^{\mathrm{C}}=\{\begin{array}{cc}-a_{ij}(1-\frac{r_{ij}}{r_{\mathrm{c}}}){\mathit{e}}_{ij}\hfill & (r_{ij}\le r_{\mathrm{c}})\hfill \\ 0\hfill & (r_{ij}>r_{\mathrm{c}})\hfill \end{array}\hfill \end{array}$$(4)

のように提案されている．$r_{\mathrm{c}}$はカットオフ半径を表しており，DPD法では長さの単位とし，$r_{\mathrm{c}}\equiv 1$とする．また，$a_{ij}$は粒子$i$，$j$間の最大斥力，そして${\mathit{r}}_{ij}={\mathit{r}}_i-{\mathit{r}}_j$，$r_{ij}=|{\mathit{r}}_{ij}|$，${\mathit{e}}_{ij}={\mathit{r}}_{ij}/r_{ij}$とする．散逸力と，ランダム力は，$\gamma$を摩擦係数，$\sigma$を揺らぎの大きさ，${\omega }^{\mathrm{D}}(r_{ij})$，${\omega }^{\mathrm{R}}(r_{ij})$を重み関数，また，${\mathit{v}}_{ij}={\mathit{v}}_i-{\mathit{v}}_j$として，保存力と同様$r_{ij}\ge r_{\mathrm{c}}=1$で0になるという条件のもと，

  

$$\begin{array}{c}\hfill {\mathit{F}}_{ij}^{\mathrm{D}}=-\gamma {\omega }^{\mathrm{D}}\left(r_{ij}\right)({\mathit{e}}_{ij}\cdot {\mathit{v}}_{ij}){\mathit{e}}_{ij}\end{array}$$(5)

  

$$\begin{array}{c}\hfill {\mathit{F}}_{ij}^{\mathrm{R}}=\sigma {\omega }^{\mathrm{R}}\left(r_{ij}\right){\theta }_{ij}{\mathit{e}}_{ij}\end{array}$$(6)

と表される．ここで${\theta }_{ij}$は以下の特性を有するガウス性白色雑音である．

  

$$\begin{array}{ccc}\hfill {\theta }_{ij}(t)& =& {\theta }_{ji}(t)\hfill \\ \hfill \mathcal{⟨}{\theta }_{ij}(t)\mathcal{⟩}& =& 0\hfill \\ \hfill \mathcal{⟨}{\theta }_{ij}(t){\theta }_{kl}(t^{\prime })\mathcal{⟩}& =& ({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})\delta (t-t^{\prime })\hfill \end{array}$$(7)

ここで${\delta }_{ij}$はKroneckerのデルタであり，$i=j$のとき${\delta }_{ij}=1$となり，それ以外の場合は${\delta }_{ij}=0$となる．これは各粒子ペアごとに${\theta }_{ij}={\theta }_{ji}$を満たすので，運動量保存が成立する．さらに，Eq. (5)，(6)中の${\omega }_{\mathrm{D}}(r_{ij}),{\omega }_{\mathrm{R}}(r_{ij})$は

  

$$\begin{array}{c}\hfill {\omega }^{\mathrm{D}}(r_{ij})=[{\omega }^{\mathrm{R}}(r_{ij}){]}^2=\{\begin{array}{cc}\hfill {(1-\frac{r_{ij}}{r_{\mathrm{c}}})}^2\hfill & \hfill (r_{ij}\le r_{\mathrm{c}})\hfill \\ 0\hfill & \hfill (r_{ij}>r_{\mathrm{c}})\hfill \end{array}\end{array}$$(8)

の関係式を満たす．

2.2　 PRNGの開発

DPD法の計算に含まれる乱雑性をNIST SP 800-22に準拠する乱数検定（以降では単にNIST乱数検定と呼ぶ）を用いて調べる．その後，粒子座標のランダムな性質を利用することで，従来の方法のように暗号化に計算コストを割くことなく，高速に乱数を生成できることを示す．この新手法の開発にあたり，以下の手順で検証した．(i) MT による乱数を用いて実行した DPDシミュレーションの粒子座標の乱雑性を NIST乱数検定を用いて各タイムステップ，粒子毎に評価した．(ii) NIST乱数検定の結果をもとに，粒子座標を元に乱数を生成し，DPD シミュレーションを実行した．(iii) DPD法では一般に一定の時間ステップ幅が用いられるが，より短い時間ステップ幅をテストし，新手法を用いた場合に特定の粒子番号に対応する座標から抽出された乱数が，連続する時間方向に品質の良い乱数を生成できているかどうかを調べた．

3　 結果

3.1　 粒子座標の乱雑性

粒子座標から取り出した数値の乱雑性を評価するために，NISTの乱数検定を用いた．これは15種類の異なる検定方法で構成されており，それぞれが1,048,576ビット単位の系列に対して行われる．さらに，この検定を1,000系列に対して実施した結果を総合的に判断する．乱数系列に対する統計的な検定のため，検定結果は各テストやサブテストごとにp値が示される．本研究では有意水準を0.01とした．ここで用いたDPDシミュレーションは倍精度の計算精度で実行した．

図1に示すように，粒子の位置を表現する浮動小数点数は，符号部（1ビット），指数部（11ビット），仮数部（52ビット）の合計64ビットで表現される．各タイムステップ間で粒子座標が変化すると，64ビットのうちどのビットも更新されるが，符号部と指数部は連続するタイムステップ間で強い相関を持つ．一方，仮数部の一部は数値が大きく変化し，前のタイムステップとの比較では純粋な乱数のように振る舞う．この乱雑性をより詳しく検証するため，64ビットの座標値から特定の32ビット（以下，ビットストリームと呼ぶ）を抽出し，それを前のタイムステップでの同じ部分と比較した．さらに，抽出するビットストリームの位置を1ビットずつ上にずらしていく操作を繰り返した．

図1  64ビット浮動小数点数で表される1粒子の粒子座標（x成分）から32ビットのビットストリームを抽出する模式図．各タイムステップにおいて，特定の粒子，座標，特定のビットストリームを時系列ごとに結合し，NISTの乱数検定セットを用いて乱数性を検証する．

ここで強調すべきは，すべて同じ粒子の座標を連続するタイムステップで評価している点である．通常，タイムステップが小さい場合，あるタイムステップから次のタイムステップにかけて粒子座標は強い相関を持つと考えられる．

図2は，このように抽出したビットストリームをNISTの乱数検定で評価した結果を示す．横軸は粒子座標の対応するビットストリームを，縦軸にNIST乱数検定の点数を15満点で示した．すでに述べたように，符号部や指数部は強い相関を持つと期待される．仮数部のうち下位側を取り出した場合は隣り合うタイムステップでも相関が無く，乱数性が高い場合が多い．しかし，仮数部の最下位ビット付近では桁上がり処理の影響を受けるため，やや乱数性が悪化する傾向にある．結果として，タイムステップ間の相関が無く乱雑性の高いビットストリームを得るには，仮数部の中位程度のビットを取り出すのがもっとも有効であることがわかる．

図2  連続するタイムステップで得られる64ビット浮動小数点数の粒子座標から，それぞれ32ビットのビットストリームを抽出してNISTの乱数検定を適用した結果．横軸は抽出した32ビットのビットストリームを選ぶ始点を示す．たとえば “First bit = 1” は仮数部の最下位から32ビットを取り出したことを意味する．縦軸はNISTの各検定において合格したテスト数（合計15）を示す．

3.2　 新たなPRNGの開発

図2の結果から，浮動小数点数の仮数部における特定範囲が前のタイムステップとの相関を抑制し，乱数のシードとして適切であることが確認された．ここでは，この知見を並列DPD計算へ応用する手法を示す．並列実行されるDPDシミュレーションでは，どの計算ノードにおいても，同一の粒子対に対して同一の乱数を用いる必要がある．従来の手法では，計算ノード間で共有される2つのシード値（タイムステップ数や粒子番号など）を用いて同じ乱数を生成する方式が一般的であった．本研究では，2つの異なる粒子の$x$座標をシードとして用いる方法を提案する．得られた2つのシード値を演算によって1つの64ビット数値に統合することで，各計算ノードにおいて同一の粒子対に対して同一の値を再現できるうえ，ノード間で追加の情報を交換する必要がなくなる．このとき，どの四則演算が連続するタイムステップ間の相関を小さくできるかを検討するために，加算および乗算を用いた場合についてNISTの乱数検定で比較を行った（なお，減算は$a-b$と$b-a$で結果が異なるため検討から除外している）．得られた乱数を前のタイムステップのものと比較し，NISTの乱数検定で評価した結果を図3に示す．乗算を用いた場合は，すべてのビットストリームにおいて良好な乱数性を示す一方，加算を用いた場合はやや性能が劣る傾向が見られた．これは，乗算によって仮数部の下位に由来する乱雑性が上位ビットにまで伝播しやすくなるためである．このように，単純な乗算演算を組み込むだけで，粒子座標をシードとする乱数生成に十分な乱雑性を得ることが可能である．したがって，計算コストの高い暗号アルゴリズムを導入する必要はないと考えられる．

図3  (a) 加算と (b) 乗算を用いて生成した乱数列をNISTの乱数検定で評価した結果を示す．

提案する$[0,1]$区間に乱数を生成する手法を，擬似コードとして表1に示す．本手法をまとめると以下の6つのステップで乱数生成を行う．

表1 提案するPRNGのソースコード例

Example of source code
union Uni
double PLAIN1
unsigned long long PLAIN2
end union
Uni uni
uni.PLAIN1 = position[i] × position[j]
uni.PLAIN2 = uni.PLAIN2 << 32
double rand = (double)uni.PLAIN2
rand = rand / ULLONG_MAX

(1) 2つの粒子のx座標を乗算し，乱数のシードを生成する．(2) 生成されたシードをunionを使って整数型として扱う．(3) 整数型として扱ったビット列を32ビット分シフトする．(4) シフトしたビット列をunsigned long long型の整数値として読み出す．(5) 読み出した整数値をdouble型にキャストする．(6) unsigned long long型の最大値で割り，$[0,1]$の範囲に正規化する．

初期条件を設定する際には注意が必要である．たとえば，粒子を完全に格子状に並べ，初期速度を0とすると，隣接粒子同士の距離が完全に同一になる場合がある．このときはいかに高度な暗号アルゴリズムを用いても，シードが同じになるため，生成される乱数も同じになり，ランダム力の計算が破綻する恐れがある．したがって，結晶構造などの規則的な初期配置が必要である場合，初期速度をMT法などに由来する乱数を用いて設定するか，初期座標をわずかにランダムにずらすなどの対策を講じる必要がある．なお，これは粒子座標をシードにする他の暗号アルゴリズムを用いる場合でも同様である．

3.3　 生成した乱数の品質評価

次に，提案するPRNGが汎用的に利用できるかどうかを検証する．DPD法ではさまざまな物質を取り扱うが，ポリマーなどのように大きな粒子が結合した系では，連続するタイムステップ間で粒子座標が大きく変化しない場合がある．そこで，座標の変化量が小さい場合に本手法で生成する乱数性質が変化するかどうかをNIST検定を用いて調べた．

具体的には，通常より小さいタイムステップを用いたシミュレーションを行い，系の座標変化量が減った場合の乱数をNISTの乱数検定で評価した．タイムステップを極端に小さくすると，粒子座標の更新もわずかとなり，乱雑性が低下すると考えられる．

水系を対象とし，GrootとWarrenの研究［5］に従ってDPD法でシミュレーションを行った．標準的な無次元タイムステップである0.04を基準とし，さらにこのタイムステップを1/10倍，1/100倍，1/1000倍…最大1/10¹⁴倍として計算を実行した．そして，連続するタイムステップで生成される乱数をNISTの乱数検定セットで評価した．結果を図4に示す．通常のタイムステップより1/10000にまで小さくした場合でも，高品質な乱数が生成できることが確認できた．つまり，液体水分子の1/10000の速度で動く，個体のような物質を含む系であっても，提案手法を利用できることを示す．

図4  異なるタイムステップ（$dt$）で実行した水のDPDシミュレーションにおいて，生成した乱数をNISTの乱数検定で評価した結果．特定の粒子ペアのみから時系列順に得られる乱数を結合し，評価した．通常のタイムステップの1/10000に相当する大きさでも品質の良い乱数が得られることが確認できる．

3.4　 計算コスト

計算効率を評価するため，本手法と従来手法であるTiny Encryption Algorithm（TEA）法を比較した．計算系として，水108,000粒子系のシミュレーションを並列実行し，その中で一様分布[0,1]に従う乱数を100億個生成した．シミュレーションの実行では，並列効率を向上させるため，FDPS（Framework for Developing Particle Simulator）［12］を用い（パラメータnleaflimit = 40, ngrouplimit = 64），Intel(R) Xeon(R) CPU E5-2680 v2 @ 2.80GHzおよびIntelコンパイラver.18.0.0で実行した．

100億個の乱数生成に必要な計算時間について，提案手法はTEA法に比べ，約62倍高速に生成できることを確認した．次に，粒子間相互作用計算（粒子間距離の計算を含む）に要する時間を比較すると，提案手法を用いることで従来のTEA法より約27%計算時間を短縮できることが確認できた．なお，本研究ではシミュレーションコードの最適化を行っていないため，コードの最適化や，粒子間距離の計算対象を削減するような手法を導入すれば，計算コストをより一層削減できることが期待できる．

4　 結論

提案手法は，計算工程に内在する乱雑性を抽出することにより，効率的に乱数を生成する．本手法はDPD法のみに限定されるものではなく，計算工程で導出される物理量に乱数性質の良い成分を含むシミュレーションであれば，特に並列計算環境でシミュレーションを行う場合に，有効である．また，ブラウン動力学，モンテカルロ（MC）法など，様々な数値シミュレーションに拡張することが期待される．さらに，本手法で生成される乱数の周期は，シミュレーションで取り扱う系のサイズにスケールする．本研究で提案する乱数生成法は粒子座標を使用するため，乱数の周期は粒子座標と速度が同一の状態間の間隔として見積もることが可能である．この周期は粒子数とともに増加するが，これは系全体の可能な状態の数が増加するためである．倍精度（64ビット），3次元，$N$粒子を仮定し，粒子座標と粒子速度の両方を考慮すると，粒子系が取り得る状態の総数は$2^{64\times 3\times N\times 2}$となる．系は確率的に遷移するので，各状態がサンプリングされる確率は等しくなる．したがって，粒子数$N$が増加するにつれて，粒子座標と粒子速度が2回目に一致する確率は減少し，周期が長くなると推定できる．粒子座標を乱数シードとして用いる場合，乱数を暗号化アルゴリズムを利用して生成する場合も同様である．したがって，非常に大規模な系に対して高いスケーラビリティを持ち，将来的に計算機資源が増大しても本技術は適用可能であると考えられる．

5　 おわりに

本研究では，DPD法の並列計算を高速化するためのPRNGとしてSHIFT法を中心に扱ってきたが，博士論文ではさらに適用範囲を広げ，より一般的な分子シミュレーション手法であるMC法に応用できる「3XORs法」を提案している．MC法では，DPD法と比べて各タイムステップで必要となる乱数の総数が多いうえ，粒子座標や粒子速度などをシードとして利用できる乱雑性が限られているため，新たなアプローチとして3XORs法を開発した．

また，DPD法やMC法にとどまらず，幅広いシミュレーション手法に対応可能なデータ駆動型のPRNG（Learned Pseudo Random Number Generator: LPRNG）を提案した．これは，敵対的生成ネットワーク（Generative Adversarial Network: GAN）を用いて乱数生成アルゴリズムの一部を半自動化するものであり，既存の乱数生成法（たとえばMT法）で得られる乱数を教師データとして学習する．Wasserstein距離を用いたGANを導入することで，NISTの乱数検定をクリアする高品質な乱数を生成することに成功した．これまで深層学習を活用したPRNGは数多く提案されてきたが，統計的検定を十分に満たす水準には至っていなかった．本研究で示したLPRNGは，一般的な利用にも耐えうる品質を初めて実現した点で大きな意義があると考えられる．興味のある読者は，そちらも参照されたい．

 謝辞

本研究の実施にあたり，指導教官である慶應義塾大学理工学部の泰岡顕治教授には，多大なるご支援と数多くのご助言を賜りました．ここに厚く御礼申し上げます．また，同大学の荒井規允准教授，山本詠士准教授，そして株式会社Cygamesの倉林修一博士には，有意義な議論や多くの助言をいただきました．深く感謝いたします．さらに，JST博士後期課程学生支援プロジェクトおよび日本学術振興会若手研究者海外挑戦プログラムの財政的支援に厚く御礼申し上げます．

参考文献

• [1]   J. R.  Perilla,  B. C.  Goh,  C. K.  Cassidy,  B.  Liu,  R. C.  Bernardi,  T.  Rudack,  H.  Yu,  Z.  Wu and  K.  Schulten, Curr. Opin. Struc. Biol., 31, 64–74 (2015).
• [2]   H. I.  Ingólfsson,  C.  Arnarez,  X.  Periole and  S. J.  Marrink, J. Cell. Sci., 129(2), 257–268 (2016).
• [3]   P. J.  Hoogerbrugge and  J. M. V. A.  Koelman, Europhys. Lett., 19(3), 155–160 (1992).
• [4]   P.  Español and  P.  Warren, Europhys. Lett., 30(4), 191–196 (1995).
• [5]   R. D.  Groot and  P. B.  Warren, J. Chem. Phys., 107(11), 4423–4435 (1998).
• [6]   B. L.  Holian,  O. E.  Percus,  T. T.  Warnock and  P. A.  Whitlock, Phys. Rev. E, 50(2), 1607–1615 (1994).
• [7]   W. B.  Langdon, Proc. 2008 IEEE Congr. Evol. Comput. (CEC), 459–465 (2008).
• [8]   A.  Zhmurov,  K.  Rybnikov,  Y.  Kholodov and  V.  Barsegov, J. Phys. Chem. B, 115(18), 5278–5288 (2010).
• [9]   S.  Tzeng and  L.-Y.  Wei, ACM Symp. Interactive 3D Graph. Games (2008).
• [10]   F.  Zafar,  M.  Olano and  A.  Curtis, High Perform. Graph., 133–141 (2010).
• [11]   C. L.  Phillips,  J. A.  Anderson and  S. C.  Glotzer, J. Comput. Phys., 230(19), 7191–7201 (2011).
• [12]   M.  Iwasawa,  A.  Tanikawa,  N.  Hosono,  K.  Nitadori,  T.  Muranushi and  J.  Makino, Publ. Astron. Soc. Jpn., 68(4), 54 (2016).

著者紹介

岡田 清志郎（博士（工学））

〔経歴〕2024年3月慶應義塾大学大学院理工学研究科後期博士課程修了，2023年から現所属．〔専門〕分子シミュレーション，疑似乱数生成手法〔趣味〕カメラ．