09_solvent - YoshioNishimoto/sandbox GitHub Wiki

溶媒効果を取り入れた計算

これまでの計算は気相中あるいは真空中での計算でした。 ただ単に量子化学計算をしたいというのであれば、真空中での計算を扱っていても全く問題が無いのですが、 多くの有機化学反応や生体反応は溶媒の中で起こる反応を取り扱うため、溶媒効果を取り入れつつ、できる限り現実的な系を扱いたいところです。

とは言え、一つの溶質分子の周りには、大量の溶媒分子があることでしょう。 しかも、溶質分子と直接相互作用する距離にある溶媒分子（first coordination sphere）のみを考慮すれば良いかと言えばそうではなく、 適切に溶媒効果を取り込もうと思うと、量子化学計算のスケールからすれば到底計算不可能な数の溶媒分子を取り扱う必要があります。 とは言え、すべての溶媒分子を量子化学計算で取り扱うのではなく、溶媒分子をより計算コストの低い手法により取り扱うことで、 溶媒分子を露わに取り扱うことができます（explicit solvent）。 例えば溶媒分子を分子力場で扱う（QM/MM）、ポテンシャルで代替する（effective fragment potential）等の手法が挙げられます。 ある意味正攻法ではあるのですが、溶媒分子がきちんと平衡になっているかを考えたり、系全体の平均（統計）を取ろうと思うと、 計算コストが高くなりがちという問題があります。 そこで、多くの近似的用いつつ、より低い計算コストで溶媒効果を取り入れるアプローチとして、 連続誘電体モデルを用いて間接的に溶媒効果を取り入れます（implicit solvent）。

連続誘電体モデルの中でもっとも有名なのは、polarizable continuum model (PCM)というものです。 反応機構の量子化学計算で溶媒効果を取り入れる場合は、9割がこのモデルを用いているでしょう。 というのも、反応機構を扱う際にGaussianが使われることが多く、Gaussianでは通常PCMが用いられるためです。 一方、PySCFでは別の連続誘電体モデル、conductor-like screening model (COSMO)が利用できます。 conductor-like PCM (C-PCMとして知られる)と似た手法ではあるのですが、私もどのように違うのかまではきちんと把握していません。

PCMによる計算の詳細

COSMOはよく分かりませんが、PCMに関して概略を説明しておきます。 PCMとCOSMOは非常に良く似ており、以下で説明する点はおそらく二つの手法で共通していると思います。

まず、溶媒分子を「分極できる点電荷」で近似します。 この点電荷をtesseraと呼び（COSMOでは別の呼び方をするかもしれません）、溶質分子の分子キャビティー（溶媒分子がこれ以上溶質分子に近づけない、という領域）の表面に貼り付けておきます。 用いるtesseraの数は溶媒分子の数と等しいわけではなく、それぞれの溶質原子に大体50から200ぐらいのtesseraを設定することが多い印象です。 準備は一応これで終わりです。 溶媒毎に（比）誘電率 ((relative) dielectric constant)というのが決まっており、 tesseraの分極のしやすさはこの誘電率に依存するようにします。 比誘電率の値（実験値を用いる）を変えることにより、様々な溶媒をモデルしています。

SCFの計算中は、まず溶質分子の密度行列を得ます。 この密度行列により溶質分子の電荷密度分布（おおざっぱに言えば、原子の電荷）を計算することができ、 この電荷密度分布が溶媒分子を模した点電荷（上のtessera）を分極させます。 例えば水分子（溶質）の場合、O原子の周りの溶媒分子は正電荷を持つ原子が配向すると予想されます。 そこで、tesseraを分極させることにより、O原子の周りに正電荷が寄ってくる状況を再現します。 この分極した点電荷は、溶質分子に対して（ゼロではない）静電ポテンシャルを生み出すこととなり、 その静電ポテンシャルを考慮してFock行列を更新する必要が出てきます。 この一連の過程は、（溶媒効果とは関係なく）SCF計算そのものと似ています。 静電ポテンシャルの話を二電子積分と置き換えれば、 Hartree--Fock法でのSCF計算と同じような計算をしていることが分かるかと思います。 計算手続きとしてはFock行列を更新すれば良いだけで、ユーザー側からは真空中の計算とどのように違うかはなかなか意識しづらいです。

最終的な「エネルギー」は、自由エネルギーとして得られます。 普通に計算して得られる内部エネルギー（ゼロ点補正を加えたエネルギーのことではなく、単なる全エネルギーのこと）に、 溶質-溶媒間の相互作用エネルギーの「半分」を足すことで、自由エネルギーを定義しています。 なぜ半分かというと、残りの半分はtesseraを分極する仕事に利用されるとのことです。

ここで考慮することができるエネルギーは、静電的な効果のみです。 溶媒和自由エネルギーを計算するためには、 他にもcaviation, dispersion（密度汎関数理論で出てきた原子間の分散項とは独立）, repulsionという効果を取り入れる必要があります。 これらはまとめて非静電相互作用と言われています。 静電相互作用は電子状態に依存する一方、 非静電相互作用は分子の構造（とパラメータ）だけで決まります。 PySCFでは、この非静電相互作用は計算できなさそうです。 他の量子化学計算ソフトウェアでも、非静電項を明示的に計算するオプションを追加する必要がある場合が多く、 非静電項の存在は見落とされがちです。

PCMとCOSMOは溶媒分子を露わに扱うわけではなく、分極できる点電荷を用いて溶媒分子との相互作用を評価するため、水素結合を上手に評価できないと言われています。 また、混合溶媒・不均一溶媒を取り扱うことができない、溶媒のミクロ・マクロな構造を得られないなど、制限があります。 別の溶媒効果を扱う手法としてreference interaction site model (RISM)という手法があり、 これらの点を改善できると言われていますが、取り扱いの際には別の難しい点が多い印象です（個人の感想です）。 あと、量子化学計算ソフトウェアでなかなか出回らないです。

ちなみにSMD (solvation model based on density)というのがあります。 これはPCMの一種、integral equation formalism PCM (IEF-PCM)というもので静電相互作用を取り扱い、 非静電相互作用を独自の関数により評価していると思われます。

実際の計算

というわけで、とりあえず計算してみます。

  1 from pyscf import gto, dft, solvent
  2 
  3 mol = gto.M(verbose=4,
  4             atom= [["O",( 0.00000000, 0.000000000,-0.452707497)],
  5                    ["H",( 0.00000000, 0.866025404, 0.047292503)],
  6                    ["H",( 0.00000000,-0.866025404, 0.047929503)]],
  7             basis="6-31G(d)")
  8 mf = dft.RKS(mol)
  9 mf.xc = 'B3LYP'
 10 
 11 solvent.ddCOSMO(mf).run()

11行目は、通常であれば例えばmf.kernel()とかするところですが、solvent.ddCOSMO(mf).run()となります。 真空中での計算からの違いはこれ（+ 1行目のimport solvent）だけです。 デフォルトでの溶媒分子は水分子という事になっているようです。 また、ddCOSMOはdomain decomposition COSMOというものらしいです。詳しいことはよく分かりません。

これで計算をしてみると、次の出力が得られると思います。

110 cycle= 7 E= -76.3684731687422  delta_E= -1.94e-11  |g|= 5.77e-07  |ddm|= 1.49e-05
111   HOMO = -0.28386397758513  LUMO = 0.0665399705657929
112 <class 'pyscf.solvent.ddcosmo.DDCOSMO'> E_diel = -0.00796938008882138
113 Extra cycle  E= -76.3684731687422  delta_E= -2.84e-14  |g|= 4.04e-07  |ddm|= 9.66e-07
114 converged SCF energy = -76.3684731687422

112行目で溶媒効果を考慮するときに表示される行で、 E_dielというのが上で触れた「相互作用エネルギーの半分」（原子単位）です。 そして最終的に得られるconverged SCF energy = -76.3684731687422は、いわゆる内部エネルギーにE_dielを足した値です。 なので、エネルギーを報告する際には、このconverged SCF energyを使っておけば良いことになります。

ここで、溶媒効果がない場合（この値を得るためには、別で計算する必要があります）・ある場合での内部エネルギーと静電相互作用（原子単位）を一応比較してみると、

　 内部エネルギー 静電相互作用 全エネルギー
溶媒効果なし -76.361013 0 -76.361013
溶媒効果あり -76.360504 -0.007969 -76.368473

であり、溶媒効果が内部エネルギーに与える影響は少々（ここでは0.3 kcal/mol程度）ということ分かります。 ただ、全エネルギーで比較すると（比較して良いのか分かりませんが）、溶媒効果を入れた方がエネルギーが低くなり、 静電項的には水分子を水中に入れるのは有利と言えそうです。 もちろん現実にはエントロピーの変化も考慮する必要がありますが。

溶媒すなわち比誘電率を変えるには、例えば次のようにすればできそうでした。 ここでは比誘電率を24.852にしており、ethanolに対応します。

  1 from pyscf import gto, dft, solvent
  2 
  3 mol = gto.M(verbose=4,
  4             atom= [["O",( 0.00000000, 0.000000000,-0.452707497)],
  5                    ["H",( 0.00000000, 0.866025404, 0.047292503)],
  6                    ["H",( 0.00000000,-0.866025404, 0.047929503)]],
  7             basis="6-31G(d)")
  8 mf = dft.RKS(mol)
  9 mf.xc = 'B3LYP'
 10 
 11 mf = mf.ddCOSMO()
 12 mf.with_solvent.eps = 24.852
 13 mf.run()

というわけで、溶媒効果を追加するのは一行加えるだけで良く非常に簡単です。 一方、PCMやCOSMOの具体的な理論・計算アルゴリズムはそれだけでHartree--Fock法と同等ぐらいには複雑です。 また、SCF以上（電子相関や励起状態）の計算を行う際の溶媒効果はいろいろ難しくなります。 MP2とかでも一応使えるはずです（ここでは取り扱いませんが）。

演習問題

計算手法は適当に設定してください。

  • CH3Cl + Br- -> CH3Br + Cl-の反応エネルギーを、真空中・水中で計算してみましょう。 どちらの値が、より現実に近そうでしょうか。
  • 適当な1,3-ジカルボニル化合物（例えばアセチルアセトン）を用いて、ケト形・エノール形の存在比率が真空中・様々な溶媒中でどのように変化するか計算してみましょう。 存在比率はボルツマン分布の式を用いると良いでしょう。 温度などの条件は、適当に設定してください。
  • （某ページから拝借したテーマ）1,2-difluoroethaneのgauche conformationとanti (trans) conformationを真空中と溶媒中で計算を行い、相対エネルギーを計算しましょう。さらに、溶媒中での相対エネルギーの増減の理由を定性的に説明してください。
⚠️ **GitHub.com Fallback** ⚠️