1E01
超並列フラグメント分子軌道法プログラム OpenFMO の開発（1 九大, 2 九州先端研, 3 筑波大, 4 東工大, 5 産総研）○稲富雄一 1, 眞木淳 2,
本田宏明 1, 石元孝佳 1, 梅田宏明 3, 渡邊寿雄 4, 池上努 5, 青柳睦 1

Development of massively parallel program of fragment molecularorbital method, OpenFMO( 1 Kyushu Univ.,

2 ISIT,

3 Univ. of Tsukuba,

4 Titech,

5 AIST)○Yuichi Inadomi1, Jun Maki2, Hiroaki Honda1, Takayoshi Ishimoto1,
Hiroaki Umeda3, Toshio Watanabe4, Tsutomu Ikegami5, Mutsumi Aoyagi1【はじめに】フラグメント分子軌道(FMO)法はタンパク質や DNA などの大規模生体分子の電子状態計算を高速に行うために開発された計算手法である．巨大分子を小さなフラグメントに分割して，フラグメント（モノマー）および，フラグメントペア（ダイマー）に対する小規模な電子状態計算を行うことで，分子全体の電子状態を近似する．各モノマー（ダイマー）に対する電子状態計算を独立に計算できること，ならびに，各小規模電子状態計算自身もさらに並列処理可能であることなどから，大規模並列処理向きの計算手法であると考えられており，数 1000 並列での FMO 計算を効率よく行えることも示されている．しかしながら，数万～数 10 万といった超並列 FMO 計算を効率よく行う
ためには，更なる並列性能向上が必要となる．また，FMO 計算をより高速に実行できる環境を整えるためには，今後開発が予定されているエクサスケールスパコンで想定されるアーキテクチャ向けのプログラム開発も必須となる．我々は，効率的な超並列 FMO 計算を目指して並列 FMO プログラム
OpenFMO の開発を行っており，高速化のための様々な最適化を行ってきた．また，今後出現が予想される計算機アーキテクチャを踏まえたコード開発を計算機科学の研究者と共同して行っているところである．本発表では，これまでの成果，および，現在行っている研究内容について報告する．
【OpenFMO とは】OpenFMO は九大と九州先端研でスクラッチから開発され C 言語で記述された
並列 FMO プログラムである．Hartree-Fock 法に基づいた FMO 計算に特化しているため，ソースコードが 5 万行強と非常にコンパクトで，最適化などを比較的容易に適用することができる．また，並．．．．
列化において MPI や OpenMP という事実上の標準となっている並列化手法を用いていること，ならびに，特定の計算機ベンダーが提供する特殊な機能を使っていないことなどから，本コードは，ほとんどの並列計算機環境で利用可能である．
【最適化１ OpenMP/MPI ハイブリッド並列化】京（コンピュータ）をはじめとする最近の大型計算機は，数～数 10 プロセッサ（コア）を搭載した多数の小型計算機（ノード）を Infiniband などの高速ネットワークで相互結合したクラスタ型並列計算機アーキテクチャを持つ．そのような計算機で効率よく動作するプログラムの開発においては，同一ノード上にあるコア間の通信（ノード内通信）とノードを跨ぐ通信（ノード間通信）の通信性能差を考慮した並列化が必要となる．そのために，我々は，ノード内並列化には OpenMP によるスレッド並列化を，また，ノード間並列化には MPI によるプロセス並列化をそれぞれ用いた OpenMP/MPI ハイブリッド並列化を分子積分コードに適用した．
この並列化手法は，京における並列化手法としても推奨されている．
【最適化２ OpenMP/MPI によるグローバルカウンタの実装】FMO 計算だけでなく一般の量子化学計算におけるカーネルコード（計算量の大きい部分）の１つである分子積分を多数のプロセスを用いて効率的に並列処理するためには，グローバルカウンタ(GC)を用いた動的負荷分散の適用が有効である．GC は並列処理を行っている全プロセスが参照可能なカウンタであり，GC にアクセスしてカウンタ値を取得し，その値に基づいて積分計算を行い，それが終了したら，次のカウンタ値を取得して，...を繰り返すことで積分計算を各プロセスに分配して並列処理を行う．GAMESS や NWChemといった汎用量子化学計算プログラムで用いられている GC の実装には ARMCI ライブラリが用いられているが，我々のターゲットマシンの 1 つである京で動作する ARMCI ライブラリは京の各ノード
上の 8 コアのうち 1 コアをライブラリが占有する（計算に使えない）実装となっているため，その利用は効率的な並列処理の妨げとなる．そこで我々は，ARMCI の代わりに OpenMP と MPI を用いて計算機資源の無駄の少ない GC を実装した．
【最適化３ 中間データの分散保存・アクセス方法の改良】FMO 計算では，並列計算を行っているすべてのプロセスからモノマー密度行列データの参照，および，更新が可能である必要がある．このモ
ノマー密度行列データは入力データの規模に比例するため，巨大な入力に対する効率的な FMO 計算を可能にするためには，データを複数のノードに分散保存して，分散保存されたデータに効率よくアクセスする必要がある．この目的を達成するために，今回，（１）MPI-2 規格の MPI ライブラリで実装されている片側通信機能を使う手法，および，（２）データ保存を行う専用プロセス（データサー
バ）を用いる手法，の２つを実装して，その性能評価を行った．その結果，片側通信を用いる方法では MPI-2 での片側通信機能の実装の影響で，密度行列データへのアクセス性能が非常に悪くなることが分かった．一方で，データサーバを用いる方法では，データアクセスによる遅延もほとんどなく，非常に効率よく密度行列データにアクセスできることが分かった．
【性能評価】 最適化を適用したコードを用いて，並列性能の評価を行った．用いた計算機は九州大学情報基盤研究開発センタ
ーの富士通 PRIMERGY CX400（8-core
Xeon×2/node, 1476 nodes,
interconnect=Infiniband FDR）である．
計算に用いるコア数を 1280～20480 まで変化させてアデノウィルス KNOB ドメイン
（PDB ID=1nob, モノマー数 576）に対する FMO 計算を行い，その計算時間を測定した（右図参照）．その結果，FMO 計算全体では 2 万並列でも非常に効率よく並列処

図：大規模並列時の性能向上比（1280 並列時を１とした場合）理できることが分かった．特に，dimer SCF 計算，および，近似ダイマー(ES dimer)計算部分は，2万並列でも，ほぼ理想的な速度向上を示した．一方で，モノマーSCF(SCC)計算部分は，2 万並列時に並列化効率が低下していることが分かる．これは，各プロセスに割り当てられる計算量が小さくなったため，プロセッサ間の負荷バランスがとりにくくなったことが原因であると考えられる．更なる効率の向上のためには，SCC 部分の並列化効率を上げるための工夫が必要となる．
【検討中の課題】現在，高性能な GC，および，データサーバの実装を容易にするためのミドルウェアの開発を九大と富士通の共同研究，および，産総研で行っている．また，次世代スパコンへの搭載が確実視されている演算加速器を用いた分子積分計算プログラムについても，九大，筑波大，東工大の共同で開発している．また，計算中に計算機の一部が故障しても，計算を続行できる（耐障害性を持つ）FMO プログラムの開発を，九大と筑波大，および，理研との共同研究で行っている．
謝辞 本研究は，科学研究費補助金基盤研究（C）「超並列フラグメント分子軌道法プログラムライブ
ラリの開発」（課題番号 22550015）， JST,CREST の研究領域「ポストペタスケール高性能計算に資するシステムソフトウェア技術の創出」の研究課題「省メモリ技術と動的最適化技術によるスケーラ
ブル通信ライブラリの開発」，ならびに，研究課題「ポストペタスケールシステムのための電力マネージメントフレームワークの開発」の支援を受けている．

1E02
コレスキー分解を用いたグリッドフリー密度汎関数計算の研究(東大生研*, 東大院工**) ◯平野敏行*、樋口恒**、佐藤 文俊*
Study of grid-free density functional calculation using Cholesky decomposition(Institute of Industrial Science, the University of Tokyo*,
School of Engineering, the University of Tokyo**)◯Toshiyuki HIRANO*, Hisashi HIGUCHI**, Fumitoshi Sato*【緒言】
	  我々は逆行列計算を行わない独自のコレスキー分解(CD)法により、分散メモリ型並列計算機でも安定かつ高速に計算可能な、クーロン・交換項積分計算法を開発した[1]。
本研究では、この計算法と相性の良い、密度汎関数(DF)法に基づく新たな高精度大規
模カノニカル分子軌道計算法を提案する。DF 法における交換相関項計算では、交換相関汎関数の複雑さから、グリッドを用いた数値計算が広く採用されている。しかし、数値計算による数値不安定性は避けられず、また、大規模分子の計算を分散メモリ型
超並列計算機で実行する場合には、巨大密度行列を用いた全グリッド上の密度・密度勾配を計算するタスクの均一化が極めて困難となる。そこで、交換相関項計算にグリッドフリー法[2]を適用することで、高速かつ解析的に交換相関項を求めることを目的とした。なお、これまでのグリッドフリー法の計算精度は決して満足のいくものではなかった。本研究では、まず計算精度が良くない原因を突き止め、改善法を提案し、さらに CD 法を用いることにより、高速かつ計算精度が制御できる交換相関項計算法を開発することに成功した。
【方法】
	  解析的に交換相関項を求めるグリッドフリー法において，交換相関項𝑴 𝑓 𝜌 は密度の行列表現𝑴 𝜌 を用いて以下のように求められる。
𝑀 𝜌 !" =!" 𝑃!" 𝑝𝑞𝑟𝑠

(1)

𝑴 𝜌 = 𝑽! 𝑴 𝜌 𝑽(2)

𝑴 𝜌 = 𝑼𝝀𝑼!

(3)

= 𝑼𝑓 𝝀 𝑼!

(4)= 𝑺𝑽𝑼𝑓 𝝀 𝑺𝑽𝑼 !

(5)

𝑴𝑓 𝜌
𝑴𝑓 𝜌

ここで、𝑷は密度行列、 𝑝𝑞𝑟𝑠 は 4 中心重なり積分、𝑽は直交化行列、𝑺は重なり積分
行列、 𝑝𝑞𝑟𝑠 は 4 中心重なり積分である。この計算法を詳細に検討した結果、(1)式では𝑴 𝜌 の質が LCAO 展開に用いる基底関数依存となるため、グリッドフリー法の計算精度に影響を与えていることを明らかにした。

そこで本研究では、密度の行列表現𝑴 𝜌 に専用の基底関数セットを用い、重なり積!
分行列𝑺の代わりに𝑆!"= 𝑝𝛼 を導入することでこの問題を解決した。すなわち、(1)と(5)式をそれぞれ以下の 2 式に置き換えることで、LCAO 展開に用いる基底関数セットよりもリッチな基底関数セットを交換相関計算に用いることができる。
𝑀 𝜌 !" =𝑴𝑓 𝜌!" 𝑃!" 𝛼𝛽𝑟𝑠

(6)

= 𝑺! 𝑽𝑼𝑓 𝝀 𝑺! 𝑽𝑼 !

(7)	  次に、(1)式も(6)式も 4 中心重なり積分は、基底関数の総数に対して 4 乗のサイズ依存性を持つため、大規模分子ではボトルネックとなる。既報の方法[2]では、(1)式の加速に RI 法を利用していたが、計算精度を制御することが難しい。本研究では、2電子反発積分でも利用している CD 法を適用することで、計算精度を制御しつつ計算効率の向上を目指した。以下に示す CDAM 法[3]とピボット CD 法[4]に基づく CD 法を用いることでコレスキーベクトルのサイズを縮小し、計算コストを削減した。
𝑉!",!" = 𝛼𝑟𝛽𝑠 ≈𝑀𝜌 =!" 𝑃!" 𝛼𝑟𝛽𝑠

≈

!

(8)

𝐿!,!" 𝐿!,!"
!

!

(9)

𝑋!,!" 𝑋!,!"

𝑋! =! 𝑄!" 𝐿!,!"

(10)

𝑃!" =

! 𝑄!" 𝑄!"

(11)

本研究で提案したグリッドフリー方法は、我々の CD 法に基づくクーロン・交換項
計算法[1]と本質的に同じであり、SCF 繰り返し計算前にコレスキーベクトル𝐿を求めることで、密度の表現行列𝑴 𝜌 がシンプルな行列演算のみで得られるといった性質も共有できる。これらの計算法を組み合わせることで、DF 計算全体に大きな利点が生まれる。すなわち、SCF 計算中に 4 中心重なり積分計算を求める必要が一切無いうえ、クーロン項、Fock 交換項、ならびに DF 法の交換相関項𝑴 𝑓 𝜌 の計算((2)~(4), (7)式)は、ScaLAPACK ライブラリ[5]を利用できる。行列の疎性を仮定した近似を用いることなく、巨大分子の DF 法に基づく標準的なカノニカル分子軌道計算が、分散メ
モリ型超並列計算機システムであっても最適に実行できるということである。
	  本方法を ProteinDF に実装し、いくつかのテスト計算を行った。講演では、本方法の性能ならびに計算例を紹介する。
【参考文献】[1] 平野敏行, 佐藤文俊, 第 6 回分子科学討論会 2012, 東京, 4P096[2] Y. C. Zheng and Jan Almlöf, Chem. Phys. Lett. 214, 397 (1993).
[3] Y. Okiyama, T. Nakano, K. Yamashita, Y. Mochizuki, N. Taguchi, S. Tanaka,

Chem. Phys. Lett., 490, 84 (2010).
[4] H. Harbrecht, M. Peters, R. Schneider, Appl. Num. Math., 62, 428 (2012).
[5] http://www.netlib.org/scalapack/1E03
マテリアル量子化学に向けた周期境界条件下での ab-initio 計算手法の開発(理研 AICS) ○島崎智実、小杉太一、中嶋隆人

Development of ab initio calculation method under periodic boundary
conditions for material quantum chemistry(RIKEN-AICS) Tomomi Shimazaki, Taichi Kosugi, and Takahito Nakajima
序論
近年、分子だけでなく表面や固体といったマテリアル分野においても量子化学手法の活用が期待されている。例えば、エネルギー問題や環境問題を解決するために，太陽電池や固体触媒のための新規マテリアルの開発を加速することが強く望まれている。このような分野のシミュレーションを行うには、従来の量子化学が主な対象としてきた孤立分子系だけでなく、表面や固体、マテリアル（材料）の化学性質を電子状態の観点から詳細に調べていくことが必要である。そこで、本研究ではマテリアル研究のために有用な周期境界条件下での ab initio
計算手法の開発を行った。特に、本手法の中では Hartree 項の高速化に注目して研究を行った。
計算手法
本手法では、周期境界条件下での Hartree 項を誤差関数および補誤差関数を用いて、短距離
成分と長距離成分の寄与に分割する。長距離成分に関しては、ガウス関数を補助基底として電子密度を展開する。本手法で用いた電荷密度に対してフィッティングを行った長距離成分の Hartree 項は次のように表される。
long-range

V H artree Q      0  r1   Q  r1 

erf   r12 

Q1



Q1

r12

 r2  d r1 d r2N
erf   r12  Q 1
0Q  C s       r1     r1  s  r2  d r1 d r2 
r12
s
 Q1


ここで、  はガウス型の基底、  はガウス型の補助基底である。これらの基底と電子密度との間には次のような関係がある。
 Q  r    Q   D   Q D   Q  r  

Q Q D r    s C s sQ  r    Q  r 

上式の C s はフィッティング密度を表す補助ガウス関数の展開係数である。ユニットセルの電荷的中性条件を確実に中性に保つために、この展開係数 C s を決めるためにラグランジュの未定乗数法を用いた。また、短距離成分に関しては、次のような 4 中心積分から求めた。
short-range

VHartree Q     D   Q D    0  r1   Q  r1 
Q 1 , Q D 

erfc   r12 
r12



Q1

 r2  Q  Q  r2  d r1 d r2
1

D

計算結果
本手法を用いた計算結果を下の図に示す。図はグラフェンのモノレイヤー層に関して
B3LYP/6-21G*の電子密度汎関数および基底を用いた計算結果から得られたバンド図である。
タイトバインディング近似による理論的解析によると、グラフェンは K (Dirac) 点でゼロギャ
ップの半導体であり、ab initio の計算結果も同様にフェルミエネルギーにおいてゼロギャップの性質を示すことが確認できる。

図 1 B3LYP 法によって計算されたグラフェンのバンド図

また、表 1 には異なった計算手法で求めた v F  1

   E  k  で定義されるフェルミ速度をまとめる。計算結果より、SVWN(LDA)および BLYP(GGA)ではフェルミ速度を実験と比べて過小 評 価 す る こ と が 分 か る 。 一 方 、 HF 法 で は フ ェ ル ミ 速 度 を 過 大 評 価 す る が 、B3LYP(hybrid-DFT)法は実験値をよく再現していることが分かる。これは、B3LYP 法がその他の手法と比較して実験に対するバンドギャップの再現性が良いためであると考えられる。

表1

フェルミ速度 [m/s]

HF

SVWN

BLYP

B3LYP

実験

2.51 × 106

0.86 × 106

0.85 × 106

1.18 × 106

1.093 × 106

まとめ
本研究では周期境界条件下での量子化学計算手法について研究した。当日は、手法や計算結果の詳細に加えて、交換相関項に遮蔽を行う理論についても議論を行う予定である。

1E04
京コンピュータを用いたウイルスの全原子分子動力学シミュレーション（名大院工）岡崎進

All-Atom Molecular Dynamics Calculation Study of Viruses
Using the K-Computer(Nagoya Univ.)Susumu Okazaki【序】 コンピュータの演算性能は現在に至っても衰えることなく向上を続け、計算機技術の進歩はとどまるところを知らないかのごとくである。このことは、京コンピュータの登場をみてもよく理解できる。このような流れの中にあって、分子シミュレーション研究も研究対象を一層巨大なシステムへと展開しようとしている。
計算機性能が向上し、計算機に応じたソフト開発がなされれば、それまで扱うことのできなかった大きな系に対する研究が可能となる。大規模計算によりはじめてシミュレーション
研究が可能となる巨大物質系は、新しいサイエンスの宝庫である。京コンピュータを用いて実際に手の届きそうなものは、概ね 106 から 107 原子系で構成される系である。もしこのような系に対して MD 計算が実行できるようになれば、ウイルスカプシドの構造安定性や感染の初期過程についての分子レベルでの研究も可能となる。その他、たとえば細胞環境に近い高濃度タンパク質水溶液、リポソーム、ミセル集合体、高分子膜等々も研究対象となり得る。
【MODYLAS】

高並列計算機を使用するに際しては、高い並列化性能を持つソフトの準備が重要となる。高並列汎用分子動力学シミュレーションソフト MODYLAS1 は、京コンピュータにおける超並列計算に対応するために、並列化効率を大幅に低下させる FFT を実質的に回避しながら、周期境界条件下にある分子集合体における長距離力に対して FMM を用いて高精度に計算できるペタフロップス級の MD 計算を目指した汎用プログラムである。京コンピ
ュータ上で行った 65,536 ノード（524,288 コア）を用いた計算では、1,000 万原子系の MD 計
算１ステップの実行を 5 ms で終える性能を有している。MODYLAS は、自由エネルギー計算
をはじめとして分子シミュレーションに関わる様々な機能を備えており、京コンピュータ等を用いたナノ分野、バイオ分野における分子集団系の研究に用いることができる。本ソフト
はソースも含めて近々に公開予定である。
【ＭＤ計算】 小児マヒウイルスカプシドは、VP1 から VP4 までの 4 種類各 60 個、合わせて 240 個のカプシドタンパク質からなる 2。ここでは、電解質水溶液中のウイルスカプシド 1個に対してシミュレーションを行った。生理学的条件 pH7.4 においては、カプシドは合わせて 4,860 個の正電荷と 5,100 個の負電荷にイオン解離しており、240 個の過剰の負電荷を持つ。溶媒としての水分子は 1,884,218 個、電解質は、Na+、K+、Cl─がそれぞれ 5,310 個、136
個、5,206 個で、全系は 6,480,236 原子で構成される。
初期構造にはＸ線結晶構造解析 PDB 1HXS3 を用いたが、RNA はＸ線で観測されず構造情報がないため、本研究においては RNA を含めていない。つまり、empty capsid に対するシミュレーションとなっている。
計算開始直後の NVE、
NVT アンサンブルにおける様々な初期化過程を経て、NPT アンサンブルでの長時間シミュレーションを実施し、熱平衡状態に到達した。平衡到
達後のスナップショットを図に示す。全計算時間は 200
ns である。
【カプシドの構造安定性】 平衡状態において、単位時間当たりのカプシド内から外へ、またカプシド外から内への水分子の移動数は同数であり、水の交換平衡が成立している。水の移動速度はいずれの方向も 1 ns あたり

図 ＭＤ計算からの小児マヒウイル
スカプシドのスナップショット。

約 9 個であり、速い。つまり、約 22 マイクロ秒でカプ
シド内のすべての水分子を自発的に交換する。このことは、小児マヒウイルスは数千気圧にものぼる圧力耐性を有する一方で、乾燥に対して容易に不活化するといったことをよく説明している。一方で、電解質は透過させない。つまり、カプシドは浸透膜として機能する。生理的条件下での電解質濃度において浸透圧はさほど高くはなく、カプシドの安定性に影響するほどではないだろう。事実、真水中でウイルスは長時間活性を保つ。
カプシド内外の電解液の状態はさらに興味深い。カプシド内部と外部の電解液の圧力を密度－圧力検量線から評価すると、カプシド内の電解液圧力は負の値を示す。ビリアル定理か
ら求めた圧力も負であり、密度の結果と定量的に良い対応を示した。負圧はバルクにおいては存在し得ないが、小さな空間であるカプシド内は負圧に保たれる。負の圧力は不均一な系、
例えば植物の維管束やナノスケールの毛細管において現実に生成されることはよく知られているが、ウイルスにおいてもそれが見られることは興味深い。この圧力差によりカプシドは収縮し、これに拮抗する復元力を発生させ力学的釣り合いを保つ。小児マヒウイルスや口蹄疫ウイルスなどピコルナ科のウイルスにおいては、自然界でも多く empty capsid が生成されるが、RNA を含む完全なウイルスと比べ不安定である。負の圧力は熱力学的には不安定であることを意味しており、empty capsid の不安定性の一因となっていると考えられる。
【感染初期過程】

一方で、感染の初期過程であるレセプターとウイルスの結合において、両者の相互作用を自由エネルギーレベルで解析しつつある。ＭＤ計算からは、CD155 レセプ
ターと小児マヒウイルスとの間に働く平均力は比較的近距離においても小さいものでしかな
く、このような巨大な構造体間でもわずかに水中のメタン－メタン間の相互作用の 10 倍程度でしかない。これは、レセプターがウイルスを認識する際に、かなり確率論的な過程が含まれていることを意味する。レセプターとの結合はエンドサイトーシスや脱殻のトリガーとなり、その結合構造は大変興味深いものであるが、現在、induced fit 機構も含めた、詳細な検討を行っているところである
本研究は、理研 AICS における京コンピュータを用いて行った。
1) Y. Andoh et al. J. Chem. Theor. Comp. 9, 3201(2013).
2) Fields Virology, ed. by D. M. Knipe and P. M. Howley, Lippincott Williams & Wilkins(2007).
3) S. T. Miller, J. M. Hogle, D. J. Filman, J. Mol. Biol. 307, 499(2001).

1E0
1E06
京コンピュータを用いたウイルスの全原子シミュレーション
２．高並列汎用分子動力学シミュレーションソフト MODYLAS の開発（名大院・工 ，名大院・工 計算セ ，立命館大・薬 ，金沢大・理工 ，分子研 ，理研 ，
富士通 ）安藤 嘉倫 ，○吉井 範行 ，藤本 和士 ，水谷 圭佑 ，小嶋 秀和 ，山田 篤志 ，
岡崎 進 ，川口 一朋 ，長尾 秀実 ，岩橋 建輔 ，水谷 文保 ，南 一生 ，
市川 真一 ，小松 秀実 ，石附 茂 ，武田 康宏 ，福島 正雄
1

2

3

4

5

6

7

1

1, 2

1

3

4

7

1

4

7

1

5

7

1

5

7

6

7

LargeK--computer
Large-scale AllAll-atom Molecular Dynamics Calculation of Viruses using K
2. Development of a Highly Parallelized GeneralGeneral-Purpose Molecular Dynamics
Simulation Program, MODYLAS(Graduate School of Engineering, Nagoya Univ.1, Center for Computational Science,
Nagoya Univ.2, Department of Pharmaceutical Sciences, Ritsumeikan Univ.3, Graduate
School of Natural Science and Technology, Kanazawa Univ.4, Institute for Molecular
Science5, Advanced Institute for Computational Science, Riken6, Computational Scienceand Engineering Solution Division, Fujitsu Limited7) Yoshimichi Andoh,1 ○Noriyuki
Yoshii,1,2 Kazushi Fujimoto,3 Keisuke Mizutani,1 Hidekazu Kojima,1 Atsushi Yamada,1Susumu Okazaki,1 Kazutomo Kawaguchi,4 Hidemi Nagao,4 Kensuke Iwahashi,5Fumiyasu Mizutani,5 Kazuo Minami,6 Shin-ichi Ichikawa,7 Hidemi Komatsu,7 Shigeru
Ishizuki,7 Yasuhiro Takeda,7 and Masao Fukushima7【はじめに】
はじめに】均一な液相や固相をはじめ、低分子や高分子水溶液、脂質膜やタンパク質からなる複合体など、さまざまな分子集合体についてその物理化学的性質を明らかにするために、分子動力学（MD）
計算が広く用いられてきた。コンピュータの性能向上に伴って大規模系の MD 計算が行われるよ
うになり、近年では長距離静電相互作用を含む系についても一千万原子を超えるものが実施されつつある。しかしながら、MD ソフトウェアに目を向けると、クラスターマシン程度では高い性能を発揮するものが見受けられるものの、超高並列コンピュータの性能を十分に引き出して大規模系を効率的に実行可能なものは現状でもほとんど開発されていない。そこで我々は「京」スパ
コンのような超高並列コンピュータを用いて大規模系の MD 計算を高効率に実行することを目指して、高並列汎用分子動力学シミュレーションソフト MODYLAS[1]の開発を進めてきた。
【開発内容】を使った研究対象として、ウイルス、タンパク質複合体、ベシクルならびに高分子
材料や有機・無機材料などといった水溶液系での自己組織化構造体の物理化学を想定している。
そのためには近距離 LJ 相互作用計算だけではなく、長距離静電相互作用および各種汎用ポテン
MODYLASシャル関数に応じた分子内相互作用計算を行うことができなければならない。また研究での実用上、各種統計アンサンブル、拘束の動力学についても実装が欠かせない。
一方、対象マシンに想定した「京」の特徴として、その膨大なノード数およびノード間を繋ぐ
トーラスネットワークがある。これらハードウェア要請から「京」上にて高並列化性能を実現するためには、（１）対 MPI 並列スケーラビリティに優れた計算アルゴリズム、（２）トーラスネッ
トワークに適した隣接通信方式、を実装する必要がある。（１）について、MD 計算では特に静電相互作用が問題となる。一般的な MD ソフトウェアは高速フーリエ変換 (FFT) にもとづく
Particle Mesh Ewald(PME)法によって静電相互作用を扱っている。しかしながら、たとえ最新のFFT ルーチンを使用したとしても、一般的には 512 MPI 並列程度で FFT 並列性能が飽和するという問題がある[2]。MODYLAS ではこの点を回避するために従来の PME 法に加え、FFT を行
わない高速多重極展開法を選択できる仕様とした。（２）については、まず基本セルをトーラスメ
ッシュ位置との対応をとる形でプロセスごとのセル・ブロックに分割した。その上で、必要に応じて周辺プロセスとのみ座標、速度、力および多極子情報を隣接一対一通信する仕様とした。通信／演算両配列の共用を前提としたメタデータ構造を実装し、超高並列時に律速となりうる配列
間コピーを隠滅した。並列方式は MPI および OpenMP のハイブリッド型であり、演算量の多い部分は SIMD 化している。力場は生体高分子系で標準的な CHARMM (with CMAP) および
AMBER、OPLS に対応している。
【実行性能】

図 1 はタンパク質分子を含む 1000 万原子系をインプットとして「京」上で性能を測定した結果である。
実際的な系に対して 65,536 MPI 並列 (524,288 コア) 時においても、対 64 ノード比で高いスケーラビ
リティを達成していることがわかる。さらに 1000 万原子系について、1 MD ステップをわずか 5ms 程度で実行している。これは３日で 100 ns の MD 計算が可能であることを意味しており、長時間計算においても MODYLAS が非常に有効であるを示している。
今後さらに発展を遂げるであろう超高並列コンピ
ュータに MODYLAS を適用することにより、HIV や
インフルエンザといったウイルス、あるいは細胞小 図1. 1000万原子からなるタンパク質水溶液器官といった極めて大規模かつ重要な系についても、 系の1MDステップあたりに要する計算時
全原子モデルとして取り扱えるようになるものと期 間(赤)と、並列ノード数N に対する加速待される。
率(青、対64ノード比）．
node【参考文献】の詳細については Y. Andoh, et al., J. Chem. Theory Comput. 9, 3201−3209(2013)を参照のこと。なお、MODYLAS は近日中に公開する予定である。

[1] MODYLAS[2] C. C. Gruber, J. Pleiss, J. Compt. Chem. 32,
32 600-606 (2011).

1E07
シクロオキシゲナーゼ-2 によるプロスタグランジン G2の生合成経路に関する酵素の効果の研究(茨城大理 1、カセサート大理 2) 〇吉村
森 聖治 1

誠慶 1, Paul Gleeson2, Supa Hannongbua2,

Role of Enzyme on Biosynthesis of Prostagrandin G2 by
Cyclooxygenase-2(Ibaraki Univ., Kasetsart Univ.) 〇Takayoshi Yoshimura1, Paul Gleeson2, Supa
Hannongbua2, Seiji Mori1【緒言】シクロオキシゲナーゼ(COX)はアラキドン酸(AA)からプロスタグランジン(PG)H2 を生成するタンパク質である。COX にはシクロオキシゲナーゼ活性部位とペルオキシダーゼ活性部位が存在し、それぞれ AA から PGG2、PGG2 から PGH2 を生成する。
COX には 3 つのイソ酵素(COX-1、COX-2、COX-3)が知られているが、その中でCOX-2 はアレルギー反応、炎症作用に深く関係している。現在 COX-2 を選択的に阻害する薬剤(NSAIDs)に注目が集まっていて、その詳しい反応機構を検討することは非常に重要である。

図 1. アラキドン酸カスケード
過去の理論的研究では、シクロオキシゲナーゼ活性部位について COX や AA をフェノキシルラジカルと基質の一部を切り出したモデルを用いた研究が行われている
1,2。本研究では、COX による AA から PGG2 の生合成経路について、ONIOM 法を用いたタンパク質全体のモデルを用いて触媒反応機構を検討した。

【計算方法】X 線結晶構造解析 3 から得られた COX-2 の単量体を初期構造として 20
ns の分子動力学(MD)シミュレーションを行い、その中から反応が進行するのに適した構造を選び、ONIOM モデル(モデル A)の初期構造とした。高レベル領域の計算には B3LYP/6-31G(d)を、低レベル領域の計算には AMBER 力場を用いた。また、タン
パク質が無い状態と比較するために、タンパク質の一部と基質の一部を切り出したモデ ル ( モ デ ル B) を 用 い て B3LYP/6-31G(d) レ ベ ル で 反 応 機 構 を 検 討 し 、
B3LYP(PCM)/6-31G(d)レベルでエネルギー一点計算を行った。計算には Gaussian09
プログラムを用いた。
【結果】ONIOM モデル、DFT モデルのどちらも反応の最初の段階である Tyr385 ラ
ジカルによって基質の水素原子を引き抜くときの反応障壁が最も高いことが分かった。このことは水素原子を重水素に置換した速度同位体効果による実験結果とも一致している。4 より詳細な結果は当日報告する。

図 2. Tyr385 ラジカルによる基質の水素原子引き抜き反応の原系(A0)、遷移状態(TS-A1)、生成系(A1)それぞれの構造。
[1] Blomberg, L. M.; Blomberg, M. R. A.; Siegbahn, P. E. M., J. Phys. Chem. B
2003, 107, 3297-3308[2] Silva, P. J.; Fernandes, P. A.; Ramos, M. J. Theor. Chem. Acc. 2003, 110,345-351[3] Vecchio, A. J.; Simmons, D. M.; Malkowski, M. G. J. Biol. Chem. 2010, 285(29),21903-21907[4] Wu, G.; Lu, J.-M.; Donk, W. A.; Kulmacz, R. J.; Tsai, A.-L. J. Inorg. Biochem.
2011, 105, 382-390

1E08
鉄触媒による窒素分子変換機構に関する理論的研究（九大先導研 1・東大院工 2）○田中宏昌 1・結城雅弘 2・佐々木晃逸 2・三宅由寛 2・西林仁昭 2・吉澤一成 1

Theoretical Study on Fe-Catalyzed Conversion of Molecular Dinitrogen(IMCE, Kyushu Univ.1, School of Engineering, The Univ. of Tokyo2)○Hiromasa Tanaka1, Masahiro Yuki2, Kouitsu Sasaki2, Yoshihiro Miyake2,
Yoshiaki Nishibayashi2, Kazunari Yoshizawa1【序論】 化学的に不活性な窒素分子をアンモニアへ変換する窒素固定反応は産業的に極めて重要なプロセスある．工業的に用いられる Haber-Bosch 法が高温高圧条件を必要とするエネルギー多消費プロセスであることから，穏和な条件下での反応を可能とする触媒の開発が求められている．遷移金属錯体を触媒とした窒素固定反応はこれまで 2 例のみ知られており，いずれもモリブデン錯体が採用されている．1,2 ごく最近
西林らの実験グループは，鉄カルボニル Fe(CO)5 やフェロセンなどの安価で単純な鉄錯体を触媒とし，THF 溶媒中，還元剤 Na およびシリル化剤 Me3SiCl 共存下，常温常圧で窒素ガスからシリルアミン N(SiMe3)3 を合成する手法を開発した．アンモニ
ア等価体であるシリルアミンは水処理によりアンモニアへ定量的に変換可能である
ことから，本反応は窒素固定反応の一種であり，「鉄錯体を触媒とした窒素固定反応」の初めての例となった．一方で，反応中間体などの実験的情報は全く得られていない．
本研究では，ある鉄二価錯体を共通する反応活性種として推定し，その触媒サイクル
を理論的に検証した．3
N2 + 6 Na + 6 Me3SiCl(1 atm)

cat.:

cat.
THF, rt

H 2O

CO
OC

Fe

OC

2 N(SiMe) 3 + 6 NaCl

CO

or

2 NH3

Fe

CO

【計算手法】Gaussian 09 プログラムを用いて密度汎関数法による計算を行った．異なるスピン状態間のエネルギー差の評価が必要となったため，計算理論に B3LYP のパ
ラメータを変更した B3LYP* を，基底関数に 6-311+G** を採用した．窒素分子と反応するケイ素化学種は，類似の実験および理論研究の結果をもとに SiMe3 ラジカルとした．各反応ステップに関わる中間体および遷移状態の構造を求め，エネルギープロ
ファイルを作成した．

【計算結果】 実験結果から反応活性種が均一系の鉄錯体であることが示唆されている
こと，触媒活性を示す鉄錯体に共通する化学種などから，活性鉄化学種として低原子
価鉄(II)錯体である [Fe(SiMe3)2(thf)] (I) を採用した．計算により得られた触媒サイク
ルおよびそのエネルギープロファイルを下図に示す．I への窒素分子付加は障壁なしに進行する．最初の 3 個のケイ素ラジカルは鉄中心から遠位の窒素原子へ付加し(distal mechanism)，N–N 結合の切断を経て 1 分子目のシリルアミンを生じる．
この際，鉄(V)化学種であるニトリド錯体 (FeN) はエネルギー的に不安定であり，Fe 上の
SiMe3 基が自発的に転位して鉄(I)化学種のアミド錯体 VI となる．VI に対してさらに3 個のケイ素ラジカルが付加することで 2 分子目のシリルアミンが生成し，鉄中心から脱離して反応サイクルが完結する．ほぼ全ての反応ステップが発熱的に進行し，最大の活性化エネルギーが 19.0 kcal/mol (SiMe3 基の転位，VIII  IX) であることから，
この機構は室温で進行する反応としてエネルギー的に妥当である．

–12.1(no barrier)

N2

Si
Si
Si Si SiSiSi
Si Fe N N
Si Fe N N
Si Fe N N
Si
–36.6
–21.3
thf
thf
thf
(0.9)
(3.1)
II
III
IV

Si
Fe thf
Si
I
N(SiMe3)3
Si
Si

Si

Si
Si Fe N NSi3
thf
V

–16.6
(0.4)

N–N cleavage
Si migration

Fe NSi 3–121.9(4.1)
N(SiMe3)3

thf Fe NSi2

thf
X
–22.3
(0.6)

–1.6
(12.0)

VI
Si

Si

Si migration
Fe NSi3

thf
IX

+7.6
(19.0)Si
Si Si
Si Fe NSi2
Fe NSi2–22.3 thf
thf
(0.6)
VIII
VII

Si
–35.6
(no barrier)図．[Fe(SiMe3)2(thf)] (I) を活性種とした窒素－シリルアミン変換サイクル．Si は SiMe3 ラジカルの略．図中の数値は各反応ステップの反応熱 (kcal/mol) で，カッコ内は活性化エ
ネルギーを表す．
1) D. V. Yandulov, R. R. Schrock, Science, 301, 76 (2003).
2) K. Arashiba, Y. Miyake, Y. Nishibayashi, Nat. Chem., 3, 120 (2011).
3) M. Yuki, H. Tanaka, K. Sasaki, Y. Miyake, K. Yoshizawa, Y. Nishibayashi, Nat. Commun.,3, 1254 (2012).

1E09
スピンクロスオーバー錯体におけるポテンシャルエネルギー曲面の交差に関する理論的研究（九大先導研）○塩田淑仁, 阿部誉史, 吉澤一成

Scanning of potential energy surfaces of a spin-crossover complex(IMCE, Kyushu Univ.) ○SHIOTA, Yoshihito; ABE Takafumi; YOSHIZAWA,
Kazunari	  トリス(2−アミノメチルピリジン)鉄(II)	
 イオン錯体	 [Fe(2-pic)3]2+の磁気的性質
は、121	 K	
 を境にして変化する。この錯体は、121	 K	
 以下の温度では主配置として低スピン状態（一重項状態）をとり、それ以上の温度では高スピン状態（五重項状態）をとる。スピン転移現象は、熱の他に、光の照射や圧力の変化などによっても生じる。
特に、光の照射によって生じるスピン転移現象は、LIESST	 (Light	
 Induced	 Excited	 Spin	
 State	 Trapping)	
 と呼ばれ、Fe(II)	
 錯体や	
 Fe(III)	
 錯体に特有のものである。
LIESST	
 のメカニズムとして、次のようなモデルが提案されている。それは、光照射によって錯体を一重項基底状態から一重項励起状態へ励起させ、その後系間交差によって三重項状態、五重項状態へと緩和させるというものである。このモデルは、LIESST	 のメカニズムとして多くの研究者に支持されている。	
 	  いっぽう、熱によるスピン転移現象のメカニズムを理論的に解析している例は少ない。その理由として、従来の計算理論は、対象錯体の構造を正確に算出するが、そのエネルギーの安定性を実験的事実と逆に評価してしまう等の問題が存在することが挙げられる。近年、Reiher	
 らにより、この問題を解決するための計算理論が提案された。この理論は、B3LYP*
法

と

呼

ば

れ

、	
 

N
2.141
2
65 NH2
N .1382.1
2+ 2.1
Fe
5
85
7
2.1
NH2
2.145
N
H2
N

Hartree–Fock 交換エネル
ギーの混合率を従来の20％(B3LYP 法)から 15％

N
2.107

N
2.115
N 2.0232

1.42 NH2

Fe2+ 2.05
0
93
2.2

N
H2

2.103
N

N 2.021

NH2

8
2.2
N
H2

2.3

45 NH22+ 2.0
99
0 Fe

2.100

NH2

N

にすることで、基底状態のスピン多重度を正確に予言するように改良された

SQM

STM

TQM

6.8 kcal/mol

12.9 kcal/mol

12.8 kcal/mol

ものである。そこで本研究

図 1

MECP の安定構造と相対エネルギー. 単位は Å とでは、[Fe(2-pic)3]2+ 錯体

kcal/mol.

の熱的に生じるスピン転移現象のメカニズムを、B3LYP*法を用いて理論的に検討した。	
 	 	  次に、各ポテンシャル面の交差領域の最小値(Minimum	
 Energy	
 Crossing	 Point,	 MECP)を求めてそのエネルギーを比較した(図 1)。さらに安定構造の振動モードの解析を行い、低波数の振動モードが効率よく MECP に到達する変形を与える事があきらかとなった。さらに、いくつかの錯体についても MECP を求めたのでそれらについても紹介する。
【計算方法】計算方法は B3LYP*/6-311G**レベルで行った。ベンチマークとして、他の汎関数による高スピン(HS)-低スピン(LS)間のエネルギー差も検討した。はじめに、錯体の安定構造とエネルギーを求め、次いで、多重度を変更して再最適化を実行する
ことにより、構造変化と各多重度のネルギー変化を追跡した。最後に得られたエネル
ギー曲面から状態間の交差領域におけるエネルギーの最低点	
 (Minimum	
 Energy	 Crossing	 Point,	
 MECP)	
 を求めた。MECP は二つのポテンシャルエネルギー曲面 A と曲
面 B とが重なる領域、つまり、EA	 =	
 EB を満たしながら構造最適化することで得られる。	
 Farazdel と Dupuis が開発した方法は以下のように目的関数 f(x)、拘束関数 c(x)を定
める。	
 
	
 

cA，cB	
 は任意な定数、x	
 =	
 {x1,	 x2,	 …,	 xn}は各原子の各座標である。このような条件のもと Lagrange 未定乗数法を用いると	
 
	
 
EA	
 =	
 EB	
 のもとで

の最適化することと等価となる。特に cA	
 +	
 cB	
 =	
 1 のとき交差領域上で L	
 =	
 EA	
 =	
 EB となり L(ｙ)の最小値に対応するｙを決定することで MECPが得られる。	
 各構造最適化にはガアウシアンを、MECP 計算には GAMESS を用いた。	
 	
 【参考文献】[1] Farazdel A.; Dupuis, D. J. Comput. Chem. 1991, 12, 276.
[2] Y. Shiota; D. Sato; G. Juhász; K. Yoshizawa J. Phys. Chem. A, 2010, 114, 5862.
[3] Sato, D.; Shiota, Y.; Gergely J.; Yoshizawa, K. J. Phys. Chem. A 2010, 114, 12928.

1E10(MBH4)n および(MAlH4)n (M=Li, Na; n=1~6)クラスターの構造と電子状態に関する理論研究（千葉工大院・工 1，慶大・理工 2）○松澤	
 秀則 1，川崎	 夕美絵 1，岩田	
 末廣 2	 	
 Theoretical Study on Geometric and Electronic Structures of(MBH4) n and (MAlH4) n (M=Li, Na; n=1~6) Clusters(Chiba Institute of Technology1, Keio University2)○Hidenori Matsuzawa1, Yumie Kawasaki1,
Suehiro Iwata2【序】水素化ホウ素および水素化アルミニウム化合物は、水素吸蔵化合物として知られ、またそれらの結晶は固体電解質として高速イオン伝導（LiBH4）など特異的な性質を示すことが報告されており、固体に対する物性研究が実験、理論ともに活発に行われている。しかし、分子やイオンクラスターとしての構造や電子状態を扱った研究
例 は き わ め て 少 な い 。 そ こ で 、 本 研 究 で は [Lim(BH4)n]+m-n, [Lim(AlH4)n]+m-n,[Nam(BH4)n]+m-n, [Nam(AlH4)n]+m-n (n=1~6)クラスターを対象として、これらの構造と電子状態、クラスターへの電子付加・電子脱離による状態変化を調べたので報告する。
【方法】各クラスターの構造は、基底関数に aug-cc-pVDZ を用い、２次の Møller-Plesset
摂動(MP2)法により求め、振動解析によって安定性を評価した。プログラムは、Mac Pro
および自然科学研究機構計算科学研究センターの SGI UV1000 上で、Gaussian09 を使
用した。

-

-

2-

e

-

【結果および考察】BH4 , BH4 および AlH4 ,AlH42- はいずれも正四面体構造で、BH42- で-eは NH4 ラジカルと同様に、Rydberg 軌道に過剰の電子が収容される。中性の BH4 は、
BH2-と H2+のイオン錯体に近い構造となり、イオン間で H-B-H の３中心２電子結合を形-[Li(BH4)] (C3v)-[Li(BH4)]4 (C3v)

図１	
  過剰電子の分布	
 成し、その結合エネルギーは 15.5eV であった。一方、中性の AlH4 は AlH2…H2 と AlH3 + H の 2 種類の構造が得られ、前者の方がエネルギーは 1.01eV 低かった。[Li(BH4)]n (n=1-6)では、１〜３量体が平面構造、４量体と５量体では、平面構造と立体構造が得られた。４量体の立体構造（図１）は「３量体＋Li(BH4)」となっており、結晶中でも類似の構造が見られた。また６量体は、２組の３量体が重なった立体構造であった。[Li(AlH4)]n (n=1-6)の場合も、基本的にはホウ素系と同じ構造だが、全体的に対称性はホウ素系よりも低かった。[Na(XH4)]n(X=B, Al) (n=1-6)も、基本的に Li と類似した構造であるが、４量体の立体構造では、２量体が単位になっていた。

-

e	  それぞれの構造に１電子を付加した場合、付加の前
a. [Li2(AlH4)]

後で構造およびエネルギー+の変化はほとんどなかった。

b. [Li2(AlH4)]c. [Li2(AlH4)2]

図２	
  過剰電子の役割

これは付加した電子が構造中に取り込まれず、過剰電子として BH4-や AlH4-から離れた位置の Li+に電子が付加する（図 1）ためである。金属原子の数が多い[Lin+1(XH4)n]+ (X=B, Al) の幾何構造では、直線構造(図 2 a)が最安定となったが、1 電子付加によって、折れ曲がり構造（図 2 b）となり、２〜４量体では、図 2 c の XH4-が e-に置き換わった形になる。したがって、
e-が金属イオン同士をつなげるはたらきをしていると考えられる。
	  電子脱離では、1 電子脱離で水素原子が脱離するケースと、2 電子脱離で水素分子が脱離するケースがあった。水素原子が脱離する場合は、残った水素が H-X-H の架橋構造を作る。このケースは金属イオンと XH4-が同数のときに起こりやすい。また水素分子が脱離するタイプは金属イオンよりも XH4-の数が多いときにみられる。これはLi+ や Na+の数が少ないことで、XH4-同士が近づきやすくなるためと考えられる。そして、XH4-の数が多いほうが、少ないエネルギーで電子が脱離することがわかった。

-

e7.94eV+[Li(AlH4)]3
1.676[Li(AlH4)]3

8.20eV
2.792[Li(AlH4)]2

1.824

-

e

-

e

8.73eV
4.148

0.769

0.999
+

[Li(AlH4)]2

2+[Li(AlH4)]2

図 3	
 電子脱離による水素原子または水素分子の脱離