• トップ
  • G-DEPについて
  • ご購入ガイド
  • サポート
  • お問い合わせ

G-DEPトップ  >  G-DEPの高速演算記  >  高速演算記 第22回 日本ゼオン株式会社 本田 隆 様

高速演算記 第22回 日本ゼオン株式会社 本田 隆 様

高速演算記 第22回は日本ゼオン株式会社 総合開発センター基盤技術研究所 本田 隆様より、GPU を利用した高分子の相分離シミュレーションの高速化-OCTA/SUSHI へのGPU 計算の導入-についてご寄稿頂きました。

GPU を利用した高分子の相分離シミュレーションの高速化-OCTA/SUSHI へのGPU 計算の導入-

1 はじめに

 高分子は種々の化学構造を有するモノマーが重合することにより形成される分子であり、天然に存在し、また、化学的に合成・重合することもできる。1 高分子中のモノマー数は104 を超えるものもある。DNA やタンパク質も高分子の仲間であり、多くのプラスチック材料も高分子である。またモノマーのシーケンスを変化させるだけでも性質が変化する。一般には、このような化学的なモノマー構造やシーケンスを一次構造と呼ぶ。一次構造の特性により、ポリエチレンの結晶化やタンパク質の螺旋やシート構造のような2 次構造を発現し、さらにそれらの凝集やホールディング等により高次構造も形成する。

  よって、高分子材料が発揮するマクロな物性の要因を理解するには、ミクロな一次構造の把握のみならず、その上のメゾ・スケール構造の把握が重要となってくる。また、熱力学的には重心の並進のエントロピーは分子量に反比例して小さくなるので、2 つの異なる種類の高分子を混合しようとしても、異なるモノマー間に働くわずかな相互作用によって相分離してしまうのが一般的である。

 そこで、高分子の材料においては、高次構造や相分離現象を把握・予測するシミュレーション技術が重要となってくる。しかし、時間・空間スケールが変化すると、主要な物理現象の要因も変化する。よって、単一の理論・手法で、全てのスケールのシュミレーションをすることは無理である。例えば、1 モノマー構造のスケール、つまり、原子・分子スケールであれば量子化学が重要で、電子状態を正確に解く必要があるが、さらに大きなスケールでは、分子は仮想的な分子力場中の質点とみなし、古典力学を用いる分子動力学法が実用的なシミュレーション手法となる。さらに大きなスケールでは高分子の絡み合いが重要となるので、フルアトミスティックな描像は捨て、高分子のダイナミックな動きをシミュレーションする粗視化分子動力学が現実的な手法となる。

 高分子溶融体の相分離現象を扱うスケールにおいて用いられるのは、高分子の密度汎関数理論である。この理論は、ミクロ・スケールにおける量子化学計算で1 電子近似により電子を密度場として捉えたように、個々の高分子の認識を止め、同一構造からなる高分子は、同一の密度場による統計的な重みに変換する。この描像により、個々の分子認識にかかる計算量を減らすと共に統計的に最も存在する可能性のある高分子の分布、つまり、相分離構造をシミュレートする。

 我々は、これまで、高分子の密度汎関数理論を用いたシミュレーション・プログラムOCTA/SUSHI(Simulation Utilities for Soft and Hard Interfaces) [1] を開発してきた。しかし、リリース時は1CPU 用のコードであり、リリース後10 年が経過したので東工大のTSUBAME2.0 の力を借り、GPU を利用して並列化した。本稿ではその手法と結果について報告する。
 

2 ガウス鎖近似

 高分子を粗視化するために、いくつかのモノマーの繋がりをセグメントという1つの疑似粒子としてみなすことにする。セグメントからなる高分子の溶融体を考えよう。粗視化を進めセグメントを構成するモノマー数を増やしてゆけば、セグメントが集合するとモノマーは同一のものであるから、個々のセグメントを構成するモノマーが互いに入り込みセグメントは重なる。つまり、粒子とみなしたセグメント間には排除体積効果が働かないような状況に到達する。このように、排除体積効果を無視できるセグメントからなる高分子を理想鎖と呼ぶ。

 理想鎖を構成するセグメントの効果的な長さはb であるとする。セグメントからなる高分子を格子上に配置すると図1(a) にようになる。高分子内の1セグメントは図1(b) のように3 次元格子点の1 つの格子点を占め、このセグメントに繋がったセグメントはその周りの格子点に存在することができる。
 


図1 格子モデル(a) 1 次元モデル, (b) 3 次元モデルによる隣接セグメントの配置[6]
 

  このように配置されるセグメントの統計的重みを考えてみよう。まず、高分子鎖のある特定の部分を抜き出したものを部分鎖と呼ぶことにする。セグメントの接続順により、のセグメントのインデックスを持つ部分鎖の両末端がそれぞれ、位置 に固定されている場合の部分鎖の統計的重みをで表すとする。隣合うセグメントは図1(b) に示す配置となるから、以下の関係式が成り立つ。

(1)




 ここで、 に隣接するサイトの位置、 は配位数を表す。第(1) 式は′ に位置する 番目のセグメントの統計的重みの和になっていることを表している。  番目のセグメントの寄与は確率的に1/ になるのである。この第(1) 式を連続体表記に変更するのに、両辺から を引くと次の式となる。

(2)

 


左辺は での微分として、右辺は の空間的な発散として考えることができるので、3 次元的( = 6) には次の連続の式を考えることができる。 

(3)

 

 
 

  ここで を連続体表記したものである。 がかかるのは拡散が1セルの断面の面積に比例するからである。以降の説明では次元に関係せず = 6 を利用することとする。
 第(3) 式は を時間と考えると の3 次元の拡散方程式の形式となっている。拡散係数は である。よって、初期値を与えれば解くことが可能である。もし、 の位置で末端が固定されていれば、初期条件は次のデルタ関数で与えることができる。 

(4)

 

よって、第(3) 式の解は次のようになる[2] 。

(5)

 

 

右辺はガウス分布になっている。よって、このモデルははガウス統計に従うのでガウス鎖と呼ぶのである。ガウス鎖は理想鎖の仮定より導出された格子の制約のない連続的な記述となっているので、この後で説明するSCF 理論に利用される。
  

3 SCF(Self-Consistent Field)理論

 ガウス鎖近似は、外的影響の無い理想鎖の性質を表しているが、異なる成分からなる高分子溶融体では異なるセグメント間の反発が起き、また、壁等の障害物があれば、非圧縮条件によりガウス鎖は変形するはずである。そのような外的要因を外場として取り入れ、高分子の溶融状態を記述するのが高分子のSCF(Self-Consistent Field)理論である[3, 4, 5, 6] 。

3.1 経路積分

 ガウス鎖の統計を説明するのに、鎖の長さ を時間と見立てたセグメントの統計的な重みが従う拡散方程式である第(3)式を利用した。この式に、セグメント間相互作用と非圧縮条件の影響を外場として取り入れると次の式となる。

(6)

 


 ここで、 番目の部分鎖を構成するセグメントの種類を表し、 はセグメントの比容積、 はセグメントに働く外場である。この方程式はEdwards 方程式と呼ばれる非保存の拡散方程式である。このセグメントの統計的重み を経路積分と呼ぶ。
この方程式は、セグメント間の相互作用以外に非圧縮性を満たすように解かれなければならないので、には2 つの寄与:セグメント-セグメント間相互作用と非圧縮条件が含まれる。
セグメント-セグメント間相互作用としては、次のようなポテンシャルを考える。

(7)

 


 ここで、 種間のセグメント間に働く相互作用エネルギーである。 種のセグメント密度は、 のセグメントからなる全ての部分鎖からの寄与であるので、次の式で求められる。 

(8)

 


 ここで、 は、系中の全ての部分鎖を指すインデックスであるとする。
 さて、あるセグメント分布における外場 が求められたとしよう。その外場は を含むとともに、
非圧縮状態であるためのポテンシャル、逆に言うと、 番目の部分鎖を構成するセグメントの拡散を阻むためのポテンシャルが含まれていなくてはいけいないので、次のように表すことができるはずである。

(9)

 ここで、右辺の第2 項は 番目の部分鎖の化学ポテンシャルに負の符号をつけたものである。つまり、セグメントの拡散を阻むためのポテンシャルである。

3.2 セグメント密度

 外場の式第(9)式は、セグメント間の相互作用を第(8)式により取り入れている。よって、計算にはセグメント密度 を求めなくてくてはならない。そのためには、まず経路積分からすべてのセグメントの寄与を求める必要がある。
 部分鎖が従う経路積分の式である第(6)式は末端を固定して考えている。末端は全ての空間に存在しうる。また、部分鎖は2つの末端があるので、両末端が全ての空間に存在するとして、両末端の位置を積分した2 つの経路積分を次のように定義しよう。

(10)

  

(11)





ここで、 は第(6) 式を満足する末端の初期条件である。上添え字 はポリマーの種類を示し、 番目の部分鎖の鎖長である。

  番目の部分鎖を構成するセグメント密度は、これら積分された経路積分 の積を全ての について積分したものに比例するはずである。つまり、両末端の存在を全ての空間に許せば、部分鎖上の任意のセグメントの密度分布は、このセグメントまでの両末端からの経路積分の積に比例するはずであるから、次の式を考えることができる。

(12)

 


ここで、 は規格化定数である。
 経路積分は、非保存系の拡散方程式であるので、第(12)式の規格化定数は次のように分母にやはり経路積分を含むものでなければならない。 

(13)

 

 

ここで 種部分鎖が含まれる 種鎖の系内での数である。結局分母、つまり、右辺の 種鎖の分配関数である。この分配関数は部分鎖のインデックスに依存しないので、片末端からの経路積分のみを利用し、次のように簡単に計算することができる。
 

(14)

 


ここまで、部分鎖に特化した説明をしてきたが、任意の構造を有するポリマーの経路積分は、経路積分の順番を考え、部分鎖の経路積分が、きちんと部分鎖の接合点に集まるようにすれば計算できる。その概念図を図2 に示す。


図2 トポロジー構造を持つ高分子の経路積分の順序の例[6]



3.3 SCF 計算

 さて、既にお気づきのことかもしれないが、外場にはセグメント間相互作用が含まれているが、セグメント密度は外場を利用した経路積分が計算された後に求められるものである。よって、これまで述べてきた計算は、非線形の計算であり、解析的には解くことができない。よって、繰り返し計算により数値解析的に求めなくてはならない。その概念図を図(3)に示す。

外場→経路積分→セグメント密度→外場 のように循環する繰り返し計算により全ての物理量が求められる。つまり、高分子は自分が形成する外場の中に存在し、その外場が無矛盾となるまで計算する必要が生じるので、この計算を自己無撞着場(Self-Consistent Field)計算と呼ぶのである。

SCF計算はさらに2つに分類できる。これらは、適当な自己無撞着場の初期条件に対して相分離構造を最適化する静的な計算方法(例えば分子力場における構造最適化に相当する)と、適当な初期の相分離構図から化学ポテンシャルを駆動力とした動力学的な相分離をシミュレートする動的な計算方法(例えば分子力場を用いたMDに相当する)である。 SUSHIには両手法を実装してある。
 

4 SCF 計算のGPU 上への実装 

 OCTAの計算エンジン群は全てC++で実装されており、SUSHIもFDM(有限差分法)を用いた C++で実装されたエンジンとなっている。 SUSHIのGPU並列化において、当初OpenCLとCUDAを比較検討した後、CUDAで実装することにした。両手法において同等な実装が可能であるが、検討開始当時のドキュメントの充実度やとっつきやすさが CUDA選択の理由であったと思う(詳細は記憶していない)。 

 


 

 図3 SCF のスキーム[6] 


 CUDAはC++のサブセットとみなせる言語であるので、C++で記述されたSUSHIのコードとは馴染みやすい。また、SUSHIの場合、オブジェクト指向でのプログラミングの傾向が強く、一般的な並列化手法(例えば、最も長いfor文を探して.....)とは少し異なり、あくまでもオブジェクト指向での改良方法である”クラスのメソッドの実装の変更”と”新たなクラスの追加”で CUDAを導入した。

オブジェクト指向が利用できるC++の場合、デバイスを隠蔽したライブラリーとしての実装も簡単であり、それもCUDAの利点となる。


4.1 場の実装の改良

 SUSHIでは場はScalarFieldと呼ばれるクラスとして実装されている。このクラスは演算子を含む種々のメソッドを有しており、SCF計算の実装は、オブジェクトと化したScalarFieldを用いた演算の形式をとっている。また、系の次元にかかわらず、ScalarFieldは場の値を倍精度の 1次元アレイとして保有している。つまり、個々のScalarFieldが大粒度のデータ構造を有しており ScalarFieldのメソッド自体がいわゆる"長いfor文”を形成している。よって、これらの実装を全てCUDAのカーネル関数で置き換えるようにした。例えば

  ScalarField& ScalarField::operator+=( const ScalarField& inpField )

のような演算子の実装をCUDAを利用してGPU上で実装するのである。



4.2 ∇2 の実装

 SUSHIのプロファイリングを行うと経路積分における の計算に時間を費やしていることがわかった。そこで、この実装をCUDAが得意とするスレッド並列化とした。図 4 の中間の図のように領域を分割する。緑の柱がCUDAの一ブロック担当領域に相当する。この柱の数はシステムサイズにより変化する。 柱の断面は1断面が8×8のメッシュサイズである。右図の青・黄・ピンクの3つの層は柱の断面の小領域を意味し、3Dなので3層必要となる。具体的には、29メッシュ点(3の3乗)上のScalarFieldの値を利用して中心位置の∇^2の値を計算する。1つのメッシュ点を1スレッドが担当し、このような計算をZ方向に続ける。例えば、96×96×96のメッシュの場合
グリッドDim 144 × ブロックDim 64 = 96×96で、Z方向に96回計算すると∇^2の計算が終了する。
この8×8×3の小さい領域は、 XY(8×8)×3面としてシェアード・メモリー上にポインタとして確保できるので、右図の赤い矢印に示すように、 Z軸方向に移動するにはポインタとして確保したXY面をシフトし、 1面のみのデータの更新をすればよい。また、XY面を8×8の小領域にしたのは、全てのデータをシェアード・メモリー上の確保し、計算の高速化をはかるためである [8] 。

 

4.3 経路積分のデータのデバイス上への確保

 実装当初は経路積分のデータ(これは高分子鎖の長さに比例するScalarFieldのアレイになるので大きなメモリ空間を占める)はホスト(メインボード側)上に確保し、ホスト・デバイス間通信を行いながらGPU上で経路積分の計算を行っていた。よって、期待したよりも実行速度が上がらなかった。そこで最後には経路積分のデータもすべてデバイス上で確保し、計算をするようにした。 
 


図4 GPU 上での経路積分の計算方法 

 

5 GPU 利用による高速化の結果

 図5 に示すのが、GPU による高速化をはかったSUSHI の1CPU 時の速度に対する相対速度である。横軸は、種々のメッシュサイズの違いを表している。また、各サイズにおいて、左側のバーが経路積分をホスト側に確保し、ホスト・デバイス間通信を行いながら計算した速度、右側が、全てのScalarField のデータをデバイス側に確保し極力ホスト・デバイス間通信を減らして計算した場合の速度である。これは、あくまでも、SCF の静的計算における一例であるが、まずシステムサイズが増大する程GPU による計算速度は向上すること、さらにホスト・デバイス間通信を減らせば高速化がはかれることがわかった。全てが1 コア+1 GPUの計算ではあるが、GPU により、963 のメッシュで、18.7 倍の高速化が達成された。


図5 GPU 上でのSUSHI の速度比較:鎖長N = 20、ブロック比f = 0:35、システムサイズ17:23(323 メッシュ時)
 

 しかるに、GPU 上にデータを確保することによる高速化は、計算規模がGPU 上のメモリー量で決まることになる。よって、今後、MPI と合わせた計算の大規模化を行いたいと考える。しかし、GPU デバイスは年々改良されており、今後、デバイスメモリの増加も期待されるところである。また、コンピューター・アーキテクチャーの進歩によりホスト・デバイス間通信速度が向上すれば、より簡単にGPU による計算の高速化ができるものと期待できる。


 図6 に示すのは、CUDA の高速計算により出力した計算例である。全てジブロック・コポリマーのGyroid構造[7] に関係するミクロ相分離構造であり、マイナー成分の等値面を表示してある。図6(a) は、無秩序状態から相分離させた場合の相分離後期のスナップ・ショットである。これまで計算負荷が大きく、実行することが難しかった963 のメッシュによる動力学を行っている。図6(b) は、Gyroid 構造とHPL 構造の界面を表示したもので、このような大きなの系も計算がなかなか困難であったが、GPU を利用することにより短時間で計算できるようになった。





図6 GPU を利用したSUSHI の計算例:システムサイズ17:23 (323 メッシュ時) (a) 動力学でのスナップショット、シミュレーション時間t = 1000、鎖長N = 20、ブロック比f = 0:35、システムサイズ51:63(963 メッシュ), (b) Gyroid 構造とHPL(Hexagonally Perforated Lamellar)構造の界面、56 ×64 × 72 メッシュ
  

6 おわりに

 密度汎関数理論を用いて高分子の相分離構造をSCF計算をするOCTA/SUSHIを GPUを利用するように改良し、ある計算例では18.7倍の高速化を達成した。この並列化手法は、非常にコスト・パフォーマンスに優れているといえる。また、今後の技術発展が非常に望まれる技術であることは言うまでもない。
尚、本研究の成果はOCTA2013として公開を予定している (http://octa.jp)。

謝辞

 
この研究は、東京工業大学の先端研究施設共用促進事業 「『みんなのスパコン』TSUBAMEによるペタスケールへの飛翔 」の採択課題「メソ構造を持つ高分子材料のマルチスケール・シミュレーション」として、東京工業大学 学術国際情報センターのスパコン TSUBAME2.0 を用いて行った。
ご支援頂きました東京工業大学の青木先生、渡邊先生、東京工業大学学術国際情報センターの佐々木様、松本様には紙面を借りまして厚く御礼申し上げます。

参考文献

[1] (http://octa.jp), 2010.
[2] M. Doi and S. F. Edwards: The Theory of Polymer Dynamics; Oxford Science: Oxford, 1986.
[3] E. Helfand and Z. R. Wasserman: Macromolecules 9, 879 (1976). E. Helfand and Z. R. Wasserman:
Macromolecules 11, 960 (1978). E. Helfand and Z. R. Wasserman: Macromolecules 11, 994 (1980).
[4] G. J. Fleer, M. A. Cohen Stuart, J. M. H. M. Scheutjens, T. Cosgrove, B. Vincent: Polymers at
Interfaces; Chapman & Hall: London, 1993.
[5] T. Kawakatsu: Statistical Physics of Polymers - An Introduction; Splinger-Verlag: Berlin, 2004.
[6] T. Honda and T. Kawaktsu: Computer simulations on nano-scale phenomena based on the dy-
namic density functional theory. Applications of SUSHI in OCTA system; in "Nanostructured Soft
Matter", A.V.Zvelindovsky, ed., 461-493. Springer-Verlag: Berlin, 2007.
[7] T. Honda and T. Kawakatsu: Macromolecules 39, 2340 (2006).
[8]青木 尊之、額田 彰: はじめてのCUDAプログラミング―驚異の開発環境[GPU+CUDA]を使いこなす!(I・O BOOKS) 工学社,2009.