エネルギーの計算（プロパティの計算含む）を一点計算（single-point calculation）と呼びます。 「一点の構造での計算」という意味です。 そうでない計算は、例えば構造最適化計算や分子動力学シミュレーションのように、分子構造を変化させていく計算です。

前回の通りPySCFのインストールが完了し、さらに次の入力ファイルでの計算が正しくできたことを前提にしています。

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が得られたかと思います。 これが、H₂分子の「エネルギー」を意味しています。 単位は書かれていませんが、その場合は「原子単位 (atomic unit)」あるいは「ハートリー (hartree)」であるという暗黙の了解があります。 「エネルギー」に関しては次回一点計算（出力ファイル）のところでもう少し説明する予定です。

この入力ファイル・出力ファイルだけ見たところで、事前知識がないと分かりづらいですが、 RHF（restricted Hartree--Fock）という手法とcc-pVDZという基底関数を組み合わせて計算を行いました。 このため、ここでの計算を発表する際には、例えば「(R)HF/cc-pVDZ」の計算を行ったと表記します。 (R)が必要かは場合によるかと思います。 なので、上記エネルギーは「H₂の(R)HF/cc-pVDZにおけるエネルギーは-1.06111199785749 hartree」と表現したりします。 ただし、実際に何のエネルギーを用いて議論するかや、どの程度の桁数が必要かは、議論する性質や単位によって変わります。 Hartree--Fock法についての詳細はSzabo: Chapter 3等をご参照ください。 現代ではHartree--Fock法の結果をそのまま論文等で報告することは稀ですが、 電子相関を考える際のベースとなる手法です。 現代では電子相関理論よりも密度汎関数理論（density functional theory）を用いた計算の方が主流ですが。

ところで、

[nisimoto@yn1 test]$ python3
Python 3.6.8 (default, Aug 10 2019, 06:54:07) 
[GCC 4.8.5 20150623 (Red Hat 4.8.5-36)] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> 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
-1.0611119978574877
>>> quit()

というように、Python内でコマンドを打ち込んでいっても良いのですが、例えばsp.pyというファイルに上記入力ファイルを書き込んでおき、pythonに流し込むという方法でも計算できます。

[nisimoto@yn1 test]$ cat sp.py 
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()
[nisimoto@yn1 test]$ python3 sp.py 
converged SCF energy = -1.06111199785749

どちらでも結果は同じになります。 使いやすい方法を使ってください。 ただし、LinuxやWSL上だと後者はviかEmacsでファイル編集しないといけないです。 私は基本的に後者の方法を用いています。

最初なので、一行ずつ説明してきましょう。

ここでは、"pyscf"というライブラリ（と呼んで良いか分かりませんが）から、"gto"と"scf"というモジュール（パッケージ？）をインポートしています。 pyscfライブラリの中に、gtoやscfというモジュールが入っていて、続く行の計算でそれらのモジュールを使えるようにしているというわけです。 C++で言う#includeと同じような感じですね。 別の計算、例えばMP2の計算であれば、さらにimport mpをしなければなりません。 なので、きちんとインポートしているか、確認する必要があります。 忘れていると、NameError: name 'mp' is not definedというようにエラーメッセージが出たりします。 「~ is not defined」とか言われたら、タイプミスをしていないか確認し、さらに必要なモジュールをインポートしているかを確認しましょう。

ここでは計算する分子に含まれる原子の種類と座標（atom）、そして基底関数（basis）を設定しています。 まず左辺のmolは任意です。 ここでは変数としての"mol"なので、自由に変更することができます。 例えば二行目・三行目を次のようにしても動くはずです

my_mol = gto.M(atom='H 0 0 0; H 0 0 1.2', basis='ccpvdz')
mf = scf.RHF(my_mol)

ただし、右辺の関数gto.Mは変更できません。 しかも大文字・小文字が区別されるので、この通りに打ち込む必要があります（gto.mはダメ）。 atomやbasisも同様に、この通りに打ち込む必要があります。 ただし、シングルクォーテーションマークはダブルクォーテーションマークにしても大丈夫のようです。 このgtoに関する更なる詳細は、マニュアルをご参照ください。 ここでご紹介する書き方以外にも、いろいろ対応しています。

atomで原子の種類と座標を設定します。 ここではH₂です。 一つ目のHを(x,y,z) = (0,0,0)という位置に置き、二つ目のHを(x,y,z) = (0,0,1.2)という位置に置くことを意味しています。 というわけで、二つの水素原子は、z軸方向に1.2だけ離れているという分子を設定したわけです。 それぞれの原子は、セミコロンで区切るか、改行する必要があります（ただし、Python環境でのみ。pythonの引数とする場合は、別の記法を用いる必要があります。）。 また、"H"ではなく、原子番号である"1"（ただし整数）を用いて指定するのもOKです。 このような三次元の直交座標のことを、カーテシアン座標（Cartesian coordinate）と読んでいます。 カーテシアン座標以外にもZ-matrixがありますが、私は自分で作成したことがないのでよく分かりません。

なお、"1.2"の単位はオングストローム（Å; 1 Å = 10^-10 m）です。 経験上多くのプログラムはオングストローム単位を用いていますが、 プログラムによっては、デフォルトで原子単位（atomic unit = Bohr）を用いることがあるので、注意してください。 間違えていると、大抵は有り得ないエネルギーが得られたり、SCFが収束しなかったりするので分かると思います。 オングストロームに1.889726ぐらいをかけると原子単位になります。 ためしにprint(my.atom_coords())というコマンドを行ってみると、原子単位での座標を出力してくれます。

basisは「基底関数」を指定します。 基底関数が量子化学計算に導入される過程は一応Szabo: Sec.3.4.2に書いてあります。 簡単に言えば、狭義のHartree--Fock方程式（積分微分方程式）を数値的に解くのは難しいため、基底関数を導入してHartree--Fock--Roothaan方程式を解くようにしています。 こうすると、微分方程式から線形代数の問題へ変換でき、ご存じ（？）の通り、繰り返し計算を行うことで割と簡単に方程式を解くことができます。

基底関数は、原子の上に置く関数のことで、分子軌道（ ）は基底関数（原子軌道）（）の線形結合により表されます。 よくLCAO (linear combinations of atomic orbitals)とかいうやつですね。

(Szabo: 式(3.133))

は分子軌道係数です。詳しくはSzabo: Sec.3.4.2や3.5.1をご参照ください。 量子化学の授業でも出てきた概念だと思います。 例えばヒュッケル法を用いてベンゼンの永年方程式を解くと、6つのp軌道（これが原子軌道）の組み合わせで6つのπ軌道・π*軌道（これが分子軌道）が出てくることを確認しました。 量子化学計算でも同じようなことをしていますが、p軌道だけではなくs軌道、d軌道、f軌道･･･と高次の角運動量を持つ軌道関数を用いることも少なくありません。 また、主量子数に関しても量子化学で勉強する以上の値を用いることも少なくありません（例えばC原子でも3s, 3p, 3d軌道を用いたり）。

予想されるとおり、基底関数が大きければ大きいほど（与えられた手法の中で）厳密な結果に近づいていきます。 無限に大きな基底関数を用いることができれば良いのですが、無限に計算時間が必要なため、有限の基底関数を用いる必要があります。 現実には、手元の計算リソースで用いることができる最大の基底関数を用いて研究を行っていく必要があります。 と言っても、ある程度妥協は必要ですが。

上記の例ではbasis='ccpvdz'としており、cc-pVDZ (correlation consistent polarized valence double zeta)という基底関数を用いて計算をしました。 当面はcc-pVDZ（ccpvdz）、6-31G(d)/6-31G*（6-31g(d)）、STO-3G（sto-3g）あたりを用いれば良いと思います。 上に挙げた順に、用いる関数の数が減っていきます（STO-3Gが上記の最小）。 指定できる基底関数の数は非常に多く、到底ここでご紹介できるものではありません。 他の選択に興味があれば、EMSL等をご参照ください。 また、どの角運動量の関数を用いているか等も知っておくと良いでしょう（Frank Jensen: Introduction to Computational Chemistryの5章が分かりやすいのですが･･･）。 例えば、cc-pVDZの基底関数を用いると、書き方はあまり正確ではないかもしれませんが、水素は1s + 2s + 2pの軌道を用い、酸素だと1s + 2s + 3s + 2p + 3p + 3dの軌道を用いることになります。

時間があれば、basis='ccpvtz'（cc-pVTZ基底関数）、basis='ccpvqz'（cc-pVQZ基底関数）を試して（計算時間が増えるので注意）みましょう。 大きな基底関数ほど、高次の角運動量の関数（f関数、g関数、...）を用いる、あるいは同じ角運動量でも多くの関数を用いていく（関数の広がり方が異なる）ようになります。

ここでは、"mf"という変数（構造体？）の中に、「"scf.RHF(mol)"という計算を行う」という情報（命令？）を代入しています。

ちなみに、scfは"self-consistent field"の事です。 日本語では「自己無撞着場（じこむどうちゃくば）」とか「自己無矛盾場（じこむむじゅんば）」とか訳されたりしますが、 逆に意味が分からないのでSCFなりself-consistent fieldなり書きましょう。 また、RHFは"restricted Hartree--Fock"の事を意味しています。日本語では「制限Hartree--Fock」と呼ばれます（こちらは通じる）。 この辺りはSzabo先生の教科書で勉強することになると思います（主にChapter 3）。 最初のimport gto, scfは、このscf.RHFという関数を使うために必要だったことが分かります。

"mf"は、いろいろな情報（変数）を持っています。 マニュアルこの部分にいろいろ書かれています。 現状だとあまり用いる意味はありませんが･･･例えば、verboseという変数があります。 これは

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.verbose = 4            # この行を追加
mf.kernel()

という感じで変数を変更します。 どうなるかは、試してみてください。

以上で計算の準備が完了しました。やっと計算を行います。 本当はmf.***のアスタリスクの所に、具体的な関数（たぶんmf = scf.RHF(mol).run()）を使って計算するところですが、PySCFでは単に「kernel」という関数を用いるだけでうまいこと計算してくれます。 具体的にどのような関数を用いることができるかは、マニュアルのSCFのページをご参照ください。

思ったよりも長くなってしまったので、出力ファイルの説明に関しては別ページにします。

H₂O（H-O距離が1.0 Angstrom、H-O-Hが120度）のHF/cc-pVDZ、HF/cc-pVTZ、HF/cc-pVQZのエネルギーを計算してみましょう。 また、それぞれの「相対エネルギー」をkcal/mol単位（小数点二桁までで良い）で計算してみましょう。 計算時間を比較してみると良いかもしれません。 今回は初めての計算ということで、私が計算した結果でも。

基底関数の大きさとしては、cc-pVDZ < cc-pVTZ < cc-pVQZです。 基底関数が大きくなると、基本的にはエネルギーが減少していきます。なぜでしょうか。

他の任意の分子に関しても、時間があれば同様の計算をしてみましょう。 ただし原子の数が多いと、分子構造の座標を作成するのが面倒（実際には可視化ソフトウェアを用いて分子構造を作成します）だったり、計算時間がとても長くなります。