フラグメント分子軌道(FMO)法を用いた散逸粒子動力学シミュレーションのための有効相互作用パラメータ算出の自動化フレームワーク

Abstract

近年，用途に応じて目的の機能を有した材料開発の要求が高まっており，大規模な構造予測のためのシミュレーションが注目を浴びている．私達は，今回，散逸粒子動力学(DPD: Dissipative Particle Dynamics)シミュレーションを例に，有効相互作用パラメータ(χパラメータと呼ばれる) をフラグメント分子軌道(FMO)法から算定するためのフレームワークとワークフローシステムを開発した．FMO計算は非経験的であるため，高分子系素材に関わるほとんどの有機化合物に対して使うことが可能であり，χパラメータの算定に対して汎用性を持つと考えられる．この報告では，本システム(FCEWS: FMO-based Chi-parameter Evaluation Workflow System)の開発の目的，理論的な背景，ならびにソフトウェアについて解説する．

1 序論

近年，マテリアルサイエンスや創薬などの分野で，高度な製品開発の要求が高まっている．そのために，用途に応じて目的の機能を有した材料の開発が強く求められている．このような材料の開発には材料構造の制御が不可欠となるが，生成物の特性には数十~数百nm単位のメソスケール単位でつくられる構造が大きな影響を及ぼすことが知られており，メソ構造の予測および制御が重要な課題となっている．こういった性質や構造制御のためには，分子シミュレーションによる予測が重要な意味を持つ．目的の機能を向上させるような材料構造を予測することができれば大幅なコスト削減と効率化につながる．

メソ構造を予測するには動的平均場理論 [1,2]および散逸粒子動力学(DPD) [3,4,5,6,7,8]などの粗視化シミュレーション技術が有効に用いられる．これらのシミュレーションにはFlory-Huggins理論 [9,10]から得られる粒子間の相互作用パラメータ(χ)を用いることができる．

χの予測の要求は以前から高く，歴史的にも多くの試みがなされている．方法の代表的なものは主に二つある．一つが，溶解度パラメータSP値からの予測である．SP値はHildebrand [11]によって考案され，原子団寄与モデル [12,13]，bicerano [14]など様々な経験的な推算モデルが考案，検証されている．単一分子の分子シミュレーションからSP値を予測する方法もある [15]．二つ目は，異種分子間の相互作用からχ値を予測する方法である．代表的なものとして，Fanらによるセグメント間の接触エネルギーからの予測や [16]，凝集体モデルの凝集エネルギーの差からの予測 [17]である．我々はFanの手法 [16]に着目した．この手法は接触する粒子間の相互作用を直接分子論的に求めることができる利点がある．しかしFanらの提案 [16]では，粒子間の相互作用を原子に割り当てられた力場パラメータ(MM法: Molecular Mechanics Method)に基づいて計算していた．MM法は計算コストが安価なものの，計算結果の正確性が経験的パラメータ頼みなこと，計算対象が開発されている経験的パラメータに依存することから 汎用性，信頼性共に改善の余地があった．

このように，計算対象，信頼性が経験値を含むパラメータによって制限されることを回避するには，粒子間の相互作用計算に，分子軌道 (Molecular Orbital: MO) 計算を導入することが有効である．MO計算は基底関数や計算方法を指定することで，一定の信頼性を有することが保証されている．我々はMO計算の中でも直接的に分子間相互作用を評価できるFMO法 [18,19]を使用し， DPDシミュレーションのための粒子間相互作用パラメータ(χ)をFMO計算から算定するフレームワークシステム FCEWS (FMO-based Chi-parameter Evaluation Workflow System)を開発した．

DPDシミュレーションでは，複数の水分子が1粒子として表現されるというような荒い解像度でシミュレーションを行う． そのため，水素結合などは表現できないが，パラメータを開発する難易度を下げ，シミュレーションするための系を大きくすることが可能であり，シミュレーションを実行する時間を短縮することが可能である． また，このフレームワークでは非経験的MO法から直接DPDシミュレーションを行うパラメータを生成しているため，現実に存在している金属以外のほぼすべての分子に対して適用可能であるという利点がある．このようなソフトウェアは高い実用性が求められる領域にこそふさわしいと考えたため，フレームワークを開発し公開することにした(この手法の元となったFanのMM法による算定手法 [16]は，材料シミュレーションソフトウェアJ-OCTA [http://www.j-octa.com]に実装されている)．

この報告では，フレームワークの開発の目的，理論的な背景，実際の計算手順，FCEWSのワークフロー•フレームワークでの実装等を解説する．

2 ベースとなる理論

Flory-Hugginsの格子理論 [9,10]では，2成分混合系の自由エネルギー変化ΔGは式(1)で示される．   

$\frac{\Delta G}{RT}=\frac{{\phi }_1}{{\text{x}}_1}\text{ln}{\phi }_1+\frac{{\phi }_2}{{\text{x}}_2}\text{ln}{\phi }_2+\chi {\phi }_1{\phi }_2$(1)

この式で，φ_i, X_i (i = 1, 2) はそれぞれ成分iの体積分率，鎖長である． χパラメータは次のように定義される．   

$\chi =\frac{Z\Delta E_{AB}}{RT}$(2)

Zはモデル格子の配位数であり，接触エネルギーΔE_ABは次の式で与えられる．   

$\Delta E_{AB}={\overline{E}}_{AB}-\frac{1}{2}\left({\overline{E}}_{AA}+{\overline{E}}_{BB}\right)$(3)

ここで，${\overline{E}}_{AB}$は成分AとBとの間の平均相互作用エネルギーである．

3 計算手順

χパラメータの計算に必須である物理量は，式(2)の配位数Zと接触エネルギーΔE_ABである． 理論的にははっきりしているのであるが，それを実際の手順として再現するにはいくつかのハードルが存在する．そこで，我々がχパラメータの算出に使用した方法を説明する．

まず，式(2)であるが，実際の計算手順では，式(4)または(5)の手続きでの算定を実装している．   

$\chi =\frac{Z\left\{\left({\overline{E}}_{AB}{S}_{AB}\text{+}{\overline{E}}_{BA}{S}_{BA}\right)-\left({\overline{E}}_{AA}{S}_{AA}+{\overline{E}}_{BB}{S}_{BB}\right)\right\}}{RT}$(4)

  

$\chi =\frac{\left({\overline{E}}_{AB}{Z}_{AB}{S}_{AB}\text{+}{\overline{E}}_{BA}{Z}_{BA}{S}_{BA}\right)-\left({\overline{E}}_{AA}{Z}_{AA}{S}_{AA}+{\overline{E}}_{BB}{Z}_{BB}{S}_{BB}\right)}{RT}$(5)

式 (4) において，Zは以下の式で算出される．   

$Z=\frac{Z_{AB}+Z_{BA}+Z_{AA}+Z_{BB}}{4}$(6)

ここで，Z_ABは，分子Aの周囲に分子Bを配置できる数，${\overline{E}}_{AB}$は，分子Aと分子Bの間の平均相互作用エネルギーである．また，S_ABは分子A,B間の相互作用の異方性の影響を示す値である(後述)．

はじめに，配位数Z_ABに関して述べる．同様の手順でZ_AA, Z_BA, Z_BBも計算可能である．この物理量は，ある分子Aの周囲に分子Bがいくつ配置されているかを示す量である．そのため，我々は，単純に原子に関して半径を設定し，いくつの分子B同士が重なることなく分子Aに接することができるかを求める計算を行うこととした．その結果配置することのできた分子Bの個数の平均を配位数Z_ABとしている．

Flory-Huggins理論では各成分の体積が同一，かつ混合による体積変化がおきないという仮定がある．この定義に従うと，式(4) のように4種のZを平均して計算する必要がある．一方，式(5) では個別のZを使用している．この算出方法では混合前後でモノマーの接触点の数が変わるために，Flory-Huggins理論の示す定義とは正確には一致しない．しかし実際にサイズが違うときものを混合すると成分間の接触点数は変わると考えられるために，式(4, 5) 両方の算出を実装した．共通する注意点として，精確なパラメータ算定には可能な限り計算対象分子の体積を等しくする必要がある．

次に，分子間の平均相互作用エネルギー(${\overline{E}}_{AB},{\overline{E}}_{BA},{\overline{E}}_{AA},{\overline{E}}_{BB}$)に関して説明する．平均相互作用エネルギー${\overline{E}}_{AB}$は分子Aと分子Bのランダムな接触構造を無数に生成し，得られる分子間の相互作用エネルギーE_ABを平均したものである．我々は，E_ABの計算方法として二種類提案する．一つ目は，フラグメント間相互作用エネルギー(IFIE [20,21])を利用する方法で，二つ目はBinding Energyを利用するものである．

まず，IFIEを利用する方法を説明する．IFIEとはFMO法から得られるフラグメント間の相互作用エネルギーである．FMO法はKitauraらにより1999年に提案された手法であり [18,19]，一般的なFMO計算は以下のように要約される．まず任意の対象分子系を，モノマーと呼ばれる適切なフラグメントに分割する． Hartree-Fock (HF)の手順を使用して，モノマー段階で環境静電ポテンシャル(ESP)下の各フラグメントのMOセットを決定する．フラグメントモノマーの大反復は，self-consistent-charge (SCC)が満足されるまで繰り返される．フラグメント内処理は並列化され，モノマーリストの手続きも並列処理(または2層並列処理)が可能である．モノマーのSCCの反復が完了すると，モノマーSCCのESPの下で，ダイマー(隣接したモノマーからからなる)の一連のHF計算がモノマーの場合と同様の並列性を持って実施される．分極および電荷移動は，モノマー段階およびダイマー段階でそれぞれ考慮される．計算結果の全電子エネルギーEは，以下の簡単な式によって得られる．   

$E=\sum _{I>J}{E}_{IJ}-\left(N_{\text{f}}-2\right)\sum _I{E}_I$(7)

ここで，E_IはモノマーIのエネルギー(N_f はモノマーの数)であり，E_IJはダイマーIJのエネルギーである． 2次のMoller-Plesset摂動論(MP2)などの電子相関によるエネルギー補正を加算的に取り込むことができる．また，式(7)は，フラグメント相互作用エネルギー(IFIE)の合計で書き直すことができる．   

$E=\sum _{I>J}\Delta {\tilde{E}}_{IJ}+\sum _I{E}_I^{\text{'}}$(8)

この式で，E'_IはESPの寄与を除いた架空のモノマーエネルギーである． IFIE値は，計算対象系におけるフラグメント相互作用の性質を分析するのに有用な量である．

次に，Binding Energyを利用する方法を説明する．Binding Energyは以下の式で計算できる．   

$E_{AB}=W_{AB}-\left(W_A+W_B\right)$(9)

W_ABは，分子Aと分子Bの会合体のFMO計算結果の全エネルギーであり，W_Aは分子AのFMO (MO)計算結果のエネルギーである．同様の手続きで，相互作用エネルギーE_AA，E_BBも計算を行っている．

これら二つの方法の違いは幾つかある．まず，IFIEでは一度のFMO計算を行う過程で計算できるが，Binding Energyでは会合体のFMO計算の他にモノマーの計算が必要になる．また，物理的意味も若干異なり，IFIEでは会合体の中での相互作用が反映されたエネルギーが計算されるが，Binding Energyではモノマーから会合体へのエネルギーの差分が計算される．フレームワークの実装では，どちらでも利用できるようにしているが，以下のような使い分けが考えられる．結合した分子内をフラグメントに分けて計算することで計算コストを削減することが可能だが，そういった場合には全エネルギーに誤差が生じるため，IFIEからの直接的な算出の方が信頼性は高い．一方，IFIEは静電相互作用を過大評価するために，分極の大きな系ではBinding Energyを用いるほうがよいと考えられ，両者の使い分けも必要となる．

Fanらのアプローチ [16]では，この相互作用エネルギーを計算する部分は，MM法での計算を用いている．一方，我々はFMO計算を使用しているため，MM法では表現できない電荷移動などの相互作用を考慮することが可能である．一方で，本稿で提案している手法においても，二体間の接触構造から凝集体の相互作用を見積もっているために長距離相互作用や，多体効果を表現できていないといった問題がある．これらの解決が今後の課題である．

4 ソフトウェアの動作

上記の計算手順に必要なパラメータは，分子の座標と原子半径，分子の電荷，MO計算の基底関数，計算レベルのみである． そのため，ソフトウェアの動作は大きく分けて，4つの部分からなる． エネルギー計算を行うための入力ファイルを自動で生成する部分, 配位数を計算する部分, エネルギー計算の結果を回収する部分, χの計算を行う部分である．

4.1 入力ファイル自動生成部

まず，エネルギー計算を行うための入力ファイル自動生成部から説明する． 必要な情報は，分子の座標と原子同士の接触半径，電荷である．座標については，予め任意の方法で構造最適化を行ったものを入力することを想定している．この部分では，2つの分子の座標を読み取り，分子Aの周囲に分子Bを配置する．配置する方向はランダムである．その後，分子Aと分子Bの中心の距離を近づける． 分子Aと分子Bの原子がちょうど接するようになるまで両方の分子の中心を近づけていく(Figure 1)．接触半径についてはデフォルトでvan der Waals 半径 [22]を基にした情報が用意されているが，ユーザーが任意に変更可能である．このような方法で，分子Aと分子Bの二分子モデルの構造を作成する．

Figure 1.

 Example of generated configuration (propane and butane pair). Green and white spheres represent carbon and hydrogen atoms, respectively.

式(4, 5)に示すように，xパラメータの算出には，分子Aを中心とした分子Bの会合体，分子Bを中心とした分子Aの会合体，分子Aの会合体，分子Bの会合体といった四種の平均相互作用エネルギー${\overline{E}}_{ij}\left(i,j=A,B\right)$が必要である．また，分子Aの単量体のエネルギーW_Aや分子Bの単量体のエネルギーW_Bが必要であり，入力ファイル自動生成部でファイルを生成している．また，相互作用エネルギーE_BAの入力ファイルは，相互作用エネルギーE_ABとほぼ同じであり，座標の中心が分子Aにあるか分子Bにあるかの違いしかない．そのため，相互作用エネルギーE_BAのファイルは生成しない．

この時，四種の平均相互作用エネルギー${\overline{E}}_{ij}(i,j=A,B$)の計算に関しては，各分子の位置や配向の仕方などで計算結果が変化するため，我々は各相互作用エネルギーの計算のために，2000構造を生成し座標ファイルとして生成している．また，この生成する構造の数はユーザーが設定することが可能である．

本ワークフロー [23]では，ABINIT-MPプログラム [24]でFMO計算を行う．座標ファイル生成と合わせて，エネルギー計算で使用しているABINIT-MPプログラムのための入力ファイルを出力している．この入力ファイルには分子全体の電荷情報，FMO計算での計算レベル，基底関数が記されている．

4.2 配位数Z計算部

ここでは，式(4, 5)で必要になる分子Aを中心にした場合の分子Bの配位数Z_ABの計算を説明する．この計算に必要な情報は，2つの分子の座標と，原子同士の接触半径である． 手順は以下の通りである．

① ある分子Aを空間に配置し，ランダムな方向と適当な距離で分子Bを配置する．

② その後，分子Aと分子Bの各原子がちょうど，一点だけで接するような距離まで，各分子の中心を近づけて行く．

③ 予め分子Bが存在している場合は，分子Bと重なっているかを判定する．

④ 重なっていなければ分子Bを配置する．

⑤ もう一度手順1, に戻り分子Bを配置する．

⑥ 分子Bが配置できなくなるまでstep1からstep5を繰り返す．

⑦ 分子Bの配置されている数を配位数Yとして出力する．

このようにして得られた分子Bの個数を配位数Yとして使用する(Figure 2)． 実際には，配置できなくなるまでを判定することは難しいため，step1からstep5までの試行を1万回行っている． 更に，この試行を10回繰り返した後のYの平均値を配位数Z_ABとして計算している．これは，計算の初期条件などに依存してYが変化するためである． また，分子Aがポリマーである場合，排除体積を考慮すると高分子鎖の方向の配位が妨げられるために，モノマーで計算した配位数よりも高分子における配位数は減少すると考えられる．そこで，先行研究 [16]に倣い，分子Aがポリマーの場合は計算したYから2を引いて配位数Z_ABとして扱っている．

Figure 2.

 Example of coordination number Z calculation. Butane molecules are arranged around a propane molecule.

4.3 計算結果回収

入力ファイル自動生成部を行った後，ABINIT-MPでのFMO計算を行った上でこの処理を行う．セクション3で計算した2000対の計算結果の各々の相互作用エネルギーを取得する．この時に発生する不具合と我々が行った対応策を説明する．

MO法の計算を行っていると，SCFが収束しないために計算結果が存在しない場合や外れ値が発生する．そのため，計算結果が存在しない場合の処理，外れ値の検出処理を行っている． 外れ値として判定された計算結果は，計算結果が存在しないものとして処理している．

具体的な外れ値の検出処理は，相互作用計算を行ったエネルギーの値の平均値や分散を計算し，平均値 ± 4.0σから外れたものを除外した上で，各構造のエネルギーが連続値を持っているかを判断し，安定構造のエネルギー値の85%以内の安定化を持つ構造が存在しない場合，その構造を除外している．

4.4 χ の計算

χの計算を行うためには，式(4, 5)を計算しなければならない． ここで，二つの問題がある． 一つ目は，温度をどのように${\overline{E}}_{ij}\left(i,j=A,B\right)$に反映させるのかということである．${\overline{E}}_{ij}$は，温度によって変化するべき値であるが，FMO計算でのエネルギーは分子が固定されている場合のエネルギーである． そのため，温度の影響を反映するためには，FMO計算ではない手段で温度の影響をχパラメータに反映しなければならない． 二つ目は，平均場近似と現実との違いをどのように表現するのかということである．平均場近似での${\overline{E}}_{ij}$は，分子Aの周囲に等方的に分子Bが配置している．現実ではそうではないため，その違いを補正する必要がある．

4.5 温度影響の χ パラメータへの反映

一つ目の問題に対しては，温度に関するボルツマン重みを用いて計算した． 実際には，Metropolis Monte Carlo法に従い，2000配置の中から任意の温度Tにおける取りうる配置を決定した [25]． u_vv'を，系が状態vから状態v'に遷移する確率として定義すると，u_vv'は次のように与えられる．   

${\text{u}}_{v{v}^{\prime }}=\text{exp}\left(-\frac{\Delta {\text{E}}_{v{v}^{\prime }}}{RT}\right)$(10)

ここで，ΔE_vv'は，状態vとv'との間のエネルギー差である．モンテカルロの選択により受け入れられた構成のエネルギーの平均値を，平均相互作用エネルギー${\overline{E}}_{ij}(i,j=A,B$)として定義した．

4.6 平均場近似と溶媒分子の密度の偏りの整合性

本手法では，Metropolis法による採用した配置のエネルギーの平均値をセグメント対当たりの平均相互作用としている． 次いで，対応する相互作用エネルギーにそれぞれの配位数を乗じることによって，周囲の場との相互作用エネルギーの合計が計算される．配位数は空間充填モデルからのみ計算される．セグメント対の特定の部分が強い安定化を示す場合，メトロポリスアルゴリズムにおいてサンプリングの潜在的な偏りが生じる．そのため，周囲の場との相互作用エネルギーの合計を過大評価する傾向がある．そのような過大評価を検証するために，中心セグメントの重心から周囲の空間を6つの領域に分割し，モンテカルロ計算の結果から，各領域i ( = 1–6)の重みを以下の式で計算した．   

$d_i=\frac{A{c}_i}{A{c}_{max}}$(11)

この式で，AC_i は領域iでの採用数，AC_max は採用数が最大の領域での採用数である．採用数の差は以下のように無次元化された自由エネルギーの差に相当する．   

$\text{log}\left(\frac{A{c}_i}{A{c}_{max}}\right)\cong \Delta G$(12)

しかし，Flory-Huggins理論では，セグメント分子は等方的であり，この方程式のΔGはゼロでなければならない． そこで，系の異方性を表すパラメータとして，以下のSを採用した．   

$S=\sum _i\frac{d_i}{Na}$(13)

ここで，Naは領域の数である．もし，系が完全に等方的であれば，d_i はすべての領域において同じ数になり，Sは1である．逆に，系が異方的であれば，S値は1よりも小さい値になる．S値も空間の分割に依存するため，分割の仕方を変えて100種類S値を計算し，平均したものをスケーリングファクターとして扱った．

S_ABを算出するためには，分子Aからみた分子Bの配向情報が必要であるが，入力ファイル自動生成部のセクションで示したように，分子Aを中心に分子Bを接触させる計算しか行っていない．そこで，分子Aを中心とした座標の情報を分子Bの座標が中心になるように変換することでS_BAを算出している．

4.7 溶媒効果の導入

MOによる計算は通常真空中のため，パラメータを求めたい分子種がイオン性の場合，相互作用を大きく過大評価してしまう．そこで，Poisson-Boltzmannによる溶媒効果を取り入れた計算(FMO-PB)をオプションとして導入している． 溶媒効果計算はOkiyamaらによってABINIT-MPに導入されている [26]．これによりイオン性の分子においてもパラメータ算定が可能である.

5 まとめ

今回，我々はxパラメータ算定のための自動化フレームワークを開発した．このフレームワークでは，非経験的分子軌道法の計算結果からxパラメータを算定する．そのため，分子軌道法で計算可能ならばどのような分子にも適用することが可能である．公開については，計算工学ナビ (http://www.cenav.org/) にてアナウンスの後，望月 ([e-mail]) にコンタクトをとることで配布を行う予定である．

我々は，実際にこのフレームワークを使用し，高分子ポリマーや脂質膜分子，電解質膜のためのxパラメータを計算している． [23,27,28] また，このような使用した事例に関しても実験値等との比較を行い，様々な物性値において良好な一致を示していることを確認している．また，多數のFMO計算をcapacity computingとして捉え，ジョブの自動投入を行う外部システム(HPC/PF [29])との連携も試みており，「京」に代表されるスーパーコンピュータ環境への親和性も意識した作業も進めている．

我々のフレームワークが皆の研究活動を加速することを願っている．

Acknowledgment

本研究開発は，文部科学省ポスト「京」重点課題6「革新的クリーンエネルギーの研究開発」から支援を受けている．また，日頃からご議論をいただいている(株) JSOLの小沢拓氏に深謝する．

参考文献

• [1]   J. G. E. M. Fraaije, B. A. C. van Vlimmeren, N. M. Maurits, M. Postma, O. A. Evers, C. Hoffmann, P. Altevogt, G. Goldbeck-Wood, J. Chem. Phys., 106, 4260 (1997).
• [2]   R. Hasegawa, M. Doi, Macromolecules, 30, 5490 (1997).
• [3]   P. J. Hoogerbrugge, J. M. V. A. Koelman, Europhysics Lett. EPL, 19, 155 (1992).
• [4]   J.M.V.A. Koelman and P.J. Hoogerbrugge, EPL (Europhysics Lett. 21, 363 (1993).
• [5]   A. D. Mackie, J. Bonet Avalos, V. Navas, Phys. Chem. Chem. Phys., 1, 2039 (1999).
• [6]   R. D. Groot, P. B. Warren, J. Chem. Phys., 107, 4423 (1997).
• [7]   N. Yamazaki, M. Motoyama, M. Nonomura, T. Ohta, J. Chem. Phys., 120, 3949 (2004).
• [8]   C. M. Wijmans, B. Smit, R. D. Groot, J. Chem. Phys., 114, 7644 (2001).
• [9]   D. R. Paul, Polymer Blends (Academic Press: New York, 1978).
• [10]   P. J. Flory, Principles ofPolymer Chemistry (Conell University Press, NY, 1953).
• [11]   H. A. Benesi, J. H. Hildebrand, J. Am. Chem. Soc., 71, 2703 (1949).
• [12]   D. W. Van Krevelen, Properties of Polymers, Elsevier, Amsterdam, (1990).
• [13]   M. M. Coleman, C. J. Serman, D. E. Bhagwagar, P. C. Painter, Polymer (Guildf.), 31, 1187 (1990).
• [14]   J. Bicerano, Prediction of Polymer Properties (CRC Press, 2002).
• [15]   B. Prathab, V. Subramanian, T. M. Aminabhavi, Polymer (Guildf.), 48, 409 (2007).
• [16]   C. F. Fan, B. D. Olafson, M. Blanco, S. L. Hsu, Macromolecules, 25, 3667 (1992).
• [17]   S. S. Jawalkar, S. G. Adoor, M. Sairam, M. N. Nadagouda, T. M. Aminabhavi, J. Phys. Chem. B, 109, 15611 (2005).
• [18]   D. G. Fedorov, K. Kitaura, The FRAGMENT MOLECULAR ORBIAL METHOD: Practical Applications to Large Molecular Systems (CRC Press, Boca Raton, 2009).
• [19]   K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayasi, Chem. Phys. Lett., 313, 701 (1999).
• [20]   S. Amari, M. Aizawa, J. Zhang, K. Fukuzawa, Y. Mochizuki, Y. Iwasawa, K. Nakata, H. Chuman, T. Nakano, J. Chem. Inf. Model., 46, 221 (2006).
• [21]   Y. Mochizuki, K. Fukuzawa, A. Kato, S. Tanaka, K. Kitaura, T. Nakano, Chem. Phys. Lett., 410, 247 (2005).
• [22]   A. Bondi, J. Phys. Chem., 68, 441 (1964).
• [23]   K. Okuwaki, Y. Mochizuki, H. Doi, T. Ozawa, J. Phys. Chem. B, 122, 338 (2018).
• [24]   S. Tanaka, Y. Mochizuki, Y. Komeiji, Y. Okiyama, K. Fukuzawa, Phys. Chem. Chem. Phys., 16, 10310 (2014).
• [25]   N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, J. Chem. Phys., 21, 1087 (1953).
• [26]   H. Watanabe, Y. Okiyama, T. Nakano, S. Tanaka, Chem. Phys. Lett., 500, 116 (2010).
• [27]   H. Doi, K. Okuwaki, Y. Mochizuki, T. Ozawa, J. Comput. Chem. Jpn., 16, 28 (2017).
• [28]   H. Doi, K. Okuwaki, Y. Mochizuki, T. Ozawa, K. Yasuoka, Chem. Phys. Lett., 684, 427 (2017).
• [29]   <http://www.cenav.org/20141104_hpcpf_intro/>.