第10章分子の量子化学計算を実践する

本章の計算は，すべて次の流れの上に乗っている．まず分子構造を用意し(10.2節)，構造緩和で平衡構造を求め(10.3節)，その構造に対して電子状態を解き，そこから状態密度・分子軌道・波動関数・振動準位を取り出す(10.4〜10.8節)．図10.1にその全体像を示す．以降の節は，この図の各ブロックを順に説明していくものだと思ってよい．

10.1 まずは Molcalc で触ってみる

いきなりPythonスクリプトを走らせる前に，ブラウザだけで量子化学計算が体験できるサービスに触れておこう．コペンハーゲン大学のJensenらが公開している Molcalc(https://molcalc.org/)である．分子を画面上で組み立て，ボタンを押すだけで構造最適化・熱力学量・振動数・分子軌道が返ってくる．

10.1.1 Molcalcの中身:GAMESSとSTO-3G

Molcalcは，内部で GAMESS(General Atomic and Molecular Electronic Structure System)という量子化学計算ソフトウェアを呼び出している．基底関数は STO-3G，すなわち第5章で導出した「Slater型軌道をGauss関数3個で近似した最小基底」である．1原子あたりの基底関数の数が最小限なので計算は速いが，そのぶん精度は粗い．原子数が多すぎる分子はサーバ側で計算不可になる．

画面の左側に Optimize(構造緩和)と Calculate Properties(電子状態計算・振動準位計算)のボタンがある．まず Optimize を押して構造を緩和し，次に Calculate Properties を押す，という順番が基本である．結果は4つのタブに分かれて返ってくる．

• Free Energies:エンタルピー・エントロピー・定圧熱容量の並進/回転/振動成分と，生成熱．
• Vibrational Frequencies:基準振動の波数($\mathrm{cm}^{-1}$)の一覧と，各モードのアニメーション．
• Molecular Orbitals:分子軌道のエネルギー(eV)一覧と，その等値面表示．
• Polarity and Solvation:双極子モーメントと溶媒和自由エネルギー．

たとえばメタン CH$_4$ について Calculate Properties を押すと，Free Energiesタブに表10.1のような数値が並ぶ．

振動のエンタルピーが 119.28 kJ mol$^{-1}$ と飛び抜けて大きいことに気づいてほしい．これは室温の熱エネルギーによる寄与ではなく，その大半が絶対零度でも残る 零点エネルギー(zero-point energy)である．10.8節で見るとおり，振動の基底状態のエネルギーは $\frac12\hbar\omega$ であってゼロではない．CH$_4$の9個の振動モードの $\frac12\hbar\omega$ を足し合わせると，確かに 100 kJ mol$^{-1}$ 台になる(演習10.5)．

10.1.2 Molcalcが返す分子軌道と，その基底依存性

Molecular Orbitalsタブには，CH$_4$について次の9個の軌道エネルギー(eV)が並ぶ．

$-300.14,\ -24.76,\ -14.13,\ -14.12,\ -14.12,\ +19.47,\ +19.48,\ +19.48,\ +20.55$

軌道が9個しかないのは，STO-3Gが最小基底だからである．C原子には $1s,\,2s,\,2p_x,\,2p_y,\,2p_z$ の5個，H原子には $1s$ が1個ずつで4個，合わせて9個の基底関数しかない．基底関数を $N$ 個用意すれば分子軌道もちょうど $N$ 個できる，という対応は第11章の永年方程式で明確になる．電子は10個(C:6個，H:4個)なので，下から5個の軌道が2個ずつ電子で埋まる．最低の $-300.14\ \mathrm{eV}$ はC-$1s$の内殻軌道であり，化学結合には関与しない．

Molcalcは手軽だが，(i) 基底がSTO-3Gに固定，(ii) 原子数の上限がある，(iii) 出力をファイルとして自分の解析に流し込めない，という制約がある．そこで以降は，第6章以来使ってきた PySCF に舞台を移し，Molcalcにある機能をすべて自分の手元で再現する．まずは計算する分子の構造を手に入れるところからである．

10.2 PubChem から分子構造を取ってくる

計算を始めるには，原子の並び(座標)が要る．水素原子やHe原子のように原子1個なら座標は原点だけで済んだが，分子ではそうはいかない．世の中の分子の構造は，幸いなことにデータベースとして整備されている．ここでは米国NIHが運営する PubChem(https://pubchem.ncbi.nlm.nih.gov/)を使う．

10.2.1 組成式で探すか，化学名で探すか

検索窓に組成式 C6H12O6 を入れると，候補が大量に出てくる．C$_6$H$_{12}$O$_6$ はグルコース(ブドウ糖)だけでなく，フルクトース，ガラクトース，マンノースなど多数の異性体が同じ組成式をもつからである．目的の物質が決まっているなら，化学名で検索するほうが速い．D-グルコースの環状構造(ピラノース形)なら D-Glucopyranose と入れればよい．すると次の情報をもつ候補が一意に見つかる．

物質を一意に指定する識別子としては CID(PubChem Compound ID)と InChIKey が便利である．論文やレポートに「グルコースを計算した」と書くだけでは，どの異性体・どの立体配置を使ったのかが伝わらない．CIDを併記しておけば誤解がない．

10.2.2 3D Conformer と SDF のダウンロード

目的の化合物のページを開くと，「1.2 3D Conformer」という項目がある．ここに三次元の分子模型が表示され，右上の Download Coordinates から座標ファイルを保存できる．SDFの欄の Save をクリックすると，Conformer3D_COMPOUND_CID_5793.sdf というファイルが手に入る．

SDF(structure-data file)は，原子の元素記号と座標に加えて，どの原子とどの原子が何重結合で結ばれているかという結合表(bond block)をもつ形式である．化学情報の世界では標準的だが，可視化ソフトによっては読めないことがあり，量子化学計算の入力としても冗長である．

10.2.3 ファイル名の作法

ダウンロードしたファイルは Conformer3D_COMPOUND_CID_5793.sdf という機械的な名前をしている．このまま作業ディレクトリに放り込むと，1週間後の自分が「5793って何だったか」と悩むことになる．そこで，SDFファイルは SDF_{物質名}.sdf，xyzファイルは XYZ_{物質名}.xyz という規則で名前を付け替えておく．こうしておけば，ls を打った瞬間に「どの物質のどんなファイルがあるか」が一目で分かる．

$ cp ~/Downloads/Conformer3D_COMPOUND_CID_5793.sdf .
$ cp Conformer3D_COMPOUND_CID_5793.sdf SDF_D-glucose.sdf
$ ls
calc_pyscf.py        SDF_D-glucose.sdf    XYZ_CH4.xyz
calc_pyscf_dos.py    sdf2xyz.py           XYZ_H.xyz
calc_pyscf_relax.py  XYZ_B.xyz            XYZ_He.xyz
calc_pyscf_vib.py    XYZ_Be.xyz           XYZ_N.xyz
calc_pyscf_wf.py     XYZ_C.xyz            XYZ_O.xyz

コード10.1 ダウンロードしたSDFファイルを，内容が分かる名前に付け替える．元のファイルも消さずに残しておくと，CIDを後から確認できる．

10.2.4 SDFからxyzへの変換と，xyz形式の書式

講義で配布されている sdf2xyz.py を使うと，SDFから結合表を落として座標だけを取り出したxyzファイルが得られる．

$ python3 sdf2xyz.py --sdf SDF_D-glucose.sdf
Converted: SDF_D-glucose.sdf -> XYZ_D-glucose.xyz

コード10.2 SDF から xyz への変換．出力ファイル名は入力名の SDF_ を XYZ_ に置き換えたものになる．

xyz形式は，量子化学計算でもっともよく使われる素朴なテキスト形式である．書式は次の3つの規則だけである．

3
H2O molecule  r(OH)=0.96 A, angle(HOH)=104.5 deg
O    0.000000    0.000000    0.000000
H    0.759062    0.587729    0.000000
H   -0.759062    0.587729    0.000000

コード10.3 水分子のxyzファイル(XYZ_H2O.xyz)．原子の並び順に決まりはないが，後で出力を読むとき混乱しないよう，重い原子から書く習慣をつけるとよい．

10.2.5 VESTAで眺める

座標ができたら，計算に投入する前に必ず目で見る．結晶構造可視化ソフト VESTA はxyz形式も読めるので，次のように開けばよい(macOSの例)．

$ open -a VESTA XYZ_D-glucose.xyz

コード10.4 xyzファイルをVESTAで開く．Windowsではファイルをダブルクリック，あるいはVESTAにドラッグ&ドロップすればよい．

VESTAでは分子をぐるぐる回して眺められるだけでなく，原子を2つクリックすれば結合長が，3つクリックすれば結合角が読める．D-グルコースなら，6員環(ピラノース環)が椅子形になっていること，C–O結合が約 1.4 Å，C–H結合が約 1.1 Å であることなどをその場で確認できる．画面下部のOutput欄には原子数と結合数(この構造では24原子・32結合)が表示される．

10.3 構造緩和を実行する

構造緩和(構造最適化，geometry optimization)とは，第8章で導入したポテンシャルエネルギー曲面

の極小点を探す作業である．$\Ham_{\mathrm{el}}$ は原子核を固定したときの電子系のハミルトニアン，第2項は核間反発である．Born–Oppenheimer近似(第8章)のおかげで，原子核の位置 $\{\bm{R}_I\}$ は「動かせるパラメータ」として扱えるのだった．極小点では，すべての原子に働く力

がゼロになる．すなわち構造緩和とは，式 \eqref{eq:10-pes} の曲面の上で式 \eqref{eq:10-force} の力がすべての原子についてゼロになる点を探す作業にほかならない．この力はHellmann–Feynmanの定理(第8章)によって波動関数から直接計算できるので，「力を計算する → 力の向きに原子を少し動かす → 電子状態を解き直す」という手続きを，力が十分小さくなるまで繰り返せばよい．

10.3.1 実行と，出力の読み方

D-グルコース(24原子)を 6-31G 基底，閉殻(スピン多重度1)で緩和してみよう．

$ python3 calc_pyscf_relax.py --xyz XYZ_D-glucose.xyz --spin 0 --basis 6-31g

コード10.5 構造緩和の実行．--spin はPySCFの流儀で $2S$(不対電子の数)を指定する．閉殻分子なら 0 である．

実行すると，内部で geomeTRIC という最適化ライブラリが呼ばれ，次のようなログが1ステップごとに流れる．

Step  20 : Displace = 3.194e-03/1.351e-02 (rms/max) Trust = 3.000e-01 (=) Grad = 8.707e-05/1.712e-04 (rms/max) E (change) = -683.0329609087 (-2.811e-06) Quality = 1.350
Step  21 : Displace = 1.732e-03/6.261e-03 (rms/max) Trust = 3.000e-01 (=) Grad = 6.012e-05/1.476e-04 (rms/max) E (change) = -683.0329621255 (-1.217e-06) Quality = 1.274
Step  22 : Displace = 6.254e-04/1.609e-03 (rms/max) Trust = 3.000e-01 (=) Grad = 2.309e-05/5.726e-05 (rms/max) E (change) = -683.0329624172 (-2.917e-07) Quality = 1.204
Converged!

========================================
 PySCF Geometry Optimization
========================================
XYZ file     : XYZ_D-glucose.xyz
Basis        : 6-31g
Method       : RHF
Theory       : scf
Charge       : 0
Spin (2S)    : 0
Text output  : relax_XYZ_D-glucose.txt
----------------------------------------
Optimization finished.
Final energy : -18586.273591 eV
Output XYZ   : XYZ_D-glucose-finish.xyz

コード10.6 構造緩和の出力(末尾)．24原子の分子では，構造緩和のためのSCF計算が23回(Step 0 から Step 22 まで)繰り返され，約960秒かかった．

読むべき数字は次の4つである．

• Grad:力(エネルギー勾配)の二乗平均と最大値．これが収束判定の主役である．上の例では最終的に最大 $5.7\times10^{-5}$ まで落ちている．
• E (change):全エネルギーと，前ステップからの変化量．単位はHartree．変化量が $10^{-6}$ Hartree(約 $3\times10^{-5}$ eV)を切れば十分収束している．
• Displace:原子がどれだけ動いたか．最終ステップでは最大 $1.6\times10^{-3}$ Bohr，すなわち $10^{-3}$ Å 程度しか動いていない．
• Converged!:上の3つがすべて閾値を下回った合図．この表示が出ていない結果を使ってはいけない．

最終的な全エネルギーは $-683.0329624172$ Hartree と表示されている．eVに直すには $1\ \mathrm{Hartree}=27.211386\ \mathrm{eV}$ を掛けて

となり，出力の Final energy : -18586.273591 eV と一致する．緩和後の構造は XYZ_D-glucose-finish.xyz というファイルに書き出される．元の XYZ_D-glucose.xyz は上書きされないので，緩和前後の比較ができる．

10.3.2 緩和前後の全エネルギーを比べる

構造緩和が意味のある仕事をしたかどうかは，緩和前と緩和後で全エネルギーを比べれば分かる．単純に，同じ手法・同じ基底で1点計算(single point calculation)を2回走らせればよい．

$ python3 calc_pyscf.py --xyz XYZ_D-glucose.xyz --spin 0 --basis 6-31g
========================================
 PySCF Calculation
========================================
XYZ file     : XYZ_D-glucose.xyz
Basis        : 6-31g
Charge       : 0
Spin (2S)    : 0
Theory       : scf
----------------------------------------
Method       : RHF
Total Energy : -18585.549008 eV
<S^2>        : 0.000000
Multiplicity : 1.0

$ python3 calc_pyscf.py --xyz XYZ_D-glucose-finish.xyz --spin 0 --basis 6-31g
...
Total Energy : -18586.273591 eV

コード10.7 緩和前(PubChemの構造)と緩和後の全エネルギーを 6-31G 基底で比較する．

差を取ると

である．確かにエネルギーが下がっている．これを化学でなじみのある単位に直すと，$1\ \mathrm{eV} = 96.485\ \mathrm{kJ\,mol^{-1}}$ より

となる．水素結合1本が 20 kJ mol$^{-1}$ 程度であることを思えば，決して無視できない大きさである．PubChemの力場構造が量子化学の平衡構造からそれなりにずれていた，ということである．

10.3.3 余談:基底関数を落とすと判定が逆転する

ここからが本節の山場である．まったく同じ2つの構造を，今度は STO-3G 基底で計算してみよう．

$ python3 calc_pyscf.py --xyz XYZ_D-glucose.xyz --spin 0
Basis        : sto-3g
Total Energy : -18352.694979 eV

$ python3 calc_pyscf.py --xyz XYZ_D-glucose-finish.xyz --spin 0
Basis        : sto-3g
Total Energy : -18352.345758 eV

コード10.8 同じ2つの構造を STO-3G 基底で計算する(--basis を省略すると既定の sto-3g になる)．

今度は

と，緩和後のほうがエネルギーが高いという結果になる．6-31G では 0.72 eV 安定化していたはずの構造が，STO-3G では 0.35 eV 不安定と判定されてしまった．符号が逆である．

10.4 状態密度(DOS)を描く

平衡構造が得られたら，次はその構造での電子状態を調べる．もっとも見通しのよい表示が 状態密度(density of states, DOS)である．

固体では軌道が連続的なエネルギーバンドをつくるので $D(E)$ は文字どおり連続関数になる(第14章)．一方，分子の軌道エネルギーは飛び飛びなので，式 \eqref{eq:10-dos-delta} のままでは無限に細い針が何本も立っているだけの絵になってしまい，見づらい．そこでデルタ関数を幅 $\sigma$ のGauss関数で置き換える．

これを Gauss拡がり(Gaussian broadening)という．$\sigma\to 0$ の極限で式 \eqref{eq:10-dos-delta} に戻る．Gauss関数の全積分は1なので，$\sigma$ をどう選んでもピークの面積は変わらない．変わるのは高さと幅だけである．calc_pyscf_dos.py の既定値は $\sigma=0.300$ eV である．

10.4.1 部分状態密度(PDOS)

DOSの本当の威力は，「どの原子の，どの軌道が，どのエネルギーに寄与しているか」を分解できる点にある．分子軌道はLCAO(基底関数の線形結合)で

と書ける($\chi_\mu$ は原子に置かれた基底関数，第5章のSTO-3Gなど)．式 \eqref{eq:10-lcao} の係数 $c_{\mu i}$ は，第11章で導く永年方程式を解いて決まる量である．この $i$ 番目の分子軌道のうち，基底関数 $\mu$ がどれだけの割合を占めているかを表す重み $w_{i\mu}$ を決められれば，

として，原子(あるいは元素，あるいは $s$/$p$/$d$ 軌道)ごとの 部分状態密度(partial DOS, PDOS)が定義できる．ここで $\mu\in\alpha$ は「グループ $\alpha$ に属する基底関数についての和」を意味する．

重みの決め方には流儀がある．素朴には $|c_{\mu i}|^2$ を使いたくなるが，基底関数どうしが直交していない($S_{\mu\nu}=\int\chi_\mu^*\chi_\nu\,\dd^3r \neq \delta_{\mu\nu}$，第5章の重なり積分)ため，$\sum_\mu |c_{\mu i}|^2 \neq 1$ となってしまう．そこで重なり行列 $S$ の平方根を使って基底を直交化してから係数を取る Löwdin射影

が広く使われる($C$ は係数 $c_{\mu i}$ を並べた行列)．calc_pyscf_dos.py の出力に Projection : lowdin と書かれているのはこの意味である．$\sum_\mu w_{i\mu}=1$ が成り立つので，すべてのPDOSを足すと必ず全DOSに戻る．

10.4.2 グルコースのDOSを描く

$ python3 calc_pyscf_dos.py --xyz XYZ_D-glucose-finish.xyz --spin 0 --basis 6-31g \
      --xrange -30 10 --element-pdos
========================================
 PySCF DOS-like Calculation
========================================
XYZ file     : XYZ_D-glucose-finish.xyz
Method       : RHF
Theory       : scf
Basis        : 6-31g
Total Energy : -18586.273591 eV
Charge       : 0
Spin (2S)    : 0
Projection   : lowdin
DOS mode     : gaussian broadening
Gaussian sigma: 0.300 eV
Energy zero  : absolute MO energy
Plot output  : DOS_XYZ_D-glucose-finish.pdf
CSV output   : DOS_XYZ_D-glucose-finish.csv
Text output  : DOS_XYZ_D-glucose-finish.txt
Element PDOS : C(s), C(p), H(s), O(s), O(p)

コード10.9 元素別・軌道別の部分状態密度を計算する．--xrange -30 10 は横軸の範囲(eV)，--element-pdos は元素ごとのPDOSを出すオプション．オプションを忘れたら python3 calc_pyscf_dos.py --help と打てばよい．

出てきたPDOSを見ると，$-30$ eV から $-11$ eV あたりまでが占有された価電子帯にあたり，そのうち もっとも高い占有軌道(HOMO)付近は O-$2p$ と C-$2p$ で構成されていることが分かる．一方 $+4$ eV より上の空軌道側は H-$s$ の寄与が大きい．グルコースは C, H, O の3元素からなるが，化学的に反応しやすい「一番外側の電子」がどの原子のどの軌道に住んでいるかが，この一枚の図から読み取れるわけである．

10.4.3 CH$_4$のDOSと，そこから読めること

グルコースは24原子もあって準位が混み合っているので，もっと単純なCH$_4$で練習しよう．

$ python3 calc_pyscf_dos.py --xyz XYZ_CH4.xyz --spin 0 --basis 6-31g --xrange -30 10
 PySCF DOS-like Calculation
----------------------------------------
Total Energy : -1093.361059 eV
Gaussian sigma: 0.300 eV

MO projected weights:
spin        MO   energy(eV)     occ        s        p        d        f
-----------------------------------------------------------------------
restricted   1   -305.00055       2   1.0000   0.0000   0.0000   0.0000
restricted   2    -25.71378       2   1.0000   0.0000   0.0000   0.0000
restricted   3    -14.76147       2   0.4275   0.5725   0.0000   0.0000
restricted   4    -14.76090       2   0.4275   0.5725   0.0000   0.0000
restricted   5    -14.76075       2   0.4275   0.5725   0.0000   0.0000
restricted   6      6.91612       0   1.0000   0.0000   0.0000   0.0000
restricted   7      8.77636       0   0.6643   0.3357   0.0000   0.0000
restricted   8      8.77645       0   0.6643   0.3357   0.0000   0.0000
restricted   9      8.77674       0   0.6644   0.3356   0.0000   0.0000
restricted  10     20.26052       0   0.2080   0.7920   0.0000   0.0000
restricted  11     20.26068       0   0.2080   0.7920   0.0000   0.0000
restricted  12     20.26099       0   0.2080   0.7920   0.0000   0.0000

コード10.10 CH$_4$のDOS計算．各分子軌道のエネルギー，占有数(occ)，$s$/$p$/$d$/$f$ 成分の重み(式 \eqref{eq:10-lowdin} のLöwdin重み)が表になって出てくる．

この表とDOSの図(図10.3)を見比べると，次のことが読み取れる．

• $-305.00$ eV の軌道はC-$1s$の内殻軌道である．図の範囲($-30$ eV 以上)には現れない．
• $-25.71$ eV に非縮退の軌道が1本，占有数2．$s$ 成分100%である．
• $-14.76$ eV に3本の軌道がほぼ完全に縮退しており，占有数はいずれも2，合わせて電子6個．$s$ が 42.75%，$p$ が 57.25% である．これがHOMOである．
• $+6.92$ eV に占有数0の軌道が1本($s$ 成分100%)．これがLUMOである．
• $+8.78$ eV に3重縮退の空軌道．

10.5 CH$_4$ の分子軌道

コード10.10の表は，実はCH$_4$の化学結合の全体像をそのまま含んでいる．読み解いていこう．

10.5.1 電子はどこに入っているか

CH$_4$の電子数は $6+4\times1=10$ 個である．表の占有数(occ)を上から足すと $2+2+2+2+2=10$ となり，5本の軌道が埋まっている．内訳は

である．$-14.76$ eV の軌道に電子が6個入っている——これが本節でいちばん覚えてほしい事実である．同じエネルギーの軌道が3本あり，それぞれに2個ずつ電子が入るからである．

10.5.2 $s$ 成分と $p$ 成分が語ること

表10.3の $s$ と $p$ の欄は，単なる数字の羅列ではない．MO 2($-25.71$ eV)は $p$ 成分がぴったり 0 である．つまりこの軌道はC-$2p$をまったく含まず，C-$2s$と4個のH-$1s$だけからできている．一方 MO 3–5($-14.76$ eV)は $s$ が 0.4275，$p$ が 0.5725 で，C-$2p$ が半分以上を占める．しかもこちらには C-$2s$ が混ざっていない．

10.5.3 分子軌道ダイアグラム

以上をまとめると，図10.4のような 分子軌道ダイアグラム(molecular orbital diagram)が描ける．左に4個の水素がつくる群軌道，右に炭素の原子軌道，中央に両者が混ざってできた分子軌道を置く．

図の読み方を確認しておこう．

1. まず4個の水素の $1s$ 軌道を，CH$_4$の対称性に合った4通りの組み合わせ(群軌道，symmetry-adapted linear combination)に組み替える．全部同符号のものが1本($a_1$)，符号が入れ替わるものが3本($t_2$)できる．
2. 炭素側の $2s$ は $a_1$ 型，$2p_x,2p_y,2p_z$ は $t_2$ 型である．
3. 同じ型どうししか混ざらない．$a_1$ と $a_1$ から結合性 $a_1$ と反結合性 $a_1^{*}$ が，$t_2$ と $t_2$ から結合性 $t_2$(3本)と反結合性 $t_2^{*}$(3本)ができる．
4. 価電子8個(炭素の $2s^2 2p^2$ で4個，水素4個で4個)を下から詰めると，$a_1$ に2個，$t_2$ に6個入り，ちょうど結合性軌道だけが満たされる．反結合性は空である．だからCH$_4$は安定な分子で，結合次数は形式的に4本ぶんになる．

10.6 光電子スペクトルとの比較

計算した軌道エネルギーは，実験で確かめられるのだろうか．確かめられる．光電子分光(photoelectron spectroscopy)がそれである．

10.6.1 光電子分光とは

物質に，エネルギー $h\nu$ の光(紫外線ならUPS，X線ならXPS)を当てると，電子が飛び出す．飛び出した電子の運動エネルギー $E_{\mathrm{kin}}$ を測れば，エネルギー保存則から，その電子がもともと束縛されていたエネルギー(イオン化エネルギー $I$)が分かる．

式 \eqref{eq:10-pes-energy} の $h\nu$ は入射光のエネルギーで既知，$E_{\mathrm{kin}}$ は測定量だから，$I$ が決まる．横軸を $I$，縦軸を「飛び出してきた電子の数」にとってプロットしたものが光電子スペクトルである．ある $I$ に電子がたくさん詰まっていれば，そこから飛び出す電子も多い．つまり光電子スペクトルは占有状態の状態密度を実験的に測っているとみなせる．

10.6.2 Koopmansの定理

計算で得られた軌道エネルギー $\epsilon_i$ と，実験のイオン化エネルギー $I$ を結ぶのが Koopmansの定理(Koopmans' theorem)である．

10.6.3 CH$_4$のスペクトルと計算DOSの比較

CH$_4$の紫外光電子スペクトルには，$I\approx 14$ eV あたりに幅の広い大きなバンドと，$I\approx 23$ eV に小さなバンドの，合わせて2本が現れる．これを計算DOSと並べたのが図10.5である．

比較の結論を整理しよう．

• ピークの本数と面積比:計算は $t_2$(6状態)と $a_1$(2状態)の2本を予言し，面積比は $3:1$ である．実験でも $T_2$ バンドの強度は $A_1$ バンドのおよそ3倍であり，一致する．これがCH$_4$の価電子が「1種類の等価な結合」ではないことの，直接的な実験的証拠である．
• ピークの位置:計算(Koopmans)は $14.76$ eV と $25.71$ eV，実験は約 $14$ eV と約 $23$ eV である．$t_2$ の誤差は 1 eV 弱，$a_1$ は 3 eV 程度で，内殻に近い深い準位ほど軌道緩和の効果が大きく，Koopmansのずれも大きい．
• ピークの幅:計算では $\sigma=0.3$ eV の人為的な幅しかないが，実験の $T_2$ バンドは 4 eV 近い幅をもつ．これは(i) イオン化に伴って分子が振動励起される(電子状態と振動状態の結合)，(ii) CH$_4^{+}$ が3重縮退した $T_2$ 状態になるためJahn–Teller効果で構造が歪む，という物理による．計算の細いピークと実験の太いバンドを，幅まで比べてはいけない．

10.7 波動関数を可視化する

軌道エネルギーと成分比だけでは，「その軌道がどんな形をしているか」は分からない．分子軌道 $\phi_i(\bm{r})$ は空間の各点で値をもつ関数なので，三次元の地図として描き出せる．

10.7.1 cubeファイルを作る

$ python3 calc_pyscf_wf.py --xyz XYZ_CH4-finish.xyz --spin 0 --basis 6-31g --homo --lumo
========================================
 PySCF Wavefunction Cube Generator
========================================
XYZ file     : XYZ_CH4-finish.xyz
Method       : RHF
Basis        : 6-31g
Charge       : 0
Spin (2S)    : 0
Grid mode    : nx, ny, nz = 80, 80, 80
Margin       : 3.000 Bohr
Output dir   : WF_XYZ_CH4-finish
----------------------------------------
Generated cube files:
spin        MO   energy(Ha)   energy(eV)   occ   file
restricted   5    -0.545033     -14.8311     2   WF_XYZ_CH4-finish/MO_005_HOMO.cube
restricted   6     0.256250       6.9729     0   WF_XYZ_CH4-finish/MO_006_LUMO.cube

コード10.11 HOMOとLUMOのcubeファイルを生成する．WF_ で始まるディレクトリが作られ，その中にcubeファイルが出力される．

cubeファイル(Gaussian cube format)は，直方体の格子点上での関数値を並べたテキストファイルである．ここでは分子を囲む箱を $80\times80\times80$ の格子に切り(全部で 512000 点)，各点での $\phi_i(\bm{r})$ の値を書き出している．Margin : 3.000 Bohr は，分子の外側に 3 Bohr($\approx1.6\ \text{Å}$)ぶんの余白をとる，という意味である．軌道は原子核から離れても指数関数的にしか減衰しないので，余白が狭すぎると軌道の裾が箱からはみ出して絵が切れてしまう．

他の軌道も見たければ，番号を直接指定すればよい．

$ python3 calc_pyscf_wf.py --xyz XYZ_CH4-finish.xyz --spin 0 --basis 6-31g --mo 1 2 3 4 7
$ cd WF_XYZ_CH4-finish
$ ls
MO_001.cube  MO_002.cube  MO_003.cube  MO_004.cube  MO_005_HOMO.cube
MO_006_LUMO.cube  MO_007.cube
$ open -a VESTA MO_005_HOMO.cube
$ open -a VESTA MO_006_LUMO.cube

コード10.12 番号を指定して任意の分子軌道を出力し，VESTAで開く．第9章で議論した内殻軌道(MO 1)がいかに原子核に張り付いているかも，この方法で確認できる．

10.7.2 等値面の読み方:黄は正，青は負

VESTAでcubeファイルを開くと，等値面(isosurface)が表示される．$\phi_i(\bm{r})=\pm c$($c$ は適当な定数，たとえば $0.02$)を満たす面を描いたもので，雲のような塊として見える．VESTAの既定の色分けは

黄色 = 波動関数が正の領域，青色 = 波動関数が負の領域

である．図10.6にCH$_4$のHOMOとLUMOの模式図を示す．

形から結合性・反結合性を読むコツは，「C–H結合の途中に節面(符号が変わる面)があるか」を見ることである．

• 結合性(bonding):結合の途中に節面がない．2つの原子の間で波動関数が同符号に重なるので，そこに電子密度が積み増しされ，原子核どうしを引きつける糊の役割を果たす．CH$_4$のHOMOがこれである．
• 反結合性(antibonding):結合の途中に節面がある．電子密度が結合の中央で減るので，電子が入ると分子は不安定になる．CH$_4$のLUMOがこれである．

10.7.3 グルコースのHOMOとLUMO

同じことをD-グルコース(24原子)でもやってみよう．

$ python3 calc_pyscf_wf.py --xyz XYZ_D-glucose-finish.xyz --spin 0 --basis 6-31g \
      --homo --lumo

コード10.13 グルコースのHOMO/LUMOをcube出力する．原子数が多いので格子点数も増え，ファイルサイズは1つあたり数十MBになる．

得られた等値面を読むと，次のことが分かる．

• HOMO:C-$2p$ と C-$2p$ が同符号で並ぶ $\sigma$ 結合性の部分と，C-$2p$ と O-$2p$ が逆符号で向かい合う $\pi^{*}$ 反結合性の部分が共存している．10.4節のPDOSで「HOMOはO-$2p$とC-$2p$で構成されている」と読んだこととも整合する．反結合性成分をもつ軌道がいちばん上にあるという事実は，この分子のどこが化学的に弱いか(酸化されやすいか)を示唆する．
• LUMO:波動関数が分子全体に大きく広がってしまい，目視での判断は難しい．PDOSではこの領域はH-$s$が支配的なので，HとC，HとOの間に何らかの反結合性をもつ軌道だと推測できる．

10.8 振動準位

ここまでは電子の話だった．最後に原子核の運動，すなわち分子振動を扱う．Born–Oppenheimer近似(第8章)では，原子核はポテンシャルエネルギー曲面 $E(\{\bm{R}_I\})$ の上を運動する．平衡構造はその極小点だから，その近くでは曲面を放物線で近似できる．

10.8.1 調和近似から振動数へ

簡単のため，まず1次元(2原子分子の伸縮)で考える．平衡結合長を $r_0$，そこからのずれを $x=r-r_0$ とし，ポテンシャル $U(x)$ を $x=0$ のまわりでTaylor展開する．

平衡点では力がゼロ，すなわち1次の項の係数 $\dd U/\dd x|_0$ が消える．エネルギーの原点を $U(0)=0$ に選び，3次以上を捨てると

を得る．これが 調和近似(harmonic approximation)であり，$k$ を 力の定数(force constant)という．単位は N m$^{-1}$ である．実際の計算では，平衡構造のまわりで原子を少しずつ動かしてエネルギーを何点か求め，$\frac12kx^2$ にフィッティングして $k$ を決める(あるいは2階微分を解析的に評価する)．

質量 $m_A,m_B$ の2原子が，ばね定数 $k$ のばねでつながれているときの振動は，換算質量

をもつ1個の粒子の運動に帰着する(相対座標に移ると重心運動が分離するため)．したがって古典的な角振動数は

である．量子力学では，この調和振動子のエネルギー準位は等間隔に量子化され，

となる．基底状態 $n=0$ でもエネルギーが $\frac12\hbar\omega$ 残る．これが10.1節で出てきた零点エネルギーである．

10.8.2 波数への換算

分光の世界では，振動数は 波数(wavenumber)$\tilde\nu$ で表す習慣である．定義は「$1$ cm あたりに波が何個入るか」で，

である($c$ は光速)．エネルギーとの関係は $E=h\nu=hc\tilde\nu$ だから，$1$ eV に対応する波数は

すなわち

である．$3000\ \mathrm{cm^{-1}}$ の振動なら $3000/8065.54=0.372$ eV，室温の熱エネルギー $k_{\mathrm{B}}T\approx0.026$ eV($\approx 209\ \mathrm{cm^{-1}}$)の14倍もある．分子振動の多くは室温ではほとんど励起されていない(基底状態 $n=0$ にいる)ことがここから分かる．

10.8.3 CH$_4$の振動計算

$ python3 calc_pyscf_vib.py --xyz XYZ_CH4-finish.xyz --basis 6-31g --spin 0
========================================
 PySCF Vibrational Level Analysis
========================================
XYZ file     : XYZ_CH4-finish.xyz
Theory       : scf
Method       : rhf
Basis        : 6-31g
Charge       : 0
Spin (2S)    : 0
SCF E_tot    : -40.1805541767 Hartree (-1093.368579 eV)
----------------------------------------
Mode  Frequency(cm^-1)  Type
   1         1516.7889  vibrational
   2         1516.7896  vibrational
   3         1516.7911  vibrational
   4         1708.6872  vibrational
   5         1708.6873  vibrational
   6         3181.4044  vibrational
   7         3295.8557  vibrational
   8         3295.8593  vibrational
   9         3295.8661  vibrational

コード10.14 CH$_4$の振動準位計算．$3N-6=9$ 個の振動モードが，縮退度 3, 2, 1, 3 のグループに分かれている．

実験値と比べてみよう．

計算値はすべて実験値より高い．しかも「実験/計算」の比が $0.86$〜$0.92$ とだいたい揃っている点に注目してほしい．

10.8.4 H$_2$Oの振動モード

3原子の水分子なら，振動は $3\times3-6=3$ 個だけである．3つのモードの動き方を図10.8に示す．

$ python3 calc_pyscf_vib.py --xyz XYZ_H2O-finish.xyz --basis 6-31g --spin 0
SCF E_tot    : -75.9839219956 Hartree (-2067.627850 eV)
----------------------------------------
Mode  Frequency(cm^-1)  Type
   1         1803.9668  vibrational
   2         4065.0470  vibrational
   3         4164.2241  vibrational
----------------------------------------
Harmonic vibrational levels up to n=3
(E_n = (n + 1/2) h nu, for each normal mode)

Mode 1  nu = 1803.9668 cm^-1        Mode 2  nu = 4065.0470 cm^-1
  n     E_n (eV)                      n     E_n (eV)
  0      0.111832                      0      0.252001
  1      0.335495                      1      0.756002
  2      0.559158                      2      1.260004
  3      0.782822                      3      1.764006

コード10.15 H$_2$Oの振動準位．各モードについて式 \eqref{eq:10-levels} の $E_n$ が表示される．$n=0$ の値がそのモードの零点エネルギーである．

実験値(気相)は 1595，3657，3756 $\mathrm{cm}^{-1}$ であり，FT-IRで水蒸気を測ると確かに $1600\ \mathrm{cm^{-1}}$ 付近と $3600$〜$3800\ \mathrm{cm^{-1}}$ に吸収が現れる．計算はこの「2本立て」の構造をよく再現している．

10.8.5 計算値はなぜ実験値より高いのか:スケーリング因子

表10.4でも水でも，Hartree–Fock計算の振動数は実験値より一様に 10% 前後高かった．原因は3つある．

1. 調和近似そのもの:式 \eqref{eq:10-taylor} で捨てた3次以降の項(非調和項)が効く．実際のポテンシャルは図10.7の実線のように，伸ばす方向では放物線より緩やか(結合が切れて平坦になる)である．したがって真の $n=0\to1$ 遷移エネルギーは調和近似の $\hbar\omega$ より小さい．実験の「基本振動数」は非調和性を含んだ値なので，必然的に低く出る．
2. 基底関数の不足:柔軟性の足りない基底では，結合を伸ばしたときの電子分布の変化を十分に記述できず，ポテンシャルが硬めに出る．
3. 電子相関の欠如:Hartree–Fockは電子相関を含まないので(第9章)，結合が伸びたときのエネルギー上昇を過大評価する．これも $k$ を大きくする方向に効く．

3つとも同じ向き(振動数を高くする方向)に効き，しかもその割合が分子やモードによらずおおむね一定なので，実務では計算値に定数を掛けて補正する．この定数を スケーリング因子(scaling factor)という．

10.9 計算をするときの心構え

最後に，本章で見た具体例から一般化できる「作法」をまとめておく．計算そのものは機械にやらせるべきだが，その結果を信用してよいかを判断するのは人間の仕事である．

10.9.1 基底関数依存性を必ず確かめる

10.3節のグルコースの一件が示すとおり，基底関数の選び方は結論の符号すら変える．実務上の原則は次の3つである．

• 絶対値は比較しない．基底が違えば全エネルギーは何十 eV も変わる(グルコースで 233 eV，CH$_4$ でも STO-3G と 6-31G で 12.3 eV)．意味があるのは，同じ基底・同じ手法で計算した2つの状態の差だけである．
• 収束テストをする．STO-3G → 6-31G → 6-31G(d) → cc-pVTZ と基底を大きくしていき，興味のある量(結合長，エネルギー差，振動数)が変化しなくなることを確認する．変化が止まらないうちは，その数値を結論に使ってはいけない．
• 空軌道は特に危ない．10.1節で見たように，最小基底ではLUMOの縮退度さえ変わる．HOMO–LUMOギャップや励起エネルギーを論じるなら，十分大きな基底が要る．

10.9.2 収束条件を意識する

計算には二重のループがあった(10.3節)．それぞれに収束条件がある．

• SCFの収束:エネルギー変化や密度行列の変化が閾値(既定では $10^{-9}$ Hartree 程度)を下回るまで繰り返す．収束しない場合は，初期推定を変える，収束加速法(DIIS)の設定を変える，対称性の高すぎる初期構造を崩す，といった対処をする．「収束しませんでした」という出力を見落として結果を使うのが最悪の事故である．
• 構造緩和の収束:力の最大値と，エネルギー変化，原子の変位がすべて閾値を下回ること．Converged! の表示を必ず確認する．
• 比較する2つの計算では，収束条件も揃える．片方だけ緩い条件で止めた計算と比べると，0.1 eV 程度の差は簡単に生まれる．

10.9.3 電荷とスピン多重度は自分で指定する

計算機は「この分子は中性で閉殻だろう」と勝手に判断してくれるが，それが正しいとは限らない．第7章で学んだとおり，スピン多重度は $2S+1$ である．PySCF系のスクリプトでは --spin に $2S$(不対電子の数)を，--charge に電荷を指定する．

指定を誤ると，たとえばO$_2$を一重項として計算してしまい，実験と合わない結合長・振動数が出る．O$_2$の基底状態が三重項であることはHundの規則(第9章)から予想できる．計算を始める前に，その系の電子配置を紙の上で考えておくのが正しい順番である．

10.9.4 全エネルギーの内訳を見る

エネルギーがおかしいと感じたら，内訳を見るとよい．第9章で分解した各項が出力できる．

$ python3 calc_pyscf.py --xyz XYZ_CH4-finish.xyz --spin 0 --basis 6-31g \
      --decompose-total-energy
Method       : RHF
Total Energy : -1093.366788 eV
--- decomposition (eV) ---
Kinetic          T    :  +1094.611539
Nuclear attr.    V_ne :  -3264.363352
Coulomb (Hartree) J   :   +888.540843
Exchange         K    :   -178.760384
Nuclear repulsion     :   +366.604566
--------------------------
Sum                   :  -1093.366788

コード10.16 全エネルギーの内訳．第9章で導いた $E=T+V_{ne}+J+K+E_{\text{核間}}$ の各項が数値として現れる．

ここで交換項 $K=-178.76$ eV が負であること，すなわち交換相互作用が系を安定化していることを確認してほしい．第9章で「同じスピンの電子どうしは互いを避けるので，Coulomb反発の数え過ぎが補正される」と述べた効果が，そのまま数値として現れている．またこの分解は検算にも使える．$T$ と全ポテンシャル $V=V_{ne}+J+K+E_{\text{核間}}$ の比は

となり，Coulomb系の平衡構造で成り立つビリアル定理 $2T=-V$ をほぼ満たしている．この比が2から大きくずれていたら，構造緩和が終わっていないか，計算のどこかが壊れていると判断できる．

10.9.5 実験値と比べるときの作法

本章では2回，計算値と実験値を比べた(光電子スペクトルと振動数)．どちらの場合も，比べる前に問うべき質問は同じである．

1. 同じ量を比べているか．垂直イオン化エネルギーか断熱イオン化エネルギーか．調和振動数か非調和を含む基本振動数か．
2. 系の状態は同じか．計算は絶対零度・気相・孤立分子である．実験が水溶液中や結晶中なら，周囲の分子との相互作用(第13章のvan der Waals力や水素結合)が数百 $\mathrm{cm^{-1}}$ 単位で振動数を動かす．
3. 手法固有の系統誤差はどちらを向いているか．Hartree–Fockの振動数は高めに出る，Koopmansのイオン化エネルギーは深い準位ほど過大になる，といった系統誤差の向きを知っていれば，ずれを見て慌てずに済む．
4. 絶対値が合わなくても，傾向は合うことが多い．置換基を変えたときの振動数の増減，分子を変えたときのHOMOの上下といった相対的な傾向は，絶対値よりずっと信頼できる．マテリアル探索で計算が力を発揮するのはこの点である(文献[9])．

10.10 まとめと演習

10.10.1 この章のまとめ

• Molcalc(GAMESSベース，STO-3G固定)はブラウザだけで構造最適化・熱力学量・振動数・分子軌道を返す．出力の型を覚えるのに便利だが，基底が固定で原子数の上限もあるため，本格的な計算はPySCFで行う．
• 構造の入手:PubChemで化学名(組成式では候補が多すぎることがある)を検索し，3D ConformerからSDFを保存する．SDF_{物質名}.sdf と名前を付け替え，sdf2xyz.py でxyz形式に変換し，VESTAで目視確認する．xyz形式は「1行目=原子数，2行目=コメント，3行目以降=元素記号と座標(Å)」．
• 構造緩和:力 $\bm{F}_I=-\partial E/\partial\bm{R}_I$ がゼロになる点を探す(第8章)．D-グルコース(24原子)ではSCFが23回繰り返されて収束した．6-31G基底では緩和により $\Delta E=-0.725$ eV(式 \eqref{eq:10-dE631g})と安定化するが，同じ2構造をSTO-3Gで計算すると $\Delta E=+0.349$ eV(式 \eqref{eq:10-dEsto3g})と符号が逆転する．荒い基底では構造の良し悪しの判定そのものが誤りうる．
• 状態密度:$D(E)=\sum_i g_i\delta(E-\epsilon_i)$(式 \eqref{eq:10-dos-delta})．分子では離散準位をGauss関数で幅 $\sigma$ にぼかして描く(式 \eqref{eq:10-dos-gauss})．Löwdin射影(式 \eqref{eq:10-lowdin})により元素別・軌道別のPDOS(式 \eqref{eq:10-pdos})が得られる．グルコースのHOMOは主にO-$2p$とC-$2p$からなる．
• CH$_4$の分子軌道:$-25.71$ eV に $a_1$(電子2個，$s$ 成分100%)，$-14.76$ eV に3重縮退した $t_2$(電子6個，HOMO，C-$2p$とH-$1s$の結合性)，$+6.92$ eV に $a_1^{*}$(占有数0，LUMO)．$a_1$ と $t_2$ が混ざらないのは対称性による(第12章)．
• 光電子スペクトル:Koopmansの定理 $I_i=-\epsilon_i$(式 \eqref{eq:10-koopmans})により，軌道エネルギーがイオン化エネルギーと対応する．CH$_4$では $T_2$ バンドと $A_1$ バンドの面積比が 3:1 になり，計算と実験が一致する．横軸はイオン化エネルギーである．
• 波動関数の可視化:cubeファイルをVESTAで開き，等値面を見る(黄=正，青=負)．結合の途中に節面がなければ結合性，あれば反結合性．グルコースのHOMOはC-$2p$どうしの $\sigma$ 結合性とC-$2p$/O-$2p$の $\pi^{*}$ 反結合性が共存する．厳密な解析には重なり密度やCOHP(第1章)が要る．
• 振動準位:調和近似 $U=\frac12kx^2$(式 \eqref{eq:10-harmonic})，$\omega=\sqrt{k/\mu}$(式 \eqref{eq:10-omega})，$E_n=(n+\frac12)\hbar\omega$(式 \eqref{eq:10-levels})．$1\ \mathrm{eV}=8065.54\ \mathrm{cm^{-1}}$(式 \eqref{eq:10-ev-cm})．H$_2$Oの実験値 1595, 3657, 3756 $\mathrm{cm^{-1}}$ に対し計算値は 1804, 4065, 4164 $\mathrm{cm^{-1}}$ と約10%高く，スケーリング因子 $0.895$ で補正できる．虚振動が出たら平衡構造ではない．

10.10.2 演習問題

10.10.3 参考文献