勾配と構造最適化

勾配と構造最適化#

ここまでのあらゆる計算は固定した構造 — あなたが打った座標 — を使いました。しかし打った構造が平衡構造であることはまれです。本章は核にかかる力（エネルギー勾配）と、それを使って分子を平衡構造へ緩和することについてです。これが「推測した構造でのエネルギー」を「自然が実際に取る構造」に変えます。

理論:勾配は力である#

Born–Oppenheimer の描像では、電子エネルギーは核位置の関数 \(E(\mathbf R)\) — ポテンシャルエネルギー曲面（PES） — です。原子 \(A\) にかかる力はその曲面の勾配の符号を反転したもの、

\[ \mathbf F_A = -\frac{\partial E}{\partial \mathbf R_A}, \]

であり、平衡構造はすべての力が消える停留点 \(\partial E/\partial \mathbf R = 0\)（PES 上の極小）です。qc-rs は勾配を解析的に — 有限差分ではなく微分積分から直接 — 計算します。速く正確です。組み込みの自己検査が、力が 並進不変性に従うことを確認します:\(\sum_A \mathbf F_A \approx 0\)（分子全体を剛体的にずらしてもエネルギーは変わらない）。

なぜ勾配に軌道応答が要らないのか — しかし Pulay 項は要る#

素朴には \(\partial E/\partial R_A\) は \(\partial\mathbf C/\partial R_A\) — 核が動くとき軌道自体がどう応答するか — を必要としそうに見えます。Hellmann–Feynman 定理はそうではないと言います:収束した SCF エネルギーは軌道回転について停留している（\(\mathbf g=0\)、SCF 収束理論の条件）ので、軌道応答項はすべて 1 次微分から打ち消し合い、演算子の明示的な \(R\) 依存性だけが残ります:

\[ \frac{\partial E}{\partial R_A} = \sum_{\mu\nu} D_{\mu\nu}\,\frac{\partial h_{\mu\nu}}{\partial R_A} + \tfrac12\sum_{\mu\nu\lambda\sigma} D_{\mu\nu}D_{\lambda\sigma}\,\frac{\partial(\mu\nu|\lambda\sigma)}{\partial R_A} \;-\; \sum_{\mu\nu} W_{\mu\nu}\,\frac{\partial S_{\mu\nu}}{\partial R_A} \;+\; \frac{\partial V_{nn}}{\partial R_A}, \]

エネルギー重み付き（Lagrange）密度 \(W_{\mu\nu}=2\sum_i \varepsilon_i\,C_{\mu i}C_{\nu i}\)（RHF;UHF はスピンごとの和、ROHF は射影子形 \(W=\sum_\sigma P_\sigma F_\sigma P_\sigma\)）が、最後の一電子項 — Pulay 力 \(-\operatorname{Tr}(WS^{[1]})\) — で通常の密度の代わりに現れます。

これは純粋に幾何学的な理由で存在します:qc-rs のガウス基底関数は核に付随しているので、核が動くとその AO も一緒に動き — 軌道係数はそれ以上微分する必要がないのに、重なり行列自体が \(R_A\) に依存します。Pulay 項を省くと、勾配はまさにこの基底追従の寄与だけ間違います — 分野で最もよくある勾配のバグの 1 つです。上の演算子項には 2 種の \(R\) 依存性が現れます:AO 中心部分（微分されたガウス関数自体が動く）と Hellmann–Feynman 演算子中心部分（核引力演算子 \(-Z_B/|\mathbf r-\mathbf R_B|\) はあらゆる核 \(B\) に明示的に依存するので、微分される殻がどこにあろうと \(B\) に力を寄与する）。qc-rs の qc-grad crate はまさにこれらの項 — 運動エネルギー、核引力（両方の部分）、二電子反発、Pulay の重なり項、\(V_{nn}\)、有効なら ECP/PCM/D3-D4 — を組み立てます。これは既に知っている energy_components のエネルギー分解を映しています。二電子項と DFT-XC 勾配の縮約を含む完全な導出は解析的微分にあります。

ある構造での勾配#

SCF の後、力はアクセサ 1 つで得られます:

import qc, numpy as np
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

done = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run()
grad = np.asarray(done.scf.gradient)          # 形 (natom, 3)、原子単位（Ha/bohr）

grad.shape                 # (3, 3)
np.abs(grad).max()         # 0.014317   最大の力成分
grad.sum(axis=0)           # ~[0, 0, 0]  並進不変性の自己検査

done.scf.gradient（qc.grad(done) と等価）は \(\partial E/\partial \mathbf R\) の [natom, 3] 配列を hartree/bohr で返します。非ゼロの勾配（ここでは最大 0.0143）は分子が平衡にないことを意味します — 力はより良い構造へ下り坂を指しています。勾配は RHF/UHF/ROHF と KS-DFT（ハイブリッド含む）に対応し、有効にした ECP・PCM・DFT-D3/D4 の寄与を自動的に含みます。

構造最適化#

その力をたどって極小へ行くには、.opt() ステップを足します。これは解析的勾配を geomeTRIC の内部座標オプティマイザで駆動し、力（と変位）がしきい値を下回るまでステップを取ります:

opt = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc="b3lyp").opt().run()

opt.opt.converged     # True
opt.opt.energy        # -76.420628   最適化構造でのエネルギー
opt.opt.e_traj        # [-76.420349, ..., -76.420628]   最適化ステップごとのエネルギー
opt.opt.gmax_traj     # ステップごとの最大勾配（0 へ収束）
opt.opt.grms_traj     # ステップごとの rms 勾配
opt.coordinates()     # 最適化された座標（bohr）

.opt() は SCF に連結します:各最適化ステップは少し動いた構造での完全な SCF で、geomeTRIC が勾配から次の動きを提案します。最後のチェックポイントは最適化された構造とその電子状態を保持し、物性や振動解析の準備が整います。オプティマイザは実行中に自身のステップごとの進捗（エネルギー・勾配・信頼半径）を表示します。

# 緩和した構造を読み戻す
c = np.asarray(opt.coordinates()) * 0.52917721092   # bohr -> angstrom
# この B3LYP/cc-pVDZ の水では:O–H ≈ 0.9687 Å、H–O–H ≈ 102.7°

Tip

opt() のオプション opt(coordsys=..., maxiter=...) は geomeTRIC の一般的な制御を公開します — coordsys は座標系を選び（"tric" 並進回転内部座標が頑健な既定）、maxiter はステップ数の上限です。全体は SCF/最適化リファレンスにあります。

収束した最適化を読む#

*_traj アクセサは最適化の物語です — 挙動の確認に使います:

import numpy as np
e = np.asarray(opt.opt.e_traj)
print("steps          :", len(e))                    # 4
print("energy lowered :", round(e[0] - e[-1], 6), "Ha")   # 0.000279
print("final max-grad :", float(np.asarray(opt.opt.gmax_traj)[-1]))  # ~1e-5、しきい値以下

健全な最適化では、エネルギーが単調に減少し、最大勾配がゼロへ縮小します。geomeTRIC の収束基準（Gaussian 流）は、エネルギー変化・RMS/最大勾配・RMS/最大変位がすべてしきい値を下回ることを要求します — だから converged=True は「エネルギーが変化しなくなった」より強い主張です。

総合例:緩和でどれだけ下がるか#

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"   # わずかにずれた入力構造

single = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc="b3lyp").run()
opt    = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc="b3lyp").opt().run()

print(f"single point : {single.scf.energy:.6f}")   # -76.420349
print(f"optimized    : {opt.opt.energy:.6f}")       # -76.420628
print(f"relaxation   : {(single.scf.energy - opt.opt.energy)*627.509:.3f} kcal/mol")  # 0.175

最適化エネルギーはより低く（そうでなければならない — 緩和は PES を下るだけ）、ここでは入力が既に近かったので約 0.18 kcal/mol です。悪い初期構造ならその利得ははるかに大きく、そして構造（結合長・角度）こそが実際に欲しいものであることも多いです。

練習 6

入力構造での SCF の後、np.abs(done.scf.gradient).max() が 2e-6 でした。分子は平衡にありますか。ここから .opt() は何をしますか。
分子によらず、sum(gradient, axis=0) ≈ 0 が勾配実装の有用な自己整合チェックになるのはなぜですか。
最適化が maxiter ステップ後に converged=False を返しましたが、e_traj はまだ着実に下がっています。考えられる原因と対処は何ですか。

解答練習 6

ほぼ平衡です — 2e-6 Ha/bohr の最大力は通常の最適化しきい値を下回るので、（ある）停留点にあります。.opt() は微小なステップを 1〜2 回取り、ほぼ即座に converged=True を報告するでしょう。
分子全体の剛体並進はエネルギーを変えられないので、正味の力はゼロでなければなりません:\(\sum_A \mathbf F_A = 0\) （厳密に）。これに反する計算勾配はバグや不整合（項の欠落、符号違い）を露呈します — 無料で分子非依存の検査です。
オプティマイザが（まだ下降中の）極小に達する前にステップを使い切っただけです。maxiter を上げる（opt(maxiter=...)）か、チェックポイントが既に保持している最後の構造から最適化を再開します。

エネルギー・相関・力・最適化構造が揃ったので、残りのガイド章は環境（溶媒和と分散）・ 可視化・ログを加え、その後に大規模な分子物性スイートへ進みます。