冗長内部座標系における測地線アプローチを用いた幾何構造最適化の加速

6957 words

35 minutes

冗長内部座標系における測地線アプローチを用いた幾何構造最適化の加速

2026-01-29

Theoretical Chemistry

Computational Chemistry

Geometry Optimization

Differential Geometry

Geodesic Equation

Redundant Internal Coordinates

last_modified: 2026-01-29

生成AIによる自動生成記事に関する免責事項: 本記事は、J. Chem. Phys.に掲載された学術論文（Hermes et al., 2021）に基づき、大規模言語モデルによって作成された解説記事です。記述内容は微分幾何学および計算化学の数理的定義に基づき正確性を期していますが、厳密な証明や詳細な実装については、原著論文および関連する標準的な教科書を参照してください。

1. 序論：構造最適化における座標系の課題#

分子、固体、その他の原子系の計算機モデリングにおいて、幾何構造最適化は基礎的かつ重要な工程である。このプロセスは、ポテンシャルエネルギー曲面（PES）上の停留点（極小点または鞍点）を探索することを目的とする。

1.1 カーテシアン座標と内部座標#

分子の幾何構造を記述する最も直感的な方法は、各原子の空間位置を指定するカーテシアン座標 $\mathbf{x} \in \mathbb{R}^{3n}$ （ $n$ は原子数）を用いることである。凝縮相システムなど、原子の動きが近接原子によって空間的に拘束されている場合、このアプローチは有効である。しかし、気相分子においては、原子団が協奏的な運動によって長距離を変位することが頻繁に生じる。カーテシアン空間においてこのような変位はノルムが大きく、非効率な最適化ステップを多数必要とする傾向がある。

この問題を回避するため、化学的に意味のある特徴量（結合長、結合角、二面角など）の集合である内部座標 $\mathbf{q} \in \mathbb{R}^{m}$ が一般的に採用される。内部座標基底においては、最適化アルゴリズムがこれらの化学的特徴を直接操作するため、曲線的なステップを通じて大規模な変位を効率的に実現可能となる。

1.2 冗長性と多様体構造#

非線形分子の自由度は $3n-6$ であるが、化学的に重要な結合や角度を全て記述しようとすると、変数の数 $m$ は $3n-6$ を超過する。これを「冗長内部座標系（Redundant Internal Coordinates）」と呼ぶ。冗長内部座標系においては、任意のベクトル $\mathbf{q}$ が必ずしも物理的に実現可能な原子配置に対応するとは限らない（例：ベンゼン環の閉環条件など）。物理的に有効な内部座標の集合は、 $m$ 次元空間内に埋め込まれた $(3n-6)$ 次元の「多様体（Manifold）」を形成する。

1.3 従来法の幾何学的欠陥#

従来の最適化手法（一般にニュートン法と称されるアプローチ）では、探索方向 $\Delta \mathbf{q}$ を算出した後、多様体への射影（Projection）または逆変換を行うことで次の構造を決定する。しかし、この直線的なステップとその後の補正操作は、多様体の「曲率（Curvature）」を考慮していない。特に、縮合環を持つ分子のように内部座標間のカップリングが強い系では、この曲率が無視できず、収束性の悪化や発散を招く要因となる。

本稿で解説するHermesらの研究は、多様体の曲率に沿った「測地線（Geodesic）」を用いることで、この幾何学的問題を解決し、最適化効率を劇的に改善するものである。

2. 理論的枠組みと数理的定式化#

本節では、測地線アプローチの基礎となる微分幾何学的な定式化を行う。

2.1 座標変換とウィルソンB行列#

カーテシアン座標 $\mathbf{x}$ から内部座標 $\mathbf{q}$ への変換は非線形であり、その局所的な線形変換はヤコビアン行列、通称ウィルソンB行列（Wilson B-matrix）によって与えられる。

B_{ij} = \frac{\partial q^i}{\partial x^j} \quad (1)

また、冗長内部座標空間における変位 $\Delta \mathbf{q}$ が多様体の接空間（Tangent Space）にあることを保証するために、非局所化内部座標（Delocalized Internal Coordinates）への射影がしばしば用いられる。これはB行列の特異値分解（SVD）により構築される射影行列 $\mathbf{U}$ を用いて行われる。

\mathbf{B} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T \quad (2)

ここで $\mathbf{U}$ はB行列の左特異ベクトルであり、非ゼロ特異値に対応する成分が接空間の基底を張る。

2.2 従来のステップ（ニュートン法）の問題点#

従来のRFO（Rational Function Optimization）やBFGS法では、変位ベクトル $\Delta \mathbf{q}$ を求めた後、対応するカーテシアン座標 $\mathbf{x}_{new}$ を求めるために、以下の反復式（ニュートン・ラフソン法のような逆変換）を用いることが一般的である。

x_{(k+1)}^{i} = x_{(k)}^{i} + (B_{(k)}^{+})_{\lambda}^{i} (q_{0}^{\lambda} + \Delta q^{\lambda} - q_{(k)}^{\lambda}) \quad (3)

ここで、 $B^+$ はB行列のムーア・ペンローズ擬似逆行列（Moore-Penrose pseudo-inverse）である。また、ギリシャ文字の添字は内部座標成分、ラテン文字の添字はカーテシアン座標成分を表す（アインシュタインの縮約記法を採用）。この式は、カーテシアン空間での直線的な変位の積み重ねによって、内部座標空間上の目標点 $q_0 + \Delta q$ に最も近い多様体上の点を探す操作である。しかし、図形的には、これは接ベクトル方向に進んだ後に多様体へ「垂直に」戻る操作に近く、多様体自体が大きく曲がっている場合、目標とする極小点への経路から逸脱する可能性がある。

2.3 測地線方程式の導入#

Hermesらが提案する手法は、変位ベクトル $\Delta \mathbf{q}$ を「測地線の接ベクトル」として解釈するものである。多様体上の2点間を結ぶ最短経路である測地線 $\mathbf{q}(\tau)$ は、以下の測地線方程式（Geodesic Equation）を満たす。ここで $\tau$ はアフィンパラメータ（ $0 \le \tau \le 1$ ）である。

\frac{d^2 q^\lambda}{d\tau^2} + \Gamma^\lambda_{\mu\nu} \frac{dq^\mu}{d\tau} \frac{dq^\nu}{d\tau} = 0, \quad \lambda = 1, \dots, m \quad (4)

ここで $\Gamma^\lambda_{\mu\nu}$ は第2種クリストッフェル記号（Christoffel symbols of the second kind）であり、座標系の接続（Connection）を規定する量である。初期条件は $\mathbf{q}(0) = \mathbf{q}_0$ 、 $\dot{\mathbf{q}}(0) = \Delta \mathbf{q}$ となる。この方程式を $\tau=1$ まで積分することで、多様体の曲率に沿った自然な変位後の座標が得られる。

注：第2種クリストッフェル記号という表現について、原典では形式的に同型だが、明示的な計量に基づく Levi-Civita 接続を仮定しているわけではないようだ。

2.4 カーテシアン基底における実装#

式(4)を直接解くことは困難である。なぜなら、内部座標 $\mathbf{q}$ は独立変数ではなく、カーテシアン座標 $\mathbf{x}$ の関数として従属的に定義されているからである。そこで、著者は測地線方程式をカーテシアン座標の変数 $\mathbf{x}(\tau)$ に対する微分方程式に変換するアプローチを採用した。

連鎖律 $\dot{q}^\lambda = B^\lambda_k \dot{x}^k$ およびその時間微分を用いると、最終的に以下の常微分方程式（ODE）が得られる。

\ddot{x}^i + (B^+)^i_\lambda \frac{\partial^2 q^\lambda}{\partial x^k \partial x^l} \dot{x}^k \dot{x}^l = 0, \quad i = 1, \dots, 3n \quad (5)

ここで、初期条件は $\mathbf{x}(0) = \mathbf{x}_0$ および $\dot{\mathbf{x}}(0) = B^+ \Delta \mathbf{q}$ である。この式(5)は極めて重要である。第2項には内部座標のカーテシアン座標に関する二階微分 $\frac{\partial^2 q}{\partial x^2}$ が含まれており、これが多様体の曲率情報を担っている。この項が存在することで、原子の運動は内部座標の定義（結合の伸縮や角度の回転）に自然に従う曲線的な軌道を描くことになる。

2.5 勾配の平行移動（Parallel Transport）#

構造最適化においては、次のステップの探索方向を決定するためにヘシアン（力の定数行列）の更新が必要となる（BFGS法など）。通常、セカント条件（Secant Condition）が用いられるが、曲がった多様体上では、異なる点 $\mathbf{q}_0$ と $\mathbf{q}_1$ におけるベクトルを単純に比較・減算することは幾何学的に正しくない。ベクトルはそれぞれの点の接空間に属しているからである。

正しい更新を行うためには、点 $\mathbf{q}_0$ における勾配ベクトル $\mathbf{g}_0$ を、測地線に沿って点 $\mathbf{q}_1$ まで「平行移動（Parallel Transport）」させる必要がある。平行移動されたベクトル $\tilde{\mathbf{g}}_0$ を用いた修正セカント条件は以下のようになる。

H_{\lambda\mu} (\dot{q}(1))^\mu = (g_1 - \tilde{g}_0)_\lambda \quad (6)

平行移動もまた微分方程式として記述され、式(5)と連立させて解くことが可能である。

\dot{y}_3^i = - (B^+)_\lambda^i \frac{\partial^2 q^\lambda}{\partial x^k \partial x^l} \dot{y}_2^k \dot{y}_3^l \quad (7)

ここで $y_3$ は勾配ベクトルに対応する補助変数である。

3. アルゴリズムと計算機実装#

本手法の実装における計算コストと数値解法について詳述する。

3.1 ODEソルバーの利用#

式(5)および式(7)は、変数を $\mathbf{y} = [\mathbf{x}, \dot{\mathbf{x}}, \mathbf{g}_{trans}]^T$ と置くことで、一階の常微分方程式系に帰着される。著者の実装（Sellaパッケージ）では、SciPyライブラリの LSODA または CVODE ソルバーが用いられている。これらは硬い（stiff）方程式とそうでない方程式の両方に適応的に対応できる多段階法を採用しており、高い精度で積分を実行できる。

3.2 計算コストの評価#

測地線アプローチの追加コストは主に以下の点にある。

B行列の微分: $\frac{\partial^2 q}{\partial x^2}$ の計算が必要となる。これは解析的に導出可能であり（Hollman & Schaefer, 2012等）、多くの標準的な内部座標（距離、角度、二面角）に対して実装されている。また、このテンソルは疎（Sparse）であるため、計算量は原子数に対して過度に増大しない。
ODEの数値積分: 構造最適化の1ステップごとにODEを解く必要があるが、電子状態計算（DFT計算等）によるエネルギー・勾配評価のコストに比べれば、この幾何学的な計算のオーバーヘッドは無視できるほど小さい。

3.3 ソフトウェア実装 “Sella”#

本手法は、Pythonパッケージ Sella としてオープンソースで公開されている。このソフトウェアは、本来は鞍点探索（Saddle Point Optimization）を主眼に置いているが、本論文ではその最小化（Minimization）機能に焦点が当てられている。

4. 結果と考察#

提案手法の有効性を検証するため、BirkholzとSchlegelによるベンチマークセット（20〜95原子の20分子）を用いた比較実験が行われた。比較対象は、従来の「ニュートン法（標準的な逆変換ステップ）」と本手法「測地線法」である。電子状態計算にはDFTB+（DFTB3パラメータ）が用いられた。

4.1 収束ステップ数の大幅な削減#

論文中のTable Iに示された結果は、測地線アプローチの優位性を明確に示している。

全般的な傾向: テストされた全ての分子において、測地線法はニュートン法よりも少ないステップ数で収束に達した。
平均ステップ数:
- ニュートン法: 198.5ステップ
- 測地線法: 71.6ステップ
- 標準偏差も測地線法の方が小さく、安定した性能を示している。
顕著な例:
- Azadirachtin (95原子): 108ステップ $\to$ 45ステップ
- Sphingomyelin (84原子): 217ステップ $\to$ 164ステップ（さらに、より低いエネルギー準位へ収束）
- Mg porphyrin (37原子): 86ステップ $\to$ 15ステップ

4.2 幾何学的考察#

Fig. 2に示された最適化軌跡の解析から、以下の特性が明らかになった。

初期段階: 最適化の初期は、勾配の主要成分が結合伸縮（Bond stretching）にあることが多い。結合伸縮はカーテシアン空間でも比較的直線的な挙動を示すため（多様体の曲率が小さい）、この段階では両手法の差は小さい。
後期段階: 結合長が緩和された後、結合角や二面角といった「柔らかい」モードの最適化が支配的になる。これらの座標変位はカーテシアン空間において強い非線形性（高い曲率）を持つ。
スパイクの抑制: ニュートン法では、曲率の高い領域でエネルギーが跳ね上がる（Spike）現象が頻繁に観測された。これは、直線的なステップが多様体から大きく外れ、無理やり射影されることで結合長などが歪むためである。一方、測地線法ではこのようなエネルギー上昇が極めて少なく、滑らかに極小点へ向かう軌跡を描いている。

4.3 物理的解釈#

測地線法が優れている本質的な理由は、ステップの「大きさ」ではなく「軌道（Trajectory）」にある。従来のニュートン法では、角度を変える際に原子を直線的に動かすため、必然的に結合距離が縮んだり伸びたりする「寄生的な」変形を伴う。測地線法では、式(5)の第二項が一種の「遠心力」や「拘束力」として働き、結合距離を保ったまま角度を変えるような、分子の自然な動き（回転運動など）を再現する。これにより、高エネルギー領域への意図しない侵入を防ぎ、効率的な探索が可能となる。

5. 結論と展望#

Hermesらの研究は、冗長内部座標系を用いた幾何構造最適化において、多様体の微分幾何学的性質を陽に取り入れることの重要性を実証した。

5.1 本研究の成果#

測地線ステップの定式化: 内部座標多様体上の測地線をカーテシアン空間のODEとして定式化し、実用的な計算コストで解く手法を確立した。
平行移動の実装: ヘシアン更新における幾何学的な不整合を解消するため、勾配ベクトルの平行移動を導入した。
実証的な加速: 複雑な分子系において、最適化ステップ数を平均で60%以上削減するという劇的な効率化を達成した。

5.2 今後の展望#

本手法は、単なる極小点探索だけでなく、鞍点探索（遷移状態探索）においても同様、あるいはそれ以上の効果を発揮することが示唆されている。反応経路探索においては、谷底の形状に沿って進むことが極めて重要であるため、曲率を考慮した測地線アプローチは標準的な手法となるポテンシャルを秘めている。

本稿で解説した技術は、計算化学における「座標系の選択」という古典的な問題に対し、現代的な幾何学の視点から与えられた強力な解答であると言える。

参考文献#

Primary Source: Hermes, E. D., Sargsyan, K., Najm, H. N., & Zádor, J. (2021). Geometry optimization speedup through a geodesic approach to internal coordinates. The Journal of Chemical Physics, 155(9), 094105.
Wilson B-Matrix: Wilson, E. B., Decius, J. C., & Cross, P. C. (1955). Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra. Dover Publications.
Standard Optimization: Pulay, P., & Fogarasi, G. (1992). Geometry optimization in redundant internal coordinates. The Journal of Chemical Physics, 96(4), 2856-2860.
B-Matrix Derivatives: Hollman, D. S., & Schaefer III, H. F. (2012). Efficient evaluation of the second derivative of the Wilson B matrix. The Journal of Chemical Physics, 137(16), 164103.
Benchmark Set: Birkholz, A. B., & Schlegel, H. B. (2016). Geometry optimization in delocalized internal coordinates: an efficient solution for large molecules. Theoretical Chemistry Accounts, 135, 84.
Sella Software: Hermes, E. D. (2021). Sella: A geometry optimization package. Zenodo. http://doi.org/10.5281/zenodo.4747052

補足:実装する際の手順を考えてみた#

1. 概要#

本稿は、冗長内部座標系における測地線アプローチの実装仕様を定義するものである。前版からの主な変更点は、常微分方程式（ODE）の右辺計算における線形代数演算の効率化指針と、入力物理量の定義の明確化にある。

2. 数理モデルと定式化#

2.1 状態ベクトル#

拡張状態ベクトル $\mathbf{Y}(\tau) \in \mathbb{R}^{9n}$ を以下のように定義する。

\mathbf{Y} = \begin{pmatrix} \mathbf{x} \\ \dot{\mathbf{x}} \\ \mathbf{p} \end{pmatrix}

2.2 支配方程式#

に基づく支配方程式は以下の通りである。

座標の時間発展: $\dot{y}_1 = y_2$
測地線方程式: $\dot{y}_2 = -(B^+) (\nabla_\mathbf{x}^2 \mathbf{q}) : (y_2 \otimes y_2)$
平行移動: $\dot{y}_3 = -(B^+) (\nabla_\mathbf{x}^2 \mathbf{q}) : (y_2 \otimes y_3)$

ここで、 $(B^+)$ はウィルソンB行列のムーア・ペンローズ擬似逆行列を表す。

3. 数値計算上の重要課題と対策#

3.1 擬似逆行列の計算コストと安定性#

ODEソルバーの各ステップ（および内部的な試行ステップ）において $B(\mathbf{x})$ は変化するため、その都度 pinv(B) を計算することは $O(m^2 n)$ のコストを要し、数値的にも不安定になりやすい。

推奨される対策:

最小二乗解法の利用: $B^+ \mathbf{v}$ を明示的な逆行列計算によって求めるのではなく、線形方程式 $B \mathbf{u} = \mathbf{v}$ の最小ノルム解として求める（例: lstsq）。
SVDの更新/再利用: 積分ステップ幅が小さい場合、前ステップの特異値分解（SVD）結果を初期値として利用する、あるいはヤコビアンの更新公式を適用することでコストを削減する。
正則化: $B B^T$ が特異に近い場合、微小な正則化項（Damping factor）を導入して solve((B @ B.T + epsilon * I), ...) の形式で解くアプローチも検討に値する。

3.2 勾配ベクトルの定義#

初期条件として与える勾配ベクトル $\mathbf{g}_0$ は、内部座標系における勾配 $\nabla_\mathbf{q} V$ でなければならない。カーテシアン勾配 $\mathbf{g}_x = \nabla_\mathbf{x} V$ との関係は $\mathbf{g}_x = B^T \mathbf{g}_q$ であるため、入力値の空間を取り違えないよう厳密に区別する必要がある。

4. Pythonによる実装（改訂版）#

以下に、数値安定性と定義の明確化を反映したコードを示す。

import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import pinv, lstsq

def apply_pseudo_inverse(B, v):
    """
    Bの擬似逆行列とベクトルvの積 (B^+ v) を計算する。
    
    Note:
        実運用では pinv(B) @ v よりも lstsq(B, v) や
        SVD分解済みオブジェクトを用いた計算が推奨される。
        ここでは可読性のため明示的な pinv を残すが、
        大規模系では計算ボトルネックとなる点に注意。
    """
    # 簡易実装: B_pinv = pinv(B); return B_pinv @ v
    
    # より安定な実装例 (最小二乗解):
    # B.T @ (B @ B.T)^-1 @ v  (ただし冗長座標でB B.Tが特異な場合は注意)
    # 一般的には numpy.linalg.lstsq を B.T に対して解く等の方法があるが
    # B^+ = B.T (B B.T)^-1 なので、ここでは簡易的に pinv を使用する。
    # 厳密には: B^+ v -> 解 x 使得 ||Bx - v|| 最小 かつ ||x|| 最小
    
    return pinv(B) @ v

def geodesic_rhs(tau, y_flat, n_atoms, calc_b_matrix, calc_b_prime_contracted):
    """
    測地線方程式の右辺ベクトル計算。
    
    Parameters:
        tau (float): 積分変数
        y_flat (ndarray): 状態ベクトル [x, x_dot, p]
        ...
    """
    n_dim = 3 * n_atoms
    
    # 状態変数の展開
    x = y_flat[0:n_dim]
    v = y_flat[n_dim:2*n_dim]
    p = y_flat[2*n_dim:3*n_dim]

    # B行列の計算
    B = calc_b_matrix(x)

    # テンソル縮約項の計算 (d2q/dx2 : u : w)
    # 計算コスト高のため、この関数の内部実装も最適化が必要
    dq2_vv = calc_b_prime_contracted(x, v, v)
    dq2_vp = calc_b_prime_contracted(x, v, p)

    # 加速度項および平行移動項の算出
    # Eq. (C2): y2_dot = - B^+ @ (d2q : v : v) [cite: 330]
    accel = -apply_pseudo_inverse(B, dq2_vv)

    # Eq. (C3): y3_dot = - B^+ @ (d2q : v : p) [cite: 337]
    p_dot = -apply_pseudo_inverse(B, dq2_vp)

    return np.concatenate([v, accel, p_dot])

def execute_geodesic_step(x0, delta_q, g0_internal, calc_b, calc_b_prime_cont):
    """
    1ステップの測地線最適化を実行する。

    Parameters:
        x0 (ndarray): 初期カーテシアン座標
        delta_q (ndarray): 目標変位ベクトル（内部座標）
        g0_internal (ndarray): 初期勾配ベクトル（内部座標）
        calc_b (callable): B行列計算関数
        calc_b_prime_cont (callable): テンソル縮約計算関数

    Returns:
        tuple: (新しいカーテシアン座標 x1, 平行移動された勾配 g_tilde)
    """
    n_dim = len(x0)
    
    # 1. 初期条件の設定
    B0 = calc_b(x0)
    
    # 初期速度 v0 = B^+ * delta_q [cite: 183]
    v0 = apply_pseudo_inverse(B0, delta_q)
    
    # 平行移動用ベクトルの初期値 p0 = B^+ * g0 [cite: 340]
    p0 = apply_pseudo_inverse(B0, g0_internal)
    
    y0 = np.concatenate([x0, v0, p0])
    
    # 2. ODE積分の実行
    # 論文に従い LSODA (または同等のstiff solver) を指定 
    sol = solve_ivp(
        fun=lambda t, y: geodesic_rhs(t, y, n_dim//3, calc_b, calc_b_prime_cont),
        t_span=(0.0, 1.0),
        y0=y0,
        method='LSODA',
        rtol=1e-6,
        atol=1e-8
    )
    
    if not sol.success:
        raise RuntimeError(f"ODE integration failed: {sol.message}")
        
    # 3. 最終状態 (tau=1) の取得
    y_final = sol.y[:, -1]
    x1 = y_final[0:n_dim]
    p1 = y_final[2*n_dim:3*n_dim]
    
    # 4. 平行移動された勾配の復元
    # Eq. (C4): \tilde{g}_0 = B(x1) @ p(1) [cite: 341]
    B1 = calc_b(x1)
    g_tilde_0 = B1 @ p1
    
    return x1, g_tilde_0

5. 結論と実装上の留意点#

本仕様書で定義されたアルゴリズムは、内部座標多様体の曲率を考慮し、最適化ステップ数を大幅に削減する。実装においては、以下の点に特に留意する必要がある。

計算コスト: $B$ 行列の微分 $\frac{\partial^2 q}{\partial x^2}$ の計算は最も高コストな処理の一つである。実用的なコードでは、自動微分ライブラリの利用や、疎行列演算による最適化が不可欠である。
数値安定性: 冗長座標系において $B$ 行列は正方行列ではなく、条件数が悪化する場合がある。擬似逆行列の計算にはSVDに基づく堅牢なソルバーを用い、必要に応じて正則化を行うことが推奨される。
ソルバーの選択: 化学系によっては方程式が硬い（stiff）性質を持つため、LSODA や CVODE 等の陰解法を含むソルバーの選択が重要である。

参考文献#

Hermes, E. D., Sargsyan, K., Najm, H. N., & Zádor, J. (2021). Geometry optimization speedup through a geodesic approach to internal coordinates. J. Chem. Phys. 155, 094105.
Hollman, D. S., & Schaefer III, H. F. (2012). Efficient evaluation of the second derivative of the Wilson B matrix. J. Chem. Phys. 137, 164103.

補遺：実装上の注意#

1. テンソル縮約におけるアンチパターン#

1.1 全テンソルおよび局所行列の構築回避#

論文中の項 $-(B^+)_\lambda^i \frac{\partial^2 q^\lambda}{\partial x^k \partial x^l} v^k v^l$ を実装する際、最大のボトルネックは二階微分テンソル $\nabla^2 \mathbf{q}$ の扱いにある。

Anti-pattern 1: 全原子に対する $3N \times 3N \times M$ の密テンソルを構築する（メモリ・計算量ともに破綻）。
Anti-pattern 2: 内部座標 $\lambda$ ごとに局所的な原子座標のみを取り出し、局所Hessian行列 $\mathbf{H}^\lambda$ （サイズ $12 \times 12$ 等）を構築して $\mathbf{v}_{loc}^T \mathbf{H}^\lambda \mathbf{v}_{loc}$ を計算する。

1.2 ベストプラクティス：二重縮約の直接評価#

実用的な実装において、行列 $\mathbf{H}^\lambda$ を明示的にメモリ上に確保する必要はない。必要なのは「縮約されたスカラー値」のみである。結合距離、結合角、二面角などの標準的な内部座標については、その二階微分とベクトルの縮約結果を解析的に直接記述できる。

\text{Result}^\lambda = \sum_{k,l \in \text{local}} \frac{\partial^2 q^\lambda}{\partial x^k \partial x^l} v^k v^l

したがって、実装は行列積ではなく、解析的なスカラー関数の総和として記述されるべきである。これは自動微分（JAX/PyTorch）を用いる場合でも同様であり、vmap 等を用いて縮約操作そのものをコンパイルすることが推奨される。

2. 擬似逆行列 $B^+$ の概念転換#

2.1 行列としての $B^+$ の排除#

数式上は $B^+$ が頻出するが、これを「行列」として定義・保持することは避けるべきである。

Anti-pattern: pinv(B) を呼び出し、その結果を行列として保持してベクトルと掛け合わせる。
Best Practice: 「 $B^+$ を掛ける」操作を、「線形方程式 $B \mathbf{x} = \mathbf{y}$ の最小ノルム解 $\mathbf{x}$ を求める」操作（Solver）として再定義する。

2.2 実装上の等価変換#

コード上では以下のように概念を置換する。

\mathbf{a} = -B^+ (\dots) \quad \Longrightarrow \quad \text{Solve } \min_{\mathbf{a}} \| B \mathbf{a} - (\dots) \|

これにより、SVDの打ち切り（cutoff）処理や、反復解法（LSQR等）の適用が容易になり、数値的安定性が向上する。

3. 平行移動の幾何学的解釈の厳密化#

本手法における平行移動は、カーテシアン空間（ユークリッド空間）に埋め込まれた部分多様体上の接続である。これは、埋め込み空間から**誘導された計量（Induced Metric）**に基づくLevi-Civita接続であり、内部座標空間に対して人為的な計量を導入したものではない。

物理的には、勾配ベクトルを移動させる際、多様体の「法線方向成分」を常に除去し続け、接空間（Tangent Space）内に拘束する操作に対応する。

4. Python実装プロトタイプ（Refactored）#

上記指針に基づき、関数名と処理ロジックを刷新した実装例を示す。

import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import lstsq

def contract_d2q_with_vectors(x, v1, v2, internal_coords_def):
    """
    二階微分テンソルと2つのベクトルの二重縮約を計算する。
    
    Design Note:
        局所Hessian行列すら構築せず、直接スカラー値を評価することを意図している。
        各内部座標ごとの縮約結果をベクトルとして返す。
    
    Returns:
        ndarray: shape (n_internal,)
                 Output[i] = sum_{kl} (d2q_i / dx_k dx_l) * v1_k * v2_l
    """
    n_internal = len(internal_coords_def)
    contraction_result = np.zeros(n_internal)
    
    # ※ 実際にはここを JAX.vmap や Numba で高速化する
    for i, coords_info in enumerate(internal_coords_def):
        # 局所的な原子インデックスとベクトル成分の抽出
        # atom_indices = coords_info.indices
        # v1_loc = v1[atom_indices] ...
        
        # 行列を作らず、解析的な縮約式を直接評価する関数を呼ぶ
        # val = eval_contraction_analytic(coords_info.type, x_loc, v1_loc, v2_loc)
        pass 
        
    return contraction_result

def solve_min_norm(B, target_vector):
    """
    線形方程式 B x = target の最小ノルム解を求める。
    
    Note:
        数式上の B^+ @ target に相当するが、行列 B^+ は構築しない。
        大規模系や特異性が問題になる場合は、SVD分解済みのオブジェクトを
        再利用する設計に切り替えること。
    """
    # lstsq は ||Bx - y||_2 を最小化する x を返す
    # rcond は計算精度に応じて調整 (e.g., 1e-6)
    x, _, _, _ = lstsq(B, target_vector, cond=1e-5)
    return x

def geodesic_rhs_refactored(tau, y_flat, internal_def, calc_b):
    """
    測地線方程式の右辺計算（Anti-pattern回避版）
    """
    n_dim = len(y_flat) // 3
    x = y_flat[0:n_dim]
    v = y_flat[n_dim:2*n_dim]
    p = y_flat[2*n_dim:3*n_dim]
    
    # 1. B行列の評価（または更新）
    B = calc_b(x)
    
    # 2. テンソル縮約（スカラー直接評価）
    #    行列生成コストを回避
    dq2_vv = contract_d2q_with_vectors(x, v, v, internal_def)
    dq2_vp = contract_d2q_with_vectors(x, v, p, internal_def)
    
    # 3. 線形ソルバーによる解の導出
    #    B^+ の構築コストを回避
    accel = -solve_min_norm(B, dq2_vv)
    p_dot = -solve_min_norm(B, dq2_vp)
    
    return np.concatenate([v, accel, p_dot])

5. 結論#

本アルゴリズムの実装においては、数式の「形」をそのままコードに写すのではなく、その「演算の意味」を計算機科学的に再解釈することが不可欠である。特に、**「行列・テンソルの構築回避（Matrix-free approach）」と「逆行列計算のソルバー化（Solver approach）」**の2点を徹底することで、論文で示された性能上の利点を、大規模系においても享受することが可能となる。