第8回 配信講義 計算科学技術特論B(2024)
2K Views
May 30, 24
スライド概要
第１回 ６月６日 ABINIT-MPプログラムによるフラグメント分子軌道(FMO)計算１
フラグメント分子軌道(FMO)法の概要を基本的な式と処理のフローを交えて解説します。具体的には、ハートリーフォック(HF)、2次および3次のメラープレセット摂動論(MP2,MP3)を取り上げ、スパコンを使ったベンチマーク計算の事例をご紹介します。また、最近のGPU対応についても触れます。
R-CCS 計算科学研究推進室
関連スライド
各ページのテキスト
配信講義 計算科学技術特論B（2024） - 2024/6/6 - #1 ABINIT-MPプログラムによる フラグメント分子軌道(FMO)計算１ 望月祐志 （立教大学理学部化学科） [email protected] 2024/5/29 1
謝辞 ◇ABINIT-MPプログラムのコード/ツールの主な開発関係者（敬称略） 望月祐志、中野達也（RIST）、坂倉耕太（FOCUS）、加藤季広（NEC）、 佐藤伸哉&山本純一（NES）、石川岳志（鹿児島大）、沖山佳生（神戸大）、 山下勝美（元NES）、奥脇弘次（JSOL/立教大）、土居英男（立教大）、 渡邊啓正（HPCS）、大島聡史（九州大）、片桐孝洋（名古屋大） ◇研究開発支援 CISS/HPCIの4プロジェクト； 東大生研&文科省 / FY2002-2015 FS2020(ポスト「京」)プロジェクト； 東大&文科省 / FY2014-2019 CRESTプロジェクト(“田中FMO”)； JST / FY2004-2009 科研費(“榊-特定領域”)； 文科省 / FY2008-2009 科研費(基盤B：代表)； 文科省 / FY2016-2018 SFR； 立教大 / FY2006-2007, 2010-2014, 2019-2020, 2022NEC様； SX-AT高速化&機能強化の共同開発 / FY2020その他企業様からの立教大宛のご寄付 JHPCN課題； jh210036-NAH, jh220010, jh230001, jh240001 「富岳」課題； hp210026, hp210261, hp220025, hp220352, hp230016, 2 hp230017, hp230375, hp240013, hp240030 2024/5/29
内容と流れ ・第一回（今回） フラグメント分子軌道(FMO)法の概要を基本的な式と処理の フローを交えて解説します。具体的には、ハートリーフォック(HF)、 2次および3次のメラープレセット摂動論(MP2,MP3)を取り上げ、 スパコンを使ったベンチマーク計算の事例も紹介します。また、 進行中の話題として積分計算のGPU対応についても触れます。 ・第二回（次回） 先ず、テンソル縮約処理が支配的な高次相関計算の扱いを紹介 します。次に、2020年度の試行的利用段階の「富岳」を使った 大規模計算の事例を新型コロナウイルスの関連タンパク質を例 にお示しします。後半は、高速化と超大規模系対応のプログラム 改修、「富岳」を使った応用計算やデータ解析例などの最近の トピック、粗視化シミュレーションへの接続についてお話します。 2024/5/29 3
https://www.springer.com/gp/book/9789811592348 / ABINIT-MPはChap. 4 - Y. Mochizuki et al., “The ABINIT-MP Program”に記載 pp. 53-67. FMO計算の発展をまとめた本 2024/5/29 2021年1月刊行 4
ABINIT-MPの全体的開発に関する論文リスト （赤字: HPC分野の先生方 / 紫字: HPC分野の技術者方） ■英文 (1) "Electron-correlated fragment-molecular-orbital calculations for biomolecular and nano systems", S. Tanaka*, Y. Mochizuki*, Y. Komeiji, Y. Okiyama, K. Fukuzawa, Phys. Chem. Chem. Phys., 16 (2014) 10310-10344. (2) "The ABINIT-MP Program", Y. Mochizuki*, T. Nakano, K. Sakakura, Y. Okiyama, H. Watanabe, K. Kato, Y. Akinaga, S. Sato, J. Yamamoto, K. Yamashita, T. Murase, T. Ishikawa, Y. Komeiji, Y. Kato, N. Watanabe, T. Tsukamoto, H. Mori, K. Okuwaki, S. Tanaka, A. Kato, C. Watanabe, K. Fukuzawa (pp. 53-67) in Recent Advances of the Fragment Molecular Orbital Method - Enhanced Performance and Applicability, ed. Y. Mochizuki, S. Tanaka, K. Fukuzawa (January 2021, Springer). ■邦文 (1) "FMOプログラムABINIT-MPの開発状況と機械学習との連携", 望月祐志*, 坂倉耕太, 秋永宜伸, 加藤幸一郎, 渡邊啓正, 沖山佳生, 中野達也, 古明地勇人,奥沢明, 福澤薫, 田中成典, J. Comp. Chem. Jpn., 16 (2017) 119-122. (2) "FMOプログラムABINIT-MPのOakForest-PACS上での多層並列化と性能評価", 渡邊啓正*, 佐藤伸哉, 坂倉耕太, 齊藤天菜, 望月祐志, J. Comp. Chem. Jpn. 17 (2018) 147-149. (3) "ABINIT-MP Openシリーズの最新の開発状況について", 望月祐志*, 秋永宜伸, 坂倉耕太, 渡邊啓正, 加藤幸一郎, 渡辺尚貴, 奥脇弘次, 中野達也, 福澤薫, J. Comp. Chem. Jpn., 18 (2019) 129-131. (4) "FMOプログラムABINIT-MPの整備状況2020", 望月祐志*, 坂倉耕太, 渡邊啓正, 奥脇弘次, 加藤幸一郎, 渡辺尚貴, 沖山佳生, 福澤薫, 中野達也, J. Comp. Chem. Jpn., 19 (2020) 142-145. (5) "FMOプログラムABINIT-MPの整備状況2021", 望月祐志*, 中野達也, 佐藤伸哉, 坂倉耕太, 渡邊啓正, 奥脇弘次, 大島聡史, 片桐孝洋, J. Comp. Chem. Jpn., 20 (2021) 132-136. (6) "FMOプログラムABINIT-MPの整備状況2022", 望月祐志*, 中野達也, 坂倉耕太, 渡邊啓正, 佐藤伸哉, 奥脇弘次, 秋澤和輝, 土居英男, 大島聡史，片桐孝洋, J. Comp. Chem. Jpn., 21 (2022) 106-110. (7) "FMOプログラムABINIT-MPの整備状況2023", 望月祐志*, 中野達也, 坂倉耕太, 奥脇弘次, 土居英男, 加藤季広, 滝沢寛之, 成瀬彰, 大島聡史, 星野哲也, 片桐孝洋, J. Comp. Chem. Jpn., 23 (2024) 4-8. 注記： (5)以降が、HPC分野の先生方とのコラボレーションによるVer. 2系の報告 2024/5/29 5
FMO法とABINIT-MPプログラム 2024/5/29 6
FMO-DFT Ref.; Y. Nishimoto et al., Phys. Chem. Chem. Phys. 18 (2016) 22047. 分子軌道（MO：量子化学計算）法の分類 ◇Hartree-Fock（HF） 単一の電子配置、平均場近似によって分子軌道を最適化 ψi＝ΣpCpiφp （基底関数による展開） → ｛Cpi｝ N4コスト ⇒ 準定量性、大規模系への適用性 ◇Post-HF 多数の電子配置による展開、電子相関の導入 Ψ＝ΣITIΦI → ｛TI｝ （展開長は～1億） N5からN7コスト 摂動論（PT）、配置間相互作用（CI）、結合クラスター（CC） ⇒ HFに対し定量性・信頼性を向上、励起状態の記述 ◇密度汎関数理論（DFT） 単一配置のまま電子相関を実効的に導入、コストはHFと同様 ⇒ 多種多様な汎関数(B3LYP等)、固体系・金属系に強み ⇒ FMO法では（DFTBを除けば）ほとんど利用されない 2024/5/29 7
「分子理論の展開」永瀬・平尾 共著（岩波書店）. / 「すぐできる量子化学計算」平尾 監修・武次 編著（講談社）. 分子軌道計算から得られる諸量 ◇プログラム（MO系） GAUSSIAN, MOLPRO, MOLCAS, GAMESS-US, NWCHEM・・・ ⇒ 各々に特徴あり、ただ実際にはGAUSSIANが寡占 ◇構造とエネルギー 安定構造や遷移状態の最適化、振動数計算 ⇒ 反応、ポテンシャルエネルギー面 ⇒ 自由エネルギー変化 ⇒ 分子間の相互作用エネルギー ◇物性値 【H のエネルギー面】 軌道(HOMO, LUMOなど)の分布 電荷（自然密度、静電ポテンシャル系）、多重極子能率 励起/発光エネルギー、(超)分極率、NMR、EPR・・・ 3 2024/5/29 8
FMO Original; K. Kitaura et al., Chem. Phys Lett. 313 (1999) 701. / General Ref.; M. Gordon et al., Chem. Rev. 112 (2012) 632. フラグメント分子軌道（FMO）法 ◇巨大分子系 生体高分子や凝集系では一般的 ⇒ タンパク質、DNA （水和状態） 数千～数万原子、数千～数十万軌道 【HIVプロテアーゼとロピナビル】 ◇分割&統合系のアプローチの一つ 北浦らが1999年に2体展開で提案 ⇒ フラグメントとその対で系のエネルギーを評価 （FMO2） ⇒ 環境静電ポテンシャル （ESP）、直接結合切断 （BDA） ⇒ 階層的な並列処理 （フラグメントリスト&内部処理） ⇒ 定量性を高める電子相関の導入も直截 フラグメント間の相互作用エネルギー （IFIE） ⇒ 計算対象の解析ツール J I I ⇒ 生物物理や創薬に向く ⇒ 材料系にも適用可能 2024/5/29 モノマー （アミノ酸単位など） ダイマー 9
S. Amari et al., J. Chem. Inf. Model. 46 (2006) 221. / GAMESS-USではPIEと呼ばれる. フラグメント間相互作用エネルギー（IFIE） FMO-HF法の全エネルギーEと全電子密度rの計算式 E EIJ ( N f 2) EI I J I r (r ) r IJ (r ) ( N f 2) r I (r ) I J I Nf: フラグメント（モノマー）の個数 EI: モノマーのエネルギー, EIJ: ダイマーのエネルギー rI: モノマーの電子密度, rIJ: ダイマーの電子密度 FMO-HF法の全エネルギーEを書き換えると、全エネルギーをフ ~ ラグメント間相互作用エネルギーDEIJ、モノマーのエネルギーEI から周囲のフラグメントとの静電相互作用エネルギーを除いた エネルギーE’Iの和で表せる。 電子相関エネルギーの補正は加成的に行えばよい。 ~ E DE IJ E ' I I J 2024/5/29 I IFIE解析 ⇒ 対象系内部の詳細情報 (Inter-Fragment Interaction Energy) 10
FMO計算のためのプログラム ◇GAMESS-US [米国Gordonグループ]； Fedorov、Gordon、北浦ら GAUSSIANに抗し得る有力なフリーソフト、世界規模 (Fortran) ⇒ 様々な機能をFMO化、多彩な計算、GDDI並列 ◇ABINIT-MP； 望月、中野ら 実用機能は十分、東大系PJ・CREST-PJなどで開発 (Fortran) ⇒ IOレス、MPI、OpenMP/MPI混成並列、スパコンと好相性 2階層の並列処理で計算資源を有効利用 ◇PAICS； 石川 FMO-MP2(RI)に特化 (C) Target system ⇒ MPI並列 Fragment momer Parallelized by fragment indices ◇OpenFMO； 稲富、鬼頭ら FMO-HF (C) ⇒ 超並列指向 Parallelized by integral 2024/5/29 indices (or dimer) Processor group 11
Ref.; “The ABINIT-MP program”, Y. Mochizuki et al. in “Recent Advances of the Fragment Molecular Method” (Springer, 2021). ABINIT-MPの主な機能 （オープンシリーズとして整備中） ・エネルギー → FMO4： HF, MP2 → FMO2： HF～CCSD(T), LRD → FMO2: CIS/CIS(D) ・エネルギー微分 → FMO4： HF, MP2 → FMO2： MP2構造最適化, MD ・その他機能 → SCIFIE, PB, sp2-BDA, α(ω) → 電子密度生成, CAFI, FILM ・並列化環境（PC~スパコン） → MPI, OpenMP/MPI混成 → 最深部はBLAS処理 2024/5/29 http://www.cenav.org/abinit-mp-open_ver-2-rev-4/ Open Ver. 1 （ポスト「京」のPJで整備） ・Rev. 5 (2016年12月) ・Rev. 10 (2018年2月) ・Rev. 15 (2019年3月) ・Rev. 22 (2020年6月)；当面は併存 Open Ver. 2 （「富岳」の時代に移行） ・Rev. 4 (2021年9月) ・Rev. 8 (2023年8月) ・Rev. 12 (2024年12月予定) （注記： Ver. 2系では、BioStation Viewer へのデータファイルの書出しを廃止した） 12
BioStation Viewerへの接続はCPFを書き出すABINIT-MPのVer. 1系まで対応. FMO計算のワークフロー （タンパク質の例） ① PDBからタンパク質構造をダウンロード 利用ソフトウェア ② 構造の補正（水素付加・力場による最適化など） 分子設計支援 ソフト等 ③ FMO計算の入力ファイルajfの作成 BioStation Viewer 入力ファイルを計算サーバにアップロード ④ FMO計算の実行（並列計算） ABINIT-MP 出力ファイルをクライアントにダウンロード ⑤ 計算結果（log、cpf ファイル）の解析 BioStation Viewer Ala3の分割例 2024/5/29 13
IFIEの可視化例 （BioStation Viewer使用） 女性ホルモン受容体と女性ホルモン間のIFIE表示 IFIE表示設定画面 IJ-PAIR DIST DIMER-ES HF-IFIE MP2-IFIE / A APPROX. / Hartree / Hartree ------------------------------------------------------2 1 0.000000 F -15.127711 -0.043867 3 1 3.533609 F 0.011185 -0.001489 3 2 0.000000 F -15.089418 -0.050422 2024/5/29 14
PIEDA Ref.; D. G. Fedorov et al., J. Comp. Chem. 28 (2007) 222. / 塚本ら, J. Comp. Chem. Jpn. 14 (2015) 1. 相互作用の成分分析 （PIEDA) 女性ホルモン受容体のファーマコフォアの解析 hydrophobic residues Phe404 Glu353 O HO C HN N O His524 H2O H 2N hydrophobic residues N H2 Arg394 C O NH ES dominates 1st and 2 nd ⇔ Due to hydrogen-bond EST HO C NH Leu387 IFIE Lig. Chrg. -107.74 -0.13 ES EX CT+mix DI -73.47 55.49 -29.61 -60.16 ・ 静電相互作用と分散相互作用がオーバーオールの安定化に寄与 ・ 電荷移動も水素結合の安定化に寄与 （Glu353など） 2024/5/29 15
Ref.; T. Nakano et al., Chem. Phys. Lett. 528 (2012) 128. FMO4展開 （ABINIT-MPのみの機能） E FMO2 EIJ ( N 2) EI DEIJ EI I J I I J E FMO4 I J K L I E ( N 2)( N 3) EIJK ( N 3) EIJ EI 2 I J K I J I DEIJK DEIJ DEIK DE JK I J K DEIJ EI I J I J K ( N 3)( N 4) EIJ 2 I J ( N 2)( N 3)( N 4) EI 6 I DEIJKL DEIJ DEIK DEIL DE JK DE JL DEKL DE IJ E IJ E I E J FMO3 EIJKL ( N 4) EIJK I J K L DEIJK DEIJ DEIK DE JK DEIJL DEIJ DEIL DE JL DEIKL DEIK DEIL DEKL DE JKL DE JK DE JL DEKL DEIJK DEIJ DEIK DE JK I J K DEIJ EI I I J DE IJK E IJK E I E J E K I DE IJKL E IJKL E I E J E K E L 3体、4体の計算は関連のモノマーが全て最近接となっているものだけを対象とする J J I I J I K Monomer 2024/5/29 Dimer L I Trimer K Tetramer 16
FMO4法の精度 （通常のMO計算に匹敵） Reference (au) FMO2 (mili-au) FMO3 (mili-au) FMO4 (mili-au) HF/6-31G* (H2O)6 (H2O)16 (Gly)5 -456.142593 -1216.421080 -1110.092757 -2.598 -10.529 -3.009 0.243 1.634 -0.060 -0.029 0.067 0.019 (Ala)5a -1305.282976 -3.510 -0.043 0.018 (Ala)5b -1305.282976 -3799.528980 880.068 -3.452 0.023 18.538 0.608 0.033 -3799.528980 1564.477 -4.895 0.019 -1.137183 -3.047880 -3.120818 -0.365 0.605 -0.284 0.053 0.868 -0.023 -0.003 0.511 0.055 (Ala)5a -3.797451 0.128 -0.096 0.036 (Ala)5b -3.797451 29.202 -0.496 0.173 a Chignolin -11.194920 3.354 0.182 0.521 Chignolinb -11.194920 48.869 -1.418 1.028 Chignolina Chignolinb MP2/6-31G* (H2O)6 (H2O)16 (Gly)5 2024/5/29 FMO4は誤差が少ない Chignolin a: これまでの主鎖による分割 b: 主鎖/側鎖による分割 ・FMO4で主鎖/側鎖分割が可能 （リガンド側も同様に小分割） ・水素結合系の精度も向上 17
GAMESS-USでは固体系はAFOによる扱いになる. FMO4法の精度 （固体の扱い可能） Reference (au) FMO2 (mili-au) FMO3 FMO4 (mili-au) (mili-au) HF/6-31G C22 H 28 -849.139220 618.483 37.943 -1.075 C35 H 36 -1345.957600 1980.173 94.602 0.246 Si22 H 28 HF/6-31G* C22 H 28 -6371.923912 203.289 77.192 3.668 -849.493710 585.694 41.553 -1.910 C35 H 36 -1346.519492 2042.928 111.708 0.633 Si22 H 28 MP2/6-31G C22 H 28 -6372.670343 150.599 31.530 0.478 -1.995930 7.573 5.658 0.589 C35 H 36 -3.182338 12.314 1.746 2.606 Si22 H 28 MP2/6-31G* C22 H 28 -1.201304 -44.628 40.083 8.050 -2.894793 6.800 5.440 0.510 C35 H 36 -4.617236 18.599 6.582 2.376 Si22 H 28 -1.815212 -11.393 17.356 2.310 2024/5/29 a: tetramantane C22H28 b: superadamantane C35H36 ・FMO2はエラーが多くて不可 ・FMO3でも十分とは言えず ・FMO4なら精度保持が可能 ・ナノ/バイオ境界の問題へ道 18
リリース状況のまとめ（2023年秋時点） ◇Ver. 1系の最終リリース版はRev. 22 （2020年６月） f関数のサポート （直交10成分のみ） MFMOによる指定領域のみでのMP2&MP3計算 加速されたPIEDA （「富岳」の新型コロナPJの中で改良） ◇「富岳」時代でVer. 1からVer. 2へ移行 - Open Ver. 2 Rev. 4 （2021年9月） 単一構造サンプルでの可視化的解析の「限界」 統計的な相互作用解析の重要性 （MDベースで百～千サンプル構造を対象） 大規模系の扱いの必要性 （タンパク質の液滴モデルは数万フラグメントに） 可視化ツールによる解析を「諦め」 （ダンプ用データがメモリを圧迫） HPC分野の専門家とコラボしての高速化 （特に「富岳」などのA64FX向け） 基本のMP2計算ではVer. 1 Rev. 22比で1.2～1.4倍の加速 ◇最新リリース版 - Open Ver. 2 Rev. 8 （2023年8月） PIEDA機能の強化 （分散力寄与の区別、静電相互作用のRESP評価） MFMOでの励起エネルギーとイオン化エネルギーの算定 MP2計算はVer. 1 Rev. 22比で1.5～2倍の加速 2万フラグメントの液滴モデルの扱い （「富岳」では多サンプル処理が可能） 2024/5/29 19
B. Drobot et al., Phys. Chem. Chem. Phys. 21 (2019) 21213. / MDサンプル構造群に対してFMO計算 / Euの相対論効果はMCPで考慮 統計的相互作用解析 （カルモジュリンとEu(III)の例） ・Ca(II)に比してEu(III)では顕著に安定化 ・部位3が相対的には安定化が小さい 部位1 （100サンプルの平均値） アミノ酸 -386.7 水 -25.9 計 -412.6 ± 15.3 kcal/mol (-245.0±8.3) 部位2 アミノ酸 -376.9 水 -29.3 計 -406.2 ± 14.0 kcal/mol 部位3 (-249.0±8.4) アミノ酸 -307.7 水 -46.8 計 -354.5 ± 14.0 kcal/mol 部位4 (-213.6.±8.0) アミノ酸 -384.2 水 -32.3 計 -416.4 ± 12.9 kcal/mol 2024/5/29 (-250.6.±7.7) 20
A. Rossberg et al., Chem. Comm. 55 (2019) 2015. / MDサンプル構造に対してFMO計算 / Uの相対論効果はMCPで考慮（5fはDZタイプ） 統計的相互作用解析 （DNAとUO22+の例） ウラニルイオンの結合で安定化が減じる 50 45 Frequency 40 35 DNA only 30 DNA+UO2 25 20 15 10 5 0 -760 -740 -720 -700 -680 -660 -640 Sum of stacking and H-bond energy (kcal/mol) -620 （IFIEの総和として評価） 2024/5/29 21
HPCI拠点でのライブラリ整備（2024年3月時点） 赤:A64FX, 緑:SX-Aurora TSUBASA, 紫:Xeon ■登録サイトと版 ・北大「Grand Chariot」： Ver. 1 Rev. 22 & Ver. 2 Rev. 4 ・東北大「AOBA-A/S」： Ver. 1 Rev. 22 & Ver. 2 Rev. 4 (vectorized version) ・東大「Wisteria / Odyssey & Aquarius」： Ver. 1 Rev. 22 & Ver. 2 Rev. 4 ・分子研「RCCS」： Ver. 1 Rev. 22 & Ver. 2 Rev. 4 ・名大「不老 / Type I & II」： Ver. 1 Rev. 22 & Ver. 2 Rev. 4 & Ver. 2 Rev. 8 ・京大「Camphor3」: Ver. 1 Rev. 22 & Ver. 2 Rev. 8 ・阪大「SQUID」： Ver. 1 Rev. 22 (vectorized version) ・R-CCS「富岳」： Ver. 1 Rev. 22 & Ver. 2 Rev. 4 & Ver. 2 Rev. 8 ・計算科学振興財団「FOCUSスパコン」： Ver. 1 Rev. 22 & Ver. 2 Rev. 4 ・九大「ITO Subsystem-A」： Ver. 1 Rev. 22 & Ver. 2 Rev. 4 (to be replaced) ・ Ver. 1 Rev. 22はBioStation Viewerとの関係で当面は併存 ・ Ver. 2 Rev. 8も2024年度内に随時整備の予定 ・ HPCI拠点のスパコン機種更新に応じた対応も適宜行う方針 （北大、東北大、東大、東工大、分子研、阪大、九大） 2024/5/29 22
「富岳」のOnDemand(Beta)に登録済 https://www.hpci-office.jp/for_users/appli_software#openondemand 2024/5/29 23
FMO-HF計算 2024/5/29 24
Ref.; T. Nakano et al., Chem. Phys. Lett. 318 (2000) 684. FMO-HF計算とBDAによる結合の切断 F x C x S x C xε x ⇔ HFの一般化固有値問題 （要反復計算） Fx H x Gx x core x H H Vx Bk k k ⇔ 1電子部分の修飾 k K u Z A r A AK K K K v PK Vx u v K ⇔ 環境静電ポテンシャル K x 1 x P 2 C* iC i G Px 2 i 1 x occ Bond Detached Atom (BDA) I R1 H N O I R2 H R3 N H sp3混成の炭素 （メタンより作成 の局在化軌道） 2024/5/29 ⇔ Hキャップをする必要なく切断可能 J O ⇔ 2電子部分 （N4） （並列処理） R1 I J O H N R3 N H R1 J O H N R3 N H O R2 H フラグメントIの変分空間 O R2 H フラグメントJの変分空間 フラグメントIの形式電荷: +1 BDAの炭素原子の核荷電: 5 フラグメントJの形式電荷: -1 BDAの炭素原子の核荷電: 1 ペプチド結合で切断されていないことには注意が要る 25
ESP approximation Ref.; T. Nakano et al., Chem. Phys. Lett. 351 (2002) 475. & CMM Ref.; 中野ら, J. Comp. Aid. Chem. 13 (2012) 44. FMO-HF計算での高速化のための近似 環境静電ポテンシャル計算の高速化 esp-aoc近似 (Mulliken AO population) L v (P L S L ) ( | ) L ・演算量がO(N3)からO(N2)に減少 ・vdW接触距離の和が指標 ・設定閾値はソフトによって異なり得る for Rmin ( X , L) Laoc Elec. Repulsion esp-ptc近似 (Mulliken atomic charge) v QA r A L AL for Rmin ( X , L) L ptc R=Ldimer-es× r1 Nuc.-Elec. Attraction r2 (vdW(r1)+vdW(r2)) QA (P LS L ) rr rA R Nuc. Repulsion Fragment I Fragment J ダイマー計算の高速化（Dimer-ES近似） I J E ' IJ E ' I E ' J Tr P I u J Tr P J u I P P ( | ) I J Dimer-ES近似の連続多重極展開（CMM）による近似 遠方の対を処理、フラグメント数が千超になるとコスト削減に有効 2024/5/29 26
FMO-HF計算の流れ 分子をフラグメントに分割し、電子を割り当てる モノマーの段階で 自己無撞着電荷 を課すのが特徴 モノマーの初期電子密度を計算する 与えられた電子密度を用いて、周囲のモノマーによる環境静 電ポテンシャル下での、モノマーのエネルギーと電子密度を 計算する No 与えられた電子密度と 得られた電子密度の差が 閾値よりも小さい モノマー Yes 収束したモノマーの電子密度(Self-Consistent Charge; SCC) を用いて、周囲のモノマーによる環境静電ポテンシャル下で の、ダイマーのエネルギーと電子密度を計算する 分子の全エネルギーおよび全電子密度を計算 ダイマー 2024/5/29 27
Ref.; S. Obara et al., J. Chem. Phys. 84 (1986) 3963. 2電子積分の生成について#1 ・ 小原のVertical Recurrence Relation (VRR)がベース ・ ジェネレータでF90コード群を自動生成 （微分も） 中野氏 ・ Gauss型関数の角運動量の昇降を利用 ・ 並列化は短縮シェルの対のループで Gijk : (r ) Nx i y k z l exp(r 2 ) 2024/5/29 28
2電子積分の生成について#2 ・ 軌道タイプの組み合わせに応じて個別のルーチンに ・ 高い軌道角運動量のルーチンの最深部ループは長い ・ 微分も同様に組み合わせ毎 （_gradが付く） sub_dddd.F90 sub_dddf.F90 sub_dddp.F90 sub_ddds.F90 sub_ddfd.F90 sub_ddff.F90 sub_ddfp.F90 sub_ddfs.F90 sub_ddpd.F90 sub_ddpf.F90 sub_ddpp.F90 sub_ddps.F90 sub_ddsd.F90 sub_ddsf.F90 sub_ddsp.F90 sub_ddss.F90 sub_dfdd.F90 sub_dfdf.F90 sub_dfdp.F90 sub_dfds.F90 sub_dffd.F90 sub_dfff.F90 sub_dffp.F90 sub_dffs.F90 sub_dfpd.F90 sub_dfpf.F90 sub_dfpp.F90 sub_dfps.F90 sub_dfsd.F90 2024/5/29 sub_dfsf.F90 sub_dfsp.F90 sub_dfss.F90 sub_dpdd.F90 sub_dpdf.F90 sub_dpdp.F90 sub_dpds.F90 sub_dpfd.F90 sub_dpff.F90 sub_dpfp.F90 sub_dpfs.F90 sub_dppd.F90 sub_dppf.F90 sub_dppp.F90 sub_dpps.F90 sub_dpsd.F90 sub_dpsf.F90 sub_dpsp.F90 sub_dpss.F90 sub_dsdd.F90 sub_dsdf.F90 sub_dsdp.F90 sub_dsds.F90 sub_dsfd.F90 sub_dsff.F90 sub_dsfp.F90 sub_dsfs.F90 sub_dspd.F90 sub_dspf.F90 sub_dspp.F90 sub_dsps.F90 sub_dssd.F90 sub_dssf.F90 sub_dssp.F90 sub_dsss.F90 sub_fddd.F90 sub_fddf.F90 sub_fddp.F90 sub_fdds.F90 sub_fdfd.F90 sub_fdff.F90 sub_fdfp.F90 sub_fdfs.F90 sub_fdpd.F90 sub_fdpf.F90 sub_fdpp.F90 sub_fdps.F90 sub_fdsd.F90 sub_fdsf.F90 sub_fdsp.F90 sub_fdss.F90 sub_ffdd.F90 sub_ffdf.F90 sub_ffdp.F90 sub_ffds.F90 sub_fffd.F90 sub_ffff.F90 sub_fffp.F90 sub_fffs.F90 sub_ffpd.F90 sub_ffpf.F90 sub_ffpp.F90 sub_ffps.F90 sub_ffsd.F90 sub_ffsf.F90 sub_ffsp.F90 sub_ffss.F90 sub_fpdd.F90 sub_fpdf.F90 sub_fpdp.F90 sub_fpds.F90 sub_fpfd.F90 sub_fpff.F90 sub_fpfp.F90 sub_fpfs.F90 sub_fppd.F90 sub_fppf.F90 sub_fppp.F90 sub_fpps.F90 sub_fpsd.F90 sub_fpsf.F90 sub_fpsp.F90 sub_fpss.F90 sub_fsdd.F90 sub_fsdf.F90 sub_fsdp.F90 sub_fsds.F90 sub_fsfd.F90 sub_fsff.F90 sub_fsfp.F90 sub_fsfs.F90 sub_fspd.F90 sub_fspf.F90 sub_fspp.F90 sub_fsps.F90 sub_fssd.F90 sub_fssf.F90 sub_fssp.F90 sub_fsss.F90 sub_pddd.F90 sub_pddf.F90 sub_pddp.F90 sub_pdds.F90 sub_pdfd.F90 sub_pdff.F90 sub_pdfp.F90 sub_pdfs.F90 sub_pdpd.F90 sub_pdpf.F90 sub_pdpp.F90 sub_pdps.F90 sub_pdsd.F90 sub_pdsf.F90 sub_pdsp.F90 sub_pdss.F90 sub_pfdd.F90 sub_pfdf.F90 sub_pfdp.F90 sub_pfds.F90 sub_pffd.F90 sub_pfff.F90 sub_pffp.F90 sub_pffs.F90 sub_pfpd.F90 sub_pfpf.F90 sub_pfpp.F90 sub_pfps.F90 sub_pfsd.F90 sub_pfsf.F90 sub_pfsp.F90 sub_pfss.F90 sub_ppdd.F90 sub_ppdf.F90 sub_ppdp.F90 sub_ppds.F90 sub_ppfd.F90 sub_ppff.F90 sub_ppfp.F90 sub_ppfs.F90 sub_pppd.F90 sub_pppf.F90 sub_pppp.F90 sub_ppps.F90 sub_ppsd.F90 sub_ppsf.F90 sub_ppsp.F90 sub_ppss.F90 sub_psdd.F90 sub_psdf.F90 sub_psdp.F90 sub_psds.F90 sub_psfd.F90 sub_psff.F90 sub_psfp.F90 sub_psfs.F90 sub_pspd.F90 sub_pspf.F90 sub_pspp.F90 sub_psps.F90 sub_pssd.F90 sub_pssf.F90 sub_pssp.F90 sub_psss.F90 sub_sddd.F90 sub_sddf.F90 sub_sddp.F90 sub_sdds.F90 sub_sdfd.F90 sub_sdff.F90 sub_sdfp.F90 sub_sdfs.F90 sub_sdpd.F90 sub_sdpf.F90 sub_sdpp.F90 sub_sdps.F90 sub_sdsd.F90 sub_sdsf.F90 sub_sdsp.F90 sub_sdss.F90 sub_sfdd.F90 sub_sfdf.F90 sub_sfdp.F90 sub_sfds.F90 sub_sffd.F90 sub_sfff.F90 sub_sffp.F90 sub_sffs.F90 sub_sfpd.F90 sub_sfpf.F90 sub_sfpp.F90 sub_sfps.F90 sub_sfsd.F90 sub_sfsf.F90 sub_sfsp.F90 sub_sfss.F90 sub_spdd.F90 sub_spdf.F90 sub_spdp.F90 sub_spds.F90 sub_spfd.F90 sub_spff.F90 sub_spfp.F90 sub_spfs.F90 sub_sppd.F90 sub_sppf.F90 sub_sppp.F90 sub_spps.F90 sub_spsd.F90 sub_spsf.F90 sub_spsp.F90 sub_spss.F90 sub_ssdd.F90 sub_ssdf.F90 sub_ssdp.F90 sub_ssds.F90 sub_ssfd.F90 sub_ssff.F90 sub_ssfp.F90 sub_ssfs.F90 sub_sspd.F90 sub_sspf.F90 sub_sspp.F90 sub_ssps.F90 sub_sssd.F90 sub_sssf.F90 sub_sssp.F90 sub_ssss.F90 sub_xxxx.F90 suberi.F90 swz_init.F90 29
C2-DIIS Ref.; H. Sellers, Intern. J. Quant. Chem. 45 (1993) 31. & SOSCF Ref.; G. B. Bacskay, Chem. Phys. 61 (1981) 385. ODA Ref.; E. Cances et al., Intern. J. Quant. Chem. 79 (2000) 82. & Anderson Expol. Ref.; D. G. Anderson, J. Ass. Comp. Machin. 12 (1965) 547. HF-SCF、monomer-SCCについて モノマーSCCの収束の例 ■SCF ・ ODA/C2-DIISがデフォルト、EDIISも可 ・ full-SOSCF、近似aug.-Hessianも可 ・ 差分Fock行列の構築は未実装 ■SCC ・ Anderson法で加速 （Fockベース） ・ PAICSの”山彦アルゴリズム”は未対応 ■SCFのトラブル ・ 構造が悪いと収束難 （特にダイマー） ・ Hの付け方の悪さが顕在化 ・ 対イオンの扱いにも依存 （MD由来） ■SCCのトラブル ・ 系全体の電荷条件が極端だと収束難 ・ 構造由来のトラブルもアリ ・ 超並列実行時は相対コストが上昇 2024/5/29 Iteration Electronic energy delta E NSCF Threshold theta ---------------------------------------------------------------------------1 -147077.1337872949 -147077.1337872949 1264 0.0100000000 0.000 2 -148886.4265209468 -1809.2927336519 584 0.0100000000 0.000 3 -148267.4538700682 618.9726508786 406 0.0100000000 0.000 4 -148431.2708061265 -163.8169360583 406 0.0100000000 0.051 5 -148369.4298011306 61.8410049959 406 0.0100000000 0.242 6 -148382.9245709342 -13.4947698036 406 0.0100000000 0.143 7 -148378.6162133594 4.3083575748 406 0.0100000000 0.403 8 -148377.5251745595 1.0910387999 406 0.0100000000 0.125 9 -148378.7268076855 -1.2016331260 406 0.0100000000 0.501 10 -148377.2939791832 1.4328285023 406 0.0100000000 0.092 11 -148378.3143993220 -1.0204201388 406 0.0100000000 0.520 12 -148377.5799346490 0.7344646730 406 0.0100000000 0.113 13 -148378.1111756845 -0.5312410356 406 0.0010000000 0.434 14 -148377.8052994407 0.3058762439 406 0.0010000000 0.196 15 -148378.0381476319 -0.2328481912 406 0.0010000000 0.223 16 -148377.9382764302 0.0998712017 406 0.0010000000 0.303 17 -148378.0168179379 -0.0785415077 406 0.0001000000 0.000 18 -148377.9999593449 0.0168585930 406 0.0001000000 0.190 19 -148378.0211514597 -0.0211921148 406 0.0001000000 0.000 20 -148378.0218082173 -0.0006567576 406 0.0001000000 0.000 21 -148378.0302730768 -0.0084648596 584 0.0000000100 0.000 22 -148378.0501479235 -0.0198748467 589 0.0000000100 1.627 23 -148378.0354601964 0.0146877271 424 0.0000000100 2.305 24 -148378.0416888229 -0.0062286265 406 0.0000000100 0.017 25 -148378.0370212285 0.0046675944 406 0.0000000100 0.000 26 -148378.0385753561 -0.0015541276 406 0.0000000100 0.000 27 -148378.0379157953 0.0006595608 406 0.0000000100 0.000 28 -148378.0382292477 -0.0003134524 406 0.0000000100 0.000 29 -148378.0381183651 0.0001108827 406 0.0000000100 0.000 30 -148378.0381876900 -0.0000693249 406 0.0000000100 0.000 31 -148378.0381702714 0.0000174186 406 0.0000000100 0.000 32 -148378.0381862515 -0.0000159801 406 0.0000000100 0.000 33 -148378.0381839521 0.0000022994 406 0.0000000100 0.000 34 -148378.0381881352 -0.0000041831 406 0.0000000100 0.000 35 -148378.0381888413 -0.0000007061 406 0.0000000100 0.000 36 -148378.0381892572 -0.0000004159 406 0.0000000100 0.000 SELF CONSISTENT CHARGE COMPLETED 30
FMO-MP2計算 2024/5/29 31
Ref.; R. J. Bartlett, Ann. Rev. Phys. Chem. 32 (1981) 359. 4次までの摂動論の図形表現 3rd-order 2nd-order MP4はMP3より処理は格段に複雑になる （1電子、3電子、4電子励起の寄与） 2024/5/29 32
摂動論による電子相関の導入 揺動ポテンシャルが摂動 2次摂動(MP2) の相関： 占有軌道iとjの電子対が“衝突”して非占有軌道aとbに“散乱” a b ⇔ i ⇔ E MP2 j ij || ab ab || ij 1 4 ijab i j a b （分子は電子間の反発積分の積、分母は軌道エネルギーの差） （電子対の相互作用） 3次摂動(MP3) の相関： 電子対間の相互作用も3通りで追加的に考慮して補正 ⇔ E E ・ MP3の方がMP2よりも縮約が複雑 ・ 空間軌道の式に落として実装 ・ DGEMM処理を上手く利用する MP3(Term 1) MP3(Term 2 ) ij || ab ab || cd cd || ij 1 8 ijabcd ( i j a b )( i j c d ) ij || ab ab || kl kl || ij 1 8 ijklab( i j a b )( k l a b ) E MP3(Term 3) ij || ab kb || cj ac || ik ( i j a b )( i k a c ) ijkabc 2024/5/29 33
Ref.; Y. Mochizuki et al., Theor. Chem. Acc. 112 (2004) 442. 積分直接利用-並列化MP2の実装 2024/5/29 ・ ファイルIO無し ・ 通信量を最小化するアルゴリズム ・ Level-1 BLASを使用 （演算数低減を重視） ・ 近似電子密度の計算も実装 34
Ref.; Y. Mochizuki et al., Chem. Phys. Lett. 396 (2004) 473. FMO-MP2計算の実用化 ・ HIV-protease+LopinavirのFMO-MP2/6-31Gジョブが14.3時間 （当時のXeon 64コア） ・ これ以後、MP2計算が実在系のタンパク質系に対して常用されることに 2024/5/29 35
1st Refs.; Y. Mochizuki et al., Chem. Phys. Lett. 396 (2004) 473. & Y. Mochizuki et al., Theor. Chem. Acc. 112 (2004) 442. 2nd Ref.; Y. Mochizuki et al., Chem. Phys. Lett. 457 (2008) 396. MP2計算の基本アルゴリズム2種： 形式的にN5 MP2相関エネルギー補正 Loop over i-batch [parallelizable when needed] Loop over [to be parallelized for worker processes] Loop over Generate (,[]) list [canonical relation for ] Do 1/4 transformation of →i [screening & DAXPY] Do 2/4 transformation of →a [DDOT] Do 3/4 transformation of →j [screening & DAXPY] End of loop over Do 4/4 transformation for →b [screening & DAXPY] End of loop over [“all-reduce” must be done for (ia|jb)] Calculate partial MP2 energy End of loop over i-batch (ia | jb)2(ia | jb) (ib | ja ) i j a b ijab EMP2 添字の変換： 4N5コスト (ia | jb) Cb Cj Ca Ci ( | ) Loop over ij-batch ! size depending on available memory Loop over ! to be parallelized Loop over Preparing (|) ! for canonical -pair Forming (i|) ! DGEMM, fixed , running over Forming (ia|) ! DGEMM, fixed , running over Forming (ia|j) ! DGEMM, fixed , direct-product for fixed End of loop over Forming (ia|jb) ! DGEMM, direct-product for fixed End of loop over ! all-reduce operation as barrier Calculate partial MP2 energy with respect to ij-batch End of loop over ij-batch 2024/5/29 【第一版】 ・ DAXPYとDDOTを使い、 閾値判断を優先 ・ 実効的演算数を下げる方針、 15年程前のチップでは有効 【第二版】 ・ DGEMMの高性能に期待 【現在：混成型がデフォルト】 ・ スカラ型のCPUでもステップ2 とステップ4はDGEMMで 36
Ref.; Y. Mochizuki et al., Chem. Phys. Lett. 457 (2008) 396. FMO-MP2計算のベクトル化 2024/5/29 37
Ref.; Y. Mochizuki et al., Chem. Phys. Lett. 457 (2008) 396. 初代地球シミュレータ(ES)でのベンチマーク Timing data of influenza HA antigen-antibody @ FMO-MP2/6-31G VPUs @ ES FMO-MP2 1024 2048 3072 4096 FMO-HF 1024 2048 3072 4096 Time (s) 10084.6 5486.2 3927.1 3204.4 4164.9 2131.7 1533.5 1205.9 Accel. TFLOPS Effic. (%) Rel. cost 1.84 2.57 3.15 1.19 2.19 3.06 3.75 14.5 13.3 12.4 11.5 1.95 2.72 3.45 1.36 2.67 3.72 4.75 16.6 16.3 15.2 14.5 2.42 2.57 2.56 2.66 Hemagglutinin (HA) Neuraminidase (NA) Membrane Protein 2 (M2) Matrix Protein (M1) Lipid Bilayer * PDB-id 1EO8 をベースにしてFMO計算用の構造を調製 * 計921残基（計911フラグメント：Cys-Cysマージ）, 計14,086原子 * 6-31G基底, 関数総数78,390 （最大級の計算：当時） * 4,096VPUで僅か53.4分でFMO-MP2ジョブが完了可能 * FMO-HFに比べて僅か2.6倍のコスト * 数千VPUを有効に利用出来る高い並列効率を実証 2024/5/29 1 Nucleoprotein (NP) Polymerase 2 3 (PB1,PB2,PA) 4 5 6 RNA 7 8 PB2 PB1 PA HA NP NA M1 NS1 M2 NEP NEP 38
Ref.; K. Takematsu et al., J. Phys. Chem. 113 (2009) 4991. インフルエンザHAとFabの解析例 現代の計算機なら1/10以下の時間で可能 インフルエンザHA の変異領域の解析 Red: Prohibited, Green: Allowed, Blue: Not measured Opteron 16 cores (2 GHz) * 6-31G （78390関数） ⇒ 8.3日で完了 * 6-31G* （121314関数） ⇒ 22.7 日で完了 Xeon 128 cores (2.66 GHz) * 6-31G* ⇒ 2.2 日で完了 ⇒ ルーチン実行可能 A: Charged, B: Hydrophobic, C: Polar 2024/5/29 39
Ref.; Y. Mochizuki et al., J. Nucl. Sci. Tech. 39 (2002) 195. / DIRACプログラムのベクトル化の際に導入. Fock行列構築のベクトル化 「粒子推進（止まり木）アルゴリズム」 ・ 積分の生成もベクトル化 （_vecで区別） ・ 配列を3次元化してアドレス衝突を回避 ・ ベクトル長は256をデフォルト 2024/5/29 40
Ref.; Y. Okiyama et al., Chem. Phys. Lett. 490 (2010) 84. General CD Ref.; T. B. Pedersen et al., Theor. Chem. Acc. 124 (2009) 1. コレスキー分解（CD)による近似： テンソル次数の低減 MP2 correlation energy (𝐸 MP2 ) Loop over I # Parallelized Coulomb matrix (J) 𝐽ሚ𝑝𝑞 = 0 Loop over I # Parallelized 𝑄𝐼 = 𝐿𝐼,𝑟𝑠 𝑃𝑟𝑠 読み替え：pqrs → μνλσ # DGEMM 𝑋𝐼,𝑝𝑖 = 𝐿𝐼,𝑝𝑟 𝐶𝑟𝑖 𝑟 # DGEMM 𝐵𝐼,𝑖𝑎 = 𝑋𝐼,𝑞𝑖 𝐶𝑞𝑎 # DDOT 𝑞 𝑟𝑠 End of loop over I 𝐽ሚ𝑝𝑞 = 𝐽ሚ𝑝𝑞 + 𝐿𝐼,𝑝𝑞 𝑄𝐼 # DAXPY End of loop over I ALLREDUCE 𝐽ሚ𝑝𝑞 𝐸 MP2 = 0 Loop over ij # Canonical ij pair (𝑖𝑎|𝑗𝑏) = 𝐵𝐼,𝑖𝑎 𝐵𝐼,𝑗𝑏 # DGEMM (partial sum of I ) I Exchange matrix (K) 𝐾෨𝑝𝑞 = 0 ALLREDUCE (𝑖𝑎, 𝑗𝑏) 2 − 𝛿ij 2 𝑖𝑎|𝑗𝑏 − 𝑖𝑏|𝑗𝑎 𝐸 MP2 = 𝐸 MP2 + 𝜀𝑖 − 𝜀𝑎 + 𝜀𝑗 − 𝜀𝑏 End of loop over ij Loop over I # Parallelized 𝑋𝐼,𝑞𝑖 = 𝐿𝐼,𝑞𝑠 𝐶𝑠𝑖 # DGEMM 𝑖𝑎|𝑗𝑏 𝑠 𝐾෨𝑝𝑞 = 𝐾෨𝑝𝑞 + 𝑋𝐼,𝑝𝑖 𝑋𝐼,𝑞𝑖 𝑖 End of loop over I ALLREDUCE 𝐾෨𝑝𝑞 2024/5/29 # DGEMM ( pq | rs ) LI , pq LI ,rs I ・ CDでは分解基底は軌道積に取る： 展開長MはN2オーダー ・ CD 1中心(1C)化は同一原子内に限る： MはNに近くなる ・ RI近似では最適化基底を別途調製： 展開長はほぼN 41
ABINIT-MPはVer.6準拠 / Regularが通常の2電子積分を使うアルゴリズム. / 1Cは1中心の制約付き. 小型タンパク質での「京」での性能評価の例#1 NLYIQWLKDGGPSSGRPPPS FMO2-MP2ジョブ @ 24ノード Basis set Time (Min) Eff. (%) Time (Min) Eff. (%) TrpCage (20 frag.) Trp23His 6-31G Regular 12.6 2.3 36.8 2.53 CD/MP2 8.2 4.4 19.5 7.13 CD/MP2&HF 5.7 15.6 14.6 18.35 1C-CD/MP2 7.6 2.5 16.7 2.97 1C-CD/MP2&HF 3.0 8.8 8.2 9.66 6-31G* Regular 38.2 3.0 136.9 3.43 CD/MP2 28.0 6.6 96.0 8.87 CD/MP2&HF 18.9 29.0 44.6 31.89 1C-CD/MP2 25.9 3.4 84.7 3.88 1C-CD/MP2&HF 8.3 18.8 22.0 21.06 cc-pVDZ Regular 76.7 2.9 267.8 3.52 CD/MP2 55.8 6.1 188.6 8.06 CD/MP2&HF 34.7 26.3 81.0 30.43 1C-CD/MP2 52.1 3.2 170.1 3.77 1C-CD/MP2&HF 16.6 14.9 39.1 18.86 2024/5/29 TEI treatment 【TrpCage】 𝑋𝐼,𝑞𝑖 = 𝐿𝐼,𝑞𝑠 𝐶𝑠𝑖 𝑠 𝐾𝑝𝑞 = 𝐾𝑝𝑞 + 𝑋𝐼,𝑝𝑖 𝑋𝐼,𝑞𝑖 ・ CDではHF計算の交換項(K)をDGEMMで処理 ・ モデル系では最高で30%を超える効率を達成 ・ トータルでの時間短縮は1Cでは最高6倍程度 𝑖 42
小型タンパク質での「京」での性能評価の例#2 FMO4-MP2ジョブ FMO2-MP2ジョブ Basis set TEI treatment Time (Min) Eff. (%) HIV-1 Protease + Lopinavir (204 frag.) @ 204 nodes 6-31G Regular 31.9 1.85 CD/MP2&HF 18.8 7.78 1C-CD/MP2&HF 14.8 3.74 6-31G* Regular 74.5 2.67 CD/MP2&HF 44.0 17.54 1C-CD/MP2&HF 29.4 8.30 cc-pVDZ Regular 161.3 2.69 CD/MP2&HF 80.4 17.47 1C-CD/MP2&HF 58.3 8.02 Infulenza HA1 + Fab (911 frag.) @ 912 nodes 6-31G Regular 123.2 1.01 CD/MP2&HF 110.7 2.27 1C-CD/MP2&HF 106.1 1.14 6-31G* Regular 212.1 1.62 CD/MP2&HF 166.5 6.47 1C-CD/MP2&HF 146.5 2.66 cc-pVDZ Regular 400.5 2.32 CD/MP2&HF 296.6 8.14 1C-CD/MP2&HF 269.9 3.37 2024/5/29 Basis set TEI treatment Time (Min) Eff. (%) Crambin (94 frag.) @ 960 nodes 6-31G Regular 23.7 2.06 CD/MP2&HF 9.1 12.06 1C-CD/MP2&HF 6.7 6.38 HIV-1 Protease + Lopinavir (363 frag.) @ 6000 nodes 6-31G Regular 144.3 2.79 CD/MP2&HF 128.0 11.85 1C-CD/MP2&HF 49.5 6.88 FMO2計算 → 数百残基の計算は容易 → 6-31G*基底の実用利用 → 実行効率は及第レベル → 1C-CDで時間短縮 （IFIEの議論は問題無） FMO4計算 → 多ノード利用のメリット → 1中心化の効果は顕著 → 詳細解析の用途を推奨 → HFのコスト低減は課題 【インフルHA+Fab】 43
「京」での大規模テスト計算の例 − − − − インフルエンザノイラミニダーゼ（386残基） FMO4-MP2/6-31G（1C-CD） フラグメント数 706 （主鎖/側鎖分割、リガンド分割） OpenMP版 THREADS=8 ・Kfast+部分的O2指定 (*OpenMP版 THREADS=8 ・Kfast指定) 20000 18000 5.4時間* + O NH O Tetramer MP2 Trimer MP2 14000 Time (sec) H 3N O 1.5時間 16000 Dimer MP2 12000 73%減 10000 8000 Monomer MP2 Tetramer SCF Trimer SCF 6000 Dimer SCF 4000 Dimer ES 2000 Monomer SCF 0 706 nodes * 2024/5/29 O- 12504 nodes （計算時間は各々5.0時間と1.1時間） 効率は6% @ 100032コア 44
CMM Ref.; 中野ら, J. Comp. Aided Chem. 13 (2012) 44.. / インフルエンザ HA3+Fab抗体 / 2325フラグメント / cc-pVDZでは36.4万関数. Dimer-ESのCMM近似 （旧FX-100@名大） (Sec.) 4096コア 40000 ・ PDB-id 1KENのインフルHA３+2Fab複合体 ・ 16スレッド – 256プロセス (128ノード) ・ FMO2-MP2、16スレッド/フラグメント ・ 連続多重極(CMM)をLdimer=5に適用 ・ 6-31Gでは当該部12.4時間が0.5時間に （全エネルギーの差は0.0002au） ・ 6-31G*とcc-pVDZのタイミングは右図 ・ 数千フラグメントの計算が容易に 35000 30000 25000 20000 Dimer Monomer 15000 10000 モノマーの コストが大 5000 0 Regular 2024/5/29 Full CD 6-31G* 1C CD Regular Full CD cc-pVDZ 2C CD 45
OpenMP/MPIの2階層並列 / Open Ver. 1 Rev. 5を使用. OakForest-PACS（JCAHPC）でのベンチマーク Xeon Phi KNL /compact/cacheモードのノード群を使用 AVX-512 - HIV-1 Protease + Lovinavir (203 fragments) FMO2-MP2/6-31G* level Time (m) #Nodes #MPI/node #Thread/MPI #Threads Regular Full CD 1C-CD 1 2 64 128 1922.8 1256.2 (x1.54) 907.4 (x2.12) 1 4 64 256 Acc 2.06 935.7 616.3 (x1.52) 432.7 (x2.17) 2 2 64 256 Acc 1.99 967.2 637.9 (x1.52) 452.6 (x2.14) 2 4 64 512 Acc 4.01 479.0 310.1 (x1.55) 216.8 (x2.21) 4 4 64 1024 Acc 7.86 244.6 157.7 (x1.56) 110.3 (x2.22) 4 16 4 256 562.0 326.6 (x1.73) 220.3 (x2.56) 4 16 8 512 293.7 175.4 (x1.68) 117.2 (x2.51) 4 16 16 1024 168.9 94.3 (x1.80) 62.5 (x2.71) 4 8 32 1024 192.0 111.3 (x1.73) 78.4 (x2.45) 8 16 16 2048 107.4 62.5 (x1.72) 40.0 (x2.69) 8 8 32 2048 100.6 56.6 (x1.78) 39.7 (x2.54) 16 8 32 4096 64.3 38.1 (x1.69) 24.1 (x2.67) 2024/5/29 ・compactモードでは物理コアに対してHT優先で割当られる ・スレッド単位は64よりも32ないし16の方がよさそう 46
OpenMP/MPI混成並列 （32スレッド/フラグメント、ノード内のMPIプロセス数は8）. Oakｆorest-PACSでのABINIT-MPのベンチ例#1 Influenza NA & Tamiflu - FMO2-MP2/6-31G* @ OFP (AVX512 on) / in Minutes 600.0 4 nodes 8 nodes 16 nodes CD不使用 32 nodes 500.0 400.0 300.0 200.0 100.0 0.0 4w/oCMM Dimer others 2024/5/29 4w/CMM Dimer MP2 8w/oCMM Dimer SCF 8w/CMM 16w/oCMM Dimer ES Monomer others 16w/CMM 32w/oCMM Monomer MP2 ・16ノードまでは旨味のある加速が得られる ・Dimer-ESを多重極展開(CMM)の効果は有り ・「その他処理」のオーバーヘッドは少ない 32w/CMM Monomer SCF フラグメント数が400弱 47
OpenMP/MPI混成並列 （32スレッド/フラグメント、ノード内のMPIプロセス数は8）. Oakｆorest-PACSでのABINIT-MPのベンチ例#2 (min.) Influenza HA monomer + Fab - FMO2-MP2/6-31G* @ OFP (AVX512 on) 500.0 16 nodes 32 nodes 64 nodes 450.0 400.0 350.0 300.0 250.0 200.0 150.0 100.0 50.0 0.0 1 2 Dimer others 2024/5/29 Dimer MP2 3 Dimer SCF Dimer ES 4 Monomer others 5 Monomer MP2 6 Monomer SCF ・911フラグメント、分子形状が細長いためDimer-ES近似が重くなる ・5Å以上のES連続多重極展開(CMM)で処理、誤差は僅か0.0001au ・ “Others”の処理が顕在化している ⇒ 対策が必要 ・ Monomer SCFのコストも目立つ ⇒ SCC条件のボトルネック 48
KMP_AFFINITY=scatter & IntelMPI の各MPIプロセスの割当先論理コアを指定 （HPCシステムズの渡邊氏による測定）. SMASH Benchmark Ref.; 齊藤ら, J. Comp. Chem. Jpn. 15 (2016) 92. OakForest-PACSでのHTの有効性チェック HT倍率2に対して1.4倍、4に対して1.6倍 2階層版 HIV-1 OFP 8ノード 2階層版 HIV-1 OFP 8ノード HT使用時 経過時間(秒) HT使用時 速度向上率 MP2-regular MP2-regular 12000.0 1.600 10000.0 1.500 8000.0 1.400 6000.0 1.300 4000.0 1.200 2000.0 1.100 0.0 1.000 16MPI×1 16MPI×2 16MPI×4 16MPI×1 16MPI×2 16MPI×4 ・分子軌道ソフトSMASHをKNCでベンチした際と同程度の結果 （倍率2で1.4倍、4で1.7倍： 計算の中身が似ている?!） ・「演算器が余っている」とも思わせる （積分計算などで改良が必要なことを示唆） 2024/5/29 49
Ref.; 渡邊ら, J. Comp. Chem. Jpn. 15 (2016) 92. FMO-MP2計算の3階層並列処理 フラグメントリストのタスクをMPIで並列処理 (1層目) Frag. Task 1 Frag. Task 2 Frag. Task 3 subroutine fmo_layer1(…) subroutine monomer_mp2(…) ・ ・ ・ テストの結果は必ずしも良くなかった… Frag. Task N 効率的な並列数を 選択可能となる メモリ増大量を 最小限におさえつつ 並列数を倍増できる subroutine MP2DRIVE_IJ(…) subroutine DMP_ENE_IJ(…) : DO IXSCS=1,IXNCS DO IXRCS=1,IXNCS : call get_tei_rs_fix(…) : 2ループを複合し MPI/OpenMP並列実行 (2層目) subroutine get_tei_rs_fix(…) : ! kth contracted shell : ! lth contracted shell : ! ith contracted shell do ixics=1,ixncs : ! jth contracted shell do ixjcs=1,ixics : call suberi(…) : 重い内部2ループを複合し OpenMP並列実行 (新設3層目) 2024/5/29 フラグメント内処理 subroutine suberi(…) : index=… select case(index) case( 0) !(ssss) call sub_ssss(…) : caee(255) !(ffff) call sub_xxxx(…) end select 50
FMO-MP3計算 2024/5/29 51
Ref.; R. J. Bartlett, Ann. Rev. Phys. Chem. 32 (1981) 359. 高次相関の寄与 3次は「逆方向の寄与」となる場合も 2024/5/29 52
MP3 Ref.; Y. Mochizuki et al., Chem. Phys. Lett. 493 (2010) 346. / MP2.5 Ref.; M. Pitonak et al., Chem. Phys. Chem. 10 (2009) 282. 3次摂動論（MP3）による相関エネルギー補正 ・ MP2に無い電子対間の相互作用を考慮 ・ OpenMP/MPI混成並列を前提に実装 ・ メモリ要求はMP2より増加、O(N6)コスト ・ 縮約/変換処理にはDGEMMを多用 ・ MP2.5スケーリングを推奨 (ia | jb)(ik | jl )[2(ka | lb) (kb | la)] ijklab ( i j a b )( k l a b ) E MP3( 2 p4 h ) (ia | jb)(ac | bd )[2(ic | jd ) (id | jc)] ijabcd ( i j a b )( i j c d ) E MP3( 4 p2 h ) [2(ia | jb) (ij | ab)][2(kc | ia) (ka | ic)][2(kc | jb) (kb | jc)] ( i j a b )( i k a c ) ijkabc E MP3(3 p3h )1 E 2024/5/29 MP 3( 3 p 3 h ) 2 (ij | ab)(ka | ic )(kb | jc ) 3 ijkabc ( i j a b )( i k a c ) 和の添字は6つ 53
Ref.; Y. Mochizuki et al., Chem. Phys. Lett. 493 (2010) 346. / 2代目の地球シミュレータ(ES2)を利用. FMO-MP3の実装 2024/5/29 54
MP3計算のループ （ijバッチとkバッチの2重構造） メモリが少ないとバッチ数が増えることに Loop over ij-batch ! size depending on available memory Loop over ! parallelized Form {(ia|jb), (ij|ab), (ik|jl)} ! 1st-4th quarter transformations End of loop over Calculate partial MP2 energy ! parallelized for ij Prepare half-backtransformed 1st-order amplitude ! parallelized for ij Perform EEO processing for “2h-4p” contribution ! parallelized for Calculate partial “2h-4p” energy ! parallelized for ij Loop over k-batch ! size depending on available memory Loop over ! parallelized Form (kc|ld) ! 1st-4th quarter transformations End of loop over Calculate partial “4h-2p” and “3h-3p” energies ! parallelized for ij End of loop over k-batch End of loop over ij-batch Sum up “2h-4p”, “4h-2p” and “3h-3p” energies for final MP3 energy 2(ic | jd ) (id | jc ) (ia | jb) MP 3( 4 p 2 h ) E (ac | bd ) ab i j a b cd ij i j c d 2024/5/29 tijcd t ij Fock行列的処理で (ac|bd)構築を回避 55
FMO-MP2とMP3のタイミング (ES2) ---------------------------------------------------------------------------------- Method Nodes Time (hours) Acc. TFLOPS ---------------------------------------------------------------------------------- HA1 HA3 NA FMO-MP2 FMO-MP3 FMO-MP2 FMO-MP3 FMO-MP2 FMO-MP3 FMO-MP2 FMO-MP3 FMO-MP3 FMO-MP3* 64 64 128 128 64 64 128 128 64 64 1.7 2.7 [x1.6] 0.8 1.3 [x1.6] 9.4 11.9 [x1.3] 4.3 5.8 [x1.3] 1.0 4.4 2.1 2.1 2.2 2.1 0.97 2.27 2.06 4.67 0.83 1.66 1.83 3.44 3.04 3.09 ---------------------------------------------------------------------------------- ・8 VPUs per node, 64 nodes = 512 VPUs, 6-31G basis without asterisk ・HA1; 14086 atoms, 921 res., 911 frag., 78390 AOs ・HA3; 36160 atoms, 2351 res., 2325 frag., 201276 AOs ・NA; 5792 atoms, 386 res., 378 frag., 32549 AOs (50447: 6-31G*) MP3/MP2の相対コストは驚くほど小さい 2024/5/29 56
ノイラミニダーゼとタミフルの相互作用解析#1 MP2 in parenthesis 黄色の分子：タミフル（黄色）と NAの各アミノ酸残基とのIFIE unit : kcal/mol タミフルと強く相互作用する残基 2024/5/29 57
ノイラミニダーゼとタミフルの相互作用解析#2 基底関数の依存性をチェック 2 2 不 安 定 化 IFIE(kcal/mol) 1.5 1 0.5 1.5 1 0.5 0 0 安 定 -0.5 化 -1 -0.5 -1 Fragment number -1.5 Fragment number -1.5 MP2法：基底関数のエネルギーの差 MP3法：基底関数のエネルギーの差 ±0.5kcal/mol 以上をピックアップ MP2・MP3どちらもd関数の 追加の効果は大きくはない 2 GLU277 1.5 GLU227 IFIE(kcal/mol) 1 ARG152 ASP151 SER179 0.5 この系では6-31G基底でも 準定量的な議論が可能 0 -0.5 -1 -1.5 2024/5/29 ARG118 GLU119 ARG224 ASN213 VAL349 ARG292 TYR406 MP2 MP3 静電相互作用と 水素結合が主 58
Ref.; A. Yoshioka et al., J. Mol. Graph. Mod. 30 (2011) 110. / この論文ではMP2.5スケーリングの値を利用. ヘマグルチニン3量体とFab２つの相互作用解析 FMO-MP3/6-31Gレベル Interacting domain HF MP2 MP3 Interacting domain HF MP2 MP3 Fab(I):HA(I) Fab(I):HA(II) Fab(I):HA(III) Fab(II):HA(I) Fab(II):HA(II) Fab(II):HA(III) -288.8 177.5 134.3 137.0 -292.7 170.8 -367.0 155.6 134.2 137.0 -380.4 157.0 -352.8 158.7 134.3 137.0 -363.7 159.5 HA(I):HA(II) HA(II):HA(III) HA(I):HA(III) Fab(I):Fab(II) HA total:Fab total -1022.4 -981.7 -1189.0 210.8 38.1 -1280.3 -1245.7 -1469.8 197.8 -163.6 -1237.1 -1200.6 -1421.3 199.5 -127.0 Dimer-ES領域では相関計算は行われない 2024/5/29 IFIE sum (kcal/mol) 59
その他の計算機能 2024/5/29 60
Ref.; T. Ishikawa et al., Theor. Chem. Acc. 118 (2007) 937. 局在化MP2を使った詳細解析 （FILM） TrpCageの籠内のCH/p相互作用の解析の例 TRP6 CH ⇔ TYR3 π (-0.1513) TRP6 π ⇔ GLY11 CH (-0.4903) TRP6 π ⇔ PRO18 CH (-0.4209) 2024/5/29 TRP6 π ⇔ PRO12 CH (-0.1788) TRP6 π ⇔ PRO18 CH (-0.3839) 61
Ref.; Y. Mochizuki, et. al., Chem. Phys. Lett. 410 (2005) 247. & Y. Mochizuki, Chem. Phys. Lett. 410 (2005) 165. 電荷移動の詳細解析 （CAFI） Vacant Space Vacant Space #1 モノマー計算は通常のFMO #2 リガンドと重要な残基を選択 #3 局在化、重み付Lowdin直交化 #4 CERFによる多体緩和の導入 （N**4、多電子論的な近似） #5 軌道毎の緩和エネルギー評価 水素結合(CT) による安定化は それぞれ -1.9kcal/mol, -1.4kcal/mol 2024/5/29 Charge Transfer Induced Polarization Occupied Space (localized) Fragment A Pre－Polarized by ESP in FMO Occupied Space (localized) Fragment B 電子の移動 方向及び電 荷の移動量 まで解析でき る。 Table CAFI energies (in hartree) for Gly3 - Gly5 hydrogen bonding in glycine pentamer Occupation No. Energy (CT) Energy (POL) Donor orbital Acceptor orbital number 1 -0.003055 0.000023 lone-pair *NH 0.00245 2 -0.002260 0.000057 p lone-pair *NH 0.00332 62
Ref.; Y. Okiyama et al.,. Chem. Phys. Lett. 509 (2011) 67. / BSSEの計算では、”幽霊軌道”を置く計算を2回追加でやることになる. BSSEの補正例 5 Thymine (1’T) 0 kcal/mol Adenine (1A) -5 6-31G 248% 211% 31% 90% 89% 1A-2A (PP) 1A-2'T 2A-2'T (HB) 1A-1'T (HB) 176% -67% -18% 2A-1'T 1'T-2'T (PP) -10 83% -15 -20 HF HF (CP) MP2 71% MP2 (CP) 83% 71% -25 5 kcal/mol Thymine (2’T) -10 81% -15 2A-2'T (HB) 1A-1'T (HB) 197% 13% 22% 2A-1'T 1'T-2'T (PP) HF (CP) MP2 74% 82% 73% HF MP2 (CP) -25 5 kcal/mol 0 -5 -10 -15 -20 2024/5/29 1A-2A 1A-2'T (PP) 46% -20 Adenine (2A) 313% 89% 89% 0 -5 6-31G* 244% -25 6-31G*(0.25) 510% -126% -242% 85% 86% 1A-2A 1A-2'T (PP) 46% 2A-2'T (HB) 1A-1'T (HB) 34% 2A-1'T 31% 1'T-2'T (PP) HF HF (CP) MP2 75% 75% 65% 67% MP2 (CP) 63
Ref.; Y. Okiyama et al., J. Phys. Chem. B 122 (2018) 4457 & 123 (2019) 957. FMO-PBジョブのタイミング例 ## ESTIMATION OF NONPOLAR SOLVATION ENERGY SASA of whole molecule / Angstrom^2 = Elapsed time = 9943.45412248 0.0 sec. ## SOLUTE TOTAL ENERGY SCF / Hartree MP2 corr. / Hartree -77621.87317209 -226.58605103 -77625.72587377 -226.07846056 -3.85270168 0.50759048 -4.60614784 ) FMO2 in vacuo in solvent difference ( ES correction ・ 系はHIV-Protease + Lopinavir ・ FMO-MP2/6-31G*レベル ・ Oakforest Pacs (OFP)の64ノードを使用 ・ フラグメント当たり32スレッド、ノードには 8プロセスをハイパースレッド下で立てた ・ 収束までの経過時間は34.3時間 ・ 3回程度で止めるのが実用的かも。。。 Total / Hartree -77848.45922312 -77851.80433433 -3.34511121 ## SOLVATION FREE ENERGY FMO2 / Hartree -3.85270168 0.11409050 -3.73861118 Electrostatic Nonpolar Total ## Time profile Number of cores (total) = Number of cores (fragment) = 512 1 THREADS (FRAGMENT) 32 Total time = 2024/5/29 = 123455.4 seconds FMO2 / kcal/mol -2417.60678689 71.59286968 -2346.01391721 Iter. Total Energy / Hartree Difference / Hartree -----------------------------------------------------1 -77625.1635597274 -3.2903876412 2 -77625.6117971116 -0.4482373842 3 -77625.6982569271 -0.0864598155 4 -77625.7183753774 -0.0201184503 5 -77625.7236744650 -0.0052990876 6 -77625.7251996845 -0.0015252195 7 -77625.7256641467 -0.0004644622 8 -77625.7258196045 -0.0001554578 9 -77625.7258510916 -0.0000314871 10 -77625.7258686437 -0.0000175522 11 -77625.7258737697 -0.0000051260 -----------------------------------------------------* Solvation calculation was successfully converged. 64
Ref.; Y. Kato et al., CBI-J. 14 (2014) 1. & J. Yamamoto et al, CBI-J. 14 (2014) 14. FMO-UHF計算 エネルギー微分（MD）も可能 開殻を持つフラグメントモノマーはUHF計算 － SCCまで 2 1 Open ・・・ 3 FMO2 I J N FMO3 開殻モノマーを含む ダイマー群 （近接） 開殻モノマーを含む トリマー群 （近接） 1-3 , 2-3 , 3-I , ・・・ 1-2-3 , 2-3-I , 3-I-J , ・・・ 閉殻のフラグメント群は通常のHF計算 相関エネルギーを求めるMP2/UMP2計算はUHF計算収束後に続けて実行する E FMO2 EIJ ( N 2) EI I J I E FMO3 EIJK ( N 3) EIJ I J K I J DE IJMP2 E IJMP2 E IMP2 E JMP2 ( N 2)( N 3) I EI 2 ~ ~ DEIJHF MP2 DEIJHF DEIJMP2 MP2 MP2 DEIJK EIJK EIMP2 E JMP2 EKMP2 ~ HF MP2 ~ HF MP2 MP2 MP2 DEIJK DEIJK DEIJK DEIJMP2 DEIK DE JK 2024/5/29 65
閉殻MP2のアルゴリズムを踏襲して実装、計3回変換を実行. / <S2> Ref.; B. Schlegel, J. Phys. Chem. 92 (1988) 3075.
FMO-UMP2計算
E
UMP 2
E
αα 2
E
ββ 2
E
1
i a , j b i a , j b ib , j a
E bb 2
2 i j ab
i j a b
S2
S2
E 2
αβ 2
0 S 2 0 S Z S Z nb S i j
UMP 2
ij
S2
UHF
2 aiajb 0 S 2 i aj b
2
1
ia , j b
E b 2
2 i ja b i j a b
2
UHF
ia, jbia, jb ib, ja
1
2 ijab
i
j
a
b
2
aiajb
ia, jb
i j a b
0 S 2 i aj b S a j Si b
i ja b
Compute -correlation energy ! ij-batch / parallel
Compute bb-correlation energy ! ij-batch / parallel
Compute b-correlation energy and <S2> ! i-batch / parallel
Accumulate UMP2 results for FMO summation
スピン汚染には注意が要る（遷移金属イオンの錯体では厳しい）
2024/5/29
66
Ref.; T. Tsukamoto et al., Chem. Phys. Lett. 535 (2012) 157. / MP2-grad. Ref.; Y. Mochizuki et al., Chem. Phys. Lett. 504 (2011) 95. MP2部分構造最適化の例 TrpCageの例 （6-31G基底） 【左 FMO-HF; 右 FMO-MP2 / 緑はNMR(MM)による構造】 ・ファーマコフォア等を選択的に最適化 ・BDAリンクした残基の微分は要計算 ・最適化法はBFGSが標準 ・HFレベルではpp相互作用は記述不可 ・MP2レベルでは実験の構造を保持 ・CH/p相互作用についても対応可能 ・速度や効率などで改良の余地は多い 2024/5/29 g (q A ) Nucl . E Total EIJ ( N 2) EI VAB q A I J q A I q A B A q A Dq B newDg DB BFGS Dg g (q new ) g (q old ) DqDqT B old DgDgT B old T Dq Dg DgT B old Dg 67
Ref.; M. Sato et al., J. Amer. Chem. Soc. 130 (2008) 2396. / FMO-MDの開発は産総研の古明地氏、鹿児島大の石川先生の貢献が大. FMO-MDの例 反応解析に適用した例 H2O + [CH3(N2)]+ 2024/5/29 68
FMO-MDの例（続き1） Tight SN2タイプでの電荷移動の様子：CAFI 5.83 ps 5.87 ps 2024/5/29 5.84 ps 5.88 ps 5.85 ps 5.86 ps 5.89 ps 5.90 ps 69
FMO-MDの例（続き2） Distance (Å) 12.0 10.0 A O-C-N 8.0 6.0 RC-O 4.0 10本の軌跡の2次元反応座標プロット Angle (degree)) Tight SN2タイプでの距離と電荷の変化 Loose type RC-N 2.0 Group charge 0.0 0.8 B qCH3 0.6 0.4 qN2 0.2 0.0 qH2O 必ずしもtight SN2の描像は成り立っていない （反応経路の多様性） 2024/5/29 70
Ref.; D. G. Fedorov et al., J. Chem. Phys. 122 (2005) 054108. 多層FMO (MFMO) ~ ~ E E ' LIi DEIJLi DEIJLi Li I J L j Li I Li J L j ILi I , J Li Layer2 Layer1 Layer2 EI : Monomer energy, EIJ : Dimer energy EI E ' I Tr D I V I , E IJ E ' IJ Tr D IJ V IJ Layer1 DE ' IJ E ' IJ E ' I E ' J , DD IJ D IJ D I D J ~ DE IJ DE ' IJ Tr DD IJ V IJ ・ 層2が高レベル領域 （例えばファーマコフォア） ・ MP2またはMP3を層2に適用 ・ 層1はモノマーSCC条件を課したHFレベルで計算 ・ 3層の扱いも可能 ・ CIS/CIS(D)による励起状態計算も今後可能に 2024/5/29 71
MFMO-MP3計算の例 ・ Chignolinを例にテスト、基底は6-31G ・ 層2には｛Tyr2, Pro4,Trp9｝をセット ・ HF/monomer-SCCは層1と層2を合わせて ・ MP3を層2に適用 ・ 最終のIFIEは層2に関してリスト ## READ NAMELIST "MFMO" : IOSTAT = MethodLow FragLow MethodMiddle FragMiddle MethodHigh FragHigh ........ ## MP3-IFIE 0 = HF = 1,2,3,4,5,6,7,8,9,10 = = = MP3 = 2,4,9 &CNTRL METHOD='MFMO' ReadGeom='Chignolin.pdb' Gradient='NO' Memory=2000 Nprint=3 / &FMOCNTRL FMO='ON' Nbody=2 Np=1 AutoFrag='ON' / &SCF / &BASIS BasisSet='6-31G' / &OPTCNTRL / &MFMO MethodLow='hf' FragLow='1,2,3,4,5,6,7,8,9,10' MethodHigh='MP2' FragHigh='2,4,9' / IJ-PAIR DIST DIMER-ES HF-IFIE MP2-IFIE GRIMME-MP2 MP3-IFIE GRIMME-MP3 PADE[2/1] / A APPROX. / Hartree / Hartree / Hartree / Hartree / Hartree / Hartree ---------------------------------------------------------------------------------------------------4 2 2.509864 F -0.000475 -0.004102 -0.003321 -0.003373 -0.003139 -0.003179 9 2 1.990745 F 0.014420 -0.018544 -0.015018 -0.015330 -0.014215 -0.015128 9 4 3.180956 F 0.003585 -0.001458 -0.001232 -0.001109 -0.001145 -0.000977 2024/5/29 72
Ref.; M. Head-Gordon et al., Chem. Phys. Lett. 219 (1994) 21. CIS(D)のスピン軌道による表式 ; CIS energy , bi a ; CIS vector U2 の寄与 （軌道緩和） ab CIS(D) ab u ij (uij ) 2 1 a a b i vi 4 ijab a b i j ia T2U1 の寄与（差分相関） ab || cj bi ab || ci b j ka || ij bk kb || ij bi c c c 1 b ca a cb b ac vi jk || bc bi a jk b j aik 2b j aik 2 jkbc aij E MP2 2024/5/29 a k a ab b ab || ij a b i j 1 ab aij ij || ab 4 ijab 計算コストを下げるために作業 配列を使って処理する ・ スピン軌道の式をスピン適合配置の式で再定式化 ・ ABINIT-MPにはMP2アルゴリズムを応用して実装 ・ DAXPYとDDOTを多用 （反省点多し） 73
Ref.; Y. Mochizuki et al., Chem. Phys. Lett. 433 (2007) 360. / MFMO-CIS Ref.; Y. Mochizuki et al., Chem. Phys. Lett. 406 (2005) 283. MFMO-CIS(D) Ref.; Y. Mochizuki et al., Theor. Chem. Acc. 117 (2007) 541. DsRedの励起/発光エネルギー算定 励起状態計算はVer. 2 Rev. 8で「復活」の予定 ◇DsRed 珊瑚由来の赤色蛍光タンパク質 （蛍光マーカー） ⇒ 吸収558nm（2.22eV）、発光583nm（2.13eV） Gln66-Tyr67-Gly68の自己反応による色素生成 Tyr67-Gly68 ⇒ Phe65とは-C=N-CO-結合 to Ser69 ◇計算と結果 階層2：“色素＋Phe65+Ser69”まで3種 Gln66 6-31G*基底 （3553原子、220フラグメント） ⇒ 実験値との差は0.1eV以内 吸収 (eV) pig pig+Phe pig+Phe+Ser 極大：実験値 CIS 3.35 3.27 3.26 CIS(D) 2.49 2.30 2.28 2.22 6-31G基底でもほぼ同じ値 2024/5/29 発光 (eV) pig pig+Phe pig+Phe+Ser 極大：実験値 with Phe65 CIS 3.25 3.20 3.18 CIS(D) 2.41 2.21 2.21 2.13 色素部はCIS最適化構造を移植 74
Ref.; N. Taguchi et al., J. Phys. Chem. B 113 (2009) 1153. PR-CIS(Ds) Ref.; Y. Mochizuki et al., Chem. Phys. Lett. mFruits（改変RFP）の励起エネルギー算定 Glu215 逆側へプロトン化 ⇒ 2.27eV -1の電荷 ⇒ 2.33eV （eV） PR-CIS(Ds)ではMP2振幅の部分再規格化と余剰緩和エネルギーが考慮される： 6-31G*基底での計算 2024/5/29 75
Ref.; N. Taguchi et al., Chem. Phys. Lett. 504 (2011) 76. Self-energy Shift Ref.; Y. Mochizuki, Chem. Phys. Lett. 453 (2008) 109 & 472 (2009) 143. BFPとYFPの励起エネルギー算定 (eV) BFP CIS(D) 3.55 PR-CIS(D) 3.43 PR-CIS(D) & SS(2) 3.36 Exptl. 3.21, 3.25 YFP 2.83 2.70 2.53 2.41 ・ 重要な周囲の3残基を含めて計算 ・ GF2補正入りCIS(D)が上手く働く ・ YFPではTyr203の寄与が大きい CIS(D)対補正による算定(6-31G*) 2024/5/29 76
Trp-Cageの中のTrpの励起エネルギー算定#1 最新のVer. 2 Rev. 8で計算可能 --------------------------------------------------------------------------PR-CIS(DS) CORRECTION WITH U2, T2*U1 AND T1*U1 CONTRIBUTIONS AND PR-SCALING --------------------------------------------------------------------------@ MP2 ENERGY OF GROUND STATE BY T2 = PARTIALLY RENORMALIZED VALUE = @ SINGLET RESULT : NUMBER OF VECTORS = -1.916060 -1.876168 NO. 1 CIS ENERGY 0.199111 ( 5.42EV ) U2 ENERGY -0.134334 ( 0.204845 5.57EV ) T2*U1 ENERGY 0.093789 T1*U1 ENERGY -0.001342 CORRECTION -0.041887 0.090764 -0.001283 -0.044853 0.129297 -0.000116 -0.013748 0.125157 -0.000107 -0.017879 T2*U1 ENERGY 0.120964 T1*U1 ENERGY -0.001211 CORRECTION 0.029488 0.117222 -0.001168 0.025789 0.127951 -0.000472 0.020604 0.123916 -0.000449 0.016592 -0.142929 PR @ TRIPLET RESULT : NUMBER OF VECTORS = NO. 1 CIS ENERGY 0.115950 ( 3.16EV ) U2 ENERGY -0.090265 PR 2 0.154984 ( 4.22EV ) PR 2024/5/29 -0.106875 ・1重項励起の実験値は4 eV程度 ・部分再規格化、誘導緩和は奏功 2 PR 2 cc-pVDZ基底 CIS(DS) ENERGY 0.157224 ( 4.28EV ) 0.154258 ( 4.20EV ) 0.191097 ( 5.20EV ) 0.186966 ( 5.09EV ) RATIO 0.79 CIS(DS) ENERGY 0.145438 ( 3.96EV ) 0.141739 ( 3.86EV ) 0.175588 ( 4.78EV ) 0.171576 ( 4.67EV ) RATIO 1.25 0.77 0.93 0.91 2 1.22 1.13 1.11 77
Koopmansの定理とグリーン関数 【Koopmansの定理】 IP t 1 G ( E ) ( E 1 ε ) ⇔ 0 HF近似は零次のGreen関数に相当 【Green関数による補正：Dyson方程式】 IP t N ECorr N 1ECorr (t ) N 1EOR (t ) 差分相関エネルギーと緩和エネルギー 【2次自己エネルギー：スピン軌道による表式】 tu( 2) ( E ) ⇔ G ( E ) G 0 ( E ) G 0 ( E )Σ( E )G ( E ) Σ Σ( 2 ) Σ(3) 自己エネルギーの摂動展開 1 ti || ab ab || ui 1 ta || ij ij || ua 2 iab E i a b 2 ija E a i j 各々、差分相関と緩和を表現 【自己エネルギーの準粒子近似と極強度】 E t tt (E ) 浅い原子価軌道では妥当 d ( E ) Pt 1 tt dE 1 妥当性の目安は0.8 （NR加速により2,3回の反復で収束する） 2024/5/29 78
Table Ref.; Y. Mochizuki., Chem. Phys. Lett. 435 (2008) 109. スケーリング補正した2次グリーン関数 【部分3次の自己エネルギー（P3）：緩和エネルギーの改善】 tt( P 3) ( E ) 1 ti || ab ab || ti 1 ta || ij ij || ta Wtaij U taij ( E ) 2 iab E i a b 2 ija E a i j （ Ortiz, JCP 104(1996)7599） 【2次の自己エネルギーのスピン成分スケーリング】 tt( 2 ) ( E ) iab (ia, tb) (it , ja) ija Molecule HF H2O Orbital 1π 3σ 1b1 3a1 1b2 2 ph S p 2h S (ia, tb) (ib, ta) (ia, tb) E i a b 2 ph T pGW2スケーリング S2 p h Sp 2 h 0.96 T2 p h Tp 2 h 0 (it, ja) (ia, jt ) (it , ja) E a i j （ Huら, J.Elec.Spec.Rel.Phen. 85(1997)39 ） KT 17.67 20.78 13.86 15.87 19.41 p2h T GF2 14.38 18.66 11.23 13.59 17.91 pGW2 15.58 19.51 12.16 14.46 18.62 P3 15.94 19.81 12.45 14.72 18.77 Expt. 16.05 20.00 12.62 14.74 18.51 2次では緩和を過剰評価 （ 単位はeV, 構造はMP2/6-31G*最適化済, IP計算では6-311++G** ） 2024/5/29 79
Ref.; Y. Mochizuki, Chem. Phys. Lett. 453 (2008) 109. CIS(D)の軌道緩和項の修正 * CIS(D)はリドベルグ状態では過小評価 * 2次グリーン関数の自己エネルギーを補正 * 高次の相関効果を実効的に繰込み * CIS(D)SSではリドベルグ状態が改善 * 原子価状態でも精度を保持 * MLFMO-CIS(D)SS として利用可能 * イオン化ポテンシャルの計算も独立に可能 6-311(2+,2+)G** (λ=0 setting) ab CIS(D)ss 1 a a b i vi 4 ijab a b ei e j ia ei i ii (ei ) pp 2024/5/29 ( 2) (uij ) 2 2次の準粒子 pi || ab ab || pi 1 pa || ij ij || pa 1 (e p ) 2 iab e p i a b 2 ija e p a i j 80
Trp-Cageの中のTrpの励起エネルギー算定#2 cc-pVDZ基底 @ LAMBDA DAMPINGS OF CIS ENERGY FOR U2 : 1.00 0.50 0.00 ・pGW2自己エネルギーシフトが有効 ・素のGF2では下がりすぎる ** SELF-ENERGY SHIFT : PGW2 ** @ SINGLET RESULT : NUMBER OF VECTORS = NO. 1 CIS ENERGY 0.199111 ( 5.42EV ) U2 ENERGY -0.140667 PR 2 T2*U1 ENERGY 0.093789 CORRECTION -0.046878 0.090764 -0.049903 -0.130630 PR -0.036841 0.090764 -0.122141 PR 2 ( 0.204845 5.57EV ) -0.148244 PR -0.028352 0.090764 -0.031377 0.129297 -0.018947 0.125157 -0.023087 -0.139067 PR -0.009770 0.125157 -0.131167 PR 2024/5/29 -0.039866 -0.013910 -0.001870 0.125157 -0.006010 CIS(D) ENERGY 0.152233 ( 4.14EV ) 0.149208 ( 4.06EV ) 0.162271 ( 4.42EV ) 0.159245 ( 4.33EV ) 0.170759 ( 4.65EV ) 0.167734 ( 4.56EV ) 0.185898 ( 5.06EV ) 0.181758 ( 4.95EV ) 0.195075 ( 5.31EV ) 0.190935 ( 5.20EV ) 0.202974 ( 5.52EV ) 0.198834 ( 5.41EV ) RATIO 0.76 0.75 0.81 0.80 0.86 0.84 0.91 0.89 0.95 0.93 0.99 0.97 81
タンパク質のIPの計算例 最新のVer. 2 Rev. 8で計算可能 水和させたUbiqutinのTyr部分 FMO-HF&GF2/6-31G*で計算 ----------------------------QP GREEN FUNCTION (IP) RESULT ----------------------------NUMBER OF TARGET OCCUPIED ORBITALS = 19 : THRESHOLD = -0.800 ITERATIVE SOLVING = YES : STATUS = COMPLETE ** KOOPMANS THEOREM ** NO. POS. -IP VALUE 1 21 -0.786556 ( 2 22 -0.703379 ( 3 23 -0.664768 ( 4 24 -0.648260 ( 5 25 -0.645418 ( 6 26 -0.619581 ( 7 27 -0.601884 ( 8 28 -0.596749 ( 9 29 -0.571902 ( 10 30 -0.547152 ( 11 31 -0.532246 ( 12 32 -0.512331 ( 13 33 -0.508945 ( 14 34 -0.500393 ( 15 35 -0.483462 ( 16 36 -0.388602 ( 17 37 -0.380923 ( 18 38 -0.347929 ( 19 39 -0.336190 ( 2024/5/29 -21.403 ) -19.140 ) -18.089 ) -17.640 ) -17.563 ) -16.860 ) -16.378 ) -16.238 ) -15.562 ) -14.889 ) -14.483 ) -13.941 ) -13.849 ) -13.616 ) -13.156 ) -10.574 ) -10.365 ) -9.468 ) -9.148 ) ** NORMAL SECOND-ORDER (GF2) RESULT / OK FLAG = 1 ** NO. POS. -IP VALUE 1 21 -0.686201 ( 2 22 -0.601090 ( 3 23 -0.574999 ( 4 24 -0.550784 ( 5 25 -0.550918 ( 6 26 -0.518404 ( 7 27 -0.507267 ( 8 28 -0.508080 ( 9 29 -0.484631 ( 10 30 -0.468994 ( 11 31 -0.453970 ( 12 32 -0.429817 ( 13 33 -0.436472 ( 14 34 -0.420589 ( 15 35 -0.406054 ( 16 36 -0.309469 ( 17 37 -0.284350 ( 18 38 -0.330536 ( 19 39 -0.320539 ( -18.673 ) -16.356 ) -15.647 ) -14.988 ) -14.991 ) -14.106 ) -13.803 ) -13.826 ) -13.188 ) -12.762 ) -12.353 ) -11.696 ) -11.877 ) -11.445 ) -11.049 ) -8.421 ) -7.738 ) -8.994 ) -8.722 ) PL.ST. OK IT. 0.8646 1 4 0.8373 1 4 0.8674 1 3 0.8635 1 3 0.8600 1 3 0.8484 1 4 0.8634 1 3 0.8651 1 3 0.8631 1 3 0.8832 1 3 0.8863 1 3 0.8708 1 3 0.8407 1 3 0.8732 1 3 0.8669 1 3 0.8780 1 3 0.8799 1 3 0.8840 1 3 0.8832 1 3 ** PARAMETERIZED GW2 RESULT / OK FLAG = 1 / OSPPH, PSPPH, OSPHH, NO. POS. -IP VALUE 1 21 -0.724863 ( 2 22 -0.642948 ( 3 23 -0.614015 ( 4 24 -0.589644 ( 5 25 -0.590154 ( 6 26 -0.563164 ( 7 27 -0.547481 ( 8 28 -0.547838 ( 9 29 -0.524784 ( 10 30 -0.504271 ( 11 31 -0.490140 ( 12 32 -0.461836 ( 13 33 -0.459240 ( 14 34 -0.460347 ( 15 35 -0.442530 ( 16 36 -0.341051 ( 17 37 -0.323869 ( 18 38 -0.338620 ( 19 39 -0.328439 ( PL.ST. OK IT. -19.725 ) 0.9044 1 3 -17.496 ) 0.8878 1 3 -16.708 ) 0.9116 1 3 -16.045 ) 0.9044 1 3 -16.059 ) 0.9040 1 3 -15.324 ) 0.8993 1 3 -14.898 ) 0.9065 1 3 -14.907 ) 0.9107 1 3 -14.280 ) 0.9105 1 3 -13.722 ) 0.9230 1 3 -13.337 ) 0.9253 1 3 -12.567 ) 0.9072 1 3 -12.497 ) 0.8680 1 3 -12.527 ) 0.9201 1 3 -12.042 ) 0.9102 1 3 -9.280 ) 0.9163 1 3 -8.813 ) 0.9193 1 3 -9.214 ) 0.9206 1 3 -8.937 ) 0.9206 1 3 82
Ref.; Y. Mochizuki et al., Chem. Phys. Lett. 418 (2006) 418. 動的分極率の計算 Dodecaneの計算例（cc-pVDZ） ・CPHF/LRによる評価 ・AO-drivenな実装、flat MPIで並列化 ・FMOの2体展開式で算定、精度は確認済み ・第二次高調波発生（SHG）の実装を検討中 ・OpenMP/MPI並列化も対応予定 V B N A 1 A 1 N B V Aia , jb 2(ia | jb) (ij | ab) ij ab ( a i ) Bia , jb 2(ia | jb) (ib | ja ) 2( N N )V XY X (M 2) XY 2024/5/29 X 83
Table Ref.; Y. Mochizuki., Chem. Phys. Lett. 435 (2008) 109. 分極率計算のA行列を軌道エネルギーシフト ABINIT-MPでの実装を検討中 pGW2 LRSS 2 LRSS Self-energy LR P3 LRSS CCSD MP2 Methane αxx = αyy = αzz 14.92 15.60 15.24 15.46 15.22 15.07 11.64 11.90 12.47 13.35 12.09 12.76 12.08 12.68 12.07 13.72 11.94 13.27 6.71 8.48 7.33 7.78 9.23 8.15 7.34 8.89 7.80 7.19 8.79 7.68 7.93 8.89 7.97 7.63 8.81 7.83 3.06 5.31 3.50 5.84 3.32 5.62 3.26 5.54 3.44 5.71 3.35 5.67 Ammonia αxx = αyy αzz Water αxx αyy αzz Hydrogen fluoride αxx = αyy αzz （ 単位はau, 構造はMP2/6-31G*最適化済, 基底は6-311+G(2d,2p)： MP2とCCSDの値はG03による ） 相関を導入したMP2やCCSDの値に近い結果が得られる!! 2024/5/29 84
Ref.; 望月ら, J. Comp. Chem. Jpn. 16 (2017) 119. 機械学習用のデータのダンプ UbiqutinのRFによるIFIE予測例 ・統計的なIFIE/PIEDA算定が今後の主流となる見込み ・機械学習処理の普及の流れ （プロキシー処理も意識） ・記述子を含む本質的なデータのみをコンパクトにダンプ （Ver. 1系までの可視化重視のCPFと異なるコンセプト） ・pythonでの処理を意識したテキストデータの並び （保存時の圧縮も容易） ・Logファイルからのデータ抽出のpythonスクリプトもアリ Gly5のFMO-MP2/cc-pVDZのデータ 2024/5/29 85
Ref.; S. Matsuoka et al., J. Comp. Chem. 45 (2024) 898. PIEDA機能の最近の強化 2024/5/29 86
LRD Ref.; Sato et al., J. Chem. Phys. 131 (2009) 224104. / 数値積分モジュールは石村氏のSMASHのモジュールを利用。 局所応答分散による補正(LRD：HF+D) 𝐴 𝐵 𝐴𝐵 𝐸𝑑𝑖𝑠𝑝 𝜌 = − 𝑚 𝑚 𝐴 𝑛𝐵 𝑉𝐴𝐵 0 𝐴 0𝐵 2 𝑛 𝑛 𝐸𝑑𝑖𝑠𝑝 𝜌 = − 𝐶𝑛𝑎𝑏 𝜌 𝑅𝑎𝑏 𝑓𝑑𝑎𝑚𝑝 𝑅𝑎𝑏 𝐴 𝜔𝑚 + 𝜔𝑛𝐵 𝑛 𝑎>𝑏 𝑛≥6 分散力の2次摂動的表現式 • • • 局所応答分散近似 ベッケの分割関数 適切なダンピング Ar-Arのポテンシャルエネルギー曲線のテスト 6-31+G* 0.2 0.2 0.15 0.15 Binding Energy (kcal/mol) Binding Energy (kcal/mol) 6-31G* 0.1 0.05 0 -0.05 -0.1 HF -0.15 HF+D(6) -0.2 HF+D(8) -0.25 HF+D(10) 0.05 0 -0.05 -0.1 HF -0.15 HF+D(6) -0.2 HF+D(8) -0.25 HF+D(10) -0.3 -0.3 3 3.5 4 4.5 5 R (Å) 2024/5/29 0.1 5.5 6 6.5 7 3 3.5 4 4.5 5 5.5 6 6.5 7 R (Å) （MP2の最小値は-0.21 kcal/mol @ 4.06Å） 87
LRDを入れたPIEDAの計算例 6-31G*基底 Pair Type Dist. IFIE ES EX CT+mix DI(LRD) DI(Rest) % LRD 1G - 2A Stacking 3.27 -10.59 -3.48 3.52 -2.08 -7.36 -1.19 86.1 1G - 1'C H-bond 1.80 -45.47 -59.84 33.75 -10.65 -3.05 -5.68 34.9 2A - 2'T H-bond 1.85 -18.53 -23.30 16.31 -5.58 -2.42 -3.54 40.6 1'C - 2'T Stacking 2.88 -4.66 -1.10 4.28 -1.75 -5.32 -0.77 87.3 Trp6 - Tyr3 p/p 2.20 -11.66 -7.56 6.58 -3.11 -5.41 -2.16 71.5 Trp6 - Pro12 CH/p 2.41 -6.49 -6.33 3.69 1.36 -4.32 -0.89 82.9 Trp6 - Pro18 CH/p 2.70 -6.84 -2.34 6.08 -3.18 -6.03 -1.37 81.5 Trp6 - Pro19 CH/p 2.59 -3.17 -2.65 2.2 1.07 -3.39 -0.4 89.4 GA:CT Trp-Cage 2024/5/29 88
ES(RESP)を使ったPIEDAの計算例 Pair IFIE Dist. IFIE(MM) ES ES(RESP) ES(MM) 1G - 2A -10.59 3.27 -10.13 -3.48 -2.14 -2.19 1G - 1'C -45.47 1.80 -24.78 -59.84 -27.94 -26.37 2A - 2'T -18.53 1.85 -12.84 -23.30 -9.46 -11.99 1'C - 2'T -4.66 2.88 -4.06 -1.10 -0.18 0.63 Trp6 - Tyr3 -11.66 2.20 -9.60 -7.56 -7.80 -5.36 Trp6 - Pro12 -6.49 2.41 -4.93 -6.33 -3.01 -1.62 Trp6 - Pro18 -6.84 2.70 -5.84 -2.34 -4.34 -1.25 Trp6 - Pro19 -3.17 2.59 -3.47 -2.65 -1.36 -0.62 GA:CT 水素結合 部分のES は過大値 Trp-Cage 𝑞𝑖 𝑞𝑗 𝑁𝐴 𝑞𝑖 𝑞𝑗 ∆𝐸𝑖𝑗 = ≈ 332.06 × 4𝜋𝜀0 𝑟𝑖𝑗 × 10−10 𝑟𝑖𝑗 (kcal/mol) ・ PIEDAのES項の過大評価問題に対する簡易対策 ・ FMO-ESPを再現するRESP電荷で古典的な式で算定 ・ あくまでも参考値としての使い方を推奨 2024/5/29 89
GPU関係（進行中） （謝意： 東京大学情報基盤センターのGPU環境移行プログラムの支援） 2022年度下期： HFミニアプリの改造はNVIDIA成瀬彰氏によります 2024/5/29 90
GPU対応の必要性 ■背景 ・GPUスパコンの遍在化 ⇒ 「Wisteria」 Aquarius&Mercury、「不老」 Type II、Oakforest PACS-II... ⇒ 使われている分野はAI、量子シミュレーション、MD… ・ネクスト「富岳」 ⇒ 「アクセラレータ」の搭載は（ほぼ）確定 ・MO系の計算がGPUに対応済 ⇒ 高名なTeraChem、GAMESS-USも対応（GPU化FMO専用コードも在り!） ⇒ 国内でも先行例在り（安田ら、梅田ら、鬼頭ら） ■対応策 ・ABINIT-MPも2電子積分生成などのGPU化が必要 ⇒ FMO-HF計算をまず改良 ⇒ MP2計算はDGEMM化で相対的には対応は容易と想定（?!） ・ABINIT-MP本体は巨大に過ぎる ⇒ HF計算のみを切り出したミニアプリを作成 ⇒ NVIDIAの成瀬彰氏のチューニングで全体で6-8倍の加速を確認 ・ABINIT-MPのGPU化を検討中 積分ルーチン群の再構成（81種→21種へ低減）から着手 2024/5/29⇒ 91
ミニアプリの作成はFOCUSの坂倉耕太氏、GPU化の改造はNVIDIAの成瀬彰氏による. ABINIT-MPのミニアプリ 2022年度下期の東大情基セの支援PJに参加 ・ FMOなしでHF計算のみを実行 ・ 2電子積分の性能評価とチューニング用 ・ ネームリストで積分タイプを使い分け可能 → ObaraはABINIT-MPと同じVRR法で 軌道の組み合わせ毎のサブルーチン → Obara-generalは汎用ループでの 計算、最適化は特にされていない → HGP-generalはHRR法での汎用ループ での計算、最適化は未 → Obara-vectorはNEC SX向けにループ 構造をベクトル計算用に変更してある ・ Obara-vectorは、GPUでの処理に向く?! ・ 演算数的にはHRRにも可能性あり ・ コンパクトなルーチンで性能が出れば最良… ・ 4タイプでの評価をお願いできれば… → 筋のよさそうなものをチューニング → 基底関数は6-31G*が基本 2024/5/29 &CNTRL ReadGeom='adamantane.pdb' eri_type=0 ! Obara ! eri_type=1 ! Obara-general ! eri_type=2 ! HGP-general ! eri_type=3 ! Obara-vector ! memory=4000 / &SCF !maxscfcyc=10 / &BASIS BasisSet='6-31G*' ! BasisSet='6-31G' ! BasisSet=‘cc-pVDZ*' / 92
Ref.; 望月ら, J. Comp. Chem. Jpn. 23 (2024) 4. ミニアプリのチューニング（NVIDIA成瀬彰氏） 成瀬氏（NVIDIA） ◇ 先ず、ABINIT-MPの通常の積分生成と同じ、小原・雑賀のVRR式のルーチン群を対象に → 元は{x, y, z, w}={s, p, d}関数の組み合わせからsub_xyzwは81種類 （自動生成） ◇ 多数の sub_xyzw をバッチ化して、GPU の多並列性を活用する → 以下の理由から単純に全部をバッチ化すればいいわけではない ・ sub_ssss：計算手順が簡単、一時配列数が少なく、そのサイズは小さい （レジスタ使用量：少）： データ並列性が高い (IJKL_ex が大きめ) ・ sub_dddd：計算手順が複雑、多数の一時配列が必要で、そのサイズは大きい （レジスタ使用量：多）： データ並列性が低い (IJKL_ex が小さい) → 計算特性の近いものをまとめる（構造の改変）、OpenACC指示文の挿入 → Fock行列への加算も、バッチ化積分ルーチンの内部で実施 （メモリアクセス量削減） ・ 81種を26種に再構成して削減 ・ 緑色が GPU 化の対象 ・ “s3d1”等は4通りを含む ・ “s1d3”、“p3d1”、”p2d2”等はレジスタ 要求量が多く、レジスタスピル発生 ・ “p1d3”、“d4” は CPU 実行 2024/5/29 93
Phe-TrpとTrpでのミニアプリのテスト（HF） 東大Wisteria-Aquariusを使用 ・ 下記の条件で再度実行(flat MPI) export NNODES=1 export NGPUS=1 export NPROCS=8 ・ Trp-Pheのテスト 6-31G (#Orb = 289, #SCF = 17) ; 27.9 s / 4.3 s → Acc = 6.5 6-31G* (#Orb = 451, #SCF = 18) ; 69.4 s / 10.8 s → Acc = 6.4 cc-pVDZ (#Orb = 520, #SCF = 19) ; 179.0 s / 21.8 s → Acc = 8.2 6-311G** (#Orb = 651, #SCF = 21) ; 257.4 s / 37.7 s → Acc = 6.8 ・ Trpのテスト 6-31G (#Orb = 159, #SCF = 17) ; 8.4 s / 1.8 s → Acc = 4.7 6-31G* (#Orb = 249, #SCF = 18) ; 21.9 s / 3.6 s → Acc = 6.1 cc-pVDZ (#Orb = 285, #SCF = 17) ; 51.4 s / 6.2 s → Acc = 8.3 6-311G** (#Orb = 357, #SCF = 21) ; 76.4 s / 10.5 s → Acc = 7.3 ・ GPUは「加速装置」としての価値アリと判断 2024/5/29 94
坂倉ら, 第29回計算工学講演会(神戸)で発表予定(2024/6/12); 要旨は2024年末に公開見込み. ABINIT-MP本体のGPU対応#0 坂倉氏（FOCUS） 積分の分類 2024/5/29 ループのスケッチ Type Unique quartet do scc-iteration [s4] (ss|ss) do 1,nf <MPI> [s3p1] (ss|sp) do scf-iteration [s3d1] (ss|sd) do one-int-loop [s2p2] (ss|pp), (sp|sp) do eri-loop with GPU [s2p1d1] (ss|pd), (sp|sd) do eri-loop (dddp,dddd) [s2d2] (ss|dd), (sd|sd) enddo [s1p3] (sp|pp) enddo [s1p2d1] (sp|pd), (sd|pp) [s1p1d2] (sp|dd), (sd|pd) [s1d3] (sd|dd) do dimer-loop <MPI> [p4] (pp|pp) do scf-iteration [p3d1] (pp|pd) do one-int-loop [p2d2] (pp|dd), (pd|pd) do eri-loop with GPU [p1d3] (pd|dd) do eri-loop (dddp,dddd) [d4] (dd|dd) enddo Count 21 enddo <mp2 processing> <mp2 processing> 95
ABINIT-MP本体のGPU対応#1 坂倉氏（FOCUS） 改造前のループ 2024/5/29 改造後のループ 96
ABINIT-MP本体のGPU対応#2 坂倉氏（FOCUS） OpenACCの 挿入例 （データ転送& 並列化） 加速例（暫定） 2024/5/29 !$acc data copyin(x,y,z,dc,pps, nsh2exp) copy(fc) do ia = 0,num_angtype-1 !$acc parallel num_workers(2) vector_length(16) !$acc loop gang worker do n4 = 1, nsize ish = nslist(1,n4) jsh = nslist(2,n4) ksh = nslist(3,n4) lsh = nslist(4,n4) 積分タイプ CPU(s) GPU+(s) Acc. (ss|ss) 59.1 0.9 65.7 (ps|ss) 148.0 3.3 44.9 (ps|ps) 112.7 4.6 25.0 (pp|ss) 58.3 2.4 24.3 (pp|ps) 120.3 3.3 36.5 (pp|pp) 42.9 1.2 36.0 Chignolin: HF/6-31G （10個のモノマーでの 積分計算部分の和） 「1CPUのみの場合と1CPU+1GPUの比較」 97
まとめ 2024/5/29 98
ABINIT-MPによるFMO計算 （その１） ◇ABINIT-MPの機能 開発から20年ほど経過、現在はOpen Ver. 2系へ、最新版はRev. 8 相互作用エネルギー解析には必要十分、FMO-MP2はルーチン実行が可能 FMO4やCIS/CIS(D)などの機能もあり ◇スパコンとの関係 OpenMP/MPI混成並列、大規模MP2計算が容易 FMO-MP3も相対コスト低く実行可 （「富岳」などの例は次回） HPCI拠点にライブラリプログラムとして整備中 GPU対応は「助走」の段階 ◇今後の方向性 高速化と超巨大系（数万フラグメント）への対応 （詳しくは次回） 機能強化、機械学習との親和性 生物物理や創薬以外の分野（マルチスケールシミュレーション等）への展開 長時間のご聴講、ありがとうございました 2024/5/29 99
ABINIT-MPによるFMO計算のロードマップ 2004年 2010年 2008年 2015年 2013年 2020年 核内受容体(ER) ~300残基 粗視化MD用 結晶－ペプチド複合系 パラメータ 新型コロナウイルス インフルエンザHA インフルエンザHA３量体 ~(SiO ) -6残基-水和 ~数万サンプル 抗原抗体系~5300残基 2 250 抗原抗体系~1000残基 抗原抗体系~2400残基 mFruits 水和DNA EGFR インフルエンザNA 12塩基対＋2500wtr チロシンキナーゼ タミフル~400残基 計算 MP2 CIS/CIS(D) 構造 PDB一点計算/ モデル埋戻し MP3 CCSD(T) FMO4 CD FMO-MD 解析 IFIE CAFI FILM 電荷、溶媒効果 2024/5/29 リガンド水和 10Å水和層 粗視化→原子復元構造 ~1万原子×サンプル数 Dimer-ES CMM MP2(p-opt) LRD FMO-DPD PIEDA BSSE FMO4-IFIE ESP/RESP 分子固体 ~千個単位 NPA SCIFIE SVD PB(SA) MD生成 多構造 統計/ML 大型液滴 100