部分ヘシアン振動解析法 (Partial Hessian Vibrational Analysis): 大規模系における振動分光の理論的枠組みと応用

最終更新：2026-01-02

注意: この記事はGemini 3.0によって自動生成されたものです。内容はJohn D. Headによる “Computation of Vibrational Frequencies for Adsorbates on Surfaces” (Int. J. Quant. Chem. 1997) および Nicholas A. Besleyらによる一連の研究論文に基づいています。正確な学術的情報は原著論文を参照してください。

序論：大規模分子系における振動解析の課題と部分ヘシアン法の位置づけ#

振動分光法（赤外分光法や高分解能電子エネルギー損失分光法：HREELSなど）は、物質の微視的な構造や化学結合の状態を同定するための強力な実験手法である。実験スペクトルの帰属や解釈において、量子化学計算による基準振動解析は不可欠な役割を果たしている。通常、調和近似（Harmonic Approximation）に基づく振動解析では、ポテンシャルエネルギー曲面の二階微分行列であるヘシアン（Hessian）行列を計算し、質量重み付き座標系において対角化することで、固有振動数と基準モードを決定する。

しかし、ヘシアン行列の計算コストは原子数 $N$ に対して $3N \times 3N$ の次元を持ち、計算量は系のサイズと共に急激に増大する。特に、固体表面上の吸着種やタンパク質のような大規模系においては、系全体のフルヘシアン（Full Hessian）を計算することは、計算資源の観点から極めて困難、あるいは非現実的である場合が多い。

この課題に対し、系の特定の部分（吸着種や特定の官能基など）に注目し、その部分空間に関連するヘシアン要素のみを計算・対角化する「部分ヘシアン法（Partial Hessian Vibrational Analysis）」が提案・発展してきた。本稿では、Head (1997) による表面吸着系への適用と収束性の議論、および Besley (2007, 2008) による生体高分子および大規模表面クラスターへの展開に基づき、本手法の数理的基礎と実用性を概観する。

1. 部分ヘシアン法の数理的定式化#

部分ヘシアン法は、全ヘシアン行列のうち、興味ある部分系（Adsorbate/Fragment）に対応するブロックのみを取り出し、残りの部分（Substrate/Environment）との結合を近似的に扱う手法である。

1.1 ヘシアン行列の分割#

$N$ 個の原子からなる全系において、振動解析の対象とする原子群を $A$ (Adsorbate/Active)、それ以外の原子群を $S$ (Substrate/Spectator) と定義する。このとき、質量重み付き力定数行列（ヘシアン） $K$ は以下のようなブロック行列として表現される。

K = \begin{pmatrix} K_{AA} & K_{AS} \\ K_{SA} & K_{SS} \end{pmatrix}

ここで、

$K_{AA}$ : 対象とする部分系内の原子座標に関する二階微分（部分ヘシアンブロック）。
$K_{SS}$ : 環境側の原子座標に関する二階微分。
$K_{AS}, K_{SA}$ : 部分系と環境間のカップリング（相互作用）を表すオフダイアゴナル項。

各要素 $K_{ij}$ は、原子質量 $m_i, m_j$ およびポテンシャルエネルギー $E$ を用いて以下のように定義される。

K_{ij} = \frac{1}{\sqrt{m_i m_j}} \frac{\partial^2 E}{\partial q_i \partial q_j}

1.2 部分ヘシアン近似 (The Partial Hessian Approximation)#

最も基本的な部分ヘシアン近似（PH近似）は、環境側の原子質量を無限大（infinite mass）と仮定することに物理的に対応する。これは数学的には、全ヘシアン行列 $K$ の対角化を行わず、部分行列 $K_{AA}$ のみを対角化して固有値を求める操作と等価である。

K_{AA} L_A = L_A \Lambda_A

ここで、 $\Lambda_A$ は部分系の近似的な固有振動数の二乗を対角成分に持つ行列、 $L_A$ は対応する固有ベクトル（基準モード）である。この近似が妥当であるためには、部分系の振動モードが環境側の振動モードと十分に分離している（デカップリングしている）必要がある。具体的には以下の条件で精度が良いとされる。

質量差: 部分系（例：水素原子や軽元素分子）の原子質量が環境（例：金属表面）よりも著しく小さい場合。
局所性: 対象とする振動モード（例：C=O伸縮）が特定の官能基に強く局在している場合。

1.3 Headによる振動数補正公式と収束性診断#

Head (1997) は、単に $K_{AA}$ を対角化するだけでなく、無視されたカップリング項 $K_{AS}$ の影響を見積もるための補正公式を導入した。これはBlackによって提案された摂動論的アプローチに基づくものである。部分ヘシアンの対角化によって得られた $i$ 番目の固有振動数を $\omega_i$ とすると、環境との相互作用を考慮した補正振動数 $\omega_i^{cor}$ は以下のように推定される。

(\omega_i^{cor})^2 = \omega_i^2 + \sum_{j \in S} \frac{(\tilde{K}_{AS})_{ij}^2}{\omega_i^2 - \Omega_j^2}

ここで、 $\Omega_j$ は環境側（基板など）のフォノン振動数である。もし、対象とする吸着種の振動数 $\omega_i$ が基板の最大フォノン振動数 $\Omega_{max}$ よりも十分に大きい場合（ $\omega_i \gg \Omega_{max}$ ）、上記の式は以下のように近似できる。

(\omega_i^{cor})^2 \approx \omega_i^2 + \frac{1}{\omega_i^2} \sum_{j \in S} (\tilde{K}_{AS})_{ij}^2

この式は、実際に $K_{SS}$ を計算・対角化することなく、エネルギー勾配（gradient）の計算過程で得られる $K_{AS}$ 要素を用いて、PH近似の妥当性を評価する指標を与える。補正項 $\Delta = \omega_i^{cor} - \omega_i$ が十分に小さければ、その振動モードは $K_{AA}$ ブロック内で収束していると判断できる。逆に補正項が大きい場合は、 $A$ の領域（バッファー領域）を拡大し、より多くの基板原子を $K_{AA}$ に含める必要があることを示唆する。

1.4 実装における計算コスト削減 (Besleyらのアプローチ)#

BesleyとMetcalf (2007) は、Q-CHEMソフトウェアパッケージ内において、Coupled Perturbed Hartree-Fock (CPHF) 方程式の解法プロセスにPH近似を組み込むことで、計算時間を大幅に短縮した。解析的振動数計算において、最もコストがかかるのは以下の2点である。

密度行列の微分の計算（CPHF方程式の反復解法）。
電子反発積分の二階微分の評価。

PH近似を用いると、摂動の対象が全原子数 $3N$ ではなく、部分系の自由度 $3N_{sub}$ に限定されるため、CPHFの求解回数および積分微分の評価回数が劇的に減少する。これにより、メモリ使用量と計算時間の両面で大幅な節約が可能となる。

2. 表面吸着系への適用事例と考察#

2.1 金属表面上の原子状吸着種 (H/Al, O/Al)#

Head (1997) は、Al(100)およびAl(111)表面上の水素(H)および酸素(O)原子吸着に対し、クラスターモデルを用いたPH法の検証を行った。

H/Al(100): H原子はAl原子に比べて質量が非常に小さいため、PH近似（H原子のみをヘシアンに含める）が極めて良好に機能する。補正公式によるシフトも僅少であり、実験値（HREELS）との良い一致を示した。
O/Al(111): O原子はH原子に比べて重く、Alとの質量比が大きくなるため、O原子のみのヘシアン（1原子ブロック）では誤差が大きいことが示された。O原子とその最近接Al原子3個を含む $Al_3O$ サブクラスターを $A$ として定義することで、振動数は劇的に改善し、補正項も小さく収束した。これは、振動モードの記述において、吸着種だけでなく、結合に関与する基板原子を明示的に含めることの重要性を示している。

2.2 半導体表面上の有機分子 (Si(100)表面)#

BesleyとBryan (2008) は、Si(100)表面上に吸着したアセチレン、エチレン、ベンゼン、1,3-ブタジエン、ナフタレンなどの有機分子に対し、PH法を適用した。

計算モデル: $Si_9H_{12}$ （1ダイマー）から $Si_{60}H_{44}$ （4ダイマー）までの大規模クラスターを使用。密度汎関数理論 (DFT/EDF1) を採用。
混合基底系: 計算コストをさらに下げるため、吸着種と結合サイトのSiには 6-31G* を、それ以外のSiには 6-31G を用いる混合基底アプローチを併用した。
結果の精度: フルヘシアン計算と比較して、PH法による振動数の誤差は典型的には $1 \sim 2$ cm $^{-1}$ 程度（C-H伸縮などの高波数領域）であり、化学的な解釈に影響を与えないレベルであることが確認された。低波数のSi-C伸縮やSi-Siモードでは多少の誤差が生じるが、定性的な傾向は再現された。
大規模吸着種の解析: ナフタレンのような大きな分子が複数のダイマー列にまたがって吸着する構造は、従来の小規模クラスターでは扱えなかった。PH法によるコスト削減により、4ダイマーを含む $Si_{60}H_{44}$ クラスターでの計算が可能となり、被覆率依存性や、隣接分子との相互作用（同じダイマー列内での吸着 vs 隣接列での吸着）によるスペクトル変化を議論できるようになった。特に1,3-ブタジエンの吸着において、[4+2]付加体のcis/trans異性体の存在が実験スペクトルの解釈に不可欠であることを理論的に示した。

3. 生体高分子（タンパク質）への展開#

タンパク質の二次構造解析において、アミドIバンド（主にC=O伸縮振動、1600-1700 cm $^{-1}$ 付近）は最も重要な指標の一つである。しかし、タンパク質全体の振動解析は巨大すぎて通常行えない。BesleyとMetcalf (2007) は、この問題に対しPH法を適用した。

3.1 アミドIバンドの局所性利用#

アミドIモードは、ペプチド結合のカルボニル基（C=O）に強く局在している。Besleyらは、以下の2種類のPHモデルを比較検討した。

C/O モデル: バックボーンのカルボニル炭素と酸素のみをヘシアンに含める。
C/O/N モデル: 上記に加え、アミド窒素もヘシアンに含める。

3.2 精度と効率のトレードオフ#

グリシン、アラニン、バリンなどのオリゴペプチド（5残基程度）を用いた検証の結果、以下の知見が得られた。

振動数: フルヘシアンと比較した平均絶対誤差（MAD）は、HF/STO-3Gレベルにおいて、C/Oモデルで約 23.5 cm $^{-1}$ 、C/O/Nモデルで約 14.5 cm $^{-1}$ であった。より高精度の DFT (EDF1/6-31G*) レベルでは、C/O/Nモデルの誤差は 10 cm $^{-1}$ 未満となった。
強度: 赤外強度は振動数よりもPH近似に敏感であり、誤差が大きくなる傾向が見られたが、スペクトルのバンドプロファイル（形状）はフルヘシアンの結果をよく再現した。
計算時間: トリプトファン5残基の系において、PH法はフルヘシアン計算の 10〜15% の時間で完了した。計算時間の短縮率は系が大きくなるほど向上する（ $O(N^3)$ vs $O(N \cdot n_{sub}^2)$ のようなスケーリングの違いに起因）。

3.3 Agitoxin 2 への適用#

Besleyらは、この手法をサソリ毒タンパク質である Agitoxin 2 (PDB: 1AGT) の二次構造要素（ $\alpha$ ヘリックス、 $\beta$ シート、ターン）に適用した。 $\beta$ シート部分は190原子を含むが、PH法を用いることでHF/6-31Gレベルでの第一原理振動解析が可能となった。計算されたスペクトルは、各二次構造に特徴的なピーク位置（ $\alpha$ ヘリックスの 1650 cm $^{-1}$ 、 $\beta$ シートの 1610/1630 cm $^{-1}$ ）を再現し、実験的に知られている経験則と一致した。これは、経験的な力場パラメータに頼らず、第一原理計算によってタンパク質の赤外スペクトルを直接予測できる可能性を示している。

4. 考察：数値微分法（Energy Fitting）との比較#

Head (1997) は、従来行われていた「エネルギー曲線の数値フィット（Energy Fitting）」による振動数算出と、本稿の「ヘシアン対角化」による手法を比較し、後者の優位性を主張している。

従来の数値フィット法: 特定の原子を予想される基準モード方向へ変位させ、そのエネルギー変化を多項式近似して振動数を求める。
- 欠点1: 基板を無限質量として扱うため、吸着種と基板の質量比が大きい場合（例：酸素原子）に誤差が生じる。
- 欠点2: 多原子吸着種の場合、真の基準モード（Normal Mode）の形状をあらかじめ知ることは困難であり、仮定したモードが不正確であれば振動数も不正確になる。
部分ヘシアン対角化法:
- 利点1: 部分系ブロックの対角化により、モードの混合（Mode Mixing）を含めた正しい基準モードが自動的に得られる。
- 利点2: 前述の補正公式を用いることで、基板とのカップリングの強さを定量的に評価でき、クラスターサイズの収束判定が容易になる。

結論#

部分ヘシアン振動解析法は、計算資源の制約が厳しい大規模分子系や表面界面系において、完全なフルヘシアン計算の代替となりうる効率的かつ実用的なアプローチである。

計算効率: 計算時間とメモリを大幅に削減し、従来不可能であったサイズの系（大規模クラスター上の吸着種、タンパク質ドメインなど）に対する第一原理振動解析を可能にする。
精度: 適切な原子群（バッファー領域を含む）をヘシアン計算の対象に含めることで、実験の解釈に十分な精度（数 cm $^{-1}$ 〜 10数 cm $^{-1}$ の誤差）を達成できる。特に、局在化した高波数モード（C-H伸縮、C=O伸縮など）に対して有効性が高い。
診断機能: Headによって示された補正公式は、計算結果が環境からの影響をどの程度受けているか、あるいは選択した部分系が適切であるかを判断する診断ツールとして機能する。

今後の展望として、並列計算機環境における効率的な実装や、溶媒効果を取り入れたQM/MM法との融合などが進むことで、より複雑な生体・材料系への適用が期待される。

付録: 部分ヘシアン対角化の概念的Pythonコード#

以下に、全ヘシアン行列から部分ブロックを抽出し、対角化して振動数を求める手順の概念的な実装を示す。なお、実際の量子化学計算では、全ヘシアンを生成せずに部分ヘシアンを直接計算するアルゴリズムが用いられることが多い。

import numpy as np
from scipy import linalg

class PartialHessianAnalysis:
    def __init__(self, full_hessian, masses, active_indices):
        """
        部分ヘシアン解析を行うクラス

        Parameters:
        -----------
        full_hessian : np.ndarray
            力定数行列 (3N x 3N, 単位: Hartree/Bohr^2 等)
        masses : np.ndarray
            全原子の質量 (N,)
        active_indices : list of int
            解析対象とする原子のインデックス (0-based)
        """
        self.H = full_hessian
        self.masses = masses
        self.active_idx = active_indices
        self.n_active = len(active_indices)
    
    def compute_frequencies(self):
        """
        部分ヘシアン行列を構築・対角化し、振動数を計算する
        
        Returns:
        --------
        freqs : np.ndarray
            振動数 (cm^-1)
        modes : np.ndarray
            質量重み付き基準モード
        """
        # 1. アクティブな原子の自由度に対応するインデックスを取得
        # (原子iに対して、3i, 3i+1, 3i+2 がx,y,z座標)
        indices = []
        for i in self.active_idx:
            indices.extend([3*i, 3*i+1, 3*i+2])
        
        # 2. 部分ヘシアンブロック K_AA の抽出
        H_sub = self.H[np.ix_(indices, indices)]
        
        # 3. 質量重み付け
        # M_sub[k] = sqrt(mass of atom corresponding to coordinate k)
        m_sub_list = []
        for i in self.active_idx:
            m_root = np.sqrt(self.masses[i])
            m_sub_list.extend([m_root, m_root, m_root])
        
        M_inv = np.diag(1.0 / np.array(m_sub_list))
        
        # 質量重み付きヘシアン K = M^(-1/2) * H_sub * M^(-1/2)
        K_AA = M_inv @ H_sub @ M_inv
        
		# 3.5. 併進と回転成分を取り除く
		K_AA = project_out_rot_and_tr(K_AA)
		
        # 4. 対角化
        eigvals, eigvecs = linalg.eigh(K_AA)
        
        # 5. 固有値から振動数への変換 (単位変換係数は系に依存)
        # 負の固有値は虚振動に対応
        freqs = []
        conversion_factor = 5140.48 # 例: Hartree/Bohr^2/AMU -> cm^-1
        
        for val in eigvals:
            sign = 1.0 if val >= 0 else -1.0
            freq = np.sqrt(abs(val)) * conversion_factor
            if sign < 0:
                freq = -freq # 虚振動は負の値で表現
            freqs.append(freq)
            
        return np.array(freqs), eigvecs

    def head_correction_diagnostic(self, sub_freqs, env_indices, coupling_block_K_AS):
        """
        Head (1997) の式 (5) に基づく補正項の推定 (簡易版)
        
        Parameters:
        -----------
        sub_freqs : np.ndarray
            K_AAから得られた振動数 omega_i
        env_indices : list
            環境側(S)のインデックス
        coupling_block_K_AS : np.ndarray
            質量重み付きオフダイアゴナル項 K_AS
            
        Returns:
        --------
        corrected_freqs_sq : np.ndarray
            補正された振動数の二乗 (推定値)
        """
        # omega_i >> Omega_j (環境側のフォノン) と仮定
        # 式: (omega_cor)^2 = omega^2 + sum(K_AS_ij^2) / omega^2
        
        n_modes = len(sub_freqs)
        corrections = np.zeros(n_modes)
        
        for i in range(n_modes):
            omega_sq = sub_freqs[i]**2
            if abs(omega_sq) < 1e-8: continue # 低波数モードはスキップ
            
            # i番目のモードに対するカップリングの二乗和
            # (固有ベクトル空間へ射影する必要があるが、ここでは概念的に記述)
            sum_coupling = np.sum(coupling_block_K_AS[i, :]**2)
            
            corrections[i] = sum_coupling / omega_sq
            
        return sub_freqs**2 + corrections

# 使用例のイメージ
# full_hess = ... (量子化学計算から取得)
# masses = ...
# analysis = PartialHessianAnalysis(full_hess, masses, active_indices=[0, 1, 2])
# freqs, modes = analysis.compute_frequencies()

参考文献#

[1] John D. Head, “Computation of Vibrational Frequencies for Adsorbates on Surfaces”, International Journal of Quantum Chemistry, Vol. 65, 827-838 (1997).

[2] Nicholas A. Besley and Katie A. Metcalf, “Computation of the amide I band of polypeptides and proteins using a partial Hessian approach”, The Journal of Chemical Physics, 126, 035101 (2007).

[3] Nicholas A. Besley and James A. Bryan, “Partial Hessian Vibrational Analysis of Organic Molecules Adsorbed on Si(100)”, The Journal of Physical Chemistry C, 112, 4308-4314 (2008).