# 勾配と構造最適化

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

## 理論:勾配は力である

[Born–Oppenheimer](../10-foundations/many-electron-and-bo.md) の描像では、電子エネルギーは核位置の関数
$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 収束理論](../10-foundations/scf-convergence-theory.md) の条件）
ので、軌道応答項はすべて 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` の[エネルギー分解](scf.md)を映しています。二電子項と DFT-XC 勾配の縮約を含む完全な導出は
[解析的微分](../10-foundations/analytic-derivatives.md) にあります。

## ある構造での勾配

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

```python
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** の内部座標
オプティマイザで駆動し、力（と変位）がしきい値を下回るまでステップを取ります:

```python
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 が勾配から次の動きを
提案します。最後のチェックポイントは**最適化された**構造とその電子状態を保持し、物性や振動解析の準備が整います。
オプティマイザは実行中に自身のステップごとの進捗（エネルギー・勾配・信頼半径）を表示します。

```python
# 緩和した構造を読み戻す
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` アクセサは最適化の物語です — 挙動の確認に使います:

```python
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` は「エネルギーが変化しなくなった」より強い主張です。

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

```python
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 です。悪い初期構造ならその利得ははるかに大きく、そして*構造*（結合長・角度）こそが実際に欲しい
ものであることも多いです。

:::{exercise}
:label: ex-grad

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

:::{solution} ex-grad
:class: dropdown

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

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