液体の積分方程式理論の解法と、電子状態計算との連成

液体の積分方程式理論の解法と
電子状態計算との連成時の問題
"アルゴリズムによる計算科学の融合と発展"
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の座標・配向を表す
 分子性の液体を扱うための様々な解法が提案され
ている

液体の構造とは？
気体と固体の場合は？
気体の構造
特定の構造を持たない
固体の構造
格子定数で一義的に決まる
（６つのパラメータ）
液体の構造
ある一定の構造を持っているが、
格子点では表現できない
ある粒子に注目したとき、その
周りに他の粒子が平均的にど
のように分布しているかを見る
液体の構造と分布関数
動径分布関数
密
密
疎
密密
疎
r
分布関数から計算出来る系の熱力学量
分子間ポテンシャル
内部エネルギー：

3
N
U  NkBT    u(r ) g (r )  4 r 2 dr
2
2 0
圧力：
1 2  du(r )
P   k BT    r
g (r )  4 r 2 dr
0
6
dr
圧縮率：
1
1
T 

 k BT k BT

2
g
(
r
)

1

4

r
dr


0
RISM, 3D-RISM, MOZ

分子性液体

2体の相互作用の記述には6次元の関数が必要
(系全体の回転・並進不変性）
ž


R 12
1
ž
2
分子１，２の配向とそれらを結ぶベクトル
液体の積分方程式理論の種類
 Molecular Ornstein-Zernike (MOZ)理論
 Reference interaction site model (RISM)理論
 Three-dimensional RISM (3D-RISM)理論
MOZ (MOLECULAR ORNSTEINZERNIKE)

分子の配向をすべて考慮（自由度は6次元）
R 12

多極子展開を用いる
ž
ž
1
 m n l 
m
n
l
g(1 2 R12 )   g (R12 )  
 R  1 R   2 R  0 Rˆ 12
mnl
           

mnl

 
多極子展開の収束の悪さから大きい分子への応用
は難しい
 近似が少ないため、物理量の再現性も良く高精度な
計算に向く
 反面、計算コストは高い


多極子展開をどこまでとるかに依存
2
RISM (REFERENCE
INTERACTION SITE MODEL)

分子を作用点（interaction site)に分け、その作用点間の
距離のみを扱う(自由度は1次元）
rab

自由度が1次元なので軽量で高速、応用範囲が広い

反面、誘電率を再現出来ない点や異方性の強い分子には
向かないなど欠点も
計算コストは
 (距離グリッド数)ｘ(分子１のサイト数)x(分子２のサイト数）

3D-RISM (THREE-DIMENSIONAL
RISM)

分子を座標に固定し、分子２の作用点の位置ベクトル
で記述（自由度は3次元）
Rb
分子１からの分子２の位置情報を得られる
 複雑な溶質（分子１）を扱うことが出来る
 計算コスト


(グリッド数)3x（分子2のサイト数）
たんぱく質
DNA
ナノチューブ
〜
3D-RISM
いらない子のように
思えるが、「理論的
厳密性」が高いため
高速に計算出来れ
ばそれなりに需要は
あるはず…
アミノ酸
核酸
〜
対象分子の規模
計算コストと対象とする系の比較
RISM
MOZ
単純な分子
(H2O, CH3CN)
計算コスト
RISM

Molecular Ornstein-Zernike

Site-Site Correlation functionの導入

RISM equation
3D-RISM

RISMと異なり、分子２についてのみ平均化


h (r)    R 2  l2  r h(1, 2)d(2)

3D-RISM equation

フーリエ変換により
%
h (k)   c%
  (k) X  (k)

X  (k)  VV (k)  hVV (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)  hVV (k)
グリッドと離散化

相関関数の離散化

g (r)  g xi , y j , zk


グリッド

３次元直交グリッド
z
y

離散フーリエ変換
x
c (si ,t j ,uk )   c (xi  , y j  , zk  )exp i(si xi   t j y j   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)  cs (r)  q f (r)
solute
f (r)   


q
erf  r  r
r  r

solute
解析的にフーリエ変換可能
f (k)   

k2
4 q  4 2 ikr
e
e
k2
LONG RANGEの問題を考慮した場合の
アルゴリズム
c (rijk )  exp    u (rijk )    (rijk )   1    (rijk )
cs (rijk )  c (rijk )  q f (rijk )
N
c (ki, j,k )  v
s

iki , j ,k ri  , j  , k 
cs (ri , j , k  )e
i , j , k 
s
%
c (kijk )  c%
 (kijk )  q f (kijk )
%
h (k ijk )   c%
  (k ijk ) X  ( k ijk )

h (ri, j,k ) 
s
1
3 k
2 
N
3

i , j , k 
  (rijk )  h (rijk )  c (rijk )
iki , j  ,k  ri , j ,k
s
h%
(k
)e

i , j , k 
3D-FFT
収束テクニック

MDIIS (Modified direct inversion in the iterative
subspace)

メリット
ロバスト
 MDIISのルーチン自体が軽量


デメリット


メモリを大量に消費
Newton-Raphson

メリット


収束性するときは速い
デメリット
Jacobianの計算が重い
 良い初期値を必要とする

MDIIS
m

(m1)
(ri )    (ri )  R (ri )   c j  ( j ) (ri )  R( j ) (ri )
(*)
(*)
j n
R( j) (ri )   ( j  ) (ri )   ( j) (ri )







1
snn
0
1
1 L
snn1 L
1 sn1n
M
1
M
smn
N
O
1   
snm   cn

  cn1

 M
smm   cm
sij    Ri (rk )Rj (rk )

m
c
i
in
1
  1 
  0 
 

  0 
 

M
 

  0 
 (m) (ri )
3D-RISM
 ( m ) (ri )
R(m) (ri )
{cn}の計算
k
 (m1) (ri )
NEWTON-RAPHSON + PICARD
HYBRID

τ(あるいはc)をcoarse partとfine
partにわけ、それぞれにNewtonRaphson、Picard法を適用する
Δτは固定
  (ri )   au, Pu (ri )    (ri )
u
coarse part
  a
(n)
1
au,(n1)

a

J


u,
v,


fine part
uv,
(n)
v,
 av, 
3D-RISM

Jacobian
の計算
aの収束判定
Δτの更新
Pはroof function, JはJacobian
coarse partとfine partを決めるため
の２重ループを行う
τの収束判定
積分方程式理論と電子状態理論の連成
積分方程式理論で電子状態理論を組み合わせること
で、溶液内分子の電子状態に溶媒効果を効率的に組
み込む
 溶媒と溶質は互いに相互作用する

溶媒(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  


 e  u (r )   (r )  1    (r)   (r)    (r)  h (r)  c (r)  h (r)
dr


   




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 
(re )
dre  
r  re

 
Z

 4   
 r  r
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
12

  




 r  r



6




グリッド数回の
グリッド数回の１電
１電子積分が必要
子積分が必要
静電ポテンシャル計算の高速化

静電ポテンシャル計算、およびSolvated Fockの計算
が3D-RISM-SCFのボトルネック


２５６３回の１電子積分･･･
高速化のアイデア

１電子積分を高速化
Pople-Hehre
 Martyna-Tuckerman
 フーリエ変換


精度を維持して高速化
１電子積分を減らす

有効静電ポテンシャル(ESP)
精度を犠牲にするが
圧倒的に高速化
有効静電ポテンシャル法(ESP)
分子の周りに分子自体が作る静電ポテンシャルを再現
するように、分子上に点電荷を配置する方法、またはそ
れによって決められた点電荷のこと
 静電ポテンシャルをフィッティングするための数千〜数
万点の１電子積分で済む
 点電荷は通常、溶質分子の原子核上におく
（かならずしもその必要はないが・・・）

π軌道などの再現は難しい（芳香族-芳香族相互作用
等）
 分子内部に埋もれた原子の点電荷の決定には難あり

ESP

位置rに溶質がつくる静電ポテンシャルは
U(r)  U N (r)  Ue (r)
ZA
U N (r)  
A r  RA
Ue (r)  tr(PA(r))
原子核の寄与
電子の寄与
(P)  P   ck c k (A)  
k

qAN,
各原子核上の部分電荷
N
q
シャルは
A
UN (r)  
A r  RA

* (r) (r)
r  r
d r
qAe ,を用いたモデルポテン
Ue (r)  
A
qAe
r  RA
原子核の分は qAN  ZAとおけるので以下では扱わない
ESP
e
 qA
は最小二乗法で決定
 n e

2
 
 l

%



U
(r
)

U
(r
)

2

q

const
 k  e k
e k 
e j
 0
qie 
 j 1

 k 1



ω、l、λはそれぞれサンプル点の重み、サンプル点、Lagrangeの未定乗数
この式から
1a 1tr(PB)  N e 1
q  a tr(PB) 
a 1
t 1
1a 1
e
1
l
(a)ij   k rki1rkj1
k 1

 1 
1   M


1


l
(B) ,i   k rki1 A (rk )
k 1
この方法によれば１電子積分はたかだかl回のみ

グリッド（サンプル点）の生成方法にもよるが、有機分子程度
で数百〜数千、アミノ酸等でも数千〜数万点
3D-RISM-SCFでのアイデア


3D-RISMの利用を考えた
場合、空間を３つの領域
に分けてそれぞれ扱いを
変えてやる
電子交換反発など
が主体になる様な
距離）ではそもそも
静電ポテンシャル
の計算の必要はな
い
問題点


領域IIとIIIの間で不連続
が発生する
IIの領域の高速化が依然
としてネック

ここは専門家に任せたい
十分離れたところ
では点電荷がつく
る静電ポテンシャ
ルで近似
分子近傍は波動関
数の広がりを考慮
した静電ポテンシャ
ルを用いる
まとめ

3D-RISM自体の高速化は以下の２点が問題
3D-FFTの高速化・高並列化
 反復回数を減らす


電子状態との連成では
静電ポテンシャル/溶媒和フォックの計算
 近似的なアプローチ（ESPなど）
 １電子積分自体の高速化


Download Report