Journal of Computer Chemistry, Japan
Online ISSN : 1347-3824
Print ISSN : 1347-1767
ISSN-L : 1347-1767
速報
調和溶媒和モデル(HSM)を用いた 凝縮系の自由エネルギー計算
中井 浩巳
著者情報
ジャーナル フリー HTML

2017 年 16 巻 4 号 p. 83-88

詳細
Abstract

We have proposed a novel quantum chemical model to evaluate condensed-phase thermodynamic properties, that is, the harmonic solvation model (HSM). This paper explains the theoretical background of HSM and the practical procedure to use for the HSM with some sample inputs. Illustrative applications demonstrate the usefulness and accuracy of HSM.

1 はじめに

今日,量子化学計算は,原子・分子,その集合体の様々な性質を評価する方法として広く用いられている.量子化学計算で扱う1個(あるいは数個)の分子のエネルギーから,温度や圧力といった集団としての振る舞いは直接的には評価できない.そこで,統計力学による分子分配関数を用いて,電子状態のみならず振動・回転・並進運動のエンタルピーおよびエントロピーを見積り,温度や圧力の効果を含んだ自由エネルギーが計算される.GAMESSなどの標準的な量子化学計算プログラムでは,理想気体の仮定に基づく自由エネルギー計算が実装されている.通常,振動数計算を行うと,自動的に熱力学計算も行われるように設定されている.高温・高圧などの状態でなければ,理想気体の取扱いは十分な精度を与える.

液相中の分子には,連続誘電体モデル(PCM)などの溶液理論により電子エネルギーに対する溶媒効果を考慮できる.さらに,PCMを用いて振動数計算を行うと,熱力学量も計算される.我々は,ここに重大な問題点があることを見出した.すなわち,液相であるにもかかわらず,理想気体による自由エネルギーの表式が用いられているのである.この影響は,エントロピーの過大評価をもたらす.そこで,我々は溶質分子とPCMのキャビティとの相互作用を調和振動子で近似する調和溶媒和モデル(HSM)を2014年に提案した [1].さらに,気液平衡 [1],有機化合物の生成熱および燃焼熱 [2],標準水素電極 [3],気体分子の溶解平衡 [4],有機電解質の酸化電位 [5]について応用し,いずれもHSMの有用性を確認してきた.

HSMは国内外の研究グループの注目を集めているが,実際の応用研究には至っていないようである.HSMを実行できるプログラムが公開されていないためかもしれない.しかし,近似的なレベルならば,既存のプログラムを用いた計算と簡単なデータ処理でHSMの計算が可能である.本稿では,HSMの理論的背景を簡単に説明し,次にプログラム開発の専門家でなくてもHSMを用いた熱力学計算を行える具体的な計算手順を解説する.最後にいくつかの応用計算の結果を示し,その重要性と必要性を述べる.

2 理論と実装

2.1 気相分子の自由エネルギー計算

気相分子の自由エネルギーGは,分子間の相互作用がないという理想気体のモデル(IGM)のもと,次式で表される.

   G = HTS  (1)   

   H = Eelec + Evib + Erot + Etrans + W  (2)   

   S = Selec + Svib + Srot + Strans  (3)   

ここで,Hはエンタルピー,Sはエントロピー,Tは絶対温度,Wは仕事である.下付のelec, vib, rot, transはそれぞれ電子,振動,回転,並進に対応する.Table 1に,各状態に対するモデル,量子的なエネルギー準位,分配関数,そして,エネルギーとエントロピーを示す.

Table 1.  Quantum energy levels, partition functions, and macroscopic energies/entropies for electronic state, vibration, rotation, and translation of molecule.

Eelecは,量子化学計算で求められるエネルギーであり,最も寄与が大きい.つまり,この項の精度が十分でないと,他の熱力学量の補正は全く意味がないといっても過言でない.通常の分子では電子基底状態と励起状態の間には数eVのエネルギー差がある.そのため,熱的な揺らぎにより励起状態に遷移する確率は非常に小さく,エントロピー項はほぼ零である.開殻系分子は基底状態が縮重しているので,状態数はスピン多自由度となる.3重項が基底状態であるO2分子では,以下のように計算される.

   Selec(O2) = R ln 3 = 9.133 J/(K·mol)   (4)   

振動状態は,通常,調和振動子を仮定して取り扱われる.基準振動数νiおよび振動量子数vに対するエネルギー準位εvから分配関数qvib,さらに,qvibからエネルギーEvibとエントロピーSvibが導かれる.Table 1に示す通り,EvibSvibには基準振動数が含まれており,いわゆる振動数計算が必要となる.電子エネルギーの核座標に対する2階微分,すなわち,ヘシアン行列を求め,原子核の質量で重みをつけて力の定数行列を計算する.これを対角化することにより,固有値として基準振動数,固有ベクトルとして基準座標が求められる.厳密には,原子数をNとすると力の定数行列の次元は3Nとなり,振動の自由度3N-6より大きい.そのため,固有値の小さな6個の状態は,並進および回転に対応すると見なし,振動項の計算には用いられない.

回転状態は,通常,剛体回転子を仮定して取り扱われる.Table 1には,回転定数B,回転量子数Jに対するエネルギー準位εJ,さらに残りの回転自由度も加えた分配関数qrotが示されている.ここで,回転対称数σは等価な原子の数に応じて値が変化する.エネルギーErotとエントロピーSrotはいずれも分配関数qrotから求められたものであるが,振動準位間のエネルギー差が小さいため,Erotはエネルギー等分配則に従う.つまり,Erotの値は分子種には依存せず,温度によってのみ変化する.一方,Srotには分子の回転定数A, B, Cが含まれており,分子に依存した値を取る.回転定数の見積りには,分子の構造情報が必要であり,いわゆる構造最適化計算の結果が用いられる.

並進状態は,通常,箱の中の粒子の取り扱いを仮定する.1次元並進運動の量子数nのエネルギー準位εnには,箱の長さLが含まれており,分配関数では3次元なので体積Vとなる.最終的には,理想気体の状態方程式により圧力の効果としてエントロピーStransに寄与する.エネルギーEtransは,回転状態と同じくエネルギー等分配則に従い,分子種には依存せず温度のみの関数となる.

仕事に関しても,理想気体の状態方程式より,分子種には依存しない温度のみの関数である.

2.2 液相分子の自由エネルギー計算

液相では分子間の相互作用が無視できない.それは,単一の分子からなる液体でも,溶媒に溶質が解けた溶液でも同様である.量子化学計算では,溶質-溶媒間相互作用はPCMなどを用いて考慮される.PCMでは,溶媒分子単体(あるいは第一溶媒和圏などのあらわな溶媒分子を含めた系)をキャビティ内に配置し,キャビティからの静電場(さらに分散力項や反撥項を含めた場)を考慮して,系の電子状態を自己無撞着場(SCF)的に求める.このようにして得られるのは,液相の電子エネルギーEelecである.振動状態はTable 1の気相分子と同じ表式を用いるが,PCMによる電子状態の変化により基準振動数が変化し,結果としてエネルギーEvibとエントロピーSvibが多少変化する.回転・並進状態を,それぞれ気相分子と同じく剛体回転子および箱の中の粒子として取り扱うと,溶質-溶媒間相互作用による運動の制限が考慮できない.

HSMでは,Figure 1のように溶質-キャビティ間相互作用を調和振動子として取り込み,回転・並進状態に対しても基準振動数を求める.厳密には,振動運動と回転・並進運動が分離できない状態なので,擬回転・擬並進状態というべきかもしれない.これらのエネルギーおよびエントロピーはそれぞれ以下のように計算される.   

E rot/trans = 1 2 R T i rot/trans h ν i k B T { 1 2 + exp ( h ν i / k B T ) 1 exp ( h ν i / k B T ) } (5)
  
S rot/trans = R i = 1 rot/trans [ ( h ν i / k B T ) exp ( h ν i / k B T ) 1 exp ( h ν i / k B T ) ln { 1 exp ( h ν i / k B T ) } ] (6)

Figure 1.

 Solvent-cavity interaction treated as harmonic oscillator in HSM.

これらの表式は,Table 1の振動状態に対応していることが理解できるであろう.ただし,(5)式において(1/2)という因子があることに注意されたい.これは,並進・回転状態に対応する基準振動数がゼロの極限で,理想気体の結果,すなわち,エネルギー等分配則の(3/2)RTに一致させるためである.言い換えれば,エネルギーは溶質と溶媒で分けるために,半分になると考えてもよい.一方,エントロピーは状態の数え上げなので,上記の因子は不要である.

HSMの振動数計算では,キャビティは無限の質量を有する壁として取り扱う.現在使用しているプログラムでは,解析的微分によるエネルギー勾配を用いた数値的微分によりヘシアン行列を用いる.すなわち,構造最適化計算までは,キャビティの位置は,通常通り原子核に完全に追随するが,振動数計算ではキャビティの位置を固定して,原子核の位置を微小変化させてエネルギー勾配を求める.そのため,前方・後方差分を用いた3点微分の場合,6N回のエネルギー勾配計算が必要となり,原子数Nの増加とともに計算コストが増大する.実際,数十原子系に対するHSM計算は困難である.この問題の解決には,HSMに対応した解析的2階微分の実装が必須である.現在,理論式の導出等を進めている段階である.

最近,この問題を近似的に回避する方法を提案した.HSMの最大の特長は,回転・並進状態に対する取り扱いである.そこで,溶媒分子を剛体(rigid body)と仮定し,回転・並進運動に対する6自由度のみについて振動解析を行う(RB-HSM).つまり,どのような溶質分子に対しても,3点微分で12 ( = 6*2)回のエネルギー勾配計算だけを行えばよい.

3 計算手順

本節では,HSMの具体的な計算手順を解説する.Table 2に,IGMを用いた気体分子とHSMおよびRB-HSMを用いた凝縮系に対する熱力学量の計算の流れを示す.Step 1では構造最適化,Step 2では振動解析,そして,Step 3で熱力学量を求める.RB-HSMでは,さらにStep 4で回転・並進状態に対する振動解析とStep 5でそれらの熱力学量の計算を行う.

Table 2.  Procedures to evaluate thermodynamic properties in gas-phase by IGM and those in condensed phase by HSM and RB-HSM.

気体分子に対する熱力学量の計算は,すべてGAMESSなどの標準的な量子化学計算プログラムを用いて実行することができる.

HSMの計算手順では,まずPCMなどの溶媒効果を考慮して,構造最適化を実行する.この際,特別な指定はせずに,PCMのキャビティは原子核に追随させる.次に,キャビティ固定のもと振動数計算を実行する.Figure 2にGAMESSを用いる場合の入力を示す.対象系は水中の水分子,計算レベルはMP2/cc-pVDZである.振動数計算のために$CONTRLネームリストでRUNTYPE = HESSIANを指定し,$FORCEネームリストでMETHOD = SEMINUMと指定して解析的なエネルギー勾配を用いた数値的なヘシアン行列の計算を実行する.$PCMネームリストでNESFP = (キャビティ数)を指定して,$PCMCAVで実際のキャビティの座標および半径,スケール因子を入力する.通常,キャビティの個数は原子数に等しく,それらの座標も原子の座標と同じものを入力する.

Figure 2.

 Sample input to perform the frequency calculations for HSM by using GAMESS program.

Figure 2の計算が正常に終了すると,Figure 3のように基準振動数が計算される.GAMESSのプログラム上の問題で振動数に対するソートは行われないが,振動数より擬回転・擬並進のモードが判別できる.Figure 3の場合,∼200 cm−1の3つのモードが擬回転,∼70 cm−1の3つのモードが擬並進である.1619.59, 3832.69, 3936.36 cm−1はそれぞれ変角,対称伸縮,逆対称伸縮に対応する振動モードである.

Figure 3.

 GAMESS output for the sample input in Figure 2.

熱力学量の計算は,Figure 3の振動数さえ得られれば,Table 1EvibSvib,(5)式のErot/trans,(6)式のSrot/transにより関数電卓やエクセルなどを用いて簡単に計算できる.我々の研究室では,独自コード [6]を用いて行っている.

RB-HSMによる熱力学量の計算手順は,もう少し複雑である.まず,標準的な量子化学計算プログラムを用いて,PCMなどの溶媒効果を考慮して,構造最適化および振動数計算,熱力学計算まで一通り実行する.そして,EvibSvibに関しては,この結果を用いる.次に,溶質分子を剛体として取り扱ってヘシアン行列を求める.そのためにまず,剛体の回転・並進の6自由度に対する前方・後方差分の12構造におけるエネルギー勾配を解析的微分により求める.ここで,前方・後方差分の12構造は,独自コード [7]を用いて発生させる.Figure 4に,解析的微分によるエネルギー勾配をGAMESSで計算する際のインプットを示す.$CONTRLネームリストでRUNTYPE = GRADIENTを指定し,$DATAネームリストで前方・後方差分の12構造を指定する.よく見ると,すべての原子のx座標が0.01 Åだけシフトしていることがわかるであろう.この計算は,12回繰り返す必要がある.キャビティに関する入力は,Figure 2のHSM計算の場合と同様である.$FORCEネームリストの入力は不要である.

Figure 4.

 Sample input to evaluate the energy gradient for RB-HSM by using GAMESS program.

12個の構造に対して,Figure 5のようなエネルギー勾配が得られる.

Figure 5.

 GAMESS output for the sample input in Figure 4.

次に,12構造のエネルギー勾配から数値的にヘシアン行列を求め,さらに力の定数行列を対角化して基準振動数を求める.最後に,(5),(6)式よりErot/transおよびSrot/transを計算する.我々の研究室では,この一連の作業も独自コード [7]により実行している.

4 応用計算

Table 3に,標準状態(298.15 K, 1 atm)における気相および液相の水分子の熱力学量を示す.液相に関しては,通常の量子化学計算で採用されているIGMの結果と我々が開発したHSMおよびRB-HSMの結果を比較している.計算レベルは,構造最適化および振動数計算がMP2/cc-pVDZ,1点計算がCCSD (T)/cc-pV6Zである.気相と液相のエンタルピーの主な差は,電子エネルギーEelecに由来している.すなわち,溶媒効果による安定化である.一方,気相と液相のエントロピーの差は,回転・並進エントロピーTSrot, TStransに由来している.気相分子の方が,これらの運動の制約が少なく,結果として大きな値となる.IGMによる取り扱いは,液相においても気相におけるTSrot, TStransとほぼ同じ値であり,HSMと比較して∼30 kJ/mol程度過大評価していることがわかる.近似的な取り扱いであるRB-HSMは,まずまずHSMの値と近い値を与えている.

Table 3.  Thermodynamic properties of H2O in gas-phase by IGM and those in condensed phase by IGM, HSM, and RB-HSM.

Figure 6は,気相および液相における水分子の自由エネルギーGの温度依存性である.絶対零度の気相分子の値を基準にしている.液相に関しては,IGMとHSMの両方の結果を示している.計算レベルは,CCSD (T)/cc-pV6Z//MP2/cc-pVTZである.(1)式からFigure 6のグラフの傾きは主にエントロピー項Sに由来している.気相分子のSTable 3に示したように値が大きく,傾きが急である.従来法のIGMでは,液相においてもSが過大評価され,急な勾配となる.結果として気相と液相のグラフが交わらない.これは0-800 Kの温度範囲で常に液相の方が安定であることを意味する.すなわち,沸点に達しないという誤った結論を導く.一方,HSMでは液相のSは気相の値より小さくなり,結果としてグラフの傾きがなだらかになる.これにより,382.8 K (109.6 °C)で気相と液相のグラフが交差し,概ね実験値の沸点を再現している.詳細な議論は,文献 [1]を参照されたい.

Figure 6.

 Temperature dependence of free energies for gaseous and liquid water.

最後に,化学反応に対する温度依存性の例を示す.エタノールアミンとCO2からプロトン化アミンとカルバメートができる反応(7)を考える.

   2RNH2(aq) + CO2(g) ↔ RNH3+(aq) + RNHCO2(aq)(R = -CH2CH2OH)   (7)   

この反応は,火力発電所や製鉄所から排出されるCO2を分離回収するCCS技術で有力とされている化学吸収法で活用されている.化学吸収法では,吸収塔において排ガスからCO2をアミン水溶液に吸収させた後,放散塔で加熱することにより気体のCO2を分離回収する.吸収塔と放散塔の通常の運転温度は,それぞれ∼40 °Cおよび∼120 °Cである.すなわち,低温の吸収塔では(7)式の正反応が,高温の放散塔では逆反応が起こる.Figure 7に,この反応の反応物と生成物の自由エネルギーの温度依存性を示す.低温では生成物の方が反応物より自由エネルギーが低く,自発的に正反応が起こることがわかる.一方,120 °C以上の高温では安定性が逆転し,逆反応が自発的に進行することがわかる.

Figure 7.

 Temperature dependence of free energies for CO2 chemical absorption method.

5 まとめ

凝縮系の熱力学計算は,広く普及している量子化学計算プログラムを用いると,誤った取り扱いのため大きな誤差を与えることを実例とともに見てきた.さらに,この点について我々が提案したHSMを紹介し,凝縮系の熱力学量に対して化学的精度での予言能を有することが示された.ここで注意すべきことは,電子エネルギーEelecを高い精度で計算することが必須であるということである.基底関数および電子相関に対する収束値を用いない場合,誤差の相殺などから誤った理解を導きかねないことに注意されたい.

我々の研究室では,HSMは凝縮系の取り扱いにおいて一つの要素技術として定着しつつある.本稿では,HSMを用いた熱力学計算の具体的な手順をできる限りわかりやすく解説したつもりである.興味をもった読者は,まずはサンプルインプットを用いて実際に計算してみて頂きたい.そして是非,それぞれの研究でHSMを活用してもらいたい.

Acknowledgment

本研究は, 独立行政法人科学技術振興機構・戦略的創造研究推進事業 (JST-CREST) 「相対論的電子論が拓く革新的機能材料設計」, 文部科学省・元素戦略研究拠点「実験と理論計算科学のインタープレイによる触媒・電池の元素戦略研究拠点 (ESICB) 」の支援を受けて実施された. また, 量子化学計算の一部は自然科学研究機構 (NINS)・計算科学研究センター (RCCS) の計算機を利用して行った.

参考文献
 
© 2017 日本コンピュータ化学会
feedback
Top