SCFが収束すると、分子軌道係数や電子密度が手に入ります。 これらの値を用いて化学的な直感を用いて理解するのに役立つ解析を行うことができます。 ここでは以下の計算を行ってみます。

• 二極子モーメント・電荷の計算
• 調和振動解析/基準振動解析
• 局在化軌道の計算

本当は「結合次数解析」というのも行いたかったです。 例えばエチレンのC-H結合とC-C結合は、結合次数解析の結果それぞれ1と2ぐらいになるような感じです。 ですが、PySCFにはたぶん実装されていないようです。 興味がある方は実装してみてください（演習問題参照）。

これらはSCFが終了した際に得られる電子密度を用いて計算できます。 二極子モーメント（dipole moment）は、エネルギーの電場による一次微分と関係しています。 電荷の計算には任意性があり様々な電荷の計算手法（Löwdin population, natural population, etc.）がありますが、ここではマリケン電荷（Mulliken charge）の計算します。

計算しますと言っても（HF/6-31G(d)レベルでの水分子の計算）、

  1 from pyscf import gto, scf
  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.047292503)]],
  7             basis="6-31g(d)")
  8 mf = scf.RHF(mol)
  9 mf.kernel()
 10 
 11 mf.analyze()

11行目にmf.analyze()と追加したように、この一行を加えるだけで二極子モーメントとマリケン電荷の計算をしてくれます。 出力として、以下が新しく追加されるかと思います（少し省略しました）。

 81 **** SCF Summaries ****
 82 Total Energy =                         -75.997330632174396
 83 Nuclear Repulsion Energy =               8.772355978237027
 84 One-electron Energy =                 -122.366294682993427
 85 Two-electron Energy =                   37.596608072582001
 86 **** MO energy ****
 87 MO #1   energy= -20.552158479917   occ= 2
 88 MO #2   energy= -1.30827382185148  occ= 2
 89 MO #3   energy= -0.706529305689577 occ= 2
...
102 MO #16  energy= 2.09223592458353   occ= 0
103 MO #17  energy= 2.63334045043919   occ= 0
104 MO #18  energy= 2.82729997931672   occ= 0
105  ** Mulliken pop on meta-lowdin orthogonal AOs  **
106  ** Mulliken pop  **
107 pop of  0 O 1s        1.99999
108 pop of  0 O 2s        1.63618
109 pop of  0 O 3s        0.00252
...
122 pop of  1 H 2s        0.01356
123 pop of  2 H 1s        0.65765
124 pop of  2 H 2s        0.01356
125  ** Mulliken atomic charges  **
126 charge of  0O =     -0.65758
127 charge of  1H =      0.32879
128 charge of  2H =      0.32879
129 Dipole moment(X, Y, Z, Debye): -0.00000, -0.00000,  1.95461

二極子モーメントとマリケン電荷以外にも多くの行が表示されており、せっかくなので簡単に説明していきます。 84行目の一電子エネルギー（電子の運動エネルギーと電子と原子核の相互作用）は

85行目の二電子エネルギー（電子と電子の相互作用）は

で計算していると思われます。 それぞれハミルトニアン中の一電子・二電子演算子の期待値に対応しています。 いわゆるelectronic energy は二項の和になりますが、 通常は全エネルギー（Total Energy; 電子エネルギー ＋ 核間反発エネルギー）を使って解析・議論をします。 電子相関を含める場合は、当然電子相関エネルギーも含める必要があります（単参照電子相関理論参照）。 ちなみにDFTの場合は

**** SCF Summaries ****
Total Energy =                         -76.361008266082152
Nuclear Repulsion Energy =               8.772355978237027
One-electron Energy =                 -122.441711847168946
Two-electron Coulomb Energy =           46.594772749151062
DFT Exchange-Correlation Energy =       -9.286425146301301

というように、DFTの交換相関エネルギー（DFT Exchange-Correlation Energy）が加えられます。 また、二電子エネルギーではなく、二電子クーロンエネルギー（Two-electron Coulomb Energy; 二電子エネルギーの式で、かっこ内の一項目）になります。 基本的にこれらのエネルギーの和はTotal Energyで表示されているので間違えることはないと思いますが、 どのようにエネルギーが計算されているか・どのエネルギーが寄与しているかを把握しておくことは重要と思います。

87から104行目で、軌道エネルギー（原子単位）と軌道の占有数を示しています。 ちなみに、軌道エネルギーの和を計算しても、電子エネルギーにはなりません（Szabo: Section 3.3.3辺りを参照？）。 単参照電子相関理論の演習問題にあるとおり、軌道エネルギーの和はMP0エネルギーに対応します。

で、マリケン電荷は126から128行目に表示されています。 酸素原子の電荷が-0.65758で、水素原子が0.32879あたりということが分かります。 酸素原子に負電荷が集まるイメージがあると思いますので、正しそうでしょう。 また、全ての原子の電荷を足すとゼロになります（系全体で中性となっている場合）。 その前のMulliken pop（Mulliken population; 107から124行目）は、 それぞれの原子軌道に何個の電子が入っているかが表示されます。 内殻に当たるO原子の1sはほとんど2個の電子に占有されていることや、 一般的にいう価電子の軌道外（と言えば良いのか分かりませんが）の3sはほとんど占有されていないことが分かります。 また、O原子は1s,2s,2p軌道だけで良いかと思いきや、やはり少しとは言えd軌道を用いて水分子の電子状態を表現していることに注意してください。 この「少し」のために、分極軌道を使っているわけです。

そして最後の行に二極子モーメントが表示されます。 単位はデバイです。 この水分子はx平面に配置しているため、x軸方向の二極子モーメントはゼロになります。 y軸方向も、対称性からゼロになります。 一方で、z軸はC₂軸に沿っているため、二極子モーメントの値は大きくなります。

･･･ですが、Mulliken populationやMulliken chargeは、GAMESS-USの値とかなり違うので、 一般的な計算アルゴリズムとは少し違う感じがします。 PySCFにはマリケン電荷以外でも計算できるので、気になる方はいろいろ試してみてください。

調和振動解析（harmonic vibrational frequency analysis?）または基準振動解析（normal mode analysis）をしてみます。 構造最適化のところで触れたとおり、 構造最適化後の構造に対して（同じ計算レベルで）これらの解析を行うことで、 その極値の性質を知ることができます。 すなわち、安定構造になっているか、鞍点（あんてん、saddle point）になっているかが分かります。 基準振動解析の結果、モードの波数が負になっている場合は、そのモードに対してポテンシャルエネルギー面は上に凸となります。 特に、波数が負になるモードの数が一つの場合、その極値のことを「遷移状態（transition state）」と言ったりします。 必ずしも毎回基準振動解析を行う必要があるわけではないですが、化学反応を扱う場合は解析を行うことが一般的です。 安定構造・鞍点になっていることを確認するのも一つの目的ですが、多くの化学反応はギブズエネルギーを用いて議論するため、ギブズエネルギーを計算するという目的もあります。

なお、一般的には構造最適化をした後で基準振動解析を行います。 極値の性質を知るという事もありますが、熱力学量の計算で一次微分がゼロになっているという近似を用いている（たしか）という理由もあります。 ですが、遷移状態の計算を始める前に、初期構造が遷移状態に十分近いかや、負にしたい振動モードが存在するかなどを確認するのに用います。 また、二次微分をプログラムに渡してやると、構造最適化がスムーズにできることもあるので、そのような使い方もあります。

とりあえず計算してみます。

  1 from pyscf import gto, scf
  2 from pyscf.hessian import thermo
  3 
  4 mol = gto.M(verbose=4,
  5             atom= [["C",( 0.0000000000, -0.0000000000, -0.6758096762)],
  6                    ["C",( 0.0000000000,  0.0000000000,  0.6415100439)],
  7                    ["H",(-0.0000000000,  0.9144074968, -1.2427442213)],
  8                    ["H",( 0.0000000000, -0.9144074968, -1.2427442213)],
  9                    ["H",( 0.0000000000,  0.9144074968,  1.2084445889)],
 10                    ["H",(-0.0000000000, -0.9144074968,  1.2084445889)]],
 11             basis="6-31G(d)")
 12 mf = scf.RHF(mol)
 13 mf.kernel()
 14 
 15 hessian = mf.Hessian().kernel()
 16 
 17 # Frequency analysis
 18 freq_info = thermo.harmonic_analysis(mf.mol, hessian)
 19 
 20 print('freq_wavenumber')
 21 print(freq_info['freq_wavenumber'])
 22 
 23 print('norm_mode')
 24 print(freq_info['norm_mode'])
 25 
 26 # Thermochemistry analysis at 298.15 K and 1 atmospheric pressure
 27 thermo_info = thermo.thermo(mf, freq_info['freq_au'], 298.15, 101325)
 28 
 29 print('Rotation constant')
 30 print(thermo_info['rot_const'])
 31 
 32 print('Zero-point energy')
 33 print(thermo_info['ZPE'   ])
 34 
 35 print('Internal energy at 0 K')
 36 print(thermo_info['E_0K'  ])
 37 
 38 print('Internal energy at 298.15 K')
 39 print(thermo_info['E_tot' ])
 40 
 41 print('Enthalpy energy at 298.15 K')
 42 print(thermo_info['H_tot' ])
 43 
 44 print('Gibbs free energy at 298.15 K')
 45 print(thermo_info['G_tot' ])
 46 
 47 print('Heat capacity at 298.15 K')
 48 print(thermo_info['Cv_tot'])

15行目でエネルギーの二次微分を計算しています。 この二次微分を行列の形にしたもの（ヘッセ行列と言います）を、 確か原子の重さでスケールして、対角化することで基準振動を得ます。 これらの操作を、18行目で行っています。 二次微分の計算にはcoupled-perturbed Hartree--Fockという方程式を解いています（Szabo: Appendix C参照。二次微分の計算には分子軌道係数の一次微分が入ってくるため）。 が、何も出力ファイルにはverboseの値に関わらず二次微分の詳細は何も表示されないようです。 かろうじてverbose=6にすると、計算ステップ毎の計算時間が表示されるようですが。

まず見ておくべきは、freq_info['freq_wavenumber']でしょう。 出力として、897.13から始まる12個の値が得られると思います。

freq_wavenumber
[ 897.13109677 1094.4184838  1098.31443855 1154.53183132 1352.73429729
 1496.50475888 1611.13761947 1856.79461874 3318.44199556 3342.02072501
 3392.84414439 3419.02324224]

これらは12個（なぜ12個でしょう？）の基準振動それぞれの、波数（単位はcm^-1）を表しています。 おそらくどこかの授業で勉強したと思いますが、3,000 cm^-1以上の振動は、水素原子（ここではC-H）の伸縮振動に対応してそうだということが分かるかと思います。 これらは、基準振動モード（print(freq_info['norm_mode'])）をMOLDENファイルとかにダンプし、 可視化ソフトウェアで見ることが出来れば良いのですが･･･PySCFではできない感じ？（いろいろ気が利かない）

27行目では、熱化学的な解析を行っています。 thermo.thermoの三つ目・四つ目の引数として指定しているとおり、ここでは298.15 K、1 atmでの結果を表示してくれます。 構造最適化計算で得られる構造・エネルギーは、 ポテンシャルエネルギー面での極小値（すなわちゼロ点エネルギーは無視されている）となります。 ですが、量子化学の授業で勉強したと思われるように、実際の分子は0 Kであってもゼロ点エネルギーを含んでいるので、その補正をしなければならない場合があります。 例えば有機化学の反応を扱う場合や実験値と比較を行う場合にはゼロ点エネルギー（とギブズエネルギーの補正）を含める事が多いです。 そのゼロ点エネルギーが、この分子では

213 Zero-point energy
214 (0.05475324671888602, 'Eh')                                                              

すなわち0.05475... hartreeとなります。 そして、0 Kの内部エネルギー（Internal energy at 0 K）は、 SCFで得られたエネルギーにゼロ点エネルギーを足したものとなっています。 さらに有限温度（ここでは298.15 Kでの）並進・回転・振動の寄与を含めた内部エネルギー（Internal energy at 298.15 K）や、 エンタルピー（Enthalpy energy...）・ギブズ自由エネルギー（Gibbs free energy...）を得ることもできます。 ちなみに、ここのHeat capacity at 298.15 Kは、モル定容（定積）熱容量のようです（GAMESS-USとの比較より）。

化学反応を議論するときは、ギブズ（自由）エネルギーを用いるのが普通です。 用いない場合もあります。 二次微分の計算はそれなりに時間がかかるため 系が大きいと諦めたりする場合もあると思われますが、 ゼロ点振動補正の影響はかなり大きい場合もあるのであまり無視したくはないです。 SCFエネルギーを使ってもギブズ（自由）エネルギーを使っても大きな差がないことも多いのですが、 分子の数が変わる場合は10 kcal/mol以上（温度次第）変化するため注意する必要があります。 論文などで報告するエネルギーがこれらの補正を含まない場合も多いのですが、 反応機構に関する論文では、ギブズ（自由）エネルギーを用いるのが基本です。

最近の熱力学・物理化学の教科書では「ギブズエネルギー」・「ヘルムホルツエネルギー」と書くようになってきていますが、 量子化学計算プログラムでは「ギブズ自由エネルギー」・「ヘルムホルツ自由エネルギー」と書く方が多いです。 このあたりの計算ルーチンは昔に書かれたプログラムをベースにしていることが多い気がするので。

通常の分子軌道（canonical orbital）は、フォック行列を対角化することにより得られる軌道です。 それぞれの分子軌道は直交化していることもあり、励起状態計算や電子相関の計算がやりやすいというメリットがありますが、 多くの場合いくつかの原子にまたがって電子密度が分布するため、視覚的にやや分かりづらいという場合があります。 一方、局在化軌道（localized orbital）は、 canonical orbitalを回転（他の軌道と混ぜると表現しても良い。また、ユニタリー回転）させることにより、例えばσ軌道一つに局在化させたような分子軌道となります。 局在化軌道は、通常占有軌道に対してしか得られません（非占有軌道でも擬似的に局在化軌道をつくる方法があります）。 また、直交しているわけではありません。 このため局在化軌道を用いたpost-HF計算は、いろいろ難しくなります。

PySCFでは、次の局在化アルゴリズムが実装されているようです（詳細はよく分かりませんが）。

• Foster--Boys
• Edmiston--Ruedenberg
• Pipek--Mezey
• meta-Löwdin
• natural atomic orbitals
• intrinsic atomic orbitals
• intrinsic bond orbitals

ここではPipek--Mezey orbital localizationを用います。 おそらく現代的には最もよく使われる方法で、CASSCFの初期軌道を作成するのに用いる場合もあります。

以下でそれっぽい計算ができます（が、GAMESSの結果と合わない。 局在化は一意に決まるわけではないのでどちらも正しいかもしれないが、 それでも以下が正しいかは怪しい。 Pipek--Mezeyはσ軌道とπ軌道が混ざらず計算できるメリットがあるはずなのに、それができていない･･･。） 。

  1 from pyscf import gto, scf
  2 
  3 mol = gto.M(verbose=4,
  4             atom= [["C",( 6.00888796e-13,-5.21236125e-14,-1.27727822e+00)],
  5                    ["C",(-5.73985844e-13, 1.59017062e-14, 1.21246132e+00)],
  6                    ["H",(-2.61455155e-13, 1.72747100e+00,-2.34915037e+00)],
  7                    ["H",(-6.70150000e-14,-1.72747100e+00,-2.34915037e+00)],
  8                    ["H",( 2.48003679e-13, 1.72747100e+00, 2.28433346e+00)],
  9                    ["H",( 5.35635242e-14,-1.72747100e+00, 2.28433346e+00)]],
 10             basis="6-31g(d)",unit="Bohr")
 11 mf = scf.RHF(mol)
 12 mf.kernel()
 13 
 14 import numpy
 15 from pyscf import lo
 16 class HubbardPM(lo.pipek.PM):
 17 # Construct the site-population tensor for each orbital-pair density.
 18 # This tensor is used in cost-function and its gradients.
 19     def atomic_pops(self, mol, mo_coeff, method=None):
 20         return numpy.einsum('pi,pj->pij', mo_coeff, mo_coeff)
 21 
 22 loc_orb_init_guess = mf.mo_coeff[:,0:8]
 23 #mol.verbose = 5
 24 locobj = HubbardPM(mol, loc_orb_init_guess)
 25 print('PM cost function  ', locobj.cost_function())
 26 loc_orb = locobj.kernel()
 27 print('PM cost function  ', locobj.cost_function())
 28 
 29 print (loc_orb)
 30 
 31 from pyscf.tools import molden
 32 with open('C2H4mo.molden', 'w') as f1:
 33     molden.header(mol, f1)
 34     molden.orbital_coeff(mol, f1, loc_orb, ene=mf.mo_energy, occ=mf.mo_occ)

これで出力されてくるC2H4mo.moldenを用いて局在化軌道を可視化してみましょう。 正しいか分かりませんが。 なお、局在化軌道に軌道の準位（軌道エネルギー）を定義することはできません。 可視化する際に軌道エネルギーが表示されるかもしれませんが、便宜上元のcanonical軌道の番号と対応して表示するようにしているだけです。

基本的な解析は以上のような感じでしょうか。 エネルギーをさらに微分するといろいろ分かります。 例えば、エネルギーを原子の座標と電場で（混合）微分して得られる値（dipole derivative）を振動モードに射影したりあれこれすると、IRスペクトルをシミュレーションすることができます。 他にもエネルギーを電場で三次微分すると、超分極率を得ることができ、振動モードに射影すると非共鳴ラマンスペクトルのシミュレーションをすることができます。 PySCFではできない解析手法も無数に存在するため、興味があればいろいろ調べて実装してみてください。

• マリケン電荷と二極子モーメントの基底関数依存性または汎関数依存性を調べてみましょう。
• 局在化軌道の計算は、占有軌道の分子軌道係数に適切なユニタリー行列をかけることで計算を行います。 通常の分子軌道を用いて計算したエネルギーと、局在化軌道を用いて計算したエネルギーが一致することを証明しましょう。 Hartree--Fockのエネルギーを使えば良いです。
• 適当な分子の基準振動解析を行い（手法は任意）、文献（実験でも計算でも良いです）を調べて波数を比較してください。
• 水分子二量体を構造最適化しましょう（手法は任意）。 基準振動解析を行い、振動の波数が実数であることを確認しましょう。 結合エネルギー、標準状態における二量化する際の内部エネルギー変化・エンタルピー変化・ギブズ自由エネルギー変化を計算してみましょう。 定温・定圧状況下で水分子の二量体は自発的に起こるでしょうか。
• HF/STO-3Gレベルで次の構造（1,3-butadiene + ethylene）の振動解析を行ってみましょう。 この構造はポテンシャルエネルギー面でどのような極値に相当するでしょうか。 freq_wavenumber中の"j"は、複素数iのことです。

C  -0.1211293313  -0.4457163427  -1.4171647487
C  -1.2739856967  -0.3091329584  -0.7145497446
C   1.0503510908   1.2853077067  -0.6792479073
H  -2.1635703465   0.0470080478  -1.2189851588
H   0.4273593840   1.9816728363  -1.2211864855
C  -1.2738735849  -0.3132249690   0.7066939692
C   1.0504420313   1.2813997884   0.6802190131
H  -2.1633816132   0.0399952783   1.2133133128
H   0.4275185790   1.9746347320   1.2262340812
C  -0.1209034547  -0.4538428678   1.4083288546
H  -0.1195898390  -0.3353062248   2.4848381172
H  -0.1199821020  -0.3209724999  -2.4929725269
H   1.9090334204   0.9219869922  -1.2249315557
H   1.9091996190   0.9149563046   1.2236916943
H   0.6961732817  -1.0393964109  -1.0372576914
H   0.6963436576  -1.0453168130   1.0248779901

結合次数解析で有名な手法は、"Wiberg bond-order analysis"と"Mayer bond-order analysis"とかでしょう。 どちらも密度行列（と、後者の場合は重なり行列）を用いるだけで、割と簡単と思います。 興味がある人はやってみてください。 Wiberg bond-order index (for restricted wavefunction)は、次の計算を行います。

AとBは原子のインデックスで、μとνは原子軌道のインデックスで、は密度行列です。 実際にプログラムする際には、係数等を若干微調整する必要があるかもしれませんが。 一方Mayer bond-order indexは、

（なぜか\mathbfが斜体になる）は、密度行列と重なり行列の積です。 すなわち

です。なので、Mayer bond-order indexの方がやや面倒ですが、おそらくこちらの方が直感に合う値になるのでしょう。 ちなみに、上の構造（1,3-butadiene + ethylene）を用いてMayer bond-order analysis (HF/STO-3G)を行うと、C1-C2間は1.633、C1-C3間は0.256、C-H間は大体0.97ぐらいになります（構造含め、GAMESS-USによる結果）。