広域分散アプリケーション特論
2005前期 月曜 ３時限
場所：情報基盤センター３Ｆ 多目的講義室
担当 青柳 睦aoyagi@cc.kyushu-u.ac.jp
5月16日（月）講義の内容，成績評価方針(再々々々々掲)

計算科学の概論
主要なシミュレーション手法の紹介
続き・・量子論的手法・・から．

講義の内容

グリッドの概要
Ｇｒｉｄコンピューティングとは ■サイエンス分野での利用■ビジネス分野での利用■

計算科学の概論
主要なシミュレーション手法■

サイエンスＧｒｉｄ ＮＡＲＥＧＩ
Globus, Unicoreの現状 ■NAREGIミドルウエア概要■連成計算とその類型化■

テーマ 考え中GT4の動向に依存・・

講義資料はWebで公開
server-500.cc.kyushu-u.ac.jp

計算科学の概論
計算機の発展と計算科学
 シミュレーション手法の特徴
 シミュレーション手法の紹介


古典系・・流体，構造解析
 量子系・・第一原理分子動力学，量子化学
 連成計算の概観量子論的手法 ・・例と課題・・
分子軌道法・・量子化学の１分野
１体近似、平均場近似、それを越えて・・
密度汎関数法エネルギーは全電子密度関数の汎関数
交換相関ポテンシャルの問題
量子モンテカルロ法
強相関量子系
格子ゲージシミュレーション

分子軌道法分子軌道：１個の電子は、原子核と他の電子がつくる平均的ポテンシャルの中を運動する・・として求めた波動関数
Hartree-Fock:一つの電子配置を仮定し、スピン軌道の反対称化された積で変分最小化．その波動関数から平均ポテンシャルを計算
２電子反発積分の演算量が膨大
Self-consistentに計算
非経験的分子軌道法ｖｓ半経験的分子軌道法
Gaussian,GAMESS,MOPAC,…

フラーレン分子の分子軌道

(参考)Hartree-Fock計算法
EHF =SCF法の流れ< Ψ0 | H | Ψ0 > ∂EHF;=0< Ψ0 | Ψ0 >∂Cμiφ = ∑C χN

Initial guess,…

Construct
Fock Matrix

λiiλλFμν = hμν + ∑Dλσ[(μν || λσ) − (μσ || λν)]
λσ

EvaluateIntegrals

FμνCν = ε ∑ SμνCν D = 2∑ C C∑ννnoiiiλσλii

New Density
Converged?
yesPopulations,…

O(N4) integral( μν || λσ ) Calcs.

vsO(N2) Fock, Density Comm.

σiTwo-Electron integrals(参考)村上先生の２電子積分ボード( μν || λσ ) の例

1000基底関数~1012 通りの演算

tei(4,4,4,4)=(((3+2*p*(4*PAx*PBx+PBx**2+PAx**2*(1+2*p*PBx**2)))*(3+2*q*(4*QCx*QDx+QDx**2+QCx**2*(1+2*q*QDx**2)))*f(0,t))/(p**2*q**2)+(4*(3+2*p*(4*PAx*PBx+PBx**2+PAx**2*(1+2*p*PBx**2)))*PQx*(QCx+QDx)*(3+2*q*QCx*QDx)*f(1,t))/(p*q*(p+q))(4*(PAx+PBx)*(3+2*p*PAx*PBx)*PQx*(3+2*q*(4*QCx*QDx+QDx**2+QCx**2*(1+2*q*QDx**2)))*f(1,t))/(p*q*(p+q))(8*(PAx+PBx)*(3+2*p*PAx*PBx)*(QCx+QDx)*(3+2*q*QCx*QDx)*(((p+q)*f(1,t))+2*p*PQx**2*q*f(2,t)))/(p*q*(p+q)**2)+(2*(3+2*p*(4*PAx*PBx+PBx**2+PAx**2*(1+2*p*PBx**2)))*(3+q*(QCx**
2+4*QCx*QDx+QDx**2))*(((p+q)*f(1,t))+2*p*PQx**2*q*f(2,t)))/(p*q**2*(p+q)**2)+(2*(3+p*(PAx**2+4*PAx*PBx+PBx**2))*(3+2*q*(4*QCx*QDx+QDx**2+QCx**2*(1+2*q*QDx**2)))*(((p+q)*f(1,t))+2*p*PQx**2*q*f(2,t)))/(p**2*q*(p+q)**2)+(4*(3+2*p*(4*PAx*PBx+PBx**2+PAx**2*(1+2*p*PBx**2)))*PQx*(QCx+QDx)*(3*(p+q)*f(2,t)+2*p*PQx**2*q*f(3,t)))/(q*(p+q)**3)¥+(8*(3+p*(PAx**2+4*PAx*PBx+PBx**2))*PQx*(QCx+QDx)*(3+2*q*QCx*QDx)*(3*(p+q)*f(2,t)+2*p*PQx**2*q*f(3,t)))/(p*(p+q)**3)(8*(PAx+PBx)*(3+2*p*PAx*PBx)*PQx*(3+q*(QCx**2+4*QCx*QDx+QDx**2))*(3*(p+q)*f(2,t)+2*p*PQx**2*q*f(3,t)))/(q*(p+q)**3)(4*(PAx+PBx)*PQx*(3+2*q*(4*QCx*QDx+QDx**2+QCx**2*(1+2*q*QDx**2)))*(3*(p+q)*f(2,t)+2*p*PQx**2*q*f(3,t)))/(p*(p+q)**3)+((3+2*p*(4*PAx*PBx+PBx**2+PAx**2*(1+2*p*PBx**2)))*(3*(p+q)**2*f(2,t)+4*p*PQx**2*q*(3*(p+q)*f(3,t)+p*PQx**2*q*f(4,t))))/(q**2*(p+q)**4)(8*(PAx+PBx)*(3+2*p*PAx*PBx)*(QCx+QDx)*(3*(p+q)**2*f(2,t)+4*p*PQx**2*q*(3*(p+q)*f(3,t)+p*PQx**2*q*f(4,t))))/(q*(p+q)**4)(8*(PAx+PBx)*(QCx+QDx)*(3+2*q*QCx*QDx)*(3*(p+q)**2*f(2,t)+4*p*PQx**2*q*(3*(p+q)*f(3,t)+p*PQx**2*q*f(4,t))))/(p*(p+q)**4)+(4*(3+p*(PAx**2+4*PAx*PBx+PBx**2))*(3+q*(QCx**2+4*QCx*QDx+QDx**2))*(3*(p+q)**2*f(2,t)+4*p*PQx**2*q*(3*(p+q)*f(3,t)+p*P
Qx**2*q*f(4,t))))/(p*q*(p+q)**4)+((3+2*q*(4*QCx*QDx+QDx**2+QCx**2*(1+2*q*QDx**2)))*(3*(p+q)**2*f(2,t)+4*p*PQx**2*q*(3*(p+q)*f(3,t)+p*PQx**2*q*f(4,t))))/(p**2*(p+q)**4)(4*p*(PAx+PBx)*(3+2*p*PAx*PBx)*PQx*(15*(p+q)**2*f(3,t)+4*p*PQx**2*q*(5*(p+q)*f(4,t)+p*PQx**2*q*f(5,t))))/(q*(p+q)**5)+(8*(3+p*(PAx**2+4*PAx*PBx+PBx**2))*PQx*(QCx+QDx)*(15*(p+q)**2*f(3,t)+4*p*PQx**2*q*(5*(p+q)*f(4,t)+p*PQx**2*q*f(5,t))))/(p+q)**5+(4*PQx*q*(QCx+QDx)*(3+2*q*QCx*QDx)*(15*(p+q)**2*f(3,t)+4*p*PQx**2*q*(5*(p+q)*f(4,t)+p*PQx**2*q*f(5,t))))/(p*(p+q)**5)(8*(PAx+PBx)*PQx*(3+q*(QCx**2+4*QCx*QDx+QDx**2))*(15*(p+q)**2*f(3,t)+4*p*PQx**2*q*(5*(p+q)*f(4,t)+p*PQx**2*q*f(5,t))))/(p+q)**5+(8*(PAx+PBx)*(QCx+QDx)*(15*(p+q)**3*f(3,t)+30*p*PQx**2*q*(p+q)*(3*(p+q)*f(4,t)+2*p*PQx**2*q*f(5,t))8*p**3*PQx**6*q**3*f(6,t)))/(p+q)**6+(2*(3+p*(PAx**2+4*PAx*PBx+PBx**2))*(15*(p+q)**3*f(3,t)30*p*PQx**2*q*(p+q)*(3*(p+q)*f(4,t)+2*p*PQx**2*q*f(5,t))+8*p**3*PQx**6*q**3*f(6,t)))/(q*(p+q)**6)+(2*(3+q*(QCx**2+4*QCx*QDx+QDx**2))*(15*(p+q)**3*f(3,t)30*p*PQx**2*q*(p+q)*(3*(p+q)*f(4,t)+2*p*PQx**2*q*f(5,t))+8*p**3*P
Qx**6*q**3*f(6,t)))/(p*(p+q)**6)(787 MUL, 261 ADD,69 FUNC)

基底関数χの組により、大小の計算粒度が混在
tei(3,1,1,1)=((PAy*(1+2*p*PAx*PBx)*(1+2*q*QCx*QDx)*f(0,t))/q+((((p+q)**4*((PAy+PQy)*q*(1+2*q*QCx*QDx)+p*(PAy+2*PAx*PAy*PQx*q+2*PAy*P
Bx*PQx*q+2*PAx*PBx*PQy*q2*PAx*PAy*q*QCx2*PAy*PBx*q*QCx2*PAy*PQx*q*QCx+2*q*((PAy*(PAx+PBx+PQx))+2*(PAy*(PAx+PBx)*PQx+PAx*PBx*PQy)*q*QCx)*QDx)2*p**2*PAx*PAy*PBx*(1+2*PQx*q*(QCx+QDx)))*f(1,t))/q)+(p+q)*((p+q)*((p+q)*(3*p*PAy+6*p**2*PAx*PAy*PQx+6*p**2*PAy*
PBx*PQx+2*p**2*PAy*PQx**2+4*p**3*PAx*PAy*PBx*PQx**2+p*PQy+2*p**2*PAx*PBx*PQy+2*p*PAy*PQx**2*q+PQy*q+2*p*PAx*PQx*PQy*q(116 MUL, 31 ADD, 2FUNC)
出展：南部＠分子研

密度汎関数理論
W. Kohn and L. J. Sham, Phys. Rev. 140, A1133(1965)

エネルギー（全電子密度の関数）＝
運動エネルギー＋Hartree項＋交換相関ポテンシャル
交換相関ポテンシャル：局所密度近似、密度勾配補正、運動エネルギー密度、…

Car-Parrinello(1985):電子状態と分子の運動を同時に計算
第一原理分子動力学法とも呼ばれる
Tight-binding法、……

量子モンテカルロ法
変分モンテカルロ法
拡散モンテカルロ法
グリーン関数モンテカルロ法
経路積分モンテカルロ法
シュレディンガー方程式の厳密解

格子ゲージシミュレーション
ゲージ理論を、ゲージ不変性を保ちつつ、４次元超立方格子上に定式化(K.Wilson)
代表的なものはQCD simulation
Gauge場(gluon)U:リンク（辺）上(SU(3))
fermion場(quark)ψ:格子点上
Quenched simulation: fermion自由度を無視。det(M(U))を１と見なす
Gluonに関するMetropolice法
Quark propagator（逆行列）を計算

格子ゲージシミュレーション
Full QCD: det(M(U))の変化を評価
確率微分方程式(Langevin方程式）
分子動力学的手法
ハイブリッド法

Quarkはハドロン（陽子、中間子など）の中に閉じこめられていることが示されたコンピュータの発展と計算科学
ENIACから地球シミュレータまで性能比は・・
109程度（２００５年 IBM Blue/Geneで１０10程度）
初期の計算機は事務用に占有され、科学技
術計算はほんの隅っこを利用しているだけであった
1970年代：スーパーコンピュータの登場
科学技術のための計算機
High performance computersの進歩

「Top500」から予想したペタフロップス
地球シミュレータ2のあるべき目標
2007年初頭に1PF

地球シミュ
レータ

IBM ASCIBlueGene/L
200TF
Cray ASCIRedStorm
25TF/60TF

IBM ASCI
Purple 60TF
バージニア工科大“Big Mac”クラスタ
10.3TF

年率186%成長
3年で6.5倍, 6年で40倍(c.f. ムーアの法則は 4倍)

並列計算機の現状
演算性能・・地球シミュレータ(NEC製)が世界第二位
Processor Node数=640
ＣＰＵ数＝5120台
640×640の単段クロスバ結合

ピーク性能＝40Tflops
主記憶容量＝10TB
外部記憶容量＝約６７０ＴＢ

TOP５００の計算機アーキテクチュア内訳
Cluster型Architectures現在の主な並列計算機
・SMP(対称型マルチプロセッサ)– メインフレーム，サーバー、科学技術計算– SMPのクラスター構成もある• PC クラスタ– 安価な PC を多数結合，スーパーコンピュータの代用となりうるのか？

• チップ内並列処理– スーパースカラ, VLIW，Dual Coreの有効活用– 動画像処理，画像生成等
SMP（16～128cpu) のクラスター構成が主流に・・

計算科学はどこへ・・？
コンピュータの能力の飛躍的発展．
Petaflops will be available soon.
多くの要素からなる系のコンピュータシミュ
レーション(particles, atoms, molecules,stars, proteins)→全体シミュレーション、連成シミュレーション
それでも現実の系のシミュレーションはできな
い(How many atoms in one gram!)Computer power is lacking!
→グリッドコンピューティング？

「Memory Wall」問題とinter-Connect
演算よりも「データの移動」が律速（高コスト）

マイクロプロセッサ
105
性能向上率

vs. DRAM
性能向上率の差 48%/年

104
103
102
101

マイクロプロセッサ
55%/年

スループットとレイテンシ問題

RISCおよびオンチップ
キャッシュの登場以降

1
1980

1990
主記憶(DRAM)
7%/年

2000

2005
Highband-Dram, Rambus,SDRAM・・
年
キャッシュ有効活用法，プリフェッチ
nonBlocking $, multiThreading, ・・【出典】J.L. Hennessy, D.A. Patterson, Computer Architecture: A Quantitative Approach 3rd Edition, Fig. 5.2

演算よりも「データの移動」が律速（高コスト）• ＣＰＵ⇔１次キャッシュ⇔2次キャッシュ⇔メモリ⇔バッファ⇔InterConnect⇔
NIC⇔Network Switch ⇔組織内ネットワーク⇔
Computing-Grid網・・・• データ移動（帯域幅と遅延）の階層性

連成計算の概観
シミュレーション対象の全体を「複数の部分系に分割」し、部分系の計算途中でお互いのデータを交換しながら全体を解く手法．
部分系Aを解く手法 ⇔ 部分系Bを解く手法(例) 分子軌道計算 / 分子動力学計算
構造力学計算 / 分子動力学計算
流体計算 / 構造力学計算
統計力学計算 /分子軌道計算 ・・・ｅｔｃ
データ交換の頻度とデータ量は、部分系間の依存度、結合の強さにより、弱連成・・強連成と区別することがある．

分散シミュレーションの形態■ シミュレータＡの出力をシミュレータＢの入力へ
シミュレーション連携、プリ・ポスト処理、可視化
標準化

ｱﾌﾟﾘA

ｱﾌﾟﾘB

ｱﾌﾟﾘA

ｱﾌﾟﾘB

■ シミュレータＡとシミュレータＢでデータ交換をしながら実行
連成シミュレーション

コント
ローラ

Mediator

ＭＰ型

RPC型
ｱﾌﾟﾘA

ｱﾌﾟﾘB

ｱﾌﾟﾘB

■ シミュレータＡがデータベースＣを参照しながら実行
標準化

ｱﾌﾟﾘA

ﾃﾞｰﾀﾍﾞｰｽＣ

ｱﾌﾟﾘA

ｱﾌﾟﾘA

ﾃﾞｰﾀﾍﾞｰｽＣ

溶媒＋溶質連成シミュレーションの例
×10～100サイクル

溶媒領域
アダティブ
メッシュ形成Mediator
離散点登録
rvdw等のデータ

溶質分子属性設定

・等間隔ﾒｯｼｭ使用
・ﾒｯｼｭ点数大

マルチスケール問題に共通に顕在化する変換技術
２体相関関数溶媒分布
計算
数GB

有効ﾎﾟﾃﾝｼｬﾙ
又は電荷ﾃﾞｰﾀ

RISM

共有メモリSMP向き, 所要メモリ量：数10GB

FMO

溶質電子
構造計算

数GB

PCクラスタ向き, 所要メモリ量： 数10GB

MPICH-G

またはGridMPI通信＠Globus

RISM-FMO連成シミュレーションにおける計算のワークフロー
RISM

Mediator 変換ﾗｲﾌﾞﾗﾘ

FMO

ｱﾀﾞﾌﾟﾃｨﾌﾞﾒｯｼｭ

等間隔ﾒｯｼｭ

離散点
位置探索

溶媒分布計算
溶媒密度
分布

異種ﾒｯｼｭ
ﾃﾞｰﾀ
変換

溶質電子状態計算
溶媒密度
分布
原子
電荷

溶媒ー溶質
ｸｰﾛﾝｴﾈﾙｷﾞｰ計算

Mediatorは、アプリケーション間での高度な意味変換機能をサポートするパッケージ群である。
計算手法の異なるアプリケーション間の物理量変換は、グリッド上でMediatorを介して自動的に実行されるため、プログラムの修正を最小限に抑えた連成シミュレーションを構築できる。
1．連成シミュレーションの構築手順

2．グリッド連成ミドルウェア Mediatorの特徴機能

Mediatorの提供するAPIを用いてレガシーコードをコンポーネント化。グリッド標準の通信フレームワーク上で連成化。
①ソースコードにMediatorのAPIを追加

Sim. AMediatorA

Sim. AMediatorA

Sim. AMediatorA

PROGRAM MAIN
CALL MPI_INIT

ナノアプリケーション
レガシーコード

CALL MPI_FINALIZE

連成シミュレーションの管理
call dex_join (color, ﾛｰｶﾙｺﾐｭﾆｹｰﾀ)
離散点情報の設定
call dex_register_quantum(位置関係, 離散化手法,離散点数,番号,座標,要素情報,属性)
データ送受信call dex_send(byte,物理量a) 又は
call dex_receive(byte,物理量b, 変換関数)

END②コンパイル及びMediatorライブラリをリンク③アプリケーションとMediatorをGridMPIジョブとして実行

MediatorB
離散点の位置関係を探索しデータを自動的に通信・変換 Mediator

Sim.
Sim.
BB

Sim. BB

・高度な意味変換機能
球内相関、矩形内相関、第一近接、最近接・離散点位置の高速・並列探索
粒子法、差分法、有限要素法
・標準通信ライブラリへの対応
MPICH、Score、MPICH-G2、GridMPIMediatorによる離散点データ交換
グリッド連成ミドルウェア（Mediator）は高度セマンティック変換をサポートすることによって、指定された相関関係にある離散点上の物理量を交換する
グリッド連成ミドルウェア library calls:
Call register (位置関係,離散手法,Call join
(ﾌﾟﾛｾｽ管理)

離散点数,番号,座標,要素情報,属性)

Sim. AMediatorA

Sim. AMediatorA

Sim. A

Sim. AMediatorB

離散点の位置関係を探索し Mediator
B
データを自動的にMediator
A
通信・変換
Mediator Middleware
A

Mediator
B

Sim.
Sim.
BB

最近接

第１近接

球内相関

矩形内相関

Sim. B

Sim. BsystemCall send(byte,物理量a)Call receive(byte,物理量b, 変換関数)

ＧｒｉｄＭＰＩによるグリッド連成ミドルウェアの応用
RISM
RISM

溶液中の電子構造

FMO
FMO

RISM

Reference Interaction
Site Model

無限系に対する溶媒分布解析プログラム
FMO

溶媒分布Mediator

電子構造
異なる離散点上のデータ変換

Site A

Fragment MolecularOrbital method

ナノ高分子の電子構造解析プログラム

Mediator

Site B
MPICH-G2, Globus

グリッド連成ミドルウェア（Mediator)を開発し、複数のナノアプリケーション（RISM，FMO)をグリッド環境下で接続し、超巨大水分子を対象に溶液中の電子構造を計算

リゾチウムよる溶菌作用を対象に、 RISM-FMO連成計算を用いて分子レベルのメカニズムを解析、連成計算の有効性を確認した。
1．リゾチウムの溶菌作用のメカニズム

2．RISM-FMO計算による水溶液中でのプロトン移動の解析

１）細胞壁を形作る糖の一種(ペプチドグリカン(図２))の結合を加水分解（切断）する
２）上の加水分解は活性部位(グルタミン酸35）からペプチドグリカンへのプロトン移動により起こる

１）分子動力学計算によりペプチドグリカンを結合したリゾチウムの水和構造を得た(図３)２） この構造に基づきRISM-FMO連成計算により水溶液中でのプロトン移動について調べた結果、プロトンの移動に伴い、ペプチドグリカンの電子構造に大きな変化が認められた(図４－Ａ，Ｂ)

リゾチウム

ペプチドグリカン
NAM
NAG(NAG)2CH3を結合した水和構造

１）ライゾザイムの２）
３）○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○図４－Ａ

グリッド上でのRISM-FMO連成計算
NAM

図１

NAM: N-アセチルムラミン酸
NAG: N-アセチルグルコサミン
図２

図３

図４－Ｂ

Structure and functionality of Lysozyme
A catalytic mechanism of hydrolysis in Lysozyme isanalyzed, in which proton transfer leads to bacteriolysis
in peptidoglycan of bacterial cell wall.
Peptidoglycan
Lysozyme
NAM
H＋Active site of Glu35

H:
C:N:O:NAG

NAM
NAM: N-Acetylmuramin
NAG: N-AcetylglucosaminLysozyme functionality in solvent
Whole electronic structure of Lysozyme in solvent isanalyzed by using RISM-FMO and electronic densitychanges according to the position of proton transfer.
Hydrate structure
binding (NAG)2CH3

Whole electronic structure

RISM-FMO
on Grids

Proton at Glu35

Solvent
distribution
Proton at Glycosidic O

複雑現象（複合，連成系）の理解へ向けた
更なる高性能のニーズ
部分 Æ 全体
個別現象 Æ 複合

金型・プレス機も含めたシミュレーション
モデルの統合

前面衝突用
側面衝突用

後面衝突用

Source: 姫野博士＠理研

全ての部品を入れてモデル化

Source: 姫野博士＠理研

更なる高性能のニーズ（続）
更なる高性能のニーズ
マルチスケール人体シミュレーション
ミクロからのアプローチ

マクロからのアプローチ
連続体シミュレーション
血管系モデリング
流体・熱・構造
骨格モデル

MD・第一原理・量子化学シミュレーション
タンパク構造解析
薬剤反応性解析

分子ネットワーク解析

タンパク・DNA細胞

1-8~-6

1-5~-4

大きさ

化学反応の連成組織

1-3~-2

全身循環器

1-110

Meso

Micro

DDS
薬剤開発
テーラーメード医療 再生医療

器官

HIFU

Macro
結石破砕
外科手術
カテーテル
マイクロマシン

治療技術

マクロな医療用生体計算での見積もり（例）•••ボクセルモデルによる人体モデリング理研東大で2009年度末に開発
生体膜の取り扱いが重要– 生体膜0.01mmをボクセル表現すると•

人体全身––––

2 m x 0.75 m x 0.5 mと仮定
平均1mm：750M点(240GW)
メモリーは十分実効１00MFLOPSで１週間、速度は速くする必要あり< 2m– 平均0.25mmだと：48G点(15TW)– 実効１００Mで６４週間（約１年）、３日/ES実効１０T– 診療用途ではこの位の計算が１時間、よって実効１P– 0.01mm：240PW、アダプティブメッシュで1/100、数PW、計算速度実効１Pで３０日< 0.5m
< 0.75m
Source: 姫野博士＠理研

もの作り分野の見積もり（例）••••

現状の計算メッシュ生成の問題
から今後はボクセルベースの解
析へと進む
問題はボクセルサイズをどこまで扱えるか
現在、ものつくりでの最大の問題
は自動車の衝突解析
スポット溶接部分を考えると最小ボクセル0.1mm– 車サイズ５ｍｘ２ｍｘ２ｍから20T
ボクセル– メモリー：１PW、現状の100億倍の計算量– アダプティブで1/100の演算量とすると実効1PFで数時間
Source: 姫野博士＠理研

今日のおわりに・・【まとめ】

第三の科学手法としてシミュレーションは研究開発現場で重要な役割を果たしている．
 構成的手法（ミクロの要素からマクロを理解） 現実の物理は連成・複合系

Multi-Scale, Multi-Physics Simulation へ．