入力ファイルからの続きです。 次の入力ファイルでの計算が正しくできたことを前提にしています。

from pyscf import gto, scf 
mol = gto.M(atom='H 0 0 0; H 0 0 1.2', basis='ccpvdz')
mf = scf.RHF(mol)
mf.kernel()

converged SCF energy = -1.06111199785749が得られているかと思います。 これが、いわゆる分子の「エネルギー」です。 全エネルギーとも言います。 この全エネルギーは核間反発、電子-核間相互作用、電子間反発、そして電子の運動エネルギーを含めたエネルギーです。 単位は原子単位（atomic unit）またはハートリー（hartree）です。 基本的にはこのエネルギー（絶対エネルギー; absolute energy）をグループでのゼミ、学会、論文においてそのまま報告することは少ないです。 いくつかの構造や状態のエネルギー差（相対エネルギー; relative energy）を計算して、さらに単位を変換（多くの場合、kcal/mol、kJ/mol、eV）してから報告することが多いです。

ですが、この計算だとエネルギーしか表示してくれないので、何をしているのかは何も分かりません。 もっとたくさんの情報を得られるようにしましょう。

from pyscf import gto, scf 
mol = gto.M(atom='H 0 0 0; H 0 0 1.2', basis='ccpvdz', verbose=4) # この行にverbose=4を追加
mf = scf.RHF(mol)
mf.kernel()

二行目に、とりあえずverbose=4を追加してみました。 これでもう一度同じ計算をしてみましょう。 すると、今回は66行ぐらいのアウトプットになるかと思います。

  1 #INFO: **** input file is /data/home/nisimoto/lib/pyscf/test/test/verbose4.py ****
  2 from pyscf import gto, scf, mp
  3 mol = gto.M(atom="H 0 0 0; H 0 0 1.2", basis="ccpvdz", verbose=4)
  4 mf = scf.RHF(mol)
  5 mf.kernel()
  6 #INFO: ******************** input file end ********************
  7 
  8 
  9 System: uname_result(system='Linux', node='yn1.fukui.kyoto-u.ac.jp', release='3.10.0-957.1.3.el7.x86_64', version='#1 SMP Thu Nov 29 14:49:43 UTC 2018',     machine='x86_64', processor='x86_64')  Threads 6
 10 Python 3.6.8 (default, Aug 10 2019, 06:54:07)
 11 [GCC 4.8.5 20150623 (Red Hat 4.8.5-36)]
 12 numpy 1.12.1  scipy 0.18.1
 13 Date: Fri Apr 10 16:31:55 2020
 14 PySCF version 1.7.1
 15 PySCF path  /usr/local/lib64/python3.6/site-packages/pyscf
 16 
 17 [CONFIG] conf_file None
 18 [INPUT] verbose = 4
 19 [INPUT] num. atoms = 2
 20 [INPUT] num. electrons = 2
 21 [INPUT] charge = 0
 22 [INPUT] spin (= nelec alpha-beta = 2S) = 0
 23 [INPUT] symmetry False subgroup None
 24 [INPUT] Mole.unit = angstrom
 25 [INPUT]  1 H      0.000000000000   0.000000000000   0.000000000000 AA    0.000000000000   0.000000000000   0.000000000000 Bohr
 26 [INPUT]  2 H      0.000000000000   0.000000000000   1.200000000000 AA    0.000000000000   0.000000000000   2.267671349478 Bohr
 27 
 28 nuclear repulsion = 0.4409810091
 29 number of shells = 6
 30 number of NR pGTOs = 14
 31 number of NR cGTOs = 10
 32 basis = ccpvdz
 33 ecp = {}
 34 CPU time:         0.41
 35 
 36 
 37 ******** <class 'pyscf.scf.hf.RHF'> ********
 38 method = RHF
 39 initial guess = minao
 40 damping factor = 0
 41 level_shift factor = 0
 42 DIIS = <class 'pyscf.scf.diis.CDIIS'>
 43 diis_start_cycle = 1
 44 diis_space = 8
 45 SCF conv_tol = 1e-09
 46 SCF conv_tol_grad = None
 47 SCF max_cycles = 50
 48 direct_scf = True
 49 direct_scf_tol = 1e-13
 50 chkfile to save SCF result = /data/home/nisimoto/lib/pyscf/test/test/tmp9t4rz7tj
 51 max_memory 4000 MB (current use 45 MB)
 52 Set gradient conv threshold to 3.16228e-05
 53 init E= -0.723311777193342
 54   HOMO = -0.384475369076883  LUMO = 0.0412923973727384
 55 cycle= 1 E= -1.06041428904839  delta_E= -0.337  |g|= 0.0437  |ddm|= 0.657
 56   HOMO = -0.489119550273032  LUMO = 0.0965414834755516
 57 cycle= 2 E= -1.06109707888331  delta_E= -0.000683  |g|= 0.00632  |ddm|= 0.027
 58   HOMO = -0.485303877789622  LUMO = 0.100214479170792
 59 cycle= 3 E= -1.06111199673962  delta_E= -1.49e-05  |g|= 4.67e-05  |ddm|= 0.0048
 60   HOMO = -0.485252323944246  LUMO = 0.100245919708501
 61 cycle= 4 E= -1.06111199785719  delta_E= -1.12e-09  |g|= 7.28e-07  |ddm|= 7.37e-05
 62   HOMO = -0.485256027853135  LUMO = 0.100245669077032
 63 cycle= 5 E= -1.06111199785749  delta_E= -2.99e-13  |g|= 4.36e-12  |ddm|= 1.23e-06
 64   HOMO = -0.485256027850804  LUMO = 0.100245669075863
 65 Extra cycle  E= -1.06111199785749  delta_E=    0  |g|= 8.01e-13  |ddm|= 5.91e-12
 66 converged SCF energy = -1.06111199785749

デフォルト（初期設定、あるいは特に何も設定しない場合）のverbose=0は出力が少なすぎるため、どんな計算を行ったのかが分からなくなってしまうので、個人的にはやめた方が良いと思います。 この引数は0（最小）から9（最大）まで変化させられます。 自分にとって都合の良い値を選んでください。

やや（書く方は）面倒ですが、出力の行を説明していきたいと思います。

これは説明の必要もないでしょう。 入力ファイルの内容をリピートしています。

これはシステムや用いているPySCFの設定等を表示しています。 まぁ、この辺りはどちらかと言えば開発者用の情報で、ここで説明する必要もないかと思います。

ここは、gto.Mを用いて入力した分子の情報を分析し、出力してくれています。 25行目・26行目でH₂を入力したことが明らかではないでしょうか。 ちなみに、AAはオングストロームを意味しています。 あとは、 verbose=4としたオプションが18行目に反映されていることが確認できるかと思います。 そして、原子の数（19行目）、電子の数（20行目）、電荷の数（21行目）、スピンの数（22行目）を内部で計算して表示してくれています。 最初はこれらの情報が自分の意図しているものと合っているか、毎回確認する方が無難です（原子の種類も含む）。 これらのセッティングが異なる場合には、エネルギーを比較することはできません。 このため、異なる組成の分子でエネルギーを比較することも、基本的にはありません。 例えば、ethaneとethyleneは原子や電子の数が異なるためエネルギーを比較することに意味はありませんが、 cis-1,3-butadieneとtrans-1,3-butadieneは配座が異なるだけなので、エネルギーの比較に意味があり、どちらが安定（存在確率が高い）かを議論することができます。

ちなみに、例えばOH^-を計算したいとかだったら、

mol = gto.M(atom="O 0 0 0; H 0 0 1.2", basis="ccpvdz", verbose=4, charge=-1)

と、露わにcharge=-1と電荷の数を指定する必要があります。 系の電子数が一つ増えます。chargeを指定しない場合は、charge=0が自動的に選ばれます。 一方OHラジカル（二重項; doublet）だったら、

mol = gto.M(atom="O 0 0 0; H 0 0 1.2", basis="ccpvdz", verbose=4, spin=1)

と、スピンの数（α電子とβ電子の差）をspin=1と入力します。 ちなみに、この場合はrestricted open-shellのHartree--Fockを計算することになりますが、詳細は開殻系の計算や適宜教科書を参照してください。

OHラジカルのように系の電子数が奇数になっている場合の計算には特に注意が必要で、 スピンの数まで正しく入力しないとエラーが出て計算ができません。 電子の数が奇数だと、軌道に電子を詰める際に二重項なのか四重項なのかが分からず、困ってしまうわけです。 同様に、電子の数が偶数の場合でも一重項なのか三重項なのかが分かりませんが、基本的には（どのプログラムでも）デフォルトとして一重項が選ばれるようになっています。 三重項の場合には、spin=2と指定する必要があるのでしょう。

核間反発（nuclear repulsion）のエネルギーを出力しています。 原子単位でのエネルギーです。 量子化学計算と言えど、この項は古典論での計算を行っています（原子核を量子論で扱う方法も存在する）。

Z_Aは原子Aの核電荷、R_ABは原子AとBの距離です。 電卓で計算してみましょう。

ここでは基底関数に関する情報を示しています。 基底関数の導入的な記述は、Szabo: Sec.3.6をご参照ください。 より詳しい情報はverbose=5にすると用いた基底関数（cc-pVDZ）についての情報を出力してくれます。

[INPUT] ---------------- BASIS SET ---------------- 
[INPUT] l, kappa, [nprim/nctr], expnt,             c_1 c_2 ... 
[INPUT] H
[INPUT] 0    0    [3    /1   ]  13.01             0.019685
                                1.962             0.137977
                                0.4446            0.478148
[INPUT] 0    0    [1    /1   ]  0.122                1
[INPUT] 1    0    [1    /1   ]  0.727                1

とは言え、基底関数の勉強をしないと理解しづらいです。 ここではs軌道が2つ、p軌道が1つ（1セット）、それぞれのH原子に乗っていると考えてください。 すると、殻の数（number of shells）は合計で6つになります。

30行目・31行目の"pGTOs"と"cGTOs"は、それぞれ"primitive Gaussian-type orbitals"と"contracted Gaussian-type orbitals"の意味です（NRはnon-relativistic?）。 基底「関数」というのは、ガウス関数というわけです。 で、pGTOsとcGTOsの計算の仕方ですが、まずcGTOsから書きます。 この"cGTOs"の"c"は"contracted"に対応しており、これは軌道をガウス関数の線形結合で表すことを意味しています。 一つ前の段落で、それぞれのH原子にはs軌道が2つ、p軌道が1つがH原子に乗っていると書きました。 s軌道は方向依存性がないので、s軌道2つは2つの原子軌道に相当します。 一方でp軌道はx、y、zの方向性を持つため、p軌道1つ（1セット）は、3つの原子軌道に相当します。 なので、一つの原子には5つの原子軌道が載っているため、合計でcGTOsは10となります。

一方でpGTOsは、用いているガウス関数の合計です。 上の表より、一つ目のs軌道（l = 0）は3つのガウス関数（expnt = 13.01, 1.962, 0.446）の線形結合で表現しており、2つめのs軌道は1つのガウス関数（expnt = 0.122）を用いています。 s軌道のpGTOsは4です。 p軌道に関しては、これも一つのガウス関数を用いていますが、やはりx、y、z方向があるため、p軌道のpGTOsは3です。 なので、一つの原子には7つのpGTOsが載っているため、合計でpGTOsは14となります。

数え方はあまり重要ではありませんが、エネルギーを比較する際は上のように原子の数や種類、電子の数だけではなく、 基底関数の数（cGTOsで数えることが多いと思います）が一致していることも確認しましょう。 前回の演習問題ではcc-pVDZ、cc-pVTZ、cc-pVQZの基底関数を用いて「相対エネルギー」を計算しましたが、基底関数の数が違う以上、大して意味のある計算ではないです。 単位の変換の練習です。

32行目は基底関数の名前を表示しています。 33行目は、ECP (effective core potential)です。 内殻電子を有効ポテンシャルで置き換えた計算を行うときに用いられます（例えば重原子･･･内殻電子は化学反応にあまり寄与しないため）。

具体的に行うHartree--Fock計算の詳細です。 全て説明できるわけではありません。 マニュアルに書いてあるものもありますので、ある程度重要そうな行・語句をご紹介します。

restricted Hartree--Fock（制限HF）を行うという宣言です。 これは、α軌道とβ軌道で同じ空間軌道を用いるという意味です。 一方で、unrestricted（非制限HF）では、別の空間軌道を用いる計算で、スピンがある場合に用いられることが多いです。 さらに、既出ですがrestricted open-shell HFという方法が存在し、スピン汚染がないメリットがありますが、SCFの収束（後述予定）がしづらいといわれています。 Szabo: Sec.2.5.3をご参照ください。

SCF（self-consistent field; 自己無撞着（むどうちゃく）・自己無矛盾？）の初期推測・初期軌道のことです。 SCF計算は、最小のエネルギーを与えるような分子軌道を計算する手続きのことです（Szabo先生の教科書で勉強することになると思います）。 これは解析的に求めることができない（一般解を数式で導出することができないという意味）ので、繰り返し計算を用いて数値的に方程式を解きます。 このような繰り返し計算では、最小値に近い解を用いて計算を始めれば、より早く（少ない繰り返し計算で）解が見つけられるような気がすると思います。

最初のサイクルはどのような軌道を用いれば良いかというと、それほど自明ではありません。 多くのプログラムでは、拡張ヒュッケル法で計算した軌道をSCFの初期軌道として用いることが多いです。 ですが、もしも過去に行った計算結果が残っており、似たような分子軌道が得られていれば、その分子軌道を用いて計算を始めるのが良いでしょう。 そこでinit_guessというオプションを使うと便利です。 例えばこんな感じ。

from pyscf import gto, scf, mp
mol = gto.M(atom="H 0 0 0; H 0 0 1.2", basis="ccpvdz", verbose=4)
mf = scf.RHF(mol)
mf.chkfile = "h2.chk" # この行で軌道などの情報をh2.chkに保存
mf.kernel()

mol = gto.M(atom="H 0 0 0; H 0 0 1.205", basis="ccpvdz", verbose=4)
mf = scf.RHF(mol)
mf.chkfile = "h2.chk" # この行で軌道などの情報をh2.chkから読み出せるようにする
mf.init_guess = "chkfile" # 初期軌道をチェックポイントファイル（h2.chk）から読み出すようにする
mf.kernel()

なんかエラーメッセージが出ましたが･･･。 initial guessを用いると、143行目のdelta_Eが、用いない場合に比べて小さくなっているはずです。 他にもいろいろ方法がありますので、PySCFのインストールディレクトリからexamples/scf/14-restart.pyをご参考ください。

direct inversion in the iterative subspaceのことです。 SCFを収束させるのは思ったよりも難しく、単純に繰り返し計算をするだけでは大抵収束しません。 そのため、DIISという方法を用いてSCFの収束をしやすく（速く）しています。

二電子積分の計算に関連しています。 direct SCFは、SCFのサイクルごとに二電子積分を計算しています。 direct SCFでない場合は、最初に二電子積分を計算し、メモリ上やディスク上に保存し、SCF中では二電子積分を計算することなくFock行列を計算していきます。 計算効率的には後者が優れていますが、中程度の分子でもメモリ・ディスクがたくさん必要になる（4乗に比例して大きくなります）ので、最近はdirect SCFが主流です。

51行目のmax_memory 4000 MB (current use 45 MB)という部分です。 量子化学計算は基底関数（や活性空間）の数が大きくなると、メモリが大量に必要になる場合があります。 ディスクに大量に書き込む手法やアルゴリズムもありますが、ここでのメモリはRAM (= random access memory; 普通のノートパソコンだと16 GB以下のアレです）のことです。 多くの量子化学計算は線型代数の手法（つまり行列演算）を用います。 行列は二次元のメモリが必要になるため、基底関数（縦や横の要素数）が2倍になると、その行列を保存するのに必要なメモリは4倍になります。 テンソル（三次元以上の行列(?)）を扱う手法もあるので、そのような場合にはメモリが足りなくなる恐れがあります。

max_memoryの方は、PySCFが使うことのできる最大値です。約4 GBです。 これを超えると、他のプログラムに影響を及ぼす恐れがあるため、自動的にエラーで強制終了することになると思います。 リミットを変更することはできる（.pyscf_conf.pyで環境変数を変更するらしい 参考）のですが、この演習ではこのリミットが問題になることはないと思います。 一方、current use 45 MBというのは、この計算で使うメモリが45 MBだと言っています。 たかが2原子なので、まだまだ余裕がありますね。

これがSCF計算で、（広義の）Hartree--Fockの計算そのものです。 詳細はSzabo: Sec. 3.4.6に書いてあります。 前述の通り、Hartree--Fock計算は例えば二次方程式みたいに解析的な答えがあるわけではありませんし、 複雑な計算を一度行えば答えが得られるわけでもありません。 Hartree--Fock方程式が非線形な方程式であるため、繰り返し計算を通じて方程式を解きます（答えが得られない場合もあります）。

• E: そのサイクルでのエネルギー
• delta_E: 前のサイクルとのエネルギー差
• |g|: 電子状態の（？）勾配ノルム
• |ddm|: 分かりません。前のサイクルとの密度行列差？
• HOMO: HOMOの軌道エネルギー（単位はおそらくeV: 27.21138 eV = 1 a.u.）
• LUMO: LUMOの軌道エネルギー

最終的に、delta_EがSCF conv_tol = 1e-09より小さくなると、「SCFが収束した」と表現し、エネルギーが決まります。 収束しない場合は、収束するまであの手この手を尽くすか、分子構造を見直すか、諦めて「SCFが収束しませんでした」と報告します。 Hartree--Fock法は一般的にHOMO--LUMOギャップが大きくなる傾向があるため、SCFの収束は容易です。 Hartree--Fock法でSCFが収束しない場合には、何かが致命的に間違っている可能性が高いです。

ここまで計算して、やっとエネルギーを得られます。 この一行･･･倍精度の数字一つを得るために、Szabo先生の教科書約230ページ分（英語版）を勉強して、膨大なプログラムを記述する必要があるというわけですね。 その辺りを全て省略して、PySCFというプログラムを用いてエネルギー計算をしました。

以上がPySCFを用いたHartree--Fockの出力ファイルの説明になります。 何か適当な分子を作成し、計算してみましょう。

• 核間反発のエネルギーを再現（6桁ぐらい）する計算式を考えてみましょう。どのような単位変換が必要でしょうか。

• 量子化学の授業では、水素原子のシュレディンガー方程式を解くと、スレーター型関数で軌道を表現すると勉強したと思います。 しかし、多くの量子化学計算プログラムではガウス型関数を用いて軌道を表現します。 スレーター型関数と比較して、ガウス型関数は距離が増えるに従って急激に減少するという違いがあるにもかかわらず、 量子化学計算にわざわざガウス型関数を用いている理由は何でしょうか。

• H₂の距離を0.5から5.0 Angstromぐらいまで変化（0.05刻み）させて計算してみましょう。 さらに、計算して得たエネルギーを、横軸を距離、縦軸をエネルギーにしてプロットしてみましょう。 ポテンシャルエネルギーカーブを描くことができるようになります。 ただし、きちんとグラフが成立するために必要な情報を書くようにしてください。 ちなみに、RHF/cc-pVDZだと0.75 Angstromに極小点が見つかります（構造最適化による値は0.748 Angstrom）。 実験値は0.74 Angstromのようです（Wikipedia）。 距離を変化させるのは手動ではなく、forループを用いることができると良いと思います。

• H₂の0.748 Angstromでの結合エネルギーを計算しましょう。 結合エネルギーは、水素原子のエネルギーの倍とH₂のエネルギーの差として計算することができます。 つまり、2E_H - E_H2で計算することができます。 このために、水素原子のエネルギーを追加で計算する必要があります（電子の数に注意）。 Szabo先生のテキスト（英語版の165ページ）によると、実験値は4.75 eVだとか。

• 水素原子間を100 Angstrom程度に離したときのH₂のエネルギーは、2つの水素原子を別々に計算したときのエネルギーの和に一致するでしょうか。

• Koopmansの定理によりH₂のイオン化エネルギーを予測してください。Szabo先生のテキスト（英語版の194ページ）によると、実験値は0.584 a.u.だとか。