フラグメント分子軌道法およびエネルギー
表示法を活用した自由エネルギー計算の
GPGPUによる高速化
Acceleration of FMO based MD/ER by GPGPU

Ryota Koga 古賀 良太
President of X-Ability Co.,Ltd. ㈱クロスアビリティ
共同研究者
古川祐貴（クロスアビリティ）、安田耕二（名古屋大学）
松林伸幸、桜庭俊（京都大学）
20 Jun 2012近畿化学協会コンピュータ化学部会＠大阪産業創造館

X-Ability Co.,Ltd. ㈱クロスアビリティ• オフィス５つ（設立：2008年1月15日）
国内：Hongo, Tokyo + The K-comp, Kobe
海外：Thailand, Indonesia, and U.S. branch• 事業２つ(1) 計算科学（化学分野が主）
GPGPU computation、分子モデリング・メタスケジューラ
・MPIデバッガ・計算用ラップトップの販売、等(2) センサネットワーク
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

計算化学で設計した薬©wikipedia
Erlotinib hydrochloride©wikipedia
Zanamivir(trade name Tarceva)

trade name Relenza
インフルエンザ阻害剤

肺がん、膵臓がん
Steve Jobs might have used
this to extend his life…

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

3

挑戦的テーマ： 量子効果を考慮したタンパク－
リガンドの結合自由エネルギー計算のプロセス• (a) FMO based MD– QM(FMO)電荷のMD• (b) MD/ER– MDトラジェクトリを使ってER計算• (a) + (b) FMO based MD/ER -> Free energy– タンパク全体計算は可能だが高精度のリガンドバ
インディングスキームは挑戦的研究テーマ

本トークは (a) GPGPU FMO, (b) MD/ER, and (a) +(b) FMO based MD/ERのためのGUI、と流れます。
4アブストラクト (a) FMO based MD• Vectorized PRISM by AVX on Sandybridge CPU• HF (Hartree-Fock) : x 3

Benchmark : multicore CPUs

• DFT by CUDA on NVIDIA GPU• J-matrix formation : x 245• Exchange correlation : x 4
NVIDIA GPU• ERI J-matrix (Coulomb potential) by CUDA on NVIDIA GPU• FMO Environmental Electrostatic Potential : x 8.8-25.0
GAMESS-US(free application)

HF, DFT, FMO ESP
with CUDA/AVX scratch code
高速化

GAMESS-US+ XA-AVX/CUDA-QM

CPU[AVX]/GPU[CUDA] ハイブリッド演算• CPUは必ず必要– CPU+GPU architectureが一般的– GPUを使わずSIMD並列化に適しているならば、SSE/AVXは良い。
• AVXは４つの倍精度演算を同時に行う• AVXは理論値ではSSEの２倍速い– AVXはIntel SandyBridge CPU以降で動作

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

AVX : Pros and Cons
基本的にベクトルプロセッサ
Pros

Cons

- SSEの倍速- SIMD並列化- FPU演算と同じメインメモリの活用-

条件分岐 and/or ランダム
メモリアクセスには適していない
自動ベクトル化が単純な
例では効率的Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

実装計画
AVX

GPU

- AVXのintrinsic関数を陽に使うのを避けるため、C++によるオーバーロード演算子
を定義- 最内部ループの分岐を避けるためのC/C++関数ポイ
ンタの利用

-

-

GPUアーキテクチャに適したベス
トチューニングのアルゴリズム
Fermiコアの共有メモリのL1キャッシュを活用
複雑なGPUスレッド制御（同期、キュー）のためのpthreadによるCPUスレッド (pthread)
CPUスレッド並列のためのOpenMPの利用
単精度・倍精度の混合演算

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

シュレディンガー方程式
H (1 N )(1 N )  E(1 N )波動関数を短縮ガウス基底を使って展開する（LCAO近似）

s

pz

dyz

 i (r )   Ck(i )  k (r )kKs ：  (r )   d k ek 1

 k ( r  A )K

2

px：  (r )   d k ( x  Ax )e  k ( r  A)

2

k 1

F (C )C  SC 固有値問題(Self consistent)
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

9

Hartree-Fock SCF (HF-SCF)の手続き
F (C )C  SC

：SCF(ab | cd )  

 a (r )  b (r )  c (r ' )  d (r ' )dr ' dr ：ERI
| r  r '|F  H    Ca C*a 2      
2/ N

core

a



1


core
 H 
  D       
2



core
 H 
 G
Density Matrix (D) : Initial guessの後の非線形方程式を解くために毎SCFサイクルでアップデートされる。
ERIs (ab|cd) : １回だけ計算してメモリに置いとけばいいかもしれないが、大量のO(N4) メモリを必要とする。よって、ディスクI/Oを減らすためにこれらは毎サイクルで再計算される(Direct SCF)。このステップがボトルネックとなる。
Density Matrix×ERI = J-matrix : クーロンポテンシャル行列
Density Matrix×ERI = K-matrix : HF交換行列。大量のレジスタを必要とするためにGPUでの加速が大変難しい。

エルミートガウス基底| p]  H t ( x  Px )e ( x  Px ) 2

 ( y factor )  ( z factor )t-th Hermite polynomial
A Product of two Gaussian functions |ab] is expanded exactly by Hermite
Gaussians |p].
(Cartecian) Gaussian

Hermite Gaussian

J ab  Dcd [ab | cd ]~Dcd | cd ]  Dq | q ]~~J p  Dq [ p | q ]~J ab | ab]  J p | p ]ERI calculation in Hermite Gaussians space is faster than that in thenormal Gaussians.
Because of the communication cost we can’t get each ERI [p|q] from GPU.
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

11

J-matrixの手続き（クーロン）
Constructshell-pair data|ab], |p], Dq
expand themin terms ofprimitive HGs

Sort |p] and |q]
for faster exitof loops.

J p   [ p | q]Dq|p]: [p|p]1/2|q]: [q|q]1/2 Dq

descending orderconverted to Jabp, q, Dqq
1. Each thread
between eachpair of [p|q].

Accumulated Jp

2. No Integral
symmetry to
avoid expensive
memory access.

Less than 10% of total cost
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

密度汎関数法(DFT)
Wave function of electron (orbital)[ 2 / 2  veff (r )] i (r )   i i (r )Kinetic
energyCoulomb potential from nuclear (Hcore) andelectron (J) + exchange-correlation potential

Hartree-Fock :coreF  H 
 J   K 

DFT :

core
F  H 
 J   v

Exchange correlation vkl    k (r )  l (r ) f (  (r ),  (r ))dr

Acceleration of J-matrix + Exchange Correlationis necessary.
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

交換相関項の手続き
E xc   wi  (  (ri ),  (ri )),  (r )   (r )   (r )i  (ri )2 (ri ) l (ri )  (ri )   l (ri )

 


 k v xc  l   wi  k (ri )
i

Quadrature point ri andweights wi are generated.

D(ri)∇(ri)

Exc, potential on
quadrature point fi, gi ( ri )f i  wiElectron density (ri),∇(ri)on quadrature points
parallelization

 (ri )   Dkl  k (ri )  l (ri )kl ( ri )g i  2 wi ( ri )fi, givkl matrix

Vxc matrix：
vkl    f i  k (ri )  g i   k (ri ) l (ri )Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

i14

DFTの加速 [GPU]Model, Basis set

Time [sec]Total energy [a.u]GAMESS

929.407-2912.2041896614

This work

305.709-2912.2041830108

GAMESS

1296.509-2927.4589680121

This work

370.807-2927.4589838167

GAMESS

2186.225-3772.6098820622

This work

651.099-3772.6098692643

GAMESS

3010.225-3792.2248655337

this work

800.743-3792.2248881864

Paclitaxel(C47H51NO14), 3-21G

Paclitaxel(C47H51NO14), 6-31G

Valinomycin(C54H90N6O18), 3-21G

Valinomycin(C54H90N6O18), 6-31GCopyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

DFT計算の詳細 [GPU]Valinomycin,BLYP/3-21G unit:sec
GAMESSThis work

speedup

HF Fock Matrix
formation

248.92776.241x 3.27DFT J matrix
formation

619.465

2.529

x 244.94(※)

DFT exchange
correlatin matrix
formation

972.364218.387

x 4.45

total

2186.225

651.099x 3.36※カーテシアンガウス基底で陽に計算されるため、オリジナルの
GAMESS J-matrixがとても遅い
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

PRISMの実装 (1) [AVX]• PRISM Algorithm
ERIを評価する最速アルゴリズムの１つOne of
the fastest algorithms to evaluate ERIs.
Gaussianが採用している
“パス”に従い漸化式を使ったERIにBoys 関数[0](m) ∝∫u^m exp(-Tu)duを変換する。異なった短縮長や軌道角運動量の最適経路を選択する。

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

PRISMの実装 (2) [AVX]
全ての漸化式(Integral A) = a * (Integral B)+ b * (Integral C)
同じ漸化式が同一シェル型（軌
道角運動量、bra/ketの短縮長）の中間積分に適用される。
SIMD並列化に向いている。

Gill and Pople(1991)

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

PRISMの実装 (3) [AVX]
class double4# operator overload{
public:double4& operator=(double d){m_d = _mm256_set1_pd(d);return *this;}double4 operator+(const double4& dd){
return double4(_mm256_add_pd( m_d, dd.m_d) );}// … define other operators__m256d m_d;}
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

弊社(XA)とGaussian社によるPRISMの違い
Gaussian

X-Ability

What is it

それぞれのシェル４つ組をbra/ket
短縮シェルの型と数によって、漸化式を動的に生成する。

全てのPRISMの“パス”のソースコー
ドを生成し、それぞれの計算対象に対して動的に適切な関数ポイントを
設定する

Pros

昔ながらのベクトルプロセッサに合うように漸化式を使うことで、
シェル４つ組を同時にループで処理することができる

８つのシェル４つ組が同時に処理され、結果的にキャッシュミスを避けることができるCons

キャッシュミスが生じるのでスカラ コードは数万行になり、バイナリはプロセッサ向きではない。コンパイ 300-400MBの容量がある
ラはプログラムを効率的に最適化できない。
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

20

Hartree-Fockの計算結果 [AVX]Time [sec]Total energy [a.u]GAMESS

184.048-2895.7814570171

This work

78.271-2895.7814570169

GAMESS

324.386-2910.6633340322

This work

153.632-2910.6633340179

GAMESS

476.829-3750.9205018138

This work

155.481-3750.9205017267

GAMESS752.839-3770.0595984968

this work

323.098-3770.0595984236

Paclitaxel, 3-21G

Paclitaxel, 6-31G

Valinomycin, 3-21G

Valinomycin, 6-31G

10-7hartree is enough accuracy. We replaced GAMESS ERI with XA PRISM.

FMO(Fragment MO)の手順

Ab initio forlarge Insulators(protein etc.)

Divide a molecule into N fragments
Generate initial density matrices for all fragments (monomers)
Bottleneck
Prepare environmental electrostatic potentials (ESP)
using previous density matrices
NO

Solve Fock equations for all monomers

FICI = SICIeI, for I=1 to N

Are all monomer
densities converged?
YES
Prepare environmental electrostatic potentials for all
fragment pairs (dimers) using monomer densities
Solve Fock equations for all monomers FICI = SICIeI, for I=2 to N, J=1 to I-1
Calculate total energy and properties of molecule

環境静電ポテンシャルの高速化(ESP)
フラグメントA• ERI J-matrixを活用

フラグメント B

J ab   ab cd Dcd

cd: basis on 隣接フラグメント

cd• タンパクをアミノ酸毎の小さい単位に分解
Conventional SCF (ERIはディスク保存) は普通のアミノ酸サイズで有効 (# of basis < 180).
FMOでERIはボトルネックではない
環境静電ポテンシャルが支配的• J-matrix 加速プログラムの活用 by
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

テストモデルインスリン (PDBID:2HIU)
Small Protein
44 amino acids

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

システム条件(1) Intel Core i7-3930K @3.2GHz (6 core) ,GTX580 x 2, 32GB, CUDA 4.0(2) Intel Core i7-3930K @3.2GHz (6 core),
GTX580 x 4, 32GB, CUDA 4.0(3) Intel Core i7-3930K @3.2GHz (6 core) ,GTX680 x 2, 32GB, CUDA 4.0+ Intel Composer XE 12.0 + MKL 10.3

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

インスリンの全体エネルギー計算
Time [sec]Total Energy [a.u.]

Original GAMESS

3279.944-21635.4488652520

Our GPGPU work

790.527-21635.4488649211(1) Intel Core i7-3930K @3.2GHz (6 core) , GTX580 x 2, 32GB, CUDA 4.0(2) Intel Core i7-3930K @3.2GHz (6 core), GTX580 x 4, 32GB, CUDA 4.0(1) + (2)-44フラグメントの全体エネルギーのエラーは小さい
各アミノ酸のエネルギーはほぼ一致している

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

FMO2-ESP and HF-SCF 全体計算
GAMESSThis work

speedup

ESP part(GPU)2571.490

170.897

x 15.0

HF-SCF part(host)708.454619.630x 1.14

Total

3279.944

790.527x 4.15(1) Intel Core i7-3930K @3.2GHz (6 core) , GTX580 x 2, 32GB, CUDA 4.0(2) Intel Core i7-3930K @3.2GHz (6 core), GTX580 x 4, 32GB, CUDA 4.0(1) + (2)

ESPの加速率は全体加速率よりダイブ良い
HF-SCFはオンメモリ計算（GPUによるDirect SCFより速い）
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

アーキテクチャの比較(FMO2 HF-SCF)
Time [sec]Total Energy [a.u.]

Original GAMESS

7026.264-21635.4488652592

Our work [GTX580 x 2](1)1670.931-21635.4488649623

Our work [GTX680 x 2](3)1708.522-21635.4488649862(1) Intel Core i7-3930K @3.2GHz (6 core) , GTX580 x 2, 32GB, CUDA 4.0(3) Intel Core i7-3930K @3.2GHz (6 core) , GTX680 x 2, 32GB, CUDA 4.0

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

FMO2 HF-SCF 6-31G* (d-orbital)
ESP [sec]

Total time [sec] Total energy [a.u.]

Original

19915.745

25360.740 -21643.9995562825

Our GPGPU
work

1138.558

5730.080 -21643.9995557664

speedupx 17.5x 4.4(3) Intel Core i7-3930K @3.2GHz (6 core) , GTX680 x 2, 32GB, CUDA 4.0

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

FMO2 HF-Gradient 6-31G
ESP [sec]

Total time [sec] Total energy [a.u.]

Original

13234.570

20222.516 -21635.4478224023

Our GPGPU
work

1497.2567917.169 -21635.4478217928

speedup

x 8.8

x 2.5

(3) Intel Core i7-3930K @3.2GHz (6 core) , GTX680 x 2, 32GB, CUDA 4.0

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

30

FMO3 HF-SCF 6-31G
ESP [sec]

Total time [sec] Total energy [a.u.]

Original

53898.693

72271.047 -21635.6075667898

Our GPGPU
work

2589.916

19302.927 -21635.6075665173

speedupx 20.8x 3.7(3) Intel Core i7-3930K @3.2GHz (6 core) , GTX680 x 2, 32GB, CUDA 4.0

FMO3は多体効果を考慮できる

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

31

FMO3 HF-SCF 6-31G*
ESP [sec]
Original
Our GPGPU
work

speedup

Total time [sec] Total energy [a.u.]

201450.867

329469.507 -21644.1937717434

8055.506

84456.268 -21644.1937714115

x 25.0x 3.9(3) Intel Core i7-3930K @3.2GHz (6 core) , GTX680 x 2, 32GB, CUDA 4.0

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

32

パフォーマンス・サマリー Results & Discussion
AVX

GPU

Summary of Hartree-Fock x 3 byperformanc Vectorized PRISM
e results
Discussion

DFT
J-matrix formation: x 245
Exchange correlations: x 4
FMO
ESP : x 8.8-20.8

予想がたたない実装 レジスタ不足で単純並列化は難しいコスト短縮ガウス基底を使ったホストとGPUの間の転送コストが低い

GPGPU is better regarding price performance ratio.
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

33

GPGPU FMOのまとめ• GTX680(kepler) はGTX580より少々遅いが、、、– レジスタ不足がボトルネックなので、十分なCUDAコアを生かせてない可能性– GTX580 -> GTX680 clock x 2/3, core x 3, FLOPS x 2• 同時に動かせるブロックが少なく、メモリレイテンシを隠ぺいできていない。コア数が過剰になっている。

– 倍精度が遅い（演算あたりクロック数は上がった）※CUDA 5.0でないとkeplerのポテンシャルは引き出せない（が、現状CUDA 5.0にバグがある）という意見もある。

• ESPの加速はFMOのエネルギーおよびGradのどちらの加速にも効いている
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

エネルギー表示法 (ER) (b) MD/ER• ２桁オーダー高速で同等の精度である自由エネル
ギー計算方法• ermod http://sourceforge.net/p/ermod/wiki/Home/– ナショプロソフト（現時点でGromacs, NAMD, Amberの出力に対応）– 超臨界流体、イオン液体、ミセル、脂質膜に適応可能– タンパク質の全原子自由エネルギー計算も可能– 高精度リガンドバインディングは理論的チャレンジが必要

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

計算フロー-45 -40 -35 -30 -25 -20 -15 -10溶液系MDsolute + solvent

分布関数
生成プログラム

MDプログラムは何でもよい

純溶媒系MDsolvent only

分布関数
生成プログラム

-5

0

5

10

 (kal/mol) e ( )自由エネルギープログラムD 0e ( )  

 0e ( , )
-20

-10

0

 (kcal/m ol)

1020

自由エネルギー計算～計算法の発展と理論の発展

計算

量子計算における
エネルギー

D = ○○ kcal/mol

現実の使用
に耐える
正確さ

近似汎関数
1000

エネルギー分布関数

熱力学積分法（厳密手法）ne̂ e       v , x i 
i

500
100
80
60
40
20
0

~

-5  (kcal/mol)0 (kcal/mol)

容易に
計算可能

・・・

純溶媒



溶液系

純溶媒系と溶液系を結ぶ（仮想的）中間状態を（多数）用意D   de    F  e  

3

計算が難しい（時間がかかる）

MDシミュレーションと溶液理論の結合手法
溶媒和自由エネルギー ()/ (Å mol/kcal)

溶液・界面・生体分子系
統計力学に基づく計算
構造を表す分布関数やダイナミクス
高速で正確な自由エネルギー

MD（典型的に）
分子内・分子間相互作用

動径分布関数とエネルギー分布関数
1000

2

親／疎

~
~
100
80

1

溶液のとらえ方005r (Å)r (Å)
原子間距離

10

原子毎

分子丸ごと

e

n

遠／近

500

3

gO-O
OO
gO-H
OH

 ()/ (Å mol/kcal)
エネルギー分布関数

動径分布関数
g(r)
g(r)

3

60
40

親溶媒

20
0
-5

0

 (kcal/mol)

分子間2体エネルギー  (kcal/mol)

1 g/cm3・25 ºCにおける
水の動径分布関数

1 g/cm3・25 ºCにおける
水のエネルギー分布関数

自由エネルギー
動径分布関数 g(r)exp(-E/kBT)

距離1次元情報の抽出
Kirkwood, 1930年代
本格的発展, 1970年代

エネルギー表示エネルギー1次元情報の抽出

構造・配置の多次元情報（アンサンブル）

エネルギー表示溶液理論の構成
スナップショット配置における2体エネルギー分布（ヒストグラム）̂       v,x i ereference solvent

10ixi : i番目の溶媒分子の座標
v(y,x) : 溶質－溶媒間の2体相互作用
: 2体エネルギー値（分布関数の横軸）

 同じ溶質－溶媒相互作用を持つ配置（構造）をグループ化（射影）する

10
10

分布関数

: 溶質の座標

溶媒和自由エネルギーD   k BT log exp    d ˆ e  
3
2

10

1

~
~

10

2.0

10

1.5
1.0

10

分布

0.5

状態密度

0

10

-5
0
5 105 1010
 どのような近似でも、同じ溶質－溶媒相互作用を2体エネルギー  (kcal/mol)持つ配置（構造）が溶媒和自由エネルギーに同じ重みで効くという統計力学の原理を破らない ⇒ 精度の高い計算

 どのようなポテンシャル関数に対しても1次元の座標 ⇒ 計算スピード向上

 溶質および溶媒分子が、共に、全体として1つのものとして扱われる⇒ 内部自由度 ⇒ 外場 ⇒ 不均一系・QM/MM系⇒ 低密度／濃度極限 ⇒ 超臨界流体・混合溶媒
 DFT（密度汎関数理論）的な理論構成 ⇒ 系統的定式化

-2
-4
-6
-8
-10

アミノ酸アナログに対するテスト計算（１）amino acidanalog soluteamino acidanalog soluteamino acidanalog soluteAla

methaneVal

propaneLeu

iso-butaneIlen-butaneSermethanolThrethanolPhetolueneTyrp-cresolCysmethanethiolMet

methyl ethyl sulfide

AsnacetamideGlnpropionamideTrp

3-methylindoleHid

4-methylimidazoleHie

4-methylimidazole

実験値

free energy D of hydration (kcal/mol)

4

平均誤差

2

現実問題への適用可能性

0.7 kcal/mol

02-

0.7 kcal/mol

46810-

実験値
自由エネルギー摂動法
エネルギー表示法

12Ala Val Leu Ile Phe Ser Thr Tyr Cys Met Asn Gln Hid Hie Trp

力場（ポテンシャ
ル関数）の性能

エネルギー表示法（近似計算、短時間）

0.5 kcal/mol
近似（汎関数）の性能

自由エネルギー摂動法（厳密計算、長時間）temperature: 300K
pressure: 1bar

タンパク質固定構造（配座）における水和自由エネルギーDの収束挙動

平衡ゆらぎにおけるタンパク質構造と水和の関係

block average over 100 ps
cumulative average

-1490
-1500
-1510
-1520
0

0.5

1.0
1.5
Calculation time (ns)2.0solvation free energy D (kcal/mol)

solute: horse heart cytochrome c (1HRC),
104 residues and heme, 1748 atoms
solvent: 20000 water molecules (TIP3P)

水和自由エネルギーD (kcal/mol)

水和自由エネルギーD (kcal/mol)D (kcal/mol)

溶媒をあらわに取扱う、タンパク質の自由エネルギー計算

-1350
-1400
-1450
-1500
-1550
-1600
-1650
-1700

Eintra + D = const

-300
-200
-100
soluteintramolecular energy Eintra (kcal/mol)
タンパク質分子内（構造） タンパク質の水和自由エネルギーのエネルギー Eintra (kcal/mol)
解析が全原子モデルで可能に 10-15 eV におよぶ、タンパク質構造エネルギーのゆらぎが、溶媒水によって、誘起・補償される

「溶質」「溶媒」「溶媒和」概念の拡張～統一的な概念構成に向けて「溶媒」＝溶液系に最初からあるもの「溶質」＝溶液系に後から入ったもの

e通常の意味での溶媒和

脂質膜・ミセルへの分子結合タンパク質への
基質結合

電子の付加
（還元）

溶液理論構成の要件
 内部自由度のある分子
 不均一系
 混合溶媒e通常の意味での溶媒和

脂質膜・ミセル タンパク質への基質結合への分子結合 気体様低密度（1対1の「結合」）から
液体様高密度領域をカバー

 量子論との結合（QM/MM法）
 非物理的挙動を導かないような数学的に堅牢な定式化

電子の付加
（還元）

Matubayasi et al.,JCP 2006, JCP 2008

ミセルとリン脂質二重膜
DMPC (1,2-dimyristoyl-sn-glycero-3SDS (sodium dodecyl sulfate)ミセル
phosphatidylcholine)二重膜
hydrophobic tailhead group
water

0.04

0.04

density (Å )

0.03

-3

-3

density (Å )

0.03

0.02

0.02

0.01

0.01

0

0
0

5

10

15
20
r (Å)

OSO3hydrophobic

headgroup

25

headgroup

glycerol

tail

30

water

0

5

10

15

20

25

30

(Å)
r z(Å)

H2O

regionaqueous

各領域における溶媒和自由エネルギーの算出

I

II

COIII

IV
V
NMe +

VI

3

O
O
CO

H2O
O

PO3-

hydrophobic glycerol headgroupaqueousDMPC二重膜20-2-4methanebenzeneethylbenzene-6-8D (kcal/mol) (kcal/mol)
溶媒和自由エネルギーD (kcal/mol) (kcal/mol)
溶媒和自由エネルギー

SDSミセル
2
0
-2

CO
CO2
benzeneethylbenzene

-4
-6
-8

I

II

hydrophobic

III

IV

OSO3

-V

headgroup

VI

H2O

bulk

aqueousI

IICO IIIIV 3+ V
NMe

O
O
CO

VI

bulk

H2O
O

PO3-

hydrophobic glycerol headgroupaqueous 疎水性領域における顕著な安定化
 疎水性・親水性領域の違いは、ミセルの場合により顕著 分子結合の強度（分配係数）とサイトを決定可能
 疎水性溶質の分布は、ミセル中でより局所的
結合寿命も評価可能
 ミセル・膜内部と水領域での、溶質の安定性の逆転

還元過程の解析－電子１個を溶質と見るー
溶質種：電子
混合溶媒種：電子が付加される分子（イオン）（普通の意味での）溶媒（水など）

溶液内

N

 QM  N 1 H QM N 1  N H QM N

酸化状態

D

e溶液内

N+1
e-

還元状態
Takahashi et al.,CPL 2008, JCP 2008

 MM  vN 1 , x   vN , x exp D   exp    QM    MM, i  
i

 N


電子分布は、溶媒配置によって揺らぐエネルギー座標を構成する相互作用は多体的
(QM)と(MM)（及び関連相関関数）から、Dを構成
FAD (Flavin Adenine Dinucleotide) の活性部位（イソアロキサジン環） の還元への適用
DQM = -39.5 kcal/mol
DMM = -40.6 kcal/mol
大きな水の効果（波動関数部分に匹敵）

ERの欠点• 高速化近似法であるため、FEPのような厳密手法に比べると、適用範囲が現時点では狭い– 現時点でリガンドバインディングはうまくいってない（サンプリングスキームの改善を検討中）• 現状、一番良い結果でR2=0.4程度• ρとρ0が重なっていないと使えない（粗視化モデルとの融合により解決可能？）• 溶媒和自由エネルギーしかやらない• MDの欠点を改善するものではなく、引き継ぐ– MDとの組み合わせによるため、力場問題は存在• 理論がちょっと難しいが、理解してなくても使える
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

47

GUI for (a) + (b)• Winmostar開発（及び販売）を徐々にクロス
アビリティが引き継いでいる– 64bit化は成功– Macintosh対応中– MD GUI対応中• 各種HPCIにジョブ投入– 商用版リリース予定時期は2013年1月以降• その前も順次アカデミックフリー・個人版はリリース
Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

48

Winmostarは計算化学のGUI
分子モデリング
初期座標作成MOPAC6、MOPAC7
CNDO/S UV-VIS計算
内蔵
分子計算ソルバー

GUI

Gaussian、GAMESS、のインターフェイス
解析結果
・分子構造
・電子状態
・諸物性

49

Winmostar

初期構造の作成
AM1 EF PRECISE GNORM=0.05 NOINTER GRAPHF

Winmostar
O 0
H 1.1
H 0.960 0
1 0
1 101.7031

計算ソルバー

0 0
0 0
1 0

0
0
0

0
1
1

0
0
2

0
0
0

分子軌道法計算AM1 EF PRECISE GNORM=0.05 NOINTER GRAPHF

Winmostar

Winmostar
ATOM
NUMBER(I)
1
2
3

CHEMICALSYMBOLO
H
H

計算結果の表示
BOND LENGTH(ANGSTROMS)
NA:I

BOND ANGLE(DEGREES)
NB:NA:I

1.10000 *.96000 *TWIST ANGLE(DEGREES)
NC:NB:NA:I

101.70310 *

NA

NB

1
1

2

NC

CARTESIAN COORDINATES
NO.

ATOM

X

Y

Z

1
O
.0000
.0000
.0000
2
H
1.1000
.0000
.0000
3
H
-.1947.9400.0000
H: (AM1): M.J.S. DEWAR ET AL, J. AM. CHEM. SOC. 107 3902-3909 (1985)
O: (AM1): M.J.S. DEWAR ET AL, J. AM. CHEM. SOC. 107 3902-3909 (1985)50

現時点のWinmostarの特長
１．軽い
２．簡単（分子モデリング）
３．Gaussian、GAMESS,MOPAC対応
４．ソルバーの起動（PC、Linux機）
５．対応が早い（質問、不具合、新機能の要望）

51

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

52

Copyright (C) 2012 X-Ability Co.,Ltd.
All rights reserved.

53

Tencube/WMからLinux機へのジョブ投入

計算ジョブ投入
SSH
計算結果
Linuxクラスター
FOCUS
TSUBAME京？
ShareTask,LSF,PBS,
SGE,NQS等
Gaussian、GAMESS
54

Winmostar

Winscpでデータ転送
Putty(Teraterm)でジョブ投入

ジョブ投入機能vi script.shbsub < script.sh

Winscpで結果ファイル受信
Winmostar

FOCUS外からは
SSL-VPN接続

55

-q s128t24h -n 12 -R "span[ptile=12]"56

(4) Animation[import](5) MO,UV,Charge,Dipole,NMR[import]

57

20年前

現在

58

クラウド
京
TSUBAME
FOCUS

並列クラスター計算機

SSH

59

ご静聴ありがとうございました