Home
3669 words
18 minutes
【計算化学】PM6法（Parametric Method 6）の数理的背景と実装論：70元素への拡張とNDDO近似の高度化
2025-12-31
Computational Chemistry
Computational Chemistry
/
PM6
/
Semi-empirical
/
NDDO
/
MOPAC
/
Stewart
/
Algorithm

最終更新：2025-12-31

注意: この記事はAIによって自動生成されたものです。正確な情報やパラメータの詳細については、必ず引用元の原著論文（Stewart, J. Mol. Model. 13, 1173 (2007)）をご確認ください。

序論：半経験的手法のルネサンス#

2007年、James J. P. Stewartによって発表された PM6 (Parametric Method 6) 法は、半経験的分子軌道法の歴史において、1989年のPM3法以来、約20年ぶりとなる実質的なメジャーアップデートであった。

1990年代から2000年代にかけて、密度汎関数法（DFT）の台頭により、半経験的手法は「過去の遺物」とみなされる傾向にあった。しかし、数千原子を超えるタンパク質やナノ材料、あるいは結晶周期境界条件（PBC）を含む固体物理の領域において、DFTの計算コスト（O(N3)O(N^3)〜O(N4)O(N^4)）は依然として大きな壁であった。PM6法は、NDDO（Neglect of Diatomic Differential Overlap）近似の枠組みを維持しつつ、遷移金属を含む周期表のほぼ全域（70元素）に対応し、かつDFTに迫る構造予測精度を実現することで、大規模系シミュレーションにおける半経験的手法の復権を果たした。

本稿では、PM6法の理論的基盤、特にPM3からの数理的な変更点である「d軌道の導入に伴う積分計算の拡張」と「コア反発関数の洗練」、そして「パラメータ最適化戦略」について、学術的かつ実装論的な観点から詳細に解説する。

1. 歴史的背景：PM3からPM6への道程#

1.1 Dewar学派とMOPACの系譜#

半経験的手法の歴史は、Michael J. S. DewarらによるMNDO (1977) とAM1 (1985) に端を発する。これらはテキサス大学オースティン校で開発され、プログラムパッケージ「MOPAC」として世界中に普及した。StewartはDewarの研究グループ出身であり、MOPACの開発と保守を主導してきた人物である。

1.2 PM3の成功と限界#

1989年にStewartが発表したPM3 (Parametric Method 3) は、AM1のパラメータを数学的な最適化アルゴリズムによって自動決定したものであり、有機化合物の生成熱予測において高い精度を誇った。しかし、PM3には明確な限界があった。

  1. d軌道の欠如: 主族元素（s, p軌道）のみを前提としていたため、遷移金属を扱えない。後にPM3(tm)として拡張されたが、パラメータは限定的であった。
  2. パラメータの物理的妥当性: 自動最適化に頼りすぎた結果、一部の原子で非物理的な電荷分布を示すことがあった。
  3. 水素結合の記述: AM1より改善されたとはいえ、依然として水二量体などの記述には誤差が含まれていた。

1.3 PM6の開発動機#

PM6の開発における最大の動機は、「周期表のほぼすべての元素（HからBa、HfからBiまで）を単一のハミルトニアンで扱えるようにすること」であった。これを実現するために、Stewartは以下の抜本的な改革を行った。

  • 参照データの拡充: 従来の気相中の生成熱や構造データに加え、X線結晶構造解析のデータ（CSD）を大量に取り入れた。
  • d軌道の完全実装: NDDOの枠組みにd軌道を正式に組み込み、すべての主族元素および遷移金属で利用可能にした。
  • コア反発関数の改良: 原子対（Diatomic）特有の補正項を導入し、系統誤差を排除した。

2. 数理的背景：NDDO近似とd軌道の導入#

PM6法は、AM1やPM3と同様にNDDO近似に基づいている。しかし、d軌道を扱うための数理的な拡張が必要となる。ここでは、Roothaan-Hall方程式の各項がPM6でどのように定式化されているかを記述する。

2.1 基底関数と密度行列#

PM6では、最小基底系（Minimal Basis Set）としてSlater型軌道（STO）を用いる。

  • H, He: 1s1s
  • Li - F: 2s,2p2s, 2p
  • Na - Cl: 3s,3p3s, 3p （※PM6では一部の超原子価原子に対しd軌道の関与を認める）
  • Sc - Zn 等: 4s,4p,3d4s, 4p, 3d

密度行列 P\mathbf{P} は通常通り定義される。 Pμν=2∑ioccCμiCνiP_{\mu\nu} = 2 \sum_{i}^{occ} C_{\mu i} C_{\nu i}

2.2 Fock行列の構成#

Fock行列の要素 FμνF_{\mu\nu} は以下のように表される。

Fμν=Hμν+∑λσPλσ[(μν∣λσ)−12(μλ∣νσ)]F_{\mu\nu} = H_{\mu\nu} + \sum_{\lambda\sigma} P_{\lambda\sigma} [(\mu\nu|\lambda\sigma) - \frac{1}{2}(\mu\lambda|\nu\sigma)]

NDDO近似により、二電子積分 (μν∣λσ)(\mu\nu|\lambda\sigma) は、μ,ν\mu, \nu が同一原子 AA 上、λ,σ\lambda, \sigma が同一原子 BB 上にある場合のみ非ゼロとなる。

2.3 d軌道を含む二中心二電子積分#

PM6の実装において最も複雑なのが、d軌道を含む二中心積分の計算である。s, p軌道のみの場合、回転操作は比較的単純だが、d軌道（5成分: dz2,dxz,dyz,dx2−y2,dxyd_{z^2}, d_{xz}, d_{yz}, d_{x^2-y^2}, d_{xy}）が含まれると、多重極展開の項数が劇的に増加する。

PM6では、二中心積分を以下の手順で計算する。

  1. 局所座標系への変換: 原子対 A−BA-B を結ぶ軸をローカルな zz 軸とする座標系へ、基底関数を回転させる。回転行列 R\mathbf{R} は、オイラー角あるいは方向余弦を用いて構築される。d軌道の変換行列は 5×55 \times 5 となる。

  2. 局所積分の計算: 局所系では、対称性により多くの積分がゼロになる。非ゼロの積分は、多重極展開（単極子、双極子、四重極子…六重極子）を用いて近似計算される。 PM6における相互作用ポテンシャルは、点電荷近似の改良版（Klopman-Ohno式など）を用いる。 (μν∣λσ)local≈1RAB2+(δμν+δλσ)2(\mu\nu|\lambda\sigma)_{local} \approx \frac{1}{\sqrt{R_{AB}^2 + (\delta_{\mu\nu} + \delta_{\lambda\sigma})^2}} ここで δ\delta は、軌道分布の大きさを表すパラメータである。d軌道が含まれる場合、高次の多重極項（例えば dd-dd 相互作用）の計算が必要となる。

  3. グローバル座標への逆変換: 計算された局所積分を、回転行列を用いてグローバル座標へ戻す。 (μν∣λσ)global=∑μ′ν′λ′σ′Rμμ′Rνν′Rλλ′Rσσ′(μ′ν′∣λ′σ′)local(\mu\nu|\lambda\sigma)_{global} = \sum_{\mu'\nu'\lambda'\sigma'} R_{\mu\mu'} R_{\nu\nu'} R_{\lambda\lambda'} R_{\sigma\sigma'} (\mu'\nu'|\lambda'\sigma')_{local}

この計算コストを削減するため、PM6の実装（MOPACなど）では、頻出する回転パターンをあらかじめテーブル化するなどの最適化が行われている。

2.4 一中心一電子エネルギー UμμU_{\mu\mu}#

PM6における UμμU_{\mu\mu} （原子価状態イオン化エネルギーに対応）は、実験値から直接取るのではなく、パラメータ最適化の対象となっている。これにより、原子内の電子環境の変化（酸化状態など）に対して柔軟に対応できる。

3. PM6の核心：改良されたコア反発関数 (CRF)#

AM1/PM3からの最大の変更点は、コア反発関数（Core Repulsion Function）の設計にある。AM1ではガウス関数を用いたが、PM6ではこれを整理し、かつ特定の原子対に対する「二原子パラメータ（Diatomic Parameters）」を導入した。

3.1 基本的な反発項#

原子 A,BA, B 間のコア反発エネルギー EcoreABE_{core}^{AB} の基本形は、MNDO/AM1を踏襲している。

EcoreAB=ZAZB(μν∣λσ)ss(1+e−αARAB+e−αBRAB)E_{core}^{AB} = Z_A Z_B (\mu\nu|\lambda\sigma)_{ss} (1 + e^{-\alpha_A R_{AB}} + e^{-\alpha_B R_{AB}})

ここで (μν∣λσ)ss(\mu\nu|\lambda\sigma)_{ss} はs軌道間の二電子積分 γAB\gamma_{AB} である。

3.2 二原子補正項 (Specific Pairwise Parameters)#

PM6では、特定の結合（例えば C-H, O-H, N-H, C-C など重要かつ誤差が出やすい結合）に対して、個別の補正項 EijaddE_{ij}^{add} を追加している。これがPM6の精度の高さを支えている。

EcoreAB,total=EcoreAB,base+∑kakexp⁡(−bk(RAB−ck)2)E_{core}^{AB, total} = E_{core}^{AB, base} + \sum_{k} a_{k} \exp(-b_{k}(R_{AB} - c_{k})^2)

この形式自体はAM1と同じだが、PM6ではパラメータの決定方法が異なる。AM1では元素ごとにガウス関数を持たせていたが、PM6では必要に応じて「特定の元素ペア (A,B)(A, B) に対してのみ」作用するガウス関数を定義している。 例えば、Si-O 結合や P-O 結合など、超原子価や固体化学で重要な結合に対して、CSDの結晶構造データを再現するように個別の補正パラメータが与えられている。

3.3 Voityukの近似について#

論文によっては「Voityukの近似」への言及があるが、これは主に二中心積分の計算効率化に関するものであり、PM6の本質的な物理モデルの変更ではない。PM6の実装においては、厳密なNDDO回転不変性を維持しつつ、d軌道積分を高速化するための数学的テクニックが駆使されている。

4. プログラム実装のためのアルゴリズム詳細#

PM6法を実装する上で、AM1/PM3の実装から拡張すべき点は「d軌道の回転行列」と「ペアワイズなコア反発補正」である。以下に擬似コードを示す。

4.1 パラメータ構造体#

PM6は元素ごとのパラメータに加え、特定の原子ペアに対するパラメータを持つ。

class PM6AtomParams:
    def __init__(self, element):
        self.U_ss, self.U_pp, self.U_dd = ... # 一電子エネルギー
        self.beta_s, self.beta_p, self.beta_d = ... # 共鳴パラメータ
        self.zeta_s, self.zeta_p, self.zeta_d = ... # 軌道指数
        self.alpha = ... # コア反発指数の基本項
        # その他のMNDOパラメータ

class PM6PairParams:
    def __init__(self):
        # 特定のペア (A, B) に対するガウス補正パラメータのリスト
        # [(a1, b1, c1), (a2, b2, c2), ...]
        self.gaussians = []

4.2 コア反発エネルギーの実装#

PM6固有のペアパラメータを処理するロジックが必要である。

import numpy as np

def calculate_pm6_core_repulsion(atom_A, atom_B, dist, pair_params_db):
    """
    PM6法におけるコア反発エネルギーを計算する。
    
    Parameters:
    atom_A, atom_B: 原子パラメータ
    dist: 核間距離 (Angstrom)
    pair_params_db: 原子ペアごとの補正パラメータデータベース {(elemA, elemB): [params]}
    """
    
    # 1. 基本項 (MNDO/AM1 like)
    # gamma_AB は s-s 二電子積分近似値
    gamma_AB = calculate_gamma_ss(atom_A, atom_B, dist)
    
    # 遮蔽項
    term_A = np.exp(-atom_A.alpha * dist)
    term_B = np.exp(-atom_B.alpha * dist)
    
    E_base = atom_A.core_charge * atom_B.core_charge * gamma_AB
    E_base *= (1.0 + term_A + term_B)
    
    # 2. PM6固有のペアワイズ補正 (Gaussian terms)
    # ペアキーの生成 (辞書検索用: アルファベット順などで正規化)
    pair_key = tuple(sorted((atom_A.element, atom_B.element)))
    
    E_correction = 0.0
    if pair_key in pair_params_db:
        gaussians = pair_params_db[pair_key]
        scaling = (atom_A.core_charge * atom_B.core_charge) / dist
        
        for (a, b, c) in gaussians:
            # PM6の論文 Eq.19 に基づく形式
            # a = 10^A_ij, b = B_ij, c = C_ij とする場合もあるが
            # ここでは一般的なガウス形式 a * exp(-b*(R-c)^2) とする
            E_correction += a * np.exp(-b * (dist - c)**2)
            
        E_correction *= scaling

    return E_base + E_correction

4.3 Fock行列構築（d軌道への拡張）#

AM1のコードに対し、d軌道ブロックのループを追加する必要がある。

def build_pm6_fock_matrix(mol, P, integrals):
    F = np.zeros_like(P)
    
    for A in mol.atoms:
        for B in mol.atoms:
            # インデックス範囲 (s, px, py, pz, dz2, dxz, dyz, dx2y2, dxy)
            idx_A = get_orbital_indices(A)
            idx_B = get_orbital_indices(B)
            
            if A == B:
                # 一中心項 (AM1と同様だがd軌道を含む)
                # U_dd や (dd|dd), (dd|sp) などの反発項を加算
                F[idx_A, idx_A] += calculate_one_center_F(A, P[idx_A, idx_A])
                
                # 他核ポテンシャル
                # V_B = Z_B * (mu nu | s_B s_B)
                # d軌道に対しても計算
            else:
                # 二中心項 (NDDO)
                # 1. 共鳴積分 H_mu_lam
                # S_mu_lam (重なり積分) に比例。
                # d軌道を含む重なり積分の計算が必要 (STOの解析解)
                S_sub = integrals.Overlap[idx_A, idx_B]
                Beta_sub = 0.5 * (A.beta + B.beta) * S_sub 
                
                # 2. 交換積分項 -0.5 * sum P * (mu nu | lam sig)
                # d軌道を含む回転行列を用いた二電子積分計算
                G_exch = calculate_two_center_exchange(A, B, P[idx_A, idx_B])
                
                F[idx_A, idx_B] += Beta_sub + G_exch
                
    return F

この calculate_two_center_exchange 内部で、d軌道を含む 5×55 \times 5 以上の回転行列演算が行われる。

5. 実利的な成果とバリデーション#

PM6法は、大規模なパラメータ最適化（数千の分子データを使用）により、劇的な精度向上を達成した。

5.1 生成熱の精度#

有機化合物（H, C, N, O, F, P, S, Cl, Brを含む系）において、生成熱の平均絶対誤差（AUE）は以下の通り改善された。

AM1: ~10.0 kcal/mol

PM3: ~6.3 kcal/mol

PM6: ~4.4 kcal/mol

特に、PM3で問題となっていた窒素化合物（ニトロ基など）や超原子価リン化合物の誤差が大幅に低減された。

5.2 結晶構造の再現性#

PM6のパラメータ決定には、気相データだけでなく、ケンブリッジ結晶構造データベース（CSD）からの固相データも用いられた。これにより、分子結晶やイオン結晶の格子定数や密度を、半経験的手法としては異例の高い精度で再現可能となった。これは、周期境界条件（PBC）を用いた固体計算においてPM6が重用される最大の理由である。

5.3 生体高分子への適用#

水素結合の記述が改善された（水二量体の結合エネルギー誤差が < 1 kcal/mol）ことに加え、計算速度が維持されているため、タンパク質の構造最適化や、QM/MM計算におけるQM領域のソルバーとして標準的な地位を確立した。特に、酵素反応中心に金属（Fe, Zn, Mg等）が含まれる場合、d軌道を扱えないPM3では不可能だった解析がPM6で可能となった。

6. 結論#

PM6法は、NDDO近似という「枯れた技術」に対し、d軌道の完全導入と、現代的なビッグデータを用いたパラメータ最適化、そして結晶構造データの活用という「新しい息吹」を吹き込むことで誕生した手法である。その数理的本質は、複雑な量子多体問題を、回転不変性を保った多重極相互作用と、巧みに調整されたペアワイズなコア反発関数へと落とし込む点にある。現在、PM6はさらに改良されたPM7法（非共有結合相互作用や分散力補正を強化）へと進化しているが、PM6が切り開いた「全元素対応・高精度・高速」というパラダイムは、計算化学の民主化（誰もが手軽に大規模計算を行える環境）に多大な貢献を果たしたと言える。

参考文献#

  • J. J. P. Stewart, “Optimization of parameters for semiempirical methods V: Modification of NDDO approximations and application to 70 elements,” J. Mol. Model. 13, 1173 (2007).
  • J. J. P. Stewart, “Application of the PM6 method to modeling proteins,” J. Mol. Model. 15, 765 (2009).
  • J. J. P. Stewart, MOPAC2012 Manual, Stewart Computational Chemistry, Colorado Springs, CO.
【計算化学】PM6法（Parametric Method 6）の数理的背景と実装論：70元素への拡張とNDDO近似の高度化
https://ss0832.github.io/posts/20251231_compchem_pm6_overview/
Author
ss0832
Published at
2025-12-31
【計算化学】NDDO近似（Neglect of Diatomic Differential Overlap）の数理的詳細と実装論：現代半経験的分子軌道法の基盤
【計算化学】Pople系基底関数系の数理と体系：Gaussian型軌道の縮約と分極・分散関数の役割