計算生命科学 基礎

平成28年11月15日

分子軌道法
北浦和夫(京都大学 福井

神戸大学計算科学教育

基礎 応用理研 AICS)

内容
1.
分子軌道法 基礎
2. FMO/MM融合法
質複合体 構造最適化
FMO/PCM法
結合自由計算3. FMO
処理支援GUI
紹介
4.
1

様大規模系
提案

量子化学計算法 Divide-and-conquer(DC), W. Yang, PRL 66, 1438 (1991)
 Elongation method(EL), A.Imamura et al., JCP 95, 5419(1991)
 Incremental correlation method(IC), H.Stoll, PRB,46, 6700(1992)
 The fragment molecular orbital (FMO) method, Kitaura et.al., CPL(1999)
 Adjustable density matrix assembler (ADMA), T.E.Exner, P.G.Mezey,
JPCA,106,11791,(2002)
 Molecular fractionation with conjugate cap approach (MFCC), D.W.Zhang,
J.Z.H. Zhang, JCP, 119, 3599 (2003).
 Molecular tailoring (MTA), K.Babu, S.R.Gadre, JCC,24, 484 (2003)
 Systematic fragmentation approach(SFA), V.J.Deev, M.A.Collins,
JCP,122,154102 (2005)
 Kernel energy method (KEM), Huang et al., IJQC, 103, 808(2005)
 Electrostatically embedded many-body expansion (EE-MB), E.E.Dahlke,
D.G.Truhlar, JCTC, 3, 46 (2007).
 Generalized X-pol(GXP), L. L. Cembran et al., JCTC, 6, 2469(2009)
For fragment-based approaches, see M.S.Gordon et al., Chem. Rev.
2112,632(2012)Outline of the Fragment MO (FMO) Method2

 A molecule is divided into N fragments
and ab initio MO calculations on thefragments (monomers), fragment pairs (dimers)
and optionally triples (trimers) are performed
under electrostatic potential from othermonomers.

 The total energy of the whole molecule (E) iscalculated using the energies of the monomer(EI), dimer (EIJ) and trimer (EIJK);
FMO2: EFMO2  E I   E IJI4113monomer SCC)
2432-1dimer

3-2I J

FMO3: E FMO3  E FMO2   E IJK  E IJ  E JK  E KI 

4-3

3-1

4-1

4-2I J K

where , E IJ  E IJ  E I  E J 

33Computational procedure of FMO
Fock equation for fragment(monomer) and fragment pair (dimer)(x=I for monomer and x=IJ for dimer)~ x x x x~ x
F C S C ε~x ~ x x
F  H G ,~xxH  H
Vx

Total energy of monomer and dimer
1~ ~Ex  Tr Dx H x F x
2
HF total energy of whole system E  EI  EIJ  EI  EJ I JICorrelation energy correction
Ecorr  EIcorr   EIJcorr  EIcorr  EJcorr

II JBasis of FMO
FMO is a many-body theory of molecular interactions,derived from the energy decomposition analysis(EDA)(IJQC,10, 325 (1976)) of many-body molecular interactions=ΔE12 =A

B++ES

EX+

PL

CT

ES electrostatic

EES  10   02 H 12 10   02  E10  E20

EX exchange repulsion

EEX   10   02 H 12  10   02  10   02 H 12 10   02

PL polarization

EPL  1   2 H12 1   2  10   02 H12 10   02

CT charge transferECT   1CT   CT
H 12  1CT   CT 10   02 H 12 10   0222











HF
HF
00Total interaction energy: E12  E12  E1  E2  EES  EEX  EPL  ECT  EMIXTaking many-body interactions intoenergies of monomer and dimer

E123  E12  E23  E31  E1  E2  E31
 E1  E2  E3  ( E12  E1  E2 )  ( E23  E2  E3 )  ( E31  E3  E1 )Total energy is written as;E   E I   E IJ  E I  E J 
II JMany-body interaction energy in FMO
Total energy
Intenal energy
Total energy interms of internalenergies

NNII JE   EI   EIJ  EI  E J E x  E x  Tr D x V x
NNII 1 J  INNII 1 J  INNNE   EI   EIJ  EI  E J    Tr D IJ V IJ
NI 1 J  I  EI   EIJ  Tr D IJ V IJ

IJIJI
 D
 D
where, Dfor   I ,  J

IJ
IJ
J
D
 D
 D

for   I ,  J


I
DIJ
J VIJ,K K

IJ
IJ
D
 D
for   I ,  J or   J ,  IN NZa
IJKV     D   |         VIJ , K ,r  Ra
 K (  I , J )k ( i , j ) aK
 , K

Pair Interaction Energy PIE)for  ,  IJ

 E IJ   E IJ  Tr  D IJ V IJFMO-MP2 calculations of water clusters

Energy of isolated moleculeEI , 2-body eIJ  EIJ  EI  EJ , and 3-body
000000000eIJK EIJK EI0  EJ0  EK0   eIJ0  eJK
 eKI
interaction energies in the series expansion
0E12N   EI0  eIJ0  eIJK
 ...  e120 NII JI J K

FMO for covalent bonded system(RHF level)
Fock equation for fragment (monomer) and fragment pair (dimer)(x=I for monomer and x=IJ for dimer)~ x x x x~ x
F C S C ε~ ~F x H x Gx ,Hybrid orbital
projector~xx
H  HVx B  ih ih  ,bond-detached atom(BDA)iElectrostatic potentialZAKV      D( )r  RAKx AKKx

Total energy of monomer/dimer 1x ~ x ~xEx  Tr D H  F2 (x  I, IJ)

Total energy of whole systemE  EI  EIJ  EI  EJ II J9Examples of fragmentation in FMO
Polypeptide
DNA/RNAS-S結合s 炭素 位置 切断 他 分子 場合 1重結合 切断
少小誤差 大最低2030原子 含 程度
安心Intra- and inter-fragment pair interaction
enegy (PIE) is obtained with FMO
NN

N

N

N

N

I

I 1 J  I

II 1 J  I

NNII 1 J  I

NNII 1 J  IN

NE   EI   EIJ  EI  EJ    EI   EIJ  EI  EJ    Tr DIJ V IJN  EI   EIJ  Tr DIJ V IJ

N  EI   EIJI 1 J  I

PIE this term can be further
decomposed into EL, EX, CT
interaction energies by PIEDA

Intramolecular interaction:
EIJ indicates interactionenergy between residues except
covalent bonded ones.
Intermolecular interaction
EIL indicates interactionenergy between a residure andligand.

O R HI O R HJ O
HNNNNNHRO RI H O RJ
EIJI OJ OO R HRHHNNNNNR H O RI H O RJ
EILLigandLEJL
11

Example of PIE analysis: IntermolecularInteraction betweenFKBP and its Ligands
Nakanishi et al., Proteins: Struct., Funct., Bioinf. 68, 145 (2007)• FK506 is an immunosuppressant.
• We calculated the following four FKBP-ligand complexes by FMO.
• The complexes have 107 amino acid residues (about 1,800 atoms).
• The ligand geometries were optimized at FMO2-RHF/3-21G level using
truncated 20 residues model complexes (about 500 atoms).
• The binding energies were calculated at the FMO2-MP2/6-31G* level.
FK506

FKBP

1fkb

1fkg
FKBP-FK506 complex (PDB:1fkf

1fkf (FK506)

1fki12

PIE between ligand and each residue
in FKBP-ligand complexes•The sum of pair interactionenergies correlates wellwith the experimental binding
affinity, 1fkb>1fkf>1fkb>1fki.
•Val55, Tyr82 and Try26 arecommon important residues
for the all ligand bindings.
•The stronger binders haveadditional interactions with
Asp37 (1fkb and 1fkf) andGlu54(fkb).

a) Empty bar:HF, filled bar: correlation energy
contribution. b)Ligand binding modes. Theproteins is shown by surface model.

•The correlation contribution
is very large: 70-80% ofbinding energy.
13

Availabel free FMO Programs
1) FMO in GAMESS, coded by D.G.Fedorov et al.
http://www.msg.ameslab.gov/GAMESS/GAMESS.html
2) ABINIT-MP, coded by T.Nakano et al.
http://www.fsis.iis.u-tokyo.ac.jp/en/result/software/3) PAICS, coded by T.Ishikawa
http://www.paics.net/get_paics.html

GUI programs for FMO pre- and post-processing
Facio for GAMESS: http://www1.bbiq.jp/zzzfelis/Facio_Jp.html
FU for GAMESS: https://sourceforge.net/projects/fusuite/
BioStation for ABINIT-MP:
http://www.advancesoft.jp/product/advance_biostation/
PAICS View: http://www.paics.net/get_paics.html
14いろんな電子状態理論のFMO法
FMO2-RHF : original FMO (2-body expansion), Kitaura et al., CPL,313,701(1999)
FMO3-RHF : 3-body expansion FMO, Fedorov et al., JCP,120, 6832 (2004)
FMO2,3-UHF: 開殻系 非制限Hartree-Fock法, Nakata et al., JCP,137,044110 (2012)
FMO2,3-ROHF: 開殻系 制限Hartree-Fock法, S.R.Pruitt et al., JCTC,6,1(2010)
FMO2,3-DFT: FMO-based density functional theory, Fedorov et al., CPL,389, 29 (2004).
FMO2,3-MP2 : FMO-based 2nd order Møller-Plesset perturbation theory, Fedorov etal., JCP, 121, 2483 (2004), Mochizuki et al.,CPL,396, 473(2004)
FMO2-MCSCF: FMO-based MCSCF, Fedorov et al.,JCP,122,54108(2005).
FMO2,3-CC: FMO-based coupled cluster theory, Fedorov et al., JCP, 123, 134103 (2005)
FMO1-CIS and CIS(D) : FMO-based configuration interaction singles, Mochizuki et al.,
CPL, 406, 283 (2005).
MFMO : FMO-based multilayer method, Fedorov et al., JPCA,109,2638 (2005).
FMO1,2-TDDFT : FMO-based time dependent DFT, Chiba et al., JCP, 127, 104108 (2007).
FMO/PCM : FMO combined with polarizable continuum model (PCM), Fedorov et al.,
JCC, 27, 976 (2006).
FMO/EFP: FMO combined with effective fragment potential method, Nagata et al.,
JCP, 131, 024101 (2009).

Capabilities of FMO in GAMESS
http://www.msg.chem.iastate.edu/gamess/capabilities.html

Recently,
Nishimoto et. al.
have developed
FMO-DFTB(JPCL 6 (2015)5034) and
implemented its
codes intoGAMESS.

1616

Useful FMO methods in drug design
分子設計結設計化合物 活性推定
水溶液中 解離OXX

NH2O

NH2O

NH 2

数万 数十万化合物
活性化合物 探索設計化合物
合構造推定

活性予測

XNH2O化合物
質

Fast geometry
optimization
calculations
Frozendomain FMO(FMO/FD)
Fedorov et al.,JPCL, 2(2011).

Computations andanalysis of molecularinteraction energy

Geometry opts
Analytic gradients

PIE/PIEDA analysis
QM:QM hybrid
method for accurate
calculations(ONIOM+FMO)

MQ/MM hybrid(FMO/MM)Asada et al., J. Phys. Chem.
Lett.,3 (2012)

Fedorov et al.,AccChem. Res(2014)

Solvation energy
PCM model(FMO/PCM)
EFP,PCM model(FMO/EFP/PCM)
Binding energy
analysis(FMO/PCM/PIEDA) Fedorov et al.,Geometry opts insolution(FMO/PCMgradients) Nagata et al.,17 JPCA, 116 (2012).
(2012)

溶媒溶媒分子考慮

溶媒 誘電体 近似
誘電率 

系全体 MD
溶質 溶媒分子 相互作用
平均値 求
熱力学積分法 TI)
計算 可能
長時間

行

自由行 必要物質誘電体 空孔) 溶質分子
質 溶媒 相互作用
分散
計算
静電他

誘電体理論
経験式 評価PBSA GBSA

PCM

溶静電Polarizable Continuum

PCM

J. Tomasi, B. Mennucci, R. Cammi, Chem. Rev., 105 (2005) 2999.

溶質 電子分布 分子軌道法 求consistent 解 分子軌道法計算 最用self-表面 誘起
溶媒溶質 溶媒相互作用(ele + rep + disp )溶媒和自由ele
rep
disp
cav
Gsol   Gin  Gsol
 Gsol Gsol Gsol G in : 溶質分子 内部変化 気相差

ele
G sol
: 溶質分子 溶媒 静電相互作用

項

rep
G sol
: 溶質分子 溶媒 交換反発
disp
G sol
: 溶質分子 溶媒 分散
cav
G sol
:

項
項

項

空孔(cav)

PCM
溶質解

静電電子分布空孔表面 誘起溶質分子 溶媒静電ele
G項
sol

電荷 self-consistent
求空孔表面 誘起電荷 qq  1 1
C V4

:溶媒 誘電率
V:表面静電C:構造行列
誘起電荷 影響下 溶質分子 電子状態 解               F  F0  W
N TSi
W   qi w
i 1iw 1
r  Ri

表面上 誘起電荷 q

PCM

項 Gcav , G rep , Gdisp

静電以外
G cav
(式Gcav(Rk)

,交換反発

分散

G rep

G disp分子間相互作用 electrostatic, exchangerepulsion, dispersion, charge-transfer)

FMO/PCM法 FMO is combined with polarizable continuum model (PCM) Solute electrostatic field V oncavity surface is calculated frommany-body series expansion;NNVi   Vi   Vi IJ  Vi I  Vi J  I 1II JZ X R   R i
 ,   x

Vi  Tr D  w i  xiw x1r  RiSeveral approximation levels are
possible depending on V and electron
density (Dx) expansion:(1) FMO/PCM[1] uses one-bodyexpansion V and one-body density(2) FMO[2] uses two-body V and twobody density(3) FMO[1(2)] uses two-body V and twobody density (no self-consistent)

FMO/PCM[1(2)] is recommended

flowchart of FMO/PCM[2]22FMO/MM: FMO-based IMOMM
D.G.Fedorov, et al., J. Phys. Chem. A 2007, 111, 2722.
IMOMM: Maseras, F.; Morokuma, K. JCC 1995,16,1170.
SIMOMM: Shoemaker J. R.; Burggraf, L. W.; Gordon, M. S. JPCA 1999,103, 3245.
The total energy E and its gradient are divided into FMO and MM contributions,E  EFMO (R1 )  EMM (R1 , R 2 ) E FMOE E MM

R 1
R 1
R 1
E
 E MMR 2
R 2

QM領域
MM領域

where R1 and R2 are a set of atomic coordinates in the FMO region andatomic coordinates in the MM region, respectively. Note that the interactionenergy between FMO and MM atoms are included in and the link atoms inFMO region do not feel MM forces.
23Geometry optimization of protein by FMO/MMD.G.Fedorov, et al., Acc. Chem. Res, Acc. , 2014, 47, 2846.
浅田直也 博士論文 京都大学 2012

Computational Details:
1. System: alpha subunit of protein kinase 2 (CK2)
and its ligand ((1-{6-[6-(cyclopentylamino)-1Hindazol-1-YL]pyrazin-2-YL}-1H-pyrrol-3-YL)aceticacid)
2. FMO2 at the RHF-D/6-31G level (D:empiricdispersion) and AMBER f99 (protein) and gaff(ligand) force fields and TIP3P (water) were
employed for MM.
3.The FMO region was defined to include the ligand
and the amino acid residues and water molecules
separated by the unitless FMO distance70 of 2.0
from the ligand.
4. We optimized all FMO atoms (667) and 1980 ofMM atoms.
5. A single point calculation took about 36 min on
the Heian cluster (112 CPU cores of Xeon 3.0 GHz)and 470 FMO steps were required for convergence.

FMO atoms are colored inpurple (ligand) and blue (thebinding pocket of the protein).
MM atoms are shown as
yellow (optimized) and green(frozen) for the protein, andcyan (optimized) and pink24(frozen) for water.

Intermolecular interaction between theprotein and ligand
The protein recognizes its ligand with various kind of nonbonding interactions at multiple interaction sites. Theseinteractions should be described in a well balanced mannar.
/CH/

NH-N
salt bridge

CH-O

S/

Dispersion interaction is indispensable toobtain reasonable binding structure♦FMO-HF/6-31G failed to reproduce the binding structure (RMSD:0.75Å)♦ FMO-MP2 reproduced experimental structure with RMSD of 0.61Å
♦ Convenient HF-D also reproduced the binding structure (RMSD:0.56Å)HF-D: RHF augmented with empirical dispersion force (S. Grimme, JCC, 25,1463(2004) ).

,,protein kinase 2CK2 complex
MP2/6-31G (0.61Å)LigandHF/6-31G(0.75Å)HF-D/6-31G (0.56Å)X-ray (cyan)

In parentheses are given RMSD between
26
optimized and experimental coordinates.

Optimized Structure of Ligand Binding Site

The left and right panels show the top and side views of the binding site, respectively.
For each distance between heavy atoms (Å), three numbers are given, corresponding
to the FMO/MM minimum, MM minimum and experimental values, respectively.

As shown in the figure, the FMO/MM and MM structures show in general
somewhat similar deviations from experiment. However, the NH...N andCH...O hydrogen bonds between the pyrazine moiety of the ligand and VAL116
are not well described by the force field. In this particular complex, multiple weakinteractions are reasonably described by both FMO/MM and MM methods.
27

Calculation of binding free energy Target system is the complexes of protein kinase 2, CK2
 CK2 (-subunit) is a potential target for nephritis therapy.
 Variety of CK2 ligands were selected for binding free energy
calculations
polar ligand

charged ligand

NO
N

N

N
N
HN

NN

N

N

H
N

N
N
H

HN

NN
HO

C

N-OOO
HNOO

B

bromo-OE

tetraquinone

BrBrNNO

ABinding siteNBrHONOOOHOOH

NNBrBrOBr

Structure of CK2 (subunitDFAffinity: A > B > C > D > E > F >GG28

Binding free energy in solution was calculated withFMO/PCM (polarizable continuum solvent model)
Binding energy of protein-ligand complex (PL) in solutioncomplex

Gbsol+proteinligandsol
Gbsol  GPL
 GPsol  GLsol

GXsol  GXinternal  GXele  GXcav  GXdisp  GXrep (X  PL, P, L) 1) FMO-MP2/PCM/6-31G*2) conductor-like PCM in GAMESS
3) cavity was created with the simplified united atomicradii, H:0.01Å, C:1.77Å, N:1.68Å, O:1.59Å, and S:2.10Å
4) 60 tesserae per sphere.
29

Entropy term:quasi-harmonic approximation
S : quasi-harmonic analysis

Levy, et. al., Macromolecules, 1984, 17, 1370-1374

covariance matrix:  ij 1

 x  x x  x ii

j1

σ '  M σM 22jquasi-mode frequenceis:   k BT  ' 2
1

QH
QHQH
 S  S PL S PL(P) S PL(L)MD simulation for entropy term30- Initial structures and restraint conditions were the same as used for thestructure refinement.
- 10 ns MD simulation : NPT ensemble (1 atm, 300K) with 1000 kcal/mol
restraint to fixed domain (production run: last 5ns)
- ff99, gaff and RESP charges were used.
30- AMBER9 programCalculated binding free energy vs ExperimentGb (calc, kcal/mol)-15

tetraquinone

-20
-25

R = 0.963-30

polar ligand-40

A±5kcal/mol

bromoC

-35

-45

DGF
Echarged ligand

B-50

Gb (expt, kcal/mol)
The calculated values are well correlated with experiments
including charged ligands
31

FMO 関包括的情報

入手先

https://staff.aist.go.jp/d.g.fedorov/
The FMO method|__ Download FMO resources (oral presentation etc)|__ Basics of FMO(
分子軌道法 基礎|__ FMO implementation in GAMESS(GAMESS FMO使用 手引

過去 開催

FMO講習会資料

http://www.cms-initiative.jp/ja/researchsupport/about-kobe-support
CMSI神戸| CMSI神戸(|__ 過去 講習会一覧 資料|__ 第23回: FMO

講習会)

FU: Open source GUI software
for GAMESS-FMO calculations

https://sourceforge.net/projects/fusuite/
Structure modelingsGAMESS Input data
generation andcomputations

Result visualisations

FU 用FMO

実習

FMO計算

実際(平成27年5月25日)

抜粋内容
1. FMO-RHF/STO-3G 1点計算
2. FMO-DFT/STO-3G点計算
3. FMO-MP2/STO-3G点計算
4. FMO-RHF/STO-3G
構造最適化計算
5. FMO-RHF-D/STO-3G 構造最適化計算
6. RHF/STO-3G(ab initio MO 1点計算
7. FMO-RHF/PCM/STO-3G計算
8. 水和
作成
FMO計算

FMO-RHF/STO-3G 1点計算...c2h5oh-h2o1) fu.exe ckick
fumodel 起動2) PDB形式 座標(名 "c2h5h2o.pdb")
読 込"AddDel/(Add)"-"Bond""Use bond length" 結合付加(fumodel
背景色 default 黒本稿白変更分子軌道法
大規模並列計算 可能
実際 京
20万用
約2万原子質複合体 FMORIMP2/6-31G*1点計算 11分 計算
力場
融合法 FMO/MM法
現実的 時間 実行
PCM溶媒計算用FMO/PCM)

構造最適化計算
溶媒和自由PIE PIEDA)
分子内 分子間単位
相互作用解析 行
巨大 複雑分子系 構造関
詳細 情報 得質 分子
認識機構 理解有用
FMO計算 容易開GUI開発公