RISTニュースNo．
６２（２０１７）

計算化学を利用した分子構造のスクリーニング
Automatedscreeni
ngo
fmo
lecul
arstructureswi
thco
mput
atio
nalc
hemistry
一般財団法人高度情報科学技術研究機構
太田 幸宏
半経験的分子軌道法計算パッケージMOPACに基づき、化学反応経路自動探索ツールGRRMの利用を可能とするインターフェイス計算コードを開発する。開発コードにより、半経験的分子軌道法による高効率な反応経路探索が可能となる。その応用として、twostepstructur
alse
archによる安定構造の網羅的探索を提案する。すなわち、半経験的分子軌道法による安定分子構造の網羅的なスクリーニング後、分子軌道法や密度汎関数法など高精度計算による再最適化により、安定構造分布を系統的に得ることができる。開発コードとtwostepstructuralsearchの適用事例として、セルロースにおけるカチオン吸着構造探索を実施する。本インターフェイス計算コードは、半経験的分子軌道法が正当な範囲で、大規模分子系における様々な化学反応経路探索に利用できると期待される。

１．はじめにン化有機物のホウ素化反応経路探索［１０］、未系統的な化学反応経路探索は、原理的な関知の炭素構造体の予測［１１］など多岐に渡る。

心に加え、多彩な応用とつながるため、計算GRRMにおいて、ポテンシャルエネルギー化学における重要な課題の一つである［１］。

面の評価が計算実行時間を支配する。特に、その実現は計算機シミュレーションによるタンパク質や生体分子など、中・大規模分子様々な材料開発、例えば、触媒反応の設計

系において、計算コストの低減化は必須とな［２］、構成物質のみによる結晶構造の生成る。GRRMの探索アルゴリズムでは、ポテン［３］、に大きく寄与する。そのためには、ポシャルエネルギー面の勾配が要求され、それテンシャルエネルギー面上において最適化問

らは外部計算パッケージ、例えばGAUSSI
AN題を解くための高効率アルゴリズム、膨大な［１２］やGAMESS［１３］などを通し、分子軌計算結果に対する効率よいデータ処理法の整道法あるいは密度汎関数法により計算され備が要請される。こうした要請を満たす手法る。よって、現実的な全自動探索の実行にについて精力的な研究が実施されているは、利用可能な計算機資源に応じて、エネル［４，５，６，７］。GRRMは化学反応経路を自ギー勾配の計算に必要となるコストを削減す動的に探索する目的で開発された計算パッることが望まれる。

ケージである［７，８］。GRRMの利用によ

本稿では、中・大規模分子系における反応り、最安定構造の周りに存在し得る局所安定経路自動探索の効率化を目指し、半経験的分構造の全探索、安定構造間をつなぐ反応経路子軌道法計算パッケージMOPAC［１４，１５］と付随する遷移状態の探索、任意の反応物かに基づいたGRRMの利用を可能とするインら生成物に至る反応経路の非経験的探索が可ターフェイス計算コード［１６］の開発につい能となる。その適用事例は、シリコン表面にて報告する。加えて、開発されたコードの応おける酸素原子の反応経路探索［９］、ハロゲ用として、twoste
pstructur
alsearchによる
－11－

RI
STニュースNo．
６２（２０１７）安定構造の網羅的探索法を提案する［１７］。これる［１８，１９］。遷移状態探索は安定構造に比の手法では、以下の手順で安定構造を探索す

べより難しいが、様々な手法で対応できる。

る。すなわち、（１）インターフェイス計算

例えば、反応物と生成物が既知の場合、そのコードを使い、半経験的分子軌道法により、間を多点補間することで、遷移状態の探索と低コストかつ低精度で安定構造候補を網羅的

反 応 経 路 を 見 出 すNudge
dEl
ast
icBand法に見つける、（２）各候補に対し、より高コス［２０］やストリング法［２１，２２］は有効である。

トかつ高精度の計算手法を用いて再構造最適

個別の安定構造や遷移状態の探索におい化を実行する。言い換えれば、計算機シミュて、上記の計算手法は効果的である。しかレーションにより、安定な分子構造の網羅的し、次に示すような大域的な全探索への適用なスクリーニングが実行される。この提案のは困難である。すなわち、ある安定構造の周適用事例としてセルロースにおけるカチオンりの異なる安定構造を網羅的に探索し、かつ吸着サイト探索の結果を報告する。

出発点となった安定構造との間に存在する遷移状態を系統的に見つけ出す課題である。ま

２．化学反応経路自動探索ツール：GRRM

た、任意の反応物から出発して、非経験的にGRRMは化学反応経路を自動的に探索す可能性のある全ての生成物とそれに至る反応ることを目指し開発されている計算パッケー経路を予測する課題も含まれる。これらは挑ジである［１，７，８］。本節では、文献［１］に戦的であるが、計算機シミュレーションによ従い、その特徴をまとめる。

る未知の反応予測を実現するため、解決すべ効率的な反応経路探索ツールを開発する動き課題である。そのため、ポテンシャルエネ機を理解するため、反応経路全探索を概観す

ルギー面の重要な構造を詳細に探索する手法

る。化学反応を理解するアプローチとして、が望まれる。特に、ポテンシャルエネルギー

BornOppenhimer近似の下、ポテンシャル面は多変数 （N原子系であれば３N６個）エネルギー面を計算することは主要な位置をの関数であることから、多次元空間における占める。それは多変数非線形関数の最適化問高効率探索アルゴリズムが課題解決に必要と題として定式化される［１６］。ポテンシャルエ

なる（素朴な格子点サンプリング法におけるネルギー面を記述する変数は、例えば結合距計算コストの見積もりは興味深い。例えば、離や結合角など、分子に含まれる内部自由度ギ酸のような５原子分子でサンプル点を１００に対応する。N個の原子からなる （直線分子点とすると、１点あたりの計算時間を１秒とでない） 分子であれば、大域的な並進と回転して、探索終了時間は「宇宙の年齢」を超えを除き、３N-６個となる。ポテンシャルエてしまう［１］）。

ネルギー面の極小値 （山の谷底） は局所安

GRRMでは、上記の大域的探索を実現すべ定構造に対応する。異なる極小値の間に存在く考案された２種のアルゴリズムを含む計算する鞍点（山の谷底と尾根に挟まれたようなパッケージである［１，７，８］。第１に、非調峠）は遷移状態に対応する。ある安定構造か

和 下 方 歪 追 跡 法 （AnharmonicDownwardら遷移状態を経て別の安定構造（あるいは別Di
stortionFollowi
ng法 、略 し て ADDF法 ）種化合物への分離に対応する解離チャネル）

は、与えられた安定構造から出発し、その周に最小のエネルギーコストで到達する経路が辺に存在し得る安定構造 （すなわち、配座異反応経路に対応すると考えらえる。安定構造性体） 、遷移状態、そして別種化合物へ分離の 探 索 の た め、最 急 降 下 法 やBroyden-する解離チャネルを全探索するアルゴリズム

FletcherGoldfarbShanno法などが利用さである。ADDF法の基本的アイディアは、分－12－

RISTニュースNo．
６２（２０１７）子系のポテンシャルエネルギー面に典型的なAFIR法のどちらも、分子系におけるエネル構造に着目したものである。すなわち、遷移ギー勾配のデータを要求する。ここでは、半状態を経て異なる安定構造をつなぐポテン

経験的分子軌道法計算パッケージMOPACに

シャルは、各安定構造周りでは、調和近似さよりエネルギー勾配を計算することでれたそれに比べ、非調和性のため小さくなるGRRMの利用を可能にするインターフェイ（非調和下方歪み）。この非調和下方歪みを追ス計算コードの開発について説明する。

跡することで、安定構造候補、およびその中まず、半経験的分子軌道法［１８］について間ある遷移状態候補が系統的に探索される。

簡潔にまとめる。半経験的分子軌道法は、典その探索はスケールされた超球面上型的な分子系の実験データを再現するよう定（ScaledHypersphere、略してSHS） で実行められた経験的なパラメータを利用することされる。SHSは、調和近似の範囲で等エネルで、分子軌道法において計算コストがかかるギー面になるよう定められる。よって、SHS手続きを簡略化した電子状態の計算手法であ上で分子のエネルギーを計算することで、非る。具体的には、内殻の電子の寄与は考え調和下方歪みが効率よく探索され得る。第２ず、価電子に関係する基底関数のみを利用すに、人 工 力 誘 起 反 応 法 （ArtificialForceる （周期表の二列目であれば、s軌道とp軌道InducedReaction法、略してAFIR法） は、のみ考慮することになる）。さらに、ハート与えられた反応物から出発し、可能な生成物リー・フォック法で最も計算コストが高い候補と付随する反応経路候補を全探索するアフォック行列の評価に関係する２電子積分にル ゴ リ ズ ム で あ る。AFIR法 の 基 本 的 ア イ

ついて簡略化を行う。すなわち、異なる原子ディアは、人工力を付与することで生成物候間で基底関数の重なりをゼロと近似すること補を探索するというものである。一般に、化で２電子積分の多くを無視し、残りの２電子学反応には共有結合の切断や再結合が含まれ積分については実験データに基づく経験的パるため、反応のポテンシャル障壁は非常に高ラメータで評価をする。こうした近似手法のいと予測される （例えば、OH結合の結合解特徴から、遷移金属を含む化合物や共有結合離エネルギーは１０１．
７６kc
al/molであり［２３］、の切断・再結合を含む化学反応の適用には注

これを温度に換算すると約５００００Kに相当す意を要する。一方、有機低分子や生体分子のる）。人工力の付与は、こうしたポテンシャ

安定構造を調べる目的には非常に適していル障壁を下げ、効率のよい生成物探索に導る。MOPAC［１４，１５］は、このような半経く。生成物候補が一旦検出されれば、その付験的分子軌道法を利用した電子状態計算パッ随する反応経路から人工力の寄与を引くことケージである。なお、遷移金属を含むようなで、反応物から生成物へ至る反応経路候補場合について半経験的な手法により電子状態（また、遷移状態候補） も見出される。人工を評価する計算パッケージとしては、DFTB

力は原子間距離について線形なポテンシャル［２４］が知られている。DFTBとGRRMを組みとして与えられ、人工力強度に対応する単一合わせた計算事例は文献［２５］で報告されてのパラメータ、modelcollisionpar
ameterをいる。
では、開発されたインターフェイス計算調整することにより、生成物候補数が系統的

コ ー ド に つ い て 説 明 す る。コ ー ド 開 発 はに制御される。

GRRM１４ ［８］とMOPAC
２０１２ ［１４，１５］

３．インターフェイス計算コードの開発の 利 用 を 念 頭 に 実 施 さ れ た。開 発 言 語 は
Fortr
an９０が 使 用 さ れ た。GRRMに お け る

GRRMの主 要ア ルゴリ ズム、ADDF法と－13－

RI
STニュースNo．
６２（２０１７）

ADDF法あるいはAFI
R法で要求されるデーである。これは、内部的にMOPACを呼び出タはエネルギーおよびエネルギー勾配であすことに相当している。この呼び出しのオーる。開発コードが呼び出されると、GRRMのバヘッドのため、今回開発されたインター要求に従い、MOPACによりこれらのデータフェイスの利用の方が効果的な探索を可能には計算される。そして、MOPACの標準的なする［１７］。例えば、セルロース単量体におけ出力データはGRRM側で利用可能なデータ

るK+吸着サイト探索では、MOPACの直接利形式に整形され、GRRM側に引き渡される仕用におけるエネルギー勾配計算１回あたりに組みである（図１）。なお、本インターフェイ対する（平均）経過時間は、GAMESS［１３］スの使用においては、エネルギー勾配の差分を経由した場合の約１/１０になることが観測からHessianを評価することを指定するオプされた［１７］。よって、インターフェイスを通ションをGRRM１４側の入力ファイルで指定すことで、GAMESSの利用と比べ、約１０倍程する必要がある。GRRM１４の実行開発され度早く探索は終了し得る。すでにGAUSSI
AN

たコードは、Git
Hub
上にGNU GPLVersionを通した半経験的分子軌道法の計算により３．
０のライセンス形態で公開されている［１６］。

GRRMを利用する研究報告はあった［２６］。

今 回 開 発 さ れ た コ ー ド は、MOPACとし か し、こ の 分 析 か ら 明 ら か な よ う に、GRRMを直接結び付けるインターフェイス、MOPACを直接利用するほうが、より効果的という特徴がある。半経験的分子軌道法は、な探索が可能となる。

GAUSSI
ANやGAMESSを通して利用が可能

図１

GRRMとMOPACのインターフェイス計算コードの概念図－14－RISTニュースNo．
６２（２０１７）E（S=１）で評価される。ここで、E（S）
ギー差E（S）-４．Twostepstr
uctur
alsearchによるセルロースへのカチオン吸着サイト探索

は識別番号sを有する安定構造のエネルギー開発コードの適用事例として、twostepである。このエネルギー差は気相における安structur
alsearchによるセルロースへのイオ定性評価を与える。水溶液中における統計的

ン吸着サイト探索の結果を示す［１７］。セル性質を記述するため、PCMでの計算データ

ロース （C6H10O5）nの単量体 （n=１） と二量を利用し、温度Tにおける分配関数を次式で体（n=２）におけるカチオン（ Li，K ，Cs）定める。

+++の吸着サイト探索を実施した。まず、半経験的分子軌道法により、網羅的なスクリーニグを実行する。エネルギー及びエネルギー勾配
はMOPAC２０１２［１４，１５］におけるPM６で計

ここで、kBはボルツマン定数である。和に算された。カチオン吸着サイト探索の問題をおけるダッシュ記号は、幾何学的に区別されGRRMにおける反応経路探索の問題へ対応る構造についてのみ和をとることを意味すさせるべく、AFIR法の発展形であるmulti-る。構造の識別は原子間距離を計算すること

componentAFI
R（MCAFIR） 法［２７，２８］で実行された。２個の安定吸着構造についが採用された。すなわち、分子およびカチオて、すべての原子間距離が０．
１Å以下のとき、ンを反応物とみなし、それらの間に反応誘起

それらは等価な幾何的構造とみなす。分配関のための人工力が付与され、生成物として分数ZPCMで定まるカノニカル分布に従い、温度子-カチオン複合体が生成される。人工力はTにおいて、最安定吸着構造が見出される確

セルロースの酸素原子とカチオンの間に加え率は次式で与えられる。

られ、modelcollisionpar
ameterは１００kJ/molに設定された。探索は、５０回連続して同一反応経路が見つけ出されるまで繰り返された。
PM６に基づくMCAFIR法によるスクリーニ分配関数を利用することで、自由エネルギーング後、GAUSSI
AN０９［１２］を利用し、得や安定構造の統計分布に関するエントロピーられた構造を初期構造として密度汎関数法にも評価することができる［１７］。

よる再構造最適化が行われた。計算レベルは

図２に、PM６によるMCAFIR法により発B
３LYP/LanL２DZとし、ノード内８プロセス見された安定吸着構造の相対エネルギー差、並列で実行された。さらに、水の溶媒効果をおよびB
３LYP/LanL２DZレベルの密度汎関数

取 り 入 れ る た め、polarizablecontinuum法により再構造最適化された安定吸着構造のmodel（PCM）［２９，３０］の範囲で密度汎関数それを示す。横軸に構造の識別番号、縦軸に法による構造最適化も実行した。

相対エネルギー差を示す。図中の差し込み図

計算データおよびその解析方法についてまには、低エネルギー側から数えて３番目までとめる。twostepstructur
alsearchの結果、のデータを示す。図２から、幅広いエネル

安定吸着構造とそのエネルギーが、PM６（低ギー領域に渡り、カチオンの安定吸着構造がレベル計算でのスクリーニング） およびB３得られたことがわかる。再最適化データ（青LYP/LanL
２DZ（高レベル計算での精密化）

三角プロット）に着目する。セルロース二量でそれぞれ求まる。エネルギーの増加に伴い体について、カチオンの原子番号が増加する

昇順で、得られた構造に識別番号を割り当てに伴い （すなわち、イオン半径が増加するにる。吸着構造の相対的な安定性は、エネル伴い）、複数の最安定吸着構造が見出される－15－

RISTニュースNo．
６２（２０１７）ことがわかる。K+については、０．
６kc
al/molの使用は、分子構造の精密な予測を与えるも（温度換算で３００Kに相当） 内に２個の異なるのではなく、高精度計算のための良い初期構配座が存在した。同様に、Csについては、３造をバイアスなく生成することを目的とす個の異なる配座が存在した。それ以外のデーる、ということである。

+

次に、図３（a）の右パネルと図３（b）にタについては、エネルギー的に最も安定な吸着目し、水溶液中における最安定構造の出現着構造がただ一つ発見された。
図３にセルロースにおけるカチオン吸着構確率について考察する。温度はT=３００Kに固造を示す。分子構造の可視化および幾何学的定して解析をした。PCMで記述される水溶解析には、VMD［３１］を利用した。図３中の液は気相中の相対構造安定性を変更する。特各パネルの下部には、カノニカル分布から定+に、Cs吸着において劇的な変化を予測した。

められた最安定構造の検出確率p（）を示し
1 T

セルロース単量体について、最安定構造を見た。まず、図３（a）に着目する。左および出す確率は０．
６以下であった。すなわち、溶中央パネルに、twostepstructuralsearchに媒効果により吸着構造分布にブロードニングおける典型的な系列、すなわちスクリーニンが起きたことを意味する。これは、気相中でグ（PM６）から再最適化（B３LYP/LanL２DZ）単一の安定構造が存在することと対照的であの系列を示す。各パネルの下部に示された番る［図２（a）］。他方、セルロース二量体に号は図２の識別番号に対応する。右パネル

ついて、溶媒効果によって最安定構造を見出に、PCM下での最安定構造を示す。すなわ

す確率は非常に高くなることが確認された

ち、PCMにおける最安定構造は、スクリーニ０．
９）。こ れ も、気 相 中 で は０．
６kcal/mol（p1=ングでは１１番目の安定構造から来ることがわ内に複数の最安定構造が存在することと対照かる。こうして、密度汎関数法による再最適的である［図２（b）］。このように、twostep化は、半経験的分子軌道法で予測された安定structur
alsearchは、スクリーニングを通し性の順序を変更する可能性がある。ここで重て多数の安定構造候補を生成できることか要なことは、MOPACと組み合わせたGRRMら、吸着構造分布を調べることに適している。

－16－

RI
STニュースNo．
６２（２０１７）

図２ 安定吸着構造におけるエネルギー差
赤丸プロットがPM６を用いたmultico
mpo
nentArtifici
alFo
rceI
nducedReaction法による結果 （スクリーニング）、青三角プロットが密度汎関数法 （B３LYP/LanL２DZ） による結果 （再最適化） に対応する。（a）セルロース単量体（monomer）。（b）セルロース二量体（dimer）。Re
printedfrom Chem.
Phys.Lett.
,６４８,Y.Ot
a,S.RuizBar
ragan,M.Mac
hida,andM.Shiga,A screenedautomatedstr
uctur
alse
arc
h wit
hsemiempi
ricalmet
hods,１１９１２３,Co
pyright(２０１６)
,wit
h per
missionfromElsevier.

－17－RISTニュースNo．
６２（２０１７）

図３ セルロースにおけるカチオン吸着構造+吸着に対する、スクリーニグから最安定構造に至るまでの系列。
（a）セルロール二量体におけるLi（b）各 種 カ チ オ ン に 対 す る 最 安 定 構 造。Reprintedfrom Che
m.Phys.Lett.
,
６４８,Y.Ot
a,S.RuizBarragan,M.Machida,andM.Shi
ga,Ascreenedautomatedstr
uctur
alsearc
hwi
t
hse
mie
mpi
ricalmet
hods,
１１９１２３,Co
pyri
ght(２０１６),wit
hper
missionfrom Elsevier.

ト探索に適用し、PCMの範囲で溶媒効果を５．まとめ
中・大規模分子系における反応経路自動探取り入れたカチオンの吸着構造分布を予測す索の効率化を目指し、MOPACに基づいたる指針を見出した。開発されたインターフェGRRMの利用を可能とするインターフェイ

イス計算コードは、半経験的分子軌道法が正ス計算コードを開発した。開発された計算当な範囲で、大規模分子系における様々な化コードにより、GAMESSを経由した半経験学反応経路探索に利用できると期待される。

的分子軌道法の利用に比べ、より高い効率で謝辞

半経験的分子軌道計算による反応経路探索が可能となる。開発コードの応用として、two-

本稿で報告した研究成果は、Ser
giRuiz-ste
pst
ructuralse
archに よ る 安定 構 造 の 網Baragan博士（RUB）、町田昌彦 博士（JAEA）、

羅的探索法を提案した［１７］。すなわち、半経志賀基之 博士 （JAEA） との共同研究 ［１７］験的分子軌道法による安定分子構造の網羅的に基づきます。本研究を進めるにあたり、前なスクリーニング後、その結果を初期構造と田理 准教授 （北大）より受けた助言に感謝した高精度計算による再最適化により、安定します。本研究における数値計算の一部は日構造分布を系統的に得ることができる。その本 原 子 力 研 究 開 発 機 構 に お け るFuj
itsu

手法をセルロースにおけるカチオン吸着サイ

PRI
MERGYBX
９００により実施されました。

－18－

RI
STニュースNo．
６２（２０１７）

参考文献

Frisc
hetal.
,Gaus
sian,Inc.
,Walling-［１］大 野 公 一，前 田 理，Mol．Sc
i．５，

for
dCT,
２０１０.

A００４２（２０１１）．

［１３］Gener
alAtomicandMolecularElectronic Structure System, M. W.

［２］K.H.Houk and P.H.
Y.Cheong,Schmidtetal.
,J.Co
mput
.Chem.
１４,Nature455,309(2008)

１３４７(１９９３).

［３］S.M.Woodley and R Cat
low,Nat.

［１４］MOPAC ２０１２, J. J. P. St
ewart,Matter.
７,９３７(２００８).

Stewart Co
mput
ational Che
mistry,［４］D.J.Wales,EnergyLands
capes２nded. (Cambridge University Press,

Version １５.
１８０L.

Cambridge,２００３); OPTIM:

http://www.
openmopac.
net/index.
ht
ml

ht
t
p:
/
/
wwwwal
es.
ch.
c
am.
ac.
uk/

［１５］J.D.C.Mal
a,G.A.U.Car
val
ho,C.

［５］A.O.Lyakho
v,A.R.Oganov,HP.Mangue
iraJr.
,S.R.Sant
ana,L.

Stokes,and Q.Zhu,Co
mput.Phys.

A.F.Cabral,and G.B.Ro
cha,J.

Co
mmun.１８４,１１７２ (２０１３);USPEX:

Chem.Theor
yCo
mput
.
８,３０７２(２０１２)
.

htt
p://
uspe
xt
eam.
or
g

［１６］gr２moc:htt
p:
//git
hub.
com/otayuki/
gr［６］S． Q． Wu， M． Ju， C． Z．

２moc

Wang， M． C． Nguyen， X．

［１７］Y. Ota, S. Ruiz
Barragan, M.

Zhao， K， Umemoto， R． M．

Machida,andM.Shiga,Chem.Phys.

Wentzcovitch， and K． M． Ho，Lett.
６４８,１１９(２０１６).

J． Phys．： Condens． Matt．２６，［１８］F.Je
nsen,Int
roductiontoCo
mput
a-０３５４０２(２０１４)

tionalChemistry ２nd ed.(Jojn&［７］S. Maeda, K. Ohno
, and K.

Wiley,WestSussex,２００７).

Morokuma, Phys
. Che
m. Chem.

［１９］J.E.De
nnis,Jr.andR.B.Sc
hnabel,Phys.
１５,３６８３(２０１３).

Numerical Methods for［８］S.Maeda,Y.Har
abuchi,Y.Osada,T.

Unconstrained Optimiz
ation andTaketsugu,K.Mo
rokuma,and K.

Nonlinear Equat
ions (SIAM,Phila-Ohno,seedel
phia,１９９６).

htt
p://grrm.
chem.
tohoku.
ac.
jp/
GRRM

［２０］G.Henkelman and H.Jó
nsson,J.

(
accept
ed１１.
０４.
１４)
.

Chem.Phys.
１１３,９９７８(２０００).

［９］S.Ohno,K.Shudo,M.Tanaka,S.

［２１］W.E.W.Re
n and E.VandenEi
jn-Maeda,andK.Ohno,J.Phys.Chem.
C１１４,１５６７１(２０１０).

den,Phys.Re
v.B６６,５２３０１(２００２).
［２２］松永康佑，統計数理６２，２８５（２０１４）．

［１０］R.Uemat
su,E.Yamamo
to,S.Maeda,［２３］S.J.Bl
anks
byandG.B.El
li
son,Acc.
Che
m.Re
s.
３６,２５５(２００３).

H.I
to,and T.Taketsugu,J.Am.
Chem.Soc.
１３７,４０９０(２０１５)［２４］M. Elstner and G. Seifert, Phil.

［１１］K． Ohno
， H． Satoh， T．

Trans.R.Soc.A.
３７２,２０１２０４８３(２０１６).

Iwamoto， H． To
koyama， And H．

［２５］H. To
koyama, H. Yamakado, S.

Yamakado， Chem． Phys．

Maeda,and K.Ohno,Bull.Chem.

Lett．６３９，１７８（２０１５）．

Soc.Jpn.
８８,１２８４(２０１５).

［１２］Gaussi
an０９,Revi
sion C.
０１,M.J.

［２６］T. Takayanagi, K. Takahas
hi, A.

－19－RISTニュースNo．
６２（２０１７）Rev.
９４,２０２７(１９９４).

Kakizaki,M.Shiga,and M.Tachi-［３０］J.To
masi,B.Mennuc
ci,R.Cammi,kawa,Che
m.Phys.
３５８,１９６(２００９).

Chem.Rev.
１０５,２９９９(２００５).

［２７］S. Maeda and K. Mo
rokuma, J.

［３１］W. Humphrey, A. Dalke, and K.

Che
m.Phys.
１３２,２４１１０２(２０１０).
［２８］S. Maeda and K. Mo
rokuma, J.
Chem.TheoryComput.
７,２３３５(２０１１).
［２９］J.Tomas
iand M.Pe
rsico,Chem.

－20－Schulten,J.Mo
l.Gr
aphi
cs１４,３３(１９９６).