Transcript ppt
液体の積分方程式理論の解法と
電子状態計算との連成時の問題
"アルゴリズムによる計算科学の融合と発展"
2009/04/22-23
吉田紀生
分子科学研究所
はじめに
液体の積分方程式理論
液体の構造を統計力学に基づいて記述する理論
連立方程式を繰り返し計算により解く
畳み込み積分を解くためにフーリエ変換を用いる
電子状態理論との連成
溶媒和分子の電子状態に溶媒効果を加えるのに効
果的な方法
溶媒和フォック行列の計算、静電ポテンシャルの計
算が律速に
有効な高速化手法を模索中・・・
液体の積分方程式理論とは
液体の積分方程式理論とは
液体の状態・構造を記述する理論
積分方程式(Ornstein-Zernike方程式)とclosure
の連立方程式
h(1,2) c(1,2) c(1,3)h(3,2)d(3)
g(1,2) exp u(1,2) h(1,2) c(1,2) b(1,2)
1,2は分子1,2の座標・配向を表す
分子性の液体を扱うための様々な解法が提案され
ている
液体の構造とは?
気体と固体の場合は?
気体の構造
特定の構造を持たない
固体の構造
格子定数で一義的に決まる
(6つのパラメータ)
液体の構造
ある一定の構造を持っているが、
格子点では表現できない
ある粒子に注目したとき、その
周りに他の粒子が平均的にど
のように分布しているかを見る
液体の構造と分布関数
動径分布関数
密
密
疎
密 密
疎
r
分布関数から計算出来る系の熱力学量
分子間ポテンシャル
内部エネルギー:
3
N
U NkBT u(r) g(r) 4 r 2 dr
2
2 0
圧力:
1 2 du(r )
P kBT r
g (r) 4 r 2 dr
0
6
dr
圧縮率:
T
1
1
kBT kBT
2
g
(
r
)
1
4
r
dr
0
RISM, 3D-RISM, MOZ
分子性液体
2体の相互作用の記述には6次元の関数が必要
(系全体の回転・並進不変性)
ž
R12
1
ž
2
分子1,2の配向とそれらを結ぶベクトル
液体の積分方程式理論の種類
Molecular Ornstein-Zernike (MOZ)理論
Reference interaction site model (RISM)理論
Three-dimensional RISM (3D-RISM)理論
MOZ (MOLECULAR ORNSTEINZERNIKE)
分子の配向をすべて考慮(自由度は6次元)
多極子展開を用いる
ž
R12
1
m n l m
n
l
g(12R12 ) g (R12 )
R 1 R 2 R 0 Rˆ 12
mnl
mnl
ž
多極子展開の収束の悪さから大きい分子への応用
は難しい
近似が少ないため、物理量の再現性も良く高精度な
計算に向く
反面、計算コストは高い
多極子展開をどこまでとるかに依存
2
RISM (REFERENCE
INTERACTION SITE MODEL)
分子を作用点(interaction site)に分け、その作用点間の
距離のみを扱う(自由度は1次元)
rab
自由度が1次元なので軽量で高速、応用範囲が広い
反面、誘電率を再現出来ない点や異方性の強い分子には
向かないなど欠点も
計算コストは
(距離グリッド数)x(分子1のサイト数)x(分子2のサイト数)
3D-RISM (THREE-DIMENSIONAL
RISM)
分子を座標に固定し、分子2の作用点の位置ベクトル
で記述(自由度は3次元)
Rb
分子1からの分子2の位置情報を得られる
複雑な溶質(分子1)を扱うことが出来る
計算コスト
(グリッド数)3x(分子2のサイト数)
たんぱく質
DNA
ナノチューブ
〜
3D-RISM
いらない子のように
思えるが、「理論的
厳密性」が高いため
高速に計算出来れ
ばそれなりに需要は
あるはず…
アミノ酸
核酸
〜
対象分子の規模
計算コストと対象とする系の比較
RISM
MOZ
単純な分子
(H2O, CH3CN)
計算コスト
RISM
Molecular Ornstein-Zernike
Site-Site Correlation functionの導入
RISM equation
3D-RISM
RISMと異なり、分子2についてのみ平均化
h (r) R2 l2 r h(1,2)d(2)
3D-RISM equation
フーリエ変換により
%
h (k) c%
(k) X (k)
X (k) VV (k) hVV (k)
Closure
g (r) exp u (r) h (r) c (r)
3D-RISMの基本的アルゴリズム
3D-RISM equation と Closureをiterativeに解く
c (r) exp u (r) (r) 1 (r)
c (k) c (r)eikr dr
%
h (k) c%
(k) X (k)
ikr
h (r) h%
dk
(k)e
(r) h (r) c (r)
solute-solvent interaction potential
u (r) 4
solute
12 6 q q
r r r r
r r
solvent-solvent correlation term
X (k) VV (k) hVV (k)
グリッドと離散化
相関関数の離散化
g (r) g xi , y j ,zk
グリッド
3次元直交グリッド
z
y
離散フーリエ変換
x
c (si ,t j ,uk ) c (xi , yj ,zk )exp i(si xi t j yj uk zk )v
i
j
k
LONG-RANGE INTERACTION
直接相関関数の長距離挙動
c (r) u (r),
u (r) 4
solute
1/r が収束が悪い
r
12 6 q q
r r r r
r r
有限のセルサイズでは、離散フーリエ変換の精度悪化
直接相関関数の分離
振る舞いの良い
関数
c (r) cs (r) q f (r)
solute
f (r)
q
erf r r
r r
4 q 2
f (k) 2 e 4 eikr
k
solute
解析的にフーリエ変換可能
k2
LONG RANGEの問題を考慮した場合の
アルゴリズム
c (rijk ) exp u (rijk ) (rijk ) 1 (rijk )
cs (rijk ) c (rijk ) q f (rijk )
c (ki, j,k ) v
s
N
i , j , k
iki, j ,k ri , j ,k
cs (ri , j , k )e
s
%
c (kijk ) c%
(kijk ) q f (kijk )
%
h (kijk ) c%
(kijk ) X ( kijk )
N
1
iki , j ,k ri, j ,k
3
s
%
h (ri, j,k )
k
h
(k
)e
i , j , k
2 3 i
, j , k
s
(rijk ) h (rijk ) c (rijk )
3D-FFT
収束テクニック
MDIIS (Modified direct inversion in the iterative
subspace)
メリット
ロバスト
MDIISのルーチン自体が軽量
デメリット
メモリを大量に消費
Newton-Raphson
メリット
収束性するときは速い
デメリット
Jacobianの計算が重い
良い初期値を必要とする
MDIIS
m
(m1)
(ri ) (ri ) R (ri ) c j ( j ) (ri ) R( j ) (ri )
(*)
(*)
j n
R( j) (ri ) ( j ) (ri ) ( j) (ri )
0 1
1
1 snn snn1
1 sn1n
M M
1 smn
N
L 1
L snm cn
cn1
O
M
smm cm
sij Ri (rk )Rj (rk )
m
i
1
0
0
M
0
3D-RISM
( m ) (ri )
R(m) (ri )
{cn}の計算
k
c 1
in
(m) (ri )
(m1) (ri )
NEWTON-RAPHSON + PICARD
HYBRID
τ(あるいはc)をcoarse partとfine
partにわけ、それぞれにNewtonRaphson、Picard法を適用する
Δτは固定
(ri ) au, Pu (ri ) (ri )
u
coarse part
a
(n)
1
au,(n1)
a
J
u,
v,
fine part
uv,
(n)
v,
av,
3D-RISM
Jacobian
の計算
aの収束判定
Δτの更新
Pはroof function, JはJacobian
coarse partとfine partを決めるため
の2重ループを行う
τの収束判定
積分方程式理論と電子状態理論の連成
積分方程式理論で電子状態理論を組み合わせること
で、溶液内分子の電子状態に溶媒効果を効率的に組
み込む
溶媒と溶質は互いに相互作用する
溶媒(g(r))
RISM-SCF (Ten-no et al.)
3D-RISM-SCF (Sato et al.)
MOZ-SCF (NY et al.)
溶質(Ψ)
溶質の電子状態(Ψ)と
溶媒の分布(g(r))を
自己無撞着に決定。
3D-RISM-SCF
Quantumに扱う溶質が液体論で古典的に扱われる溶
媒中に無限希釈で溶けている状態を考える
系のTotal Helmholtz Free energyは
A Esolute
溶質の寄与
Esolute solute Hisolated solute
溶媒の寄与
v d u (r)g (r)dr
1
0
1
1
2
h
(r)
c
(r)
h
(r)c
(r)
v 2
2
FORMALISM OF 3D-RISM-SCF
Lagrangian
L A im (Sim im )
i
m
これは相関関数(c, h, τ)および溶質の波動関数(ϕi)お
よびMO係数の汎関数と見なせるので、偏分は
L
u (r) (r)
dr
e
1 (r) (r) (r) h (r) c (r) h (r)
%
1
%
%
dk
h
(k)
c
(k)
(k)
c%
(k)
(2 )3
2 i ij h ijkl gkl ij
ij
kl
q
r r
g (r )dr ij j
SCHEME OF 3D-RISM-SCF
真空中で溶質となる分子の電子状態を計算
isolated
溶質分子が作るポテンシャルを計算
u (r) q
12
6
(re )
Z
dre
4
r r
r r
r re
r r
3D-RISMを解く
g r
溶媒分布からSolvated Fock Matrixを計算
Fij F
isolated
ij
q
i* (re ) j (re )
r re
dre g (r)dr
溶媒分布の下での溶質分子の電子状態を計算
solvated
グリッド数回の
グリッド数回の1電
1電子積分が必要
子積分が必要
静電ポテンシャル計算の高速化
静電ポテンシャル計算、およびSolvated Fockの計算
が3D-RISM-SCFのボトルネック
2563回の1電子積分・・・
高速化のアイデア
1電子積分を高速化
Pople-Hehre
Martyna-Tuckerman
フーリエ変換
精度を維持して高速化
1電子積分を減らす
有効静電ポテンシャル(ESP)
精度を犠牲にするが
圧倒的に高速化
有効静電ポテンシャル法(ESP)
分子の周りに分子自体が作る静電ポテンシャルを再現
するように、分子上に点電荷を配置する方法、またはそ
れによって決められた点電荷のこと
静電ポテンシャルをフィッティングするための数千〜数
万点の1電子積分で済む
点電荷は通常、溶質分子の原子核上におく
(かならずしもその必要はないが・・・)
π軌道などの再現は難しい(芳香族-芳香族相互作用
等)
分子内部に埋もれた原子の点電荷の決定には難あり
ESP
位置rに溶質がつくる静電ポテンシャルは
U(r) UN (r) Ue (r)
ZA
UN (r)
A r RA
Ue (r) tr(PA(r))
原子核の寄与
電子の寄与
(P) P ck c k (A)
k
qAN,
各原子核上の部分電荷
N
q
シャルは
A
UN (r)
A
* (r) (r)
r r
d r
qAe ,を用いたモデルポテン
r RA
qAe
Ue (r)
A r RA
原子核の分は qAN ZAとおけるので以下では扱わない
ESP
e
qA
は最小二乗法で決定
n e
2
l
%
U
(r
)
U
(r
)
2
q
const
0
e k
e j
e k e k
qi
j 1
k1
ω、l、λはそれぞれサンプル点の重み、サンプル点、Lagrangeの未定乗数
この式から
1a1tr(PB) Ne 1
q a tr(PB)
a 1
t 1
1a 1
e
1
l
(a)ij k rki1rkj1
k1
1
1 M
1
l
(B) ,i k rki1A (rk )
k1
この方法によれば1電子積分はたかだかl回のみ
グリッド(サンプル点)の生成方法にもよるが、有機分子程度
で数百〜数千、アミノ酸等でも数千〜数万点
3D-RISM-SCFでのアイデア
3D-RISMの利用を考えた
場合、空間を3つの領域
に分けてそれぞれ扱いを
変えてやる
電子交換反発など
が主体になる様な
距離)ではそもそも
静電ポテンシャル
の計算の必要はな
い
問題点
領域IIとIIIの間で不連続
が発生する
IIの領域の高速化が依然
としてネック
ここは専門家に任せたい
十分離れたところ
では点電荷がつく
る静電ポテンシャ
ルで近似
分子近傍は波動関
数の広がりを考慮
した静電ポテンシャ
ルを用いる
まとめ
3D-RISM自体の高速化は以下の2点が問題
3D-FFTの高速化・高並列化
反復回数を減らす
電子状態との連成では
静電ポテンシャル/溶媒和フォックの計算
近似的なアプローチ(ESPなど)
1電子積分自体の高速化