Givens回転と誤差逆伝播法を組み合わせた自己無撞着場計算の高効率収束法の開発

Abstract

We have recently proposed a highly efficient convergence method for self-consistent (SCF) calculations combining Givens rotation and error back-propagation (EBP) algorithms, referred to as direct Givens rotation (DGR) method [J. Chem. Phys. 162, 014108 (2025)]. The Givens rotation corresponds to unitary transformations that guarantee the orthogonality of molecular orbitals. Complicated gradients constructed through sequential Givens rotations were computed using the EBP technique without deriving explicit forms. This article reviews the proposed DGR method and compares it with conventional methods such as the standard SCF procedure, the second-order SCF method, and the direct inversion in iterative subspace technique for H₂O molecule. The DGR method exhibited a convergence speed comparable to that of the SOSCF method while achieving a lower and more reliable energy.

Translated Abstract

We have recently proposed a highly efficient convergence method for self-consistent (SCF) calculations combining Givens rotation and error back-propagation (EBP) algorithms, referred to as direct Givens rotation (DGR) method [J. Chem. Phys. 162, 014108 (2025)]. The Givens rotation corresponds to unitary transformations that guarantee the orthogonality of molecular orbitals. Complicated gradients constructed through sequential Givens rotations were computed using the EBP technique without deriving explicit forms. This article reviews the proposed DGR method and compares it with conventional methods such as the standard SCF procedure, the second-order SCF method, and the direct inversion in iterative subspace technique for H₂O molecule. The DGR method exhibited a convergence speed comparable to that of the SOSCF method while achieving a lower and more reliable energy.

1 はじめに

自己無撞着場 (SCF) の手続きは，Hartree–Fock (HF) およびKohn–Sham密度汎関数理論計算を解くための標準的な手法であるが，収束の理論的根拠はない．一方，拡張ラグランジュ法 (ALM) [1]や2次 (SO) SCF法 [2]などの直接最小化法では，行列の繰り返し対角化を行わず，勾配を用いてラグランジアンを直接最小化する．SOSCF法では，指数関数形式で表現されたユニタリー行列を介して分子軌道が最適化されるため，軌道の規格直交性条件が保証される．しかし，指数関数は通常一次のTaylor展開で打ち切られるため，ユニタリー性を保証するためにGram–Schmidtの直交化が必要となる．我々は最近，Givens回転の積により構築されたユニタリー行列とALMのラグランジアンを用いた直接最小化法，Direct Givens rotation (DGR) 法を提案した [3]．Givens回転は占有軌道と仮想軌道を混合するユニタリー変換であり，SOSCF法で必要なGram–Schmidtの直交化を回避できる．最小化に必要な複雑な勾配は，誤差逆伝播法 (EBP) [4]を適用することで算出した．EBPはニューラルネットワーク (NN)で用いられる技法の一つであり，勾配を連鎖律に基づいて出力層から入力層方向に計算し，明示的な定式化を回避しつつ順次求める．本稿ではDGR法について簡単に解説し，H₂O分子を対象に制限 HF (RHF) 法へ適用した結果を報告する．

2 理論と実装

SCFによる分子軌道C^SCFは，規格直交条件を満たす初期軌道C⁽⁰⁾とユニタリー行列U^SCFにより表現される．

  

$C^{\text{SCF}}=C^{\left(0\right)}{U}^{\text{SCF}}$(1)

DGR法ではU^SCFをGivens回転の積として構築する．

  

$U^{\text{SCF}}=\prod _{i,a}^{}{G}_{}^{\text{SCF}}\left(i,a\right)$(2)

ここで，添字iとaは，それぞれ占有軌道と仮想軌道を表す．Givens回転G(i,a)は，ii, ia, ai, aa成分にそれぞれcosθ, −sinθ, sinθ, cosθを持ち，それ以外の対角成分および非対角成分には，それぞれ1および0を有する．

  

$G\left(i,a\right)=\left(\begin{array}{ccccccc}1& \cdots & 0& \cdots & 0& \cdots & 0\\ ⋮& \ddots & ⋮& & ⋮& & ⋮\\ 0& \cdots & \cos {\theta }_{ia}& \cdots & -\sin {\theta }_{ia}& \cdots & 0\\ ⋮& & ⋮& \ddots & ⋮& & ⋮\\ 0& \cdots & \sin {\theta }_{ia}& \cdots & \cos {\theta }_{ia}& \cdots & 0\\ ⋮& & ⋮& & ⋮& \ddots & ⋮\\ 0& \cdots & 0& \cdots & 0& \cdots & 1\end{array}\right)$(3)

Givens回転を構成する{θ_ia}は，繰り返し計算によって最適化される．

  

${\Theta }^{\left(n+1\right)}={\Theta }^{\left(n\right)}+\Delta {\Theta }^{\left(n\right)}$(4)

ここで，上付き添え字nは繰り返し回数を示し，Θ⁽ⁿ⁾はn回目の{θ_ia}を要素に持つ行列，ΔΘ⁽ⁿ⁾は変位を表す．本研究では準Newton法により近似的なHesse行列$\tilde{B}$を構築し，ΔΘ⁽ⁿ⁾を以下の式により求めた．

  

$\Delta {\Theta }^{\left(n\right)}=-\alpha {\left({\tilde{B}}^{\left(n\right)}\right)}^{-1}{g}^{\left(n\right)}$(5)

ここで，αは更新幅，gはALMのラグランジアンL^ALMのΘ = Θ⁽ⁿ⁾における勾配である．

  

$g_{ia}^{\left(n\right)}={\left(\frac{\partial L^{\text{ALM}}}{\partial {\theta }_{ia}^{}}\right)}_{\Theta ={\Theta }^{\left(n\right)}}$(6)

L^ALMはRHF法の場合，以下の式で表される．

  

$\begin{array}{c}{L}^{\text{ALM}}=E_0^{\text{RHF}}-2\text{tr}\left(\left(C_{\text{occ}}^{†}S{C}_{\text{occ}}^{}-I_{\text{occ}}\right){\epsilon }_{\text{occ}}^{}\right)\\ +\frac{\mu }{2}{\left(\left(C_{\text{occ}}^{†}S{C}_{\text{occ}}^{}-I_{\text{occ}}\right)\right)}_{}^2\end{array}$(7)

ここで，$E_0^{\text{RHF}}$はRHFエネルギー，C_occはMO係数行列，Sは重なり行列，I_occは単位行列であり，ε_occおよびμはパラメータである．

DGR法ではgをEBPに基づいて算出する．(i, a)を新たにlとおき，初期軌道C⁽⁰⁾に対して，Givens回転をl回作用させた中間量C^(l,n)を導入する．

  

$C^{\left(l,n\right)}=C^{\left(l-1,n\right)}{G}^{\left(l,n\right)}$(8)

ここで，G^(l,n)はΘ⁽ⁿ⁾から構築されたl回目のGivens回転を表す．C^(l,n)を用いた連鎖律により，式(6)は最終的に以下の式に変形される．

  

$g_l^{\left(n\right)}=\mathrm{tr}\left(X^{\left(n\right)}{G}^{\left(N_{\text{pair}},n\right)}{}^{†}\cdots G^{\left(l+1,n\right)}{}^{†}{\left(C^{\left(l,n\right)}{Y}^{\left(l,n\right)}\right)}^{†}\right)$(9)

ここで，Yはia, ai成分にそれぞれ−1, 1, それ以外の成分に0を有する行列であり，XはFock行列Fから構築される．

  

$X^{\left(n\right)}=2F^{\left(n\right)}{C}_{\text{occ}}^{\left(n\right)}-2S{C}_{\text{occ}}^{\left(n\right)}{\epsilon }_{\text{occ}}^{\left(n\right)}$(10)

C^(l)は以下の式により逐次的に計算することが可能である．

  

$C^{\left(l,n\right)}=C^{\left(l+1,n\right)}{G}^{\left(l+1,n\right)†}$(11)

C^(l,n)をNNにおけるl層目のノードとみなすと，式(9)はXに対して出力層側から入力層方向にG^(l,n)を順次作用することで勾配を算出している．これはNNにおけるEBPに対応する．

3 結果と考察

本稿ではH₂O分子に対してDGR法を適用した結果を示す．標準的なSCFの手続きであるHF–Roothaan (HFR) 法に加え，direct inversion in iterative subspace (DIIS) 法 [5] およびSOSCF法と比較した．近似的なHesse行列はBroyden–Fletcher–Goldfarb–Shanno (BFGS) アルゴリズムにより構築し，更新幅はArmijo条件にしたがって決定した．基底関数として6-31Gを用いた．O-H結合距離を1.0, 2.5, 5.0 Åとした3種類の構造を作成した．HOH結合角はいずれも104.5°である．

Fig. 1にエネルギーの収束過程を示す．Fig1a-cのO-H結合距離はそれぞれ1.0, 2.5, 5.0 Åである．横軸は繰り返し回数，縦軸は全エネルギーである．黒色，緑色，紫色，赤色の点は，それぞれHFR，DIIS，SOSCF，DGR法を示す．Table 1は各手法の収束エネルギー，繰り返し回数，計算時間を示す．

Fig. 1.

Convergence behavior of HFR, DIIS, SOSCF and DGR methods for H₂O. O-H bond lengths were (a) 1.0, (b) 2.5, and (c) 5.0 Å.

Table 1. Comparison of convergence energy, number of iterations, and computational time for HFR, DIIS, SOSCF, and DGR methods. O-H bond lengths were (a) 1.0, (b) 2.5, and (c) 5.0 Å.

	Method	Energy[hartree]	Iteration	Wall-clockTime [s]
(a)	SCF	−75.98016	19	1.0
	DIIS	−75.98016	12	0.7
	SOSCF	−75.98015	20	44.1
	DGR	−75.98016	28	32.5
(b)	SCF	Not Conv.	-	-
	DIIS	−75.41214	15	1.2
	SOSCF	−75.44791	30	65.8
	DGR	−75.44895	35	38.2
(c)	SCF	Not Conv.	-	-
	DIIS	Not Conv.	-	-
	SOSCF	−75.38698	49	92.2
	DGR	−75.38698	48	34.9

結合距離が短いFig. 1aにおいてはHFR，DIIS, DGR法が同じエネルギーに収束した．SOSCF法はそれよりも0.00001 hartree高いエネルギーに収束した．収束までの繰り返し計算の数はDIIS, HFR, SOSCF, DGR法の順に少なかった．

結合距離を伸ばしたFig. 1bにおいては4つの手法のうちHFR法は2つの状態を振動し，収束しなかった．他の3つの手法については，収束までの繰り返し計算の数はDIIS, SOSCF, DGR法の順に少なかった．一方，収束エネルギーはDGR, SOSCF, DIIS法の順に高くなり，変分原理の観点からはDGR法が最も正確であると言える．

結合距離をさらに伸ばしたFig. 1cにおいてはHFR, DIIS法は収束せず，SOSCF, DGR法のみが収束した．収束したエネルギーはDGR法とSOSCF法の間で一致した．また，収束までの繰り返し計算の数はDGR法が48回，SOSCFが49回と同程度であった．

4 まとめ

本稿ではGivens回転とEBPを組み合わせたDGR法について解説した．H₂O分子を対象にDGR法をHFR，DIIS，SOSCF法と比較した．DGR法はSOSCF法と同程度の収束速度を示し，SOSCF法よりも低いエネルギーに収束した．

謝辞

本研究は日本科学技術振興財団 (JST) 次世代研究者挑戦的研究プログラムJPMJSP2128の支援を受けたものである．本研究で行った量子化学計算の一部は，自然科学研究機構 (NINS) 計算科学研究センター (RCCS) の計算機を利用して行った (Project: 23-IMS-C039)．

参考文献

• [1]   M. J. D. Powell, "A Method for Nonlinear Constraints in Minimization Problems," in Optimization, edited by R. Fletcher., Academic Press, 1969, pp. 283–298.
• [2]   J. Douady, Y. Ellinger, R. Subra, B. Levy, J. Chem. Phys., 72, 1452 (1980). https://doi.org/10.1063/1.439369
• [3]   R. Oshima, H. Nakai, J. Chem. Phys., 162, 014108 (2025). https://doi.org/10.1063/5.0232518 PMID:39745160
• [4]   D. E. Rumelhart, G. E. Hinton, R. J. Williams, Nature, 323, 533 (1986). https://doi.org/10.1038/323533a0
• [5]   P. Pulay, Chem. Phys. Lett., 73, 393 (1980).https://doi.org/10.1016/0009-2614(80)80396-4