last_modified: 2026-02-02

生成AIによる自動生成記事に関する免責事項: 本記事は、計算化学における数値安定化アルゴリズムを解説するために生成されたものです。記述内容は数理的な厳密性と実装の再現性を重視していますが、特定の実用環境における動作を保証するものではありません。また、検算はしておりますが、不正確な可能性がございます。ご了承ください。

1. 序論：計算化学における座標系の特異性#

分子シミュレーション、特に分子動力学（MD）法や分子軌道法における構造最適化（Geometry Optimization）において、ポテンシャルエネルギー曲面の勾配（力）を正確に評価することは、計算の信頼性を担保する上で最も基本的かつ重要な要件である。

分子のポテンシャルエネルギー関数 $V$ は、通常、結合長、結合角、二面角といった内部座標（Internal Coordinates）を用いて定義される。しかし、運動方程式の積分や力の評価は、計算効率の観点からデカルト座標（Cartesian Coordinates）系で行われることが一般的である。この「定義空間」と「演算空間」の不一致は、座標変換に伴うヤコビアン（Jacobian）の特異点問題を引き起こす。

特に、3つの原子 $i, j, k$ によって定義される結合角 $\theta$ に関する調和ポテンシャル $V(\theta) = \frac{1}{2}k(\theta - \theta_0)^2$ は、結合角が $0^\circ$（重なり）または $180^\circ$（直線構造）に接近した際、その勾配計算において数値的な特異点が生じる。これは物理的な発散ではなく、数式表現上の特異点（Removable Singularity）であるが、倍精度浮動小数点演算においては NaN (Not a Number) や極端な桁落ちを引き起こし、シミュレーションの破綻要因となる。

本稿では、この角度ポテンシャルの特異点問題に対し、変数を余弦 $u = \cos \theta$ に変換し、境界領域において多項式展開を用いることで特異点を除去する「数値的正則化（Numerical Regularization）」手法について詳述する。また、その数理的背景が量子化学計算におけるBoys関数の処理と同型であることを示し、本手法の学術的妥当性を論じる。

2. 角度勾配における特異点の数学的構造#

2.1 連鎖律による特異項の出現#

原子 $i, j, k$ のデカルト座標を $\mathbf{r}_i, \mathbf{r}_j, \mathbf{r}_k$ とし、結合ベクトルを $\mathbf{r}_{ji} = \mathbf{r}_i - \mathbf{r}_j$、$\mathbf{r}_{jk} = \mathbf{r}_k - \mathbf{r}_j$ とする。結合角 $\theta$ は内積を用いて次のように定義される。

$$\cos \theta = \frac{\mathbf{r}_{ji} \cdot \mathbf{r}_{jk}}{|\mathbf{r}_{ji}| |\mathbf{r}_{jk}|} = u$$

ポテンシャルエネルギー $V(\theta)$ の原子 $i$ に対する勾配 $\nabla_{\mathbf{r}_i} V$ は、連鎖律（Chain Rule）により以下のように展開される。

$$\nabla_{\mathbf{r}_i} V(\theta) = \frac{dV}{d\theta} \cdot \frac{d\theta}{du} \cdot \nabla_{\mathbf{r}_i} u$$

ここで、各項を解析すると問題の所在が明らかとなる。

1. エネルギー微分項: $\frac{dV}{d\theta} = k(\theta - \theta_0)$。これは全領域で有限である。
2. 幾何微分項: $\nabla_{\mathbf{r}_i} u$。これも $|\mathbf{r}| \neq 0$ である限り有限である。
3. 変換ヤコビアン項: $\frac{d\theta}{du}$。

逆三角関数の微分公式 $\frac{d}{du} \arccos u = -\frac{1}{\sqrt{1-u^2}}$ より、

$$\frac{d\theta}{du} = -\frac{1}{\sin \theta}$$

この項が、$\theta \to 0$ ($u \to 1$) および $\theta \to \pi$ ($u \to -1$) において発散する。これが特異点の正体である。

2.2 不定形の極限と数値的破綻#

平衡角 $\theta_0 = 0$ の場合、力の大きさは以下の極限を持つ。

$$\lim_{\theta \to 0} \left( \frac{dV}{d\theta} \cdot \frac{d\theta}{du} \right) = \lim_{\theta \to 0} \left( k\theta \cdot \frac{-1}{\sin \theta} \right) = -k \lim_{\theta \to 0} \frac{\theta}{\sin \theta} = -k$$

数学的には、極限値は $-k$ という有限の値に収束する。したがって、この特異点は「除去可能（Removable）」である。しかし、計算機上では以下の問題が発生する。

• ゼロ除算: $\sin \theta$ が計算機の最小表現数（underflow）を下回る場合、除算例外が発生する。
• 桁落ち: $1 - u^2$ の計算において、$u \approx 1$ のとき有効数字が著しく失われる。
• 不定形: 分子が $0$ ($\theta \to 0$)、分母も $0$ ($\sin \theta \to 0$) となる $0/0$ の不定形計算を、極限操作なしに直接実行することになる。

これにより、直線状の分子（例：アセチレン誘導体やニトリル基）や、立体障害により結合角が極端に歪んだ遷移状態において、数値的な不安定性が生じる。

3. 数値的正則化のアルゴリズム#

前述の特異点を回避するため、特異点近傍領域において超越関数（$\arccos$）の使用を避け、多項式近似による正則化を行う。

3.1 領域分割アプローチ#

余弦 $u = \cos \theta$ の値に基づき、定義域を3つの領域に分割する。ここで $\eta$ は微小な閾値（例: $10^{-4}$）である。

1. 通常領域 ($|u| \le 1 - \eta$): 特異点から十分に離れているため、通常の $\arccos$ 関数を用いた厳密計算を行う。
2. 前方特異領域 ($u > 1 - \eta$, $\theta \approx 0$): $0^\circ$ 近傍。$\arccos u$ の代わりに、$\Delta u = 1 - u$ を用いた展開を行う。
3. 後方特異領域 ($u < -1 + \eta$, $\theta \approx \pi$): $180^\circ$ 近傍。$\Delta u = 1 + u$ を用いた展開を行う。

3.2 近似式の導出と係数の決定#

$\theta \approx 0$ 近傍における $\theta^2$ の多項式展開係数は、余弦関数のマクローリン展開の逆関数として一意に決定される。

偏差 $\Delta u = 1 - \cos \theta$ を $\theta$ の偶数次冪級数として展開する。

$$\Delta u = 1 - \left( 1 - \frac{\theta^2}{2!} + \frac{\theta^4}{4!} - \frac{\theta^6}{6!} + O(\theta^8) \right) = \frac{1}{2}\theta^2 - \frac{1}{24}\theta^4 + \frac{1}{720}\theta^6 + O(\theta^8)$$

ここで、$\theta^2$ を $\Delta u$ の冪級数として表現することを仮定し、未定係数法を適用する。

$$\theta^2 = c_1 \Delta u + c_2 (\Delta u)^2 + c_3 (\Delta u)^3 + O((\Delta u)^4)$$

この式を上式に代入し、$\Delta u$ の各次数の係数を比較する。

1次項 ($\Delta u$):

$$\Delta u = \frac{1}{2}(c_1 \Delta u) \implies 1 = \frac{c_1}{2} \implies c_1 = 2$$

2次項 ($(\Delta u)^2$):$\theta^4 = (2\Delta u)^2 = 4(\Delta u)^2$ より、

$$0 = \frac{1}{2}c_2 - \frac{1}{24}(4) \implies c_2 = \frac{1}{3}$$

3次項 ($(\Delta u)^3$):$\theta^6 \approx (2\Delta u)^3 = 8(\Delta u)^3$、および交差項 $2(c_1 c_2) = 4/3$ を考慮し、

$$0 = \frac{1}{2}c_3 - \frac{1}{24}\left(2 \cdot 2 \cdot \frac{1}{3}\right) + \frac{1}{720}(8) \implies c_3 = \frac{4}{45}$$

以上の導出により、採用すべき展開式は数学的に厳密に以下と定まる。

$$\theta^2 \approx 2\Delta u + \frac{1}{3}(\Delta u)^2 + \frac{4}{45}(\Delta u)^3$$

この係数セットは、$\arccos(1-\Delta u)$ の直接展開よりも、$\theta^2$ を直接対象とすることで平方根計算を回避し、かつ数値的安定性に優れる。

3.3 一般角 $\theta_0 \neq 0$ への拡張#

平衡角 $\theta_0$ が $0$ でない場合、$\theta \approx 0$ はポテンシャルの極小点ではなく、高いエネルギーを持つ領域となる。しかし、衝突時や遷移状態においてこの領域を通過する可能性があるため、同様の処置が必要である。

この場合、$\theta \approx \sqrt{P(\Delta u)}$ として扱う。

$$V(u) = \frac{1}{2}k (\sqrt{P(\Delta u)} - \theta_0)^2$$

勾配は以下のようになる。

$$\frac{dV}{du} = k(\sqrt{P} - \theta_0) \cdot \frac{1}{2\sqrt{P}} \frac{dP}{d(\Delta u)} \cdot (-1)$$

一見すると $\frac{1}{\sqrt{P}}$ が特異点に見えるが、$P \approx 2\Delta u$ より、

$$\frac{1}{\sqrt{P}} \frac{dP}{d(\Delta u)} \approx \frac{1}{\sqrt{2\Delta u}} \cdot 2 = \sqrt{\frac{2}{\Delta u}}$$

となり、依然として $\Delta u \to 0$ で発散するように見える。しかし、実際には平衡点が $\theta_0 \neq 0$ の場合、$\theta \to 0$ への接近は物理的に極めて高いエネルギー障壁を登る過程であり、数値計算上は勾配の発散よりも、エネルギーの連続性が重視される。 なお、前述の極限の議論 $\lim_{\theta \to 0} \frac{dV}{du} = -k (1 - \frac{\theta_0}{0})$ となり、$\theta_0 \neq 0$ の場合は物理的にも力が無限大になる（原子核同士の衝突）ため、この発散は物理的に正しい挙動であるとも言える。本手法の主眼は、$\theta_0=0$ における「偽の特異点」の除去と、$\theta_0 \neq 0$ における「数値計算の堅牢性（NaN防止）」にある。

3.4 一般角 $\theta_0 \neq 0$ における特異点の物理的解釈#

平衡角 $\theta_0 \neq 0$ の場合、$\theta \to 0$ の極限においてポテンシャル勾配は物理的に発散する。$\lim_{\theta \to 0} \frac{dV}{d\theta}\frac{d\theta}{du} = \lim_{\theta \to 0} \left[ k(\theta - \theta_0) \cdot \left(-\frac{1}{\sin \theta}\right) \right] \to \infty$数値計算上では $0/0$ の不定形や桁落ちによる NaN の発生は回避せねばならない。本手法における正則化（Regularization）は、この発散を抑制するのではなく、発散に至る過程を滑らかな多項式で記述し、計算機が扱える最大値（有限の巨大な値）として確定させる ことに主眼を置く。これにより、衝突時においても計算の破綻を防ぎ、安定した勾配評価が可能となる。

4. 量子化学計算との理論的相同性#

本手法の妥当性は、量子化学の標準的な手法との対比によっても支持される。特に、Gaussian基底関数を用いた2電子積分の計算に現れる Boys関数 の処理は、本件と数学的に完全に同型である。

4.1 Szabo-OstlundにおけるBoys関数#

Hartree-Fock法等の計算において、核引力積分や電子反発積分には以下のBoys関数 $F_n(x)$ が頻出する。

$$F_n(x) = \int_0^1 t^{2n} e^{-xt^2} dt$$

この関数は解析的には誤差関数 $\text{erf}(\sqrt{x})$ を用いて表されるが、$x \to 0$ の極限において $\frac{\text{erf}(\sqrt{x})}{\sqrt{x}}$ という $0/0$ 型の不定形となる。 量子化学計算プログラムを実装する際に、以下の処理が選択肢の1つとして、数値誤差を避けるために実装される。

1. 閾値判定: 引数 $x$ がある値（例: $10^{-8}$）より小さいかを判定。
2. 多項式展開: 小さい場合、漸近展開（マクローリン展開）を使用する。 $F_n(x) \approx \frac{1}{2n+1} - \frac{x}{2n+3} + \cdots$

この「特異点近傍での多項式スイッチング」は、計算化学において確立された正統な数値正則化手法であり、本稿で提案する角度ポテンシャルの処理は、この歴史的かつ標準的なアプローチの応用例と位置づけられる。

5. 実装仕様と参照コード#

以下に、提案手法の完全なPython実装を示す。このコードは、PyTorch等の自動微分ライブラリでの使用も想定し、テンソル演算に対応した形式で記述されているが、ロジック自体は言語非依存である。

5.1 Pythonによるリファレンス実装#

import torch

def compute_angle_potential_robust(u, k, theta_0, epsilon=1e-4):
    """
    3点間角度ポテンシャルエネルギーおよび勾配を計算する（特異点除去・数値的正則化版）。
    
    Args:
        u (Tensor): cos(theta)の値。範囲 [-1, 1]。
        k (float): 力の定数。
        theta_0 (float): 平衡角度（ラジアン）。
        epsilon (float): 近似に切り替える閾値。推奨値 1e-4。
        
    Returns:
        E (Tensor): ポテンシャルエネルギー。
        grad (Tensor): uに関する勾配 dE/du。
    """
    
    # 定数の定義 (PyTorch tensorとして確保)
    PI = torch.tensor(3.141592653589793, dtype=u.dtype, device=u.device)
    
    # 多項式係数 (theta^2 approx c1*x + c2*x^2 + c3*x^3)
    coeff_1 = 2.0
    coeff_2 = 1.0 / 3.0
    coeff_3 = 4.0 / 45.0
    
    def poly_val(x):
        """P(x) = 2x + x^2/3 + 4x^3/45"""
        return x * (coeff_1 + x * (coeff_2 + x * coeff_3))
    
    def poly_deriv(x):
        """P'(x) = 2 + 2x/3 + 4x^2/15"""
        return coeff_1 + x * (2.0 * coeff_2 + x * 3.0 * coeff_3)

    # --- 領域判定と計算 ---
    
    # 1. 前方特異点 (theta ~ 0, u ~ 1)
    if u > 1.0 - epsilon:
        delta = 1.0 - u
        theta_sq = poly_val(delta)
        dp_d_delta = poly_deriv(delta)
        
        if theta_0 == 0.0:
            E = 0.5 * k * theta_sq
            dE_du = -0.5 * k * dp_d_delta
        else:
            # theta_sim = sqrt(P(delta))
            theta_sim = torch.sqrt(theta_sq)
            E = 0.5 * k * (theta_sim - theta_0)**2
            
            # dE/du = -0.5 * k * (1 - theta0/theta) * P'
            term = 1.0 - theta_0 / theta_sim
            dE_du = -0.5 * k * term * dp_d_delta

    # 2. 後方特異点 (theta ~ pi, u ~ -1)
    elif u < -1.0 + epsilon:
        delta = 1.0 + u
        # phi = pi - theta, phi^2 approx P(delta)
        phi_sq = poly_val(delta)
        dp_d_delta = poly_deriv(delta)
        
        if theta_0 == PI:
            E = 0.5 * k * phi_sq
            # dE/du = +0.5 * k * P' (符号反転に注意: d(delta)/du = +1)
            dE_du = 0.5 * k * dp_d_delta
        else:
            phi_sim = torch.sqrt(phi_sq)
            theta_sim = PI - phi_sim
            E = 0.5 * k * (theta_sim - theta_0)**2
            
            # dE/du = k(theta - theta0) * d(theta)/du
            # d(theta)/du = -P' / (2*phi)
            term = theta_sim - theta_0
            dE_du = k * term * (-1.0 * dp_d_delta / (2.0 * phi_sim))

    # 3. 通常領域
    else:
        theta = torch.acos(u)
        E = 0.5 * k * (theta - theta_0)**2
        sin_theta = torch.sqrt(1.0 - u**2)
        dE_du = -k * (theta - theta_0) / sin_theta

    return E, dE_du

6. 結論#

デカルト座標系を用いた分子シミュレーションにおける角度ポテンシャルの特異点問題に対し、領域分割と多項式展開を組み合わせた数値的正則化手法を提示した。本手法は以下の特徴を持つ。

整合性: 量子化学におけるBoys関数の処理と同等の数理的アプローチに基づいており、アドホックな修正ではない。

数値的堅牢性: ゼロ除算や桁落ちを原理的に排除し、直線構造を含むあらゆる構造に対して安定である。

高精度: 6次オーダー（$\Delta u$の3次）の近似を用いることで、接続点における $C^1$ 連続性を実効的に確保している。

このアルゴリズムの実装は、自作のMDコードや構造最適化ライブラリの信頼性を向上させるための要件であると考えられる。

参考文献#

• Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory; Dover Publications: Mineola, NY, 1996. (Discussion on Boys Function and numerical integration)
• Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, 1987. (Fundamentals of MD and potential functions)
• Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes: The Art of Scientific Computing, 3rd ed.; Cambridge University Press: New York, 2007. (Numerical regularization techniques)
• Leach, A. R. Molecular Modelling: Principles and Applications, 2nd ed.; Pearson Education: Harlow, 2001.

検証コード：#

import torch
import matplotlib.pyplot as plt
import numpy as np

def verify_continuity():
    # 倍精度での検証が必須
    torch.set_default_dtype(torch.float64)
    
    # 閾値設定 (1e-4 推奨)
    ETA = 1e-4
    k = 100.0  # 力の定数
    
    # u = 1.0 (0度) 付近のデータを生成
    # 閾値をまたぐように細かく刻む
    steps = 1000
    u_vals = torch.linspace(1.0 - 2.0*ETA, 1.0, steps, requires_grad=True)
    
    energies_exact = []
    grads_exact = []
    
    energies_taylor = []
    grads_taylor = []
    
    # 計算ループ
    for u in u_vals:
        delta_u = 1.0 - u
        
        # --- 1. 厳密解 (acos) ---
        if u < 1.0: 
            theta = torch.acos(u)
            E_exact = 0.5 * k * theta**2
            # create_graph=False で勾配計算
            g_exact = torch.autograd.grad(E_exact, u, create_graph=False)[0].item()
        else:
            
            theta = torch.tensor(0.0, dtype=torch.float64)
            E_exact = torch.tensor(0.0, dtype=torch.float64)
            g_exact = 0.0 
            
        energies_exact.append(E_exact.item())
        grads_exact.append(g_exact)
        
        # --- 2. 近似解 (Taylor) ---
        # theta^2 approx 2du + 1/3 du^2 + 4/45 du^3
        theta_sq_approx = 2.0*delta_u + (1.0/3.0)*delta_u**2 + (4.0/45.0)*delta_u**3
        E_taylor = 0.5 * k * theta_sq_approx
        
        # 近似側の勾配計算
        g_taylor = torch.autograd.grad(E_taylor, u, create_graph=False)[0].item()
        
        energies_taylor.append(E_taylor.item())
        grads_taylor.append(g_taylor)

    # --- 3. 境界での不連続性チェック ---
    boundary_u = 1.0 - ETA
    # 最も近いインデックスを取得
    idx = (torch.abs(u_vals - boundary_u)).argmin()
    
    val_exact = energies_exact[idx]
    val_taylor = energies_taylor[idx]
    
    grad_exact_val = grads_exact[idx]
    grad_taylor_val = grads_taylor[idx]
    
    diff_energy = abs(val_exact - val_taylor)
    diff_grad = abs(grad_exact_val - grad_taylor_val)
    
    print(f"=== 検証結果 (Threshold eta = {ETA}) ===")
    print(f"境界点 u ≈ {boundary_u:.8f}")
    print(f"Energy (Exact) : {val_exact:.8e}")
    print(f"Energy (Taylor): {val_taylor:.8e}")
    print(f"Energy差分     : {diff_energy:.2e}")
    print("-" * 30)
    print(f"Grad (Exact)   : {grad_exact_val:.8e}")
    print(f"Grad (Taylor)  : {grad_taylor_val:.8e}")
    print(f"Grad差分       : {diff_grad:.2e}")
    
    # 判定ロジック
    if diff_energy < 1e-15 and diff_grad < 1e-10:
        print(">> 判定: 数値的に連続 (Smooth enough)")
    else:
        print(">> 判定: 不連続性が検出されたため、パラメータの再考が必要")

if __name__ == "__main__":
    verify_continuity()

=== 検証結果 (Threshold eta = 0.0001) ===
境界点 u ≈ 0.99990000
Energy (Exact) : 1.00101770e-02
Energy (Taylor): 1.00101770e-02
Energy差分     : 1.42e-16
------------------------------
Grad (Exact)   : -1.00003337e+02
Grad (Taylor)  : -1.00003337e+02
Grad差分       : 2.71e-12
>> 判定: 数値的に連続 (Smooth enough)

※Energy差分の判定はpythonのsysモジュールから出力される浮動小数点型のマシンイプシロンの大きさの10倍程度である。検証結果ではこのマシンイプシロンよりも小さな値であるため、妥当と判断できるかと考えられる。

※Grad差分の判定がこれで妥当かどうかは検討の余地あり。

Python 3.13.5 | packaged by Anaconda, Inc. | (main, Jun 12 2025, 16:09:02) [GCC 11.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import sys;print(sys.float_info.epsilon)
2.220446049250313e-16

一般化: $\theta_0 \neq 0$ の場合#

1. 問題の所在：特異点の数学的構造#

通常の調和ポテンシャルとその勾配は以下の通りです。

$$V(\theta) = \frac{1}{2} k (\theta - \theta_0)^2$$

勾配（力）をデカルト座標 $\mathbf{r}$ で計算する際、連鎖律により以下の項が現れます。

$$\nabla_{\mathbf{r}} V = \frac{dV}{d\theta} \nabla_{\mathbf{r}} \theta = k(\theta - \theta_0) \cdot \left( -\frac{1}{\sin \theta} \nabla_{\mathbf{r}} (\cos \theta) \right)$$

ここで $\theta \to 0, \pi$ のとき、$\sin \theta \to 0$ となり、項全体が $0/0$ の不定形、または発散を引き起こします。これを回避するため、エネルギー $V$ を $\cos \theta$ の多項式として再定義します。

2. $0^\circ$ 近傍 ($u \to 1$) における導出#

変数 $u$ を以下のように定義します。$u = \cos \theta$ $\theta \approx 0$ におけるマクローリン展開を用います。

Step 1: $u$ の $\theta$ による展開#

$$u = \cos \theta = 1 - \frac{\theta^2}{2!} + \frac{\theta^4}{4!} - \frac{\theta^6}{6!} + O(\theta^8)$$

ここで、$\Delta u = 1 - u$ （$1$ からの偏差）を定義します。

$$\Delta u = \frac{\theta^2}{2} - \frac{\theta^4}{24} + \frac{\theta^6}{720} - \dots \quad \cdots (1)$$

Step 2: 逆関数の展開（係数比較法）#

$\theta^2$ を $\Delta u$ の級数として表現することを考えます。

$$\theta^2 \approx c_1 (\Delta u) + c_2 (\Delta u)^2 + c_3 (\Delta u)^3$$

この式に (1) を代入し、$\theta$ の次数の係数を比較します。

$\theta^2$ の項:

$$c_1 \left( \frac{\theta^2}{2} \right) = \theta^2 \quad \Longrightarrow \quad \frac{c_1}{2} = 1 \quad \therefore c_1 = 2$$

$\theta^4$ の項:式には $c_1$ 由来の項と $c_2$ 由来の項が現れます。

$$c_1 \left( -\frac{\theta^4}{24} \right) + c_2 \left( \frac{\theta^2}{2} \right)^2 = 0$$ $$2 \left( -\frac{1}{24} \right) + \frac{c_2}{4} = 0 \quad \Longrightarrow \quad -\frac{1}{12} + \frac{c_2}{4} = 0 \quad \therefore c_2 = \frac{1}{3}$$

$\theta^6$ の項:

$$c_1 \left( \frac{\theta^6}{720} \right) + c_2 \left( 2 \cdot \frac{\theta^2}{2} \cdot \frac{-\theta^4}{24} \right) + c_3 \left( \frac{\theta^2}{2} \right)^3 = 0$$ $$\frac{2}{720} - \frac{1}{3} \cdot \frac{1}{24} + \frac{c_3}{8} = 0$$ $$\frac{1}{360} - \frac{1}{72} + \frac{c_3}{8} = 0$$

通分して解くと:

$$\frac{1 - 5}{360} + \frac{c_3}{8} = 0 \quad \Longrightarrow \quad -\frac{4}{360} = -\frac{c_3}{8} \quad \therefore c_3 = \frac{4}{45}$$

Step 3: 近似式の確定#

以上より、$0^\circ$ 近傍における $\theta^2$ の近似式が得られます。

$$\theta^2 \approx 2(1-u) + \frac{1}{3}(1-u)^2 + \frac{4}{45}(1-u)^3$$

平衡値 $\theta_0$ に対するポテンシャルエネルギー $V_{near0}$ は、$\theta \approx \sqrt{\theta^2}$ を用いて以下のように記述されます。

$$V_{near0}(u) = \frac{1}{2} k \left( \sqrt{2(1-u) + \frac{1}{3}(1-u)^2 + \frac{4}{45}(1-u)^3} - \theta_0 \right)^2$$

3. $180^\circ$ 近傍 ($u \to -1$) における導出#

直線からの偏差角 $\delta$ を定義します。

$$\delta = \pi - \theta \quad (\delta \approx 0)$$

Step 1: $u$ の $\delta$ による展開#

$$u = \cos(\pi - \delta) = -\cos \delta$$ $$u = -\left( 1 - \frac{\delta^2}{2} + \frac{\delta^4}{24} - \dots \right) = -1 + \frac{\delta^2}{2} - \frac{\delta^4}{24} + \dots$$

ここで、$\Delta u' = 1 + u$ （$-1$ からの偏差）を定義します。

$$\Delta u' = \frac{\delta^2}{2} - \frac{\delta^4}{24} + \dots$$

この形式は $0^\circ$ 近傍の式 (1) と全く同一です。

Step 2: 近似式の確定#

係数比較の結果も同一となるため、$\delta^2 = (\pi - \theta)^2$ の近似式は以下の通りです。

$$(\pi - \theta)^2 \approx 2(1+u) + \frac{1}{3}(1+u)^2 + \frac{4}{45}(1+u)^3$$

平衡値 $\theta_0$ に対するポテンシャルエネルギー $V_{near180}$ は、$\theta = \pi - \delta$ であるため、

$$V_{near180}(u) = \frac{1}{2} k \left( \left[ \pi - \sqrt{2(1+u) + \frac{1}{3}(1+u)^2 + \frac{4}{45}(1+u)^3} \right] - \theta_0 \right)^2$$

4. 勾配の正則性（なぜ発散しないか）の証明#

ポテンシャル $V$ を $u$ の多項式 $P(u)$ として表したとき、勾配は以下のようになります。

$$\nabla_{\mathbf{r}} V = \frac{dV}{du} \nabla_{\mathbf{r}} u$$

$\nabla_{\mathbf{r}} u$ の項:$u = \frac{\mathbf{v}_1 \cdot \mathbf{v}_2}{|\mathbf{v}_1||\mathbf{v}_2|}$ の微分は、ノルムが $0$ にならない限り（原子が重ならない限り）常に有限です。$\frac{dV}{du}$ の項:$0^\circ$ 近傍の主要項（第1項のみ）を考えると、$V \approx \frac{1}{2} k (2(1-u)) = k(1-u)$ となります。

$$\frac{dV}{du} \approx \frac{d}{du} k(1-u) = -k$$

これは定数であり、分母に $0$ が含まれません。高次項を含めても $u$ の多項式であるため、微分結果も $u$ の多項式となり、全領域で有限の値が保証されます。結論以上の導出により、$\arccos$ (およびその微分の $1/\sin \theta$) を介さずとも、$u$ の3次多項式を用いることで角度の6次精度まで厳密に一致するエネルギーと力を算出可能であることが示されました。

一般角 $\theta_0 \neq 0$ における数値的安定化：線形外挿 (Linear Extrapolation) ($\theta_0 \neq 0$かつ$\theta→0$の場合1)#

平衡角 $\theta_0 \neq 0$ の系において、結合角 $\theta$ が $0$ に接近する極限（$u \to 1$）では、ポテンシャル勾配（力）が物理的に発散する。

$$\lim_{\theta \to 0} F(\theta) = \lim_{\theta \to 0} -k(\theta - \theta_0) \frac{1}{\sin \theta} \to \infty$$

この発散は、原子同士の重なりを阻止する幾何学的な障壁として機能するが、数値計算においては Inf (Infinity) やオーバーフローによるシミュレーションの停止を招く。これを回避するため、微小角 $\theta_{cut}$ 以下の領域において、ポテンシャルを接線方向へ線形外挿（Linear Extrapolation）し、力を定数化する処理（Clamping）を導入する。

線形化の定式化閾値#

$\theta_{cut}$（例: $10^{-5}$ rad）を定義し、以下の条件分岐を行う。

物理領域 ($\theta > \theta_{cut}$):通常の調和ポテンシャル（または前述のテイラー展開近似）を使用する。

$$V(\theta) = \frac{1}{2}k(\theta - \theta_0)^2$$ $$F(\theta) = -\frac{dV}{d\theta} = -k(\theta - \theta_0)$$

外挿領域 ($\theta \le \theta_{cut}$):$\theta = \theta_{cut}$ におけるエネルギーと傾き（力）を用いて、一次関数で近似する。

$$V_{lin}(\theta) = V(\theta_{cut}) + \left. \frac{dV}{d\theta} \right|_{\theta_{cut}} (\theta - \theta_{cut})$$

このとき、力は以下の定数となる。

$$F_{lin}(\theta) = -\frac{dV_{lin}}{d\theta} = -\left. \frac{dV}{d\theta} \right|_{\theta_{cut}} = \text{const.}$$

連続性の保証#

この接続手法により、以下の連続性が保証される。

エネルギー ($C^0$ 連続): 境界 $\theta_{cut}$ において値が一致する。

力 ($C^0$ 連続 / エネルギーの $C^1$ 連続): 境界において勾配が一致する。

これにより、ポテンシャル曲面上に不連続な段差（Step）や特異点（Kink）を生じさせることなく、無限大の反発力を「極めて巨大な有限の反発力」へと正則化することが可能となる。これは、MDシミュレーションにおけるエネルギー保存則の破綻を防ぐための必須処理である。

角度ポテンシャルにおける特異点問題とCosine Expansion法による正則化 ($\theta_0 \neq 0$かつ$\theta→0$の場合2)#

※調和ポテンシャルの形を崩すことになるのでこちらは非推奨。

1. 問題の所在：$\theta \to 0$ ($\theta_0 \neq 0$) における力の特異性#

調和ポテンシャル $V(\theta) = \frac{1}{2}k(\theta - \theta_0)^2$ をデカルト座標空間での勾配（力）に変換する際、連鎖律により特異点が発生する可能性があります。力の大きさ $F$ は以下の連鎖律で記述されます（$u = \cos\theta$）。

$$F \propto \left| \frac{dV}{d\theta} \frac{d\theta}{du} \right| = \left| k(\theta - \theta_0) \cdot \frac{-1}{\sin\theta} \right|$$

$\theta \to 0$ の極限において $\sin\theta \approx \theta$ であるため、力の挙動は以下のようになります。

$$F \approx \left| \frac{k(\theta - \theta_0)}{\theta} \right|$$

Case A: $\theta_0 = 0$ （直線分子）の場合

$$\lim_{\theta \to 0} \frac{k\theta}{\theta} = k \quad (\text{有限値})$$

特異点は見かけ上のものであり、数値計算上の注意は必要ですが物理的には破綻しません。

Case B: $\theta_0 \neq 0$ （一般角）の場合

$$\lim_{\theta \to 0} \frac{k(-\theta_0)}{\theta} \to \infty \quad (\text{発散})$$

これは、$\theta_0 \neq 0$ かつ $\theta \to 0$ の場合、ポテンシャル $V$ を $u=\cos\theta$ で記述しようとしてテイラー展開を行っても、その1階微分（力）が無限大となるため展開係数が発散してしまうことを意味します。すなわち、力の発散は物理的特異性ではなく、座標系の特異性です。このため、通常の調和ポテンシャル定義では、構造が直線配置（$\theta=0$）を通過することは数値的に不可能です。

2. 解決策：Cosine Expansion による正則化#

この物理的な特異点および計算上の特異点（$1/\sin\theta$）を排除するため、ポテンシャル関数自体を $u = \cos\theta$ の多項式として再定義します。

$$V_{\text{cos}}(u) = \sum_{n=0}^N C_n u^n$$

これにより、力 $\frac{dV}{du}$ は常に $u$ の多項式（有限値）となり、定義域 $u \in [-1, 1]$ 全体で微分可能（正則）となります。これは「無限のポテンシャル障壁」を「有限のエネルギー障壁」へと物理的に改変し、シミュレーションの安定性を保証する唯一の解決策です。

3. 数学的導出#

プロセス係数 $C_n$ を決定するために、ポテンシャルを $u_0 = \cos\theta_0$ 近傍でテイラー展開します。平衡角の状態に応じて2つのアプローチに分岐します。

3.1. Case 1: $\theta_0 \neq 0, \pi$ （一般角）の場合#

この場合、厳密な合成関数の微分を用いて係数を導出します。定義:

$$A(u) = -k(\arccos u - \theta_0), \quad B(u) = (1-u^2)^{-1/2}$$

1階微分は $V'(u) = A(u)B(u)$ と記述されます。評価点 $u=u_0$ において、$A(u_0)=0$ です。また、導関数は以下の通りです。

$$A'(u) = kB, \quad B'(u) = uB^3$$

各階微分の導出:1階微分 $V'(u_0)$

$$V'(u_0) = A(u_0)B(u_0) = 0 \cdot B(u_0) = 0$$

（$u_0 \neq \pm 1$ より $B(u_0)$ は有限確定値）2階微分 $V''(u_0)$

$$V''(u) = A'B + AB'$$

$u=u_0$ で $A=0$ となる項を消去：$V''(u_0) = A'B = (kB)B = kB^2 = \frac{k}{1-u_0^2}$3階微分 $V'''(u_0)$

$$V'''(u) = (A'B + AB')' = A''B + 2A'B' + AB''$$

$u=u_0$ で $A=0$ 項を消去：

$$V'''(u_0) = A''B + 2A'B'$$

ここで $A'' = (kB)' = k u B^3$ を代入：

$$V'''(u_0) = (kuB^3)B + 2(kB)(uB^3) = 3kuB^4 = \frac{3ku_0}{(1-u_0^2)^2}$$

4階微分 $V''''(u_0)$ライプニッツ則および $A=0$ となる項の消去を行い、残る項のみを計算します。

$$V''''(u_0) = A'''B + 3A''B' + 3A'B''$$

各要素：$A''' = k(B^3 + 3u^2B^5)$

$B'' = B^3 + 3u^2B^5$これらを代入し、$B^6$ で整理すると：

$$V''''(u_0) = k(4+11u_0^2)B^6 = \frac{k(4+11u_0^2)}{(1-u_0^2)^3}$$

展開係数 $K_n$:$V(u) \approx \sum K_n (u-u_0)^n$ としたとき、$K_n = \frac{1}{n!}V^{(n)}(u_0)$ より、

$$K_2 = \frac{k}{2(1-u_0^2)}, \quad K_3 = \frac{ku_0}{2(1-u_0^2)^2}, \quad K_4 = \frac{k(4+11u_0^2)}{24(1-u_0^2)^3}$$

3.2. Case 2: $\theta_0 \approx 0$ （直線分子）の場合#

上記の一般式は、分母に $(1-u_0^2)$ を含むため、$\theta_0 \to 0$ ($u_0 \to 1$) で発散し使用できません。この場合、$\theta \approx \sqrt{2(1-u)}$ の関係を用いた直接展開（マクローリン展開ベース）を使用します。

展開式:

$$V_{\text{linear}}(u) = \frac{1}{2}k \theta^2 \approx \frac{1}{2}k \left[ 2(1-u) + \frac{1}{3}(1-u)^2 + \frac{4}{45}(1-u)^3 \right]$$

これにより、特異点を回避した有限な多項式が得られます。

4. 実装アルゴリズムの提言#

以上の議論より、ロバストな力場実装には以下の条件分岐が必須となります。

入力: 平衡角 $\theta_0$, 力の定数 $k$判定: $|\sin \theta_0| < \epsilon$ (例: $10^{-4}$)

Yes ($\theta_0 \approx 0$) 2 の係数 ($2, 1/3, 4/45$) を用いて多項式 $V(u)$ を構成する。

No ($\theta_0 \neq 0$) 1 の厳密係数 $K_2, K_3, K_4$ ($u_0$ 依存) を計算し、多項式 $V(u)$ を構成する。

余談#

筆者が実装してきたコードに対する反省を兼ねてまとめた。

これまで微小な数値を足すことでごまかしていた数値計算部分を正確にするためにはどうすればいいかを調査しまとめた。