06_open shell - YoshioNishimoto/sandbox GitHub Wiki

開殻系の計算

開殻系の計算を行います。 これまでは水素分子、窒素分子、水分子、エチレン･･･というように、どれも偶数の原子を持つ系を取り扱ってきました。 量子化学の授業で勉強したように、それぞれの軌道には二電子ずつ詰めるのが都合が良さそうなので、偶数の電子を持つ系は多くの場合そのように電子を詰めていきます。 このようにすべての軌道に二電子ずつ詰まっているような系のことを「閉殻系」や「閉殻分子」と言う場合があります。 逆に開殻系というのは、すべての軌道に二電子ずつ入っていない系のこととなります（e.g. 水素原子や各種ラジカル）。

閉殻と開殻はWikipediaによると、それぞれ「原子の最外殻に最大数の電子が入っている状態」と「外殻に最大数の電子が入っていない状態」のことと書かれています。 量子化学計算の場合は外殻ではなく分子軌道と読み替えれば良いと考えられます。 例えば水素原子は一個の電子を持っていますので開殻、 水素分子は二個の電子を持ちσ軌道に二個の電子が詰まるので閉殻と言われます。 偶数原子の系が総じて閉殻になるような気がしますが、例えば酸素分子は16個の電子を持っているので偶数電子の系であるにもかかわらず、基底状態は開殻三重項です（縮重した二つのπ*軌道に上向きスピンの電子が一つずつ入っている感じ）。

ここでは開殻系の計算として、非制限（unrestricted）波動関数と制限開殻（restricted open-shell）波動関数というものを扱います。 前者は往々にしてUHF (unrestricted Hartree--Fock)、後者はROHF (restricted open-shell Hartree--Fock)と呼ばれたりします。 詳細はSzabo: Section 2.5や3.8に書いてあると思います。 簡単に書くと、UHFはα軌道とβ軌道で別々の空間軌道（=分子軌道係数）を用い、 ROHFはRHFと同様にα軌道とβ軌道で同じ空間軌道を用いますが、いくつかの分子軌道の電子占有数が1になります。 RHFやROHFに比べるとUHFは変分空間が広がるため、RHFやROHFでは記述できない電子状態（それが正しいかは置いておいて）をカバーできます。

UHFのデメリットは、波動関数が全スピン演算子の二乗の固有値が正しくならない（スピン汚染）という問題があります。 一方でROHFはスピン汚染の問題を解決しますが、SCFが収束しづらいと言われています。 UHFかROHFのどちらが良いかと言われると何とも言えませんが、収束するならROHFを選ぶべきと思います（あるいは多配置SCF計算を行う）。

ここまでHartree--Fock法をベースにして書いてきましたが、DFTの場合でも開殻系の計算を用いることができます。 UHFとROHFのことを、それぞれUKSとROKSと呼んだりたりします（HFをKSに変えただけ）。

入力ファイル

というわけで、まずはUHF波動関数を用いてOHラジカル（二重項）の計算をしてみます。 後でROHF波動関数を用いた計算もご紹介します。

  1 from pyscf import gto, scf
  2 
  3 mol = gto.M(verbose=4,
  4             atom= [["O",(0.0, 0.0, 0.0)],
  5                    ["H",(0.0, 0.0, 1.0)]],
  6             basis="cc-pVDZ",spin=1)
  7 mf = scf.UHF(mol)
  8 mf.kernel()

とりあえず、6行目でspin=1というのを追加しました。 不対電子が一つある（二重項; スピン多重度はM = 2S + 1で表す）という意味です。 7行目は、これまでmf = scf.RHF(mol)としていましたが、ここでは非制限の計算をするためmf = scf.UHF(mol)としています。 なお、ここでは実際の計算をしていませんが、DFTの場合はmf = dft.UHF(mol)とすることで非制限のDFT計算を行うことができます（最初にimport dftも必要）。 入力ファイルを変化させるのはこれぐらいです。

出力ファイル

上の入力ファイルでPySCFの計算をすると、たぶん次のような感じになります。

 41 ******** <class 'pyscf.scf.uhf.UHF'> ********
 42 method = UHF
... （省略）
 56 number electrons alpha = 5  beta = 4
 57 Set gradient conv threshold to 3.16228e-05
 58 init E= -74.9994176773109
 59   alpha nocc = 5  HOMO = -0.4275772974238  LUMO = 0.133303166176253
 60 
 61 WARN: beta  nocc = 4  HOMO -0.427577297423801 >= LUMO -0.4275772974238
 62 
 63 
 64 WARN: system HOMO -0.4275772974238 >= system LUMO -0.4275772974238
 65 
 66 cycle= 1 E= -75.3762184403704  delta_E= -0.377  |g|= 0.185  |ddm|= 0.849
 67   alpha nocc = 5  HOMO = -0.527505653630867  LUMO = 0.165409437526525
 68   beta  nocc = 4  HOMO = -0.489211980105246  LUMO = 0.14079262608657
... （省略）
 93 cycle= 10 E= -75.3923228369339  delta_E= -1.84e-10  |g|= 1.75e-06  |ddm|= 2.25e-05
 94   alpha nocc = 5  HOMO = -0.543984216401393  LUMO = 0.179602376930914
 95   beta  nocc = 4  HOMO = -0.497936040208671  LUMO = 0.138967196052139
 96 Extra cycle  E= -75.3923228369366  delta_E= -2.66e-12  |g|= 6.6e-07  |ddm|= 2.18e-06
 97 converged SCF energy = -75.3923228369366  <S^2> = 0.75504013  2S+1 = 2.0050338

42行目method = UHFになっており、UHF計算をしていることが分かります。 あとは56行目で、α電子が5個、β電子が4になっていることを確認してください。 "WARN:"は、たぶん気にしなくて良いと思います。 初期推測でHOMOのエネルギーがLUMOのエネルギーよりも高くなっている場合に表示されるようです。 SCFが収束した後にこうなっていれば問題ですが、初期推測時点では問題ありません。 94・95行目でα軌道とβ軌道のHOMOとLUMOエネルギーが表示されています。 これらのエネルギーが異なっているため、確かにαとβで別々の空間軌道を用いていることが分かります。

エネルギー（converged SCF energy = -75.3923228369366 (hartree)）の隣に、<S^2> = 0.75504013 2S+1 = 2.0050338が表示されています。 この<S^2>が、スピン汚染の尺度（α軌道とβ軌道の違いの尺度）を表すもので、 これが理論値（二重項では0.75 (=S(S + 1))、三重項では2.0）から離れるに従ってスピン汚染が強くなっていると言えます。 スピン汚染が強い場合は、出てきた結果を正しいと信じる強い心を持つか、多配置計算を行うのが良いです。 具体的な計算式はWikipediaやSzabo: Section 2.5.3等をご参照ください。 2S+1は、おそらく<S^2>からSを計算し、それで2S + 1を計算し直していると思われます（通常のソフトウェアは<S^2>を計算するだけ）。

ROHF

一方でROHF計算は、入力ファイルの7行目を変えます（UHF→ROHF）。

 7 mf = scf.ROHF(mol)

すると、出力ファイル中の42行目が変化（method = ROHF-RHF）し、 56行目も次のようになります

num. doubly occ = 4  num. singly occ = 1

二電子が詰まった軌道が4つと、一電子入った軌道が1つの計算を行っています。 <S^2>が表示されていませんが、これは理論値となることが保証されているためです。

こちらのエネルギーはconverged SCF energy = -75.3884393376133 (hartree)で、UHFの結果と比べると少しエネルギーが高いことが分かります。

系によっては、開殻一重項が安定ということが知られています（多核金属錯体はその傾向が強い）。 片方の金属にα電子スピンが集まり、もう片方の金属にβ電子スピンが集まっているような雰囲気です（M(↑↑↑↑)......M(↓↓↓↓)）。 開殻一重項の計算をしたいところですが、やり方が分かりません。 自分でプログラムを作成しないとできないとか？

演習問題

  • ROHFのエネルギーがUHFのエネルギーよりも高くなるのは正しいでしょうか。正しければROHFのエネルギーが高くなる定性的な理由を、正しくないならばなぜ正しくない結果を与えるかの理由を指摘しましょう。

  • H2をRHFとUHF（またはROHF）で計算（Hartree--FockではなくDFTを用いても良いです。その場合汎関数は任意）して、横軸を原子間距離、縦軸をエネルギー（水素原子×2をエネルギーゼロと定義するのが良い）としてプロットしてみましょう。どちらが正しそうでしょうか。

  • どのような系の計算をすると、スピン汚染が強くなってしまうでしょうか。適当な分子（と計算レベル・計算条件）を用いて計算し、できるだけ理論値(S(S + 1))から離れるようにしましょう。

  • 好きな基底三重項分子を見つけましょう。一重項が基底状態とならないことを示す必要があると思われます。

  • Fe2の最も安定な（エネルギーが低い）スピン状態を見つけてみましょう。HFよりDFTを用いる方が良いとのことです。SCFが収束しない場合は、SCFのマニュアルに書いてあるオプションをいろいろ試して、できるだけ収束させられるように頑張りましょう。手法は任意（Hartree--Fockでも良いし、任意の汎関数によるDFTでも良い）。

⚠️ **GitHub.com Fallback** ⚠️