4E18
量子モンテカルロ法による 1,4-ジヨードベンゼン分子結晶多形予測	 
（北陸先端科学技術大学院大学*）	  ○本郷研太*，前園涼*	 

Quantum Monte Carlo study of predicting 1,4-diiodobenzene polymorphism(JAIST*) ○Kenta Hongo*, Ryo Maezono*【緒言】	 
分子結晶はアスピリンなどの薬理物質や、近年では、クリーンエネルギー実現のための有機半導体機能性物質として注目されている。その特徴は、系の安定性が分子間力(水素結合・分散力)で決まるため、比較的小さいエネルギー(数〜数十 kcal/mol)で容易に構相転移を起こし、その立体配列の違い(結晶多形)から、物理化学的性質が大きく変化する点にある[1]。シミュレーションによる正確な結晶構予測が可能となれば、ドラッグ・デザインやマテリアル・デザインにおける多大な組	 合せ探索に劇的な効率化がもたらされる。しかしながら、現行の枠組みにおいては、その構造を決定づけている分散力自体の理論記述が不十分なため、この方向での実用化が阻害されている。分散力は高次の電子相関効果であり、分子間力の中でも特に、電子間相互作用の精緻な取扱いが必要であり、ドラッグデザインや物質デザインにおける挑戦課題となっている。	 	 

本研究で対象とする 1,4-ジヨードベンゼン(para-diiodobenzene; p-DIB)分子結晶は、高効率太陽電池材料の一つとして注目を集め、その特性理解のために、密度汎関数法による先行研究が行われたが、標準的な密度汎関数(DFT)法による理論予測では結晶多形の相対的安定性について、実験結果を説明できないことが報告されている[2]。実験によれば、p-DIB は 326K で α 相から β 相への構造相転移を示し、従って、絶対零度においては、α 相の方が β 相よりも安定であることが知られている(図 1)。他方、CCSD(T)法などの高精度分子軌道法を直接分子結晶系に適用しようとすると、現状の超並列計算機であっても、計算資源不足から計算が走らないといった問題に直面する。そのため、計算精度と計算コストの両面でバランスのとれた代替法が必要となる。計算精度の観点から、第一原理量子モンテカルロ(QMC)法[3]は、分散力の取扱い信頼性の高い手法として知られており、最近では、生体分子系において分散力を高精度に再現している[4]。高い並列効率(98%超)故に超並列計算機との親和性が高く、近年、実用性を伸ばしつつある。申請者らは	 最近、前述の p-DIB 結晶多形の問題に本手法を適用し、その結晶多形の理論予測を行い、p-DIB 図	 1 :	 p -DI B 分子結晶の２つの 結晶多形：α相（上図）およ びβ相（下図）	 分子結晶多形の相対的安定性を正しく再現した[5,6]。	 	 
当該先行研究で実施した QMC シミュレーションでは、しかしながら、計算機資源の制限から、1×1×1 シミュレーション・セルでの QMC 計算を行い、それに対して、Kwee らによって提案されたDFT に基づく手法（KZK 補正スキーム[7]）によって有限サイズ効果を見積もっていた。結晶の周期性を再現するためには、本来、異なるシミュレーション・セルサイズで計算を行い、無限大への外挿が必要となるが、QMC 法は、DFT 法と異なり、非常に大きなセルサイズの計算をすることは実際上不可能であり、信頼性の高い外挿値を得ることは難しく、依然として、有限サイズ効果[8]に起因する誤差の見積もりが必要となる。そこで、本研究では、先行研究よりも大きなシミュレーション・セルサイズでのQMC 計算を行って、有限サイズ効果を見積もる。	 	 
【計算方法】	 	 
本研究では、p-DIB の多形である α と β の 2 つの相の相対的安定性を調べる。これら 2 相に対して、
本研究では、QMC 法の中でも、特に分散力の扱いに高い信頼性を持つ、拡散モンテカルロ(DMC)法を適用し、2 相のエネルギー差を評価する。本研究の DMC 計算は CASINO プログラムを用いる[3]。
DMC 計算で考慮するシミュレーション・セルサイズは、1×1×1 と 1×3×3 の 2 つであり、そこでの有限サイズ効果を調べるとともに、また、1×1×1 セルサイズに対しては、節固定近似の影響を調べるた
めに、3 つの試行節を考慮する。本研究で採用した試行関数（試行節）は、スレーター・ジャストロー型関数である。スレーター部分については、1×1×1 セルサイズに対して、LDA/PBE/B3LYP 汎関数による DFT 計算を行い、1×3×3 セルサイズに対して、PBE 汎関数のみを考慮した。（これ以降、DMC計算スペックを表す略記号として「DMC/SCF/1×n×n」を用いる。ここで、SCF = LDA/PBE/B3LYP、n = 1/3 である。）結晶軌道は平面波で展開し、そのカットオフエネルギー値は 40 Hartree とした。
試行関数生成のための DFT 計算は Qunatum Esrepsso を利用したが、1×1×1 セルの LDA 節生成計算
は先行研究と同様に ABINIT を利用した[4]。DMC 計算では、この結晶軌道をスプライン(blip)基底により再度、基底展開している。ジャストロー部分については、CASINO コードに実装されている多項式表式を採用し、合計で 36 個の変分パラメータを含む。これらの変分パラメータは、分散最小化法
によって、最適化した。p-DIB を構成する炭素・ヨウ素・水素の各元素につき、内殻電子部分をTrail-Needs 型擬ポテンシャルを用いて置き換えた。局所近似に基づく擬ポテンシャル評価で生じるバイアスの軽減を目的として、T-move スキームを採用した。固体周期系での電子間相互作用の取扱いとしては、Ewald と MPC (Model Periodic Coulumb)	 相互作用の 2 つを考慮した。ただし、DMC 虚時間発展における分岐過程では、Ewald 相互作用のみを考慮して評価した、局所エネルギー値を用いている。Ewald と MPC で得られる結果の相違は、有限サイズ効果の大小を判断する上での良い指標となることが知られている[8]。DMC の計算条件は次の通りである：目標ポピュレーション数は、1280(1×1×1)、20480 (1×3×3)である。平衡化ステップ数は 1500 として、統計蓄積ステップ数は、1×105(1×1×1)、および、7×103 (1×3×3)	 とした。1×1×1 計算は PC クラスター計算機の 32 コアを利用し、1×3×3計算はスーパーコンピュータ京の2048 コアを利用して、並列計算を行った。	 	 
【結果と考察】	 	 
表	 1:	 DMC による p-DIB の 2 相 α および β のエネルギー差 ΔE (定義については本文を参照のこと)。
エネルギーの単位は、kcal/mol/cell である。各項目中の記号+/-の右側の数値はエラーバーを表す。	 	 	 LDA/1x1x1

PBE/1x1x1

B3LYP/1x1x1

PBE/1x3x3

Ewald-3.8 +/- 1.6-0.5 +/- 1.5-0.3 +/- 1.4-1.8 +/- 1.3

MPC-3.2 +/- 1.6

2.8 +/- 1.5-6.5 +/- 1.4-1.5 +/- 1.3

表 1 に、DMC 計算によって得られた、2 つの結晶多形のエネルギー差	 ΔE = E(α)- E(β)	 の評価値を示す。ここで、E(α/β)は α/β の全エネルギーであり、ΔE < 0 は、α 相の方が β 相よりも安定であることを示し、すなわち実験結果との一致を意味する。表１から分かるように、1×1×1 セルサイズの場
合、DMC/LDA、DMC/PBE、DMC/B3LYP の 3 つの結果を比較すると、相互作用にかかわらず、結果が固定節に強く依存することが分かる。特に、DMC/PBE と DMC/B3LYP では、Ewald と MPC の結果に大きな差が存在するが、それに対して、DMC/LDA では、Ewald と MPC の差は小さく、両者はエラーバーの範囲内で一致しいる先行研究では、DMC/LDA/1×1×1 計算を行って、実験と一致する結果が得られていたが、本研究の結果から、先行研究は相対的安定性の再現に偶々成功したものであると結論づけられる。次に、DMC/PBE/1×3×3 の結果に移ると、Ewald と MPC の結果が非常に良く一致しており、2 相間のエネルギー差に関する限り、有限サイズ効果の影響は極めて小さいと結論づけられる。また、DMC/PBE/1×1×1 では、ΔE の符号を正しく再現できなかったが、シミュレーション・セルを 1×3×3 まで大きく取ることで、DMC/PBE/1×3×3 では、ΔE < 0 を正しく再現していることが分かった。当該結果を得るために、いくつかの準備計算を行っており、それらについての詳細は当日に発表する。	 
【参考文献】	 	 
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]

Bernstein, J. Polymorphism in Molecular Crystals; Oxford University Press: New York, 2002.
Brillante, A.; Della Valle, R. G.; Farina, L.; Venuti, E.; Cavazzoni, C.; Emerson, A. P. J.; Syassen, K. J.
Am. Chem. Soc. 2005, 127, 3038–43.
R. J. Needs, M. D. Towler, N. D. Drummond, P. Lopez Rvos, J. Phys.: Condens. Matter 2010, 22,023201.
Hongo, K.; Cuong, N. T.; Maezono, R. J. Chem. Theory Comput. 2013, 9, 1081–1086.
Hongo, K.; Watson, M. A.; Sánchez-Carrera, R. S.; Iitaka, T.; Aspuru-Guzik, A. J. Phys. Chem. Lett.
2010, 1, 1789–1794.
Watson, M. A.; Hongo, K.; Iitaka, T.; Aspuru-Guzik, A. Advances in Quantum Monte Carlo; Chapter 10,
pp 101–117.
Kwee,H.;Zhang,S.;Krakauer,H. Phys. Rev. Lett. 2008, 100, 126404.
Drummond, N. D.; Needs, R. J.; Sorouri, A.; Foulkes, W. M. C. Phys. Rev. B 2008, 78, 125106.