内部座標多重構造近似による低周波数モード結合とねじれ非調和性の熱化学的処理：理論的枠組みと実装

7023 words

35 minutes

内部座標多重構造近似による低周波数モード結合とねじれ非調和性の熱化学的処理：理論的枠組みと実装

2026-01-21

Physical Chemistry

Computational Chemistry

Thermochemistry

Torsional Anharmonicity

Mode Coupling

Partition Function

Voronoi Tessellation

last_modified: 2026-01-21

生成AIによる自動生成記事に関する免責事項: 本記事は、論文 “Practical methods for including torsional anharmonicity in thermochemical calculations on complex molecules: The internal-coordinate multi-structural approximation” (Phys. Chem. Chem. Phys., 2011, 13, 10885-10907) の内容に基づき、その数理的背景、アルゴリズムの実装詳細、および実証結果を技術仕様書的な観点から記述したものです。記述の正確性には万全を期していますが、実計算への適用にあたっては必ず原著および関連する量子化学計算パッケージのマニュアルを参照してください。

1. 序論：調和近似の限界とモード結合の問題#

気相における複雑な分子系の熱力学的性質（エントロピー $S$ 、熱容量 $C_p$ 等）の算出において、調和振動子（Harmonic Oscillator, HO）近似は計算コストの観点から標準的に用いられている。しかし、ねじれ（Torsion）運動のような低周波数モードは、ポテンシャルエネルギー曲面（PES）上で非調和性が著しく、かつ複数の極小点（配座）を持つため、単一構造に基づくHO近似（SS-HO）はしばしば定性的に誤った結果を与える。

特に深刻な問題は、低周波数モード間の結合（Coupling） である。多くの複雑な分子において、ねじれモードは他のねじれモードや低周波数の変角（Bending）モードと強く混成する。この結果、特定の法線モード（Normal Mode）を特定の内部回転に対応させることが困難となり、Pitzer-Gwinn近似のような1次元の分離モデルを適用する際の障害となる。

Zheng, Truhlar ら（2011）によって提案された「多重構造（Multi-Structural, MS）法」は、法線座標ではなく内部座標（Internal Coordinates） を用いることでこの問題を回避する。本稿では、特にモード間の結合がいかにして理論に取り込まれているかに焦点を当て、その数理的枠組みと実装を解説する。

2. 理論的枠組みと数理的背景#

MS法は、複数の配座（Conformer）からの寄与を包括的に扱い、かつ各配座におけるねじれ非調和性とモード間結合を補正する。

2.1 全分配関数の一般形式#

電子、並進を除く、配座-回転振動（Conformational-Rovibrational）分配関数 $Q_{con-rovib}^{MS-AS}$ は、区別可能な全構造 $J$ にわたる総和として定義される。

Q_{con-rovib}^{MS-AS} = \sum_{j=1}^{J} Q_{rot,j} \exp(-\beta U_j) Q_{j}^{HO} Z_j \prod_{\tau=1}^{t} f_{j,\tau}

ここで、 $Q_{rot,j}$ は回転分配関数、 $U_j$ は相対エネルギー、 $Q_{j}^{HO}$ は通常の法線モード調和振動子分配関数である。本理論の核心は、内部座標ねじれ非調和性係数 $f_{j,\tau}$ と、モード結合および高温極限補正係数 $Z_j$ にある。

コラム：なぜ「法線モード」ではダメなのか？ —— 混成の問題#

「法線モード（Normal Mode）解析ではねじれが純粋に定義できない」という点について、直観的なイメージで解説します。

1. 法線モードは「全体運動」である#

分子の振動（法線モード）は、分子全体の原子が協奏的に動く「波」のようなものです。例えば、長い鎖状の分子がブルンと震えるとき、その動きは単純ではありません。「ある結合がねじれる動き」と、「分子全体が弓なりにたわむ動き（変角）」がセットになって起こります。理想的なねじれ: 雑巾を絞るように、軸に対して回転だけしている状態。実際の法線モード: 雑巾を絞りながら、同時に振り回している状態。これを数式的に言うと、「低振動数の法線モードベクトルは、ねじれ座標と変角座標の線形結合（混成）になっている」ということです。

2. 何が困るのか？#

従来の1次元回転モデル（Pitzer-Gwinn近似など）は、「純粋な回転運動」に対して適用する公式です。しかし、法線モードから得られる振動数 $\omega$ は、「ねじれ＋変角」の混ざったエネルギーを表しています。この「不純物が混ざった振動数」を、「純粋な回転の式」に無理やり代入してしまうと、エントロピーや分配関数の計算結果に大きな誤差（物理的な矛盾）が生じます。

3. 内部座標による解決策#

そこで本手法（MS法）では、法線モード（全体運動）の結果を一度脇に置き、**「局所的な視点」に切り替えます。アプローチ: 「分子全体がどう揺れようと知ったことではない。この4つの原子（A-B-C-D）で作る二面角 $\phi$ だけを動かしたときの、バネの強さ（力定数）はどうなっているか？」このように内部座標（二面角）**を直接指定して力定数 $k_{j,\tau}$ を計算すれば、他の変角モードが混ざる余地はなくなり、「純粋なねじれの硬さ」だけを抽出できます。この「純粋なねじれ情報」を使って補正係数 $f$ を作ることで、法線モード解析の「混ざってしまう弱点」を克服しているのです。

2.2 内部座標によるねじれ非調和性 ( $f_{j,\tau}$ ) の導出#

法線モード解析では、ねじれと変角が混成するため、純粋なねじれ振動数を定義することが困難である。そこで本手法では、内部座標（二面角） に基づく局所的な記述を採用する。

構造 $j$ におけるねじれモード $\tau$ の内部座標力定数 $k_{j,\tau}$ は、ポテンシャル $V$ の二面角 $\phi_\tau$ に対する2階微分として定義される。

k_{j,\tau} \equiv \frac{\partial^2 V}{\partial \phi_\tau^2}\bigg|_{\phi_{\tau, eq, j}}, \quad \bar{\omega}_{j,\tau} = \sqrt{\frac{k_{j,\tau}}{I_{j,\tau}}}

ここで $\bar{\omega}_{j,\tau}$ は内部座標周波数、 $I_{j,\tau}$ は Pitzer の定義による（非結合の）慣性モーメントである。

この内部座標を用いて、参照ポテンシャル（コサイン型）に対する古典的分配関数と、調和近似分配関数の比をとることで、非調和性補正係数 $f_{j,\tau}$ が解析的に導出される。

f_{j,\tau} = \frac{\sqrt{2\pi \beta k_{j,\tau}}}{M_{j,\tau}} \exp\left( -\frac{\beta k_{j,\tau}}{M_{j,\tau}^2} \right) I_0 \left( \frac{\beta k_{j,\tau}}{M_{j,\tau}^2} \right)

ここで $M_{j,\tau}$ は構造 $j$ のねじれ $\tau$ における局所的な周期性（Wellsの数に相当）、 $I_0$ は第1種変形ベッセル関数である。この式は、慣性モーメント $I_{j,\tau}$ に依存せず、力定数と周期性のみで決定される点が重要である。これにより、慣性項に由来するモード結合の問題を、次節の $Z_j$ 項へと分離することが可能となる。

2.3 モード間結合の補正 ( $Z_j$ と Kilpatrick-Pitzer 行列)#

低周波数モード間の結合は、ポテンシャル的な結合（力定数行列の非対角項）と、運動学的な結合（慣性テンソルの非対角項）に大別される。 $f_{j,\tau}$ は主に単一モード内のポテンシャル非調和性を補正するものであるため、モード間の結合は補正係数 $Z_j$ によって取り扱われる。

$Z_j$ は、低温域（調和近似が有効な領域）と高温域（自由回転および古典的結合が支配的な領域）を接続するように設計されている。

Z_j = \frac{g_j Q_{rot,j} Q_j^{HO} + (1-g_j) Q_j^{imp}}{Q_{rot,j} Q_j^{HO}}

ここで $Q_j^{imp}$ は改善された分配関数であり、以下の要素を含む。

Q_j^{imp} = Q_{rot,j} Q_j^{HO} Z_j^{int} Z_j^{coup}

2.3.1 運動学的結合の処理 ( $Z_j^{coup}$ )#

$f_{j,\tau}$ の導出において、各ねじれは独立した Pitzer 慣性モーメント $I_{j,\tau}$ を持つと仮定された。しかし実際には、あるねじれ運動は他のねじれや分子全体の回転と連動するため、運動エネルギー行列の非対角成分（結合）が生じる。

これを補正するために、Kilpatrick と Pitzer (1949) による内部回転運動エネルギー行列 $D$ を導入する。 $D$ 行列はねじれ間の運動学的結合を記述するものであり、補正係数 $Z_j^{coup}$ は以下のように定義される。

Z_j^{coup} = \left( \frac{\det D}{\prod_{\tau=1}^{t} I_{j,\tau}} \right)^{1/2}

この項は、個々の慣性モーメントの積（非結合近似）に対する、結合を含んだ行列式（厳密解）の比を表しており、高温極限において運動学的結合を分配関数に正しく反映させる役割を果たす。

2.3.2 振動空間の整合性 ( $Z_j^{int}$ )#

高温極限において、法線モード振動数 $\omega_{j,m}$ の積は、内部座標振動数 $\bar{\omega}_{j,\tau}$ と残りの振動モード $\bar{\bar{\omega}}_{j,m}$ の積へと置き換えられるべきである。 $Z_j^{int}$ はこの座標変換に伴うヤコビアン的な補正を行う。

コラム： $Z_j$ って何？ —— 「別々に回る」のと「繋がって回る」のは違う#

数式に出てくる補正係数 $Z_j$ （特に $Z_j^{coup}$ ）は、一言で言えば**「動きの連動（カップリング）による補正」**です。

1. 運動学的結合 ( $Z_j^{coup}$ )：フィギュアスケートのペア#

ねじれ運動における「運動学的結合」を、フィギュアスケートのペアに例えてみましょう。独立した動き ( $I_{j,\tau}$ の積):2人のスケーターが、リンク上で離れて、それぞれ勝手にスピンしている状態です。Aさんの回転はBさんに何の影響も及ぼしません。計算は簡単で、Aの慣性＋Bの慣性です。

結合した動き ( $D$ 行列):2人が手をつないでスピンしている状態です。この時、Aさんが回転の勢いを変えると、つないだ手を介してBさんも振り回されます。また、2人全体の重心も回転に関わります。この「手をつないでいる効果（非対角項）」を含めて厳密に計算したのが、KilpatrickとPitzerの $D$ 行列です。 $Z_j^{coup}$ の式 $\sqrt{\det D / \Pi I}$ は、「手をつないで回る現実の複雑さ」を「バラバラに回る単純な計算」で割った比率です。高温になって分子が激しく回転（自由回転）し始めると、この「手をつないでいる効果（遠心力やコリオリ力のような相互作用）」が無視できなくなるため、この係数で補正をかけるのです。

2. 振動空間の整合性 ( $Z_j^{int}$ )：通貨の両替レート#

法線モード振動数 $\omega$ と、内部座標振動数 $\bar{\omega}$ は、言ってみれば**「通貨」**が違います。

法線モード: 分子全体の協奏的な揺れ（円）内部座標: 局所的なねじれの揺れ（ドル）計算の途中で、ねじれ部分だけを「円（法線モード）」から「ドル（内部座標）」に置き換えて計算しました。しかし、そのままでは全体の勘定が合いません。 $Z_j^{int}$ は、この座標変換に伴う**「両替レート（ヤコビアン）」**のようなものです。これを掛けることで、座標系を入れ替えても、位相空間の全体積（状態の総数）が正しく保たれるように調整しています。

3. 結論#

低温では、分子はカチコチに固まっているので、単純なバネ（調和振動子）として扱えます（ $g_j \approx 1$ ）。高温では、分子はぐにゃぐにゃ動き回り、互いに引っ張り合います（ $g_j \approx 0$ ）。 $Z_j$ は、この「低温の秩序ある世界」と「高温の無秩序で絡み合った世界」を、矛盾なく滑らかに繋ぐためのスイッチの役割を果たしているのです。

2.4 強結合系へのボロノイ分割の適用#

一部の分子では、立体障害等によりねじれモード間の結合が極めて強く（Strong Coupling）、各ねじれ座標に対して独立した周期性パラメータ $M_{j,\tau}$ を割り当てることが幾何学的に不可能な場合がある。このような系（例：1-ペンチルラジカル）に対し、本手法ではボロノイ分割（Voronoi Tessellation） を用いたアプローチを採用する。

ねじれ角空間（ $\phi_1, \dots, \phi_{t_{SC}}$ ）において、各構造 $j$ を母点とするボロノイセルを生成し、そのセル体積 $\Omega_j^{SC}$ を計算する。強結合しているねじれ群（ $t_{SC}$ 個）に対する有効周期性パラメータ $M_j^{SC}$ は、以下のように定義される。

M_j^{SC} = \frac{2\pi}{(\Omega_j^{SC})^{1/t_{SC}}}

このアプローチにより、モード間結合によって歪んだポテンシャル井戸の形状（相空間上の体積）を、単一のパラメータ $M$ に縮約して評価することが可能となる。これは、分離不可能な多次元ポテンシャル面を扱うための強力な幾何学的解決策である。

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

3.1 MS-AS法 (Multi-Structural All Structures)#

全構造を考慮する最も厳密な手法である。

配座探索: PES上の全ての極小点 $j=1 \dots J$ を特定する。
局所解析: 各構造において、回転定数、法線モード振動数、および内部座標力定数 $k_{j,\tau}$ を算出する。
結合の評価: Kilpatrick-Pitzer 行列 $D$ を構築し、 $Z_j^{coup}$ を計算する。
パラメータ決定: 構造の対称性に基づき $M_{j,\tau}$ を決定する。強結合系の場合はボロノイ分割を用いる。
総和: 式(4)に基づき全分配関数を算出する。

3.2 MS-RS法 (Reference Structure)#

計算コスト削減のため、参照構造（Reference Structure）から独立ねじれ近似によって生成される構造のみを考慮する手法である。大規模分子において構造数は指数関数的に増大するため、本手法の線形スケーリング特性（ $O(t)$ ）は実用上重要である。ただし、参照構造からの摂動として扱えないほど結合が強い場合、精度は低下する可能性がある。

4. 低周波数モード結合の実証的評価#

4.1 エタノール：モード混成の解決#

エタノール（ $CH_3-CH_2-OH$ ）は、C-C結合およびC-O結合周りの2つのねじれを持つ。標準的な法線モード解析では、これら2つのねじれは低周波数領域で強く混成しており、純粋な「C-Cねじれ」や「C-Oねじれ」として分離することはできない。

MS法の成果: 内部座標を用いることで、この混成の影響を受けずに各ねじれの寄与を評価できる。
結果: 1000 K における分配関数は、単一構造HO近似（SS-HO）と比較して約3倍の値となる。これは、3つの配座（trans, gauche+, gauche-）の存在だけでなく、各配座における非調和性とモード結合が正しく考慮された結果である。

4.2 1-ペンチルラジカル：立体的強結合とボロノイ分割#

1-ペンチルラジカルは4つのねじれ自由度を持つ。単純な分離近似を仮定した場合、ラジカル中心（ $CH_2$ ）の回転対称性等を考慮すると、理論上は非常に多くの構造（分離可能であれば最大162個程度）が予想される。しかし、実際には立体障害による強いモード間結合のため、安定な構造はわずか15個しか存在しない。

強結合の処理: C1-C2ねじれとC2-C3ねじれは強く相関しており、独立したパラメータ $M_{j,\tau}$ を定義することが困難である。
ボロノイ分割の適用: 結合したねじれ空間をボロノイ分割し、体積 $\Omega_j^{SC}$ から有効パラメータを算出した結果、異なる分割スキーム（NS=2:2 や NS=4:0）間でも分配関数は良好に一致した（400 Kで比率1.30以内）。
結論: これは、ボロノイ分割を用いたアプローチが、配座空間の「体積」を正しく評価することで、強結合系の熱力学量をロバストに予測できることを示している。

4.3 1-ブタノール：浅い極小点と結合#

1-ブタノールでは、立体障害により一部の配座が極めて浅い極小点となり、実質的に遷移状態のような振る舞いをする。このような系では、低周波数モードが広範な変位を持つため、モード間の非調和結合が顕著になる。 MS-AS法による計算結果は、実験値に基づくエントロピーと高温域まで良好に一致し、多数の低エネルギー配座とそれらの結合を無視するSS-HO法（1000 Kで約31倍の過小評価）の欠陥を浮き彫りにした。

5. 結論#

Zhengらによる多重構造法は、複雑分子の熱化学計算におけるボトルネックであった「ねじれ非調和性」と「モード間結合」の問題に対し、以下の解決策を提示した。

内部座標の採用: 法線モードの混成問題を回避し、局所的な非調和性を記述可能にした。
運動学的結合の厳密な補正: Kilpatrick-Pitzer 行列式を含む $Z_j^{coup}$ 項の導入により、ねじれ間の慣性相互作用を高温極限で正しく再現した。
幾何学的アプローチによる強結合処理: ボロノイ分割を用いることで、立体障害等によりポテンシャル面が複雑に歪んだ系（Strong Coupling）に対しても、物理的に妥当な分配関数を算出する道を拓いた。

本手法は、配座探索と振動数計算のみ（遷移状態探索不要）で実行可能であり、中規模以上の柔軟な分子に対する標準的な高精度熱化学計算手法として位置づけられる。

参考文献#

Jingjing Zheng et al., Phys. Chem. Chem. Phys., 2011, 13, 10885.
B. M. Wong and W. H. Green, Mol. Phys., 2005, 103, 1027.
H. E. O’Neal and S. W. Benson, in Free Radicals, ed. J. H. Koshi, Wiley, New York, 1973, p. 275.
K. S. Pitzer, J. Chem. Phys., 1946, 14, 239.
Zheng et al., Phys. Chem. Chem. Phys., 2011, 13, p. 10888, Eq (4).
J. E. Kilpatrick and K. S. Pitzer, J. Chem. Phys., 1949, 17, 1064.
Zheng et al., Phys. Chem. Chem. Phys., 2011, 13, p. 10891.

補遺：強結合ねじれ系におけるボロノイ分割の幾何学的・物理的解釈#

1. 導入：なぜボロノイ分割が必要か#

多重構造（MS-AS）法において、ねじれ非調和性補正係数 $f_{j,\tau}$ を計算するためには、各構造 $j$ の各ねじれモード $\tau$ に対する局所的な周期性パラメータ $M_{j,\tau}$ を決定する必要があります。

弱結合（Nearly Separable）の場合: ねじれ座標間が直交に近いため、各次元を独立に分割することで $M_{j,\tau}$ を決定できます（例： $2\pi$ を3つの構造で等分するなら $M=3$ ）。
強結合（Strongly Coupled）の場合: 立体障害等により、あるねじれ角の変位が他のねじれ角の許容範囲を大きく制限します。この場合、配置空間上の安定領域（Basin）は単純な超直方体ではなく、複雑に歪んだ多面体となります [cite: 338-343]。

この歪んだ多次元領域の「体積」を定量的に評価し、等価な有効パラメータ $M$ に縮約するために、ボロノイ分割（Voronoi Tessellation）が導入されます。

2. 直観的な説明：ポテンシャル井戸の「勢力圏」#

化学反応におけるポテンシャルエネルギー曲面（PES）を、複数の谷（極小点）を持つ地形図としてイメージしてください。

母点（サイト）の配置: ねじれ角空間上に、構造最適化で見つかった各配座（極小点）をプロットします。これらがボロノイ分割における「母点」となります。
勢力圏の定義: 空間内の任意の位置（ねじれ角のセット）が、「どの母点に最も近いか」を考えます。空間上のすべての点を、最も近い母点に帰属させることで、空間は多面体のセル（小部屋）に分割されます。化学的には、これは**「ある構造から出発して、勾配に従って緩和させたときに到達する領域（Attraction Basin）」の幾何学的な近似**とみなせます。
体積とエントロピー: ある構造 $j$ のボロノイセルの体積 $\Omega_j$ が大きいということは、その構造が配置空間内で広い領域を占有していることを意味します。分配関数（状態和）の観点では、「広い井戸」は多数の状態を含むため、エントロピー的に有利になります。

MS法では、この複雑な形状をしたセルの体積 $\Omega_j$ を計算し、「もしこのセルが単純な立方体だったとしたら、一辺の長さ（有効的なポテンシャル幅）はどれくらいか？」という逆算を行うことで、有効周期性 $M_j^{SC}$ を決定します。

3. 厳密な説明：数理的定義とアルゴリズム#

線形代数の言葉を用いて、強結合系（Strongly Coupled, SC）におけるパラメータ決定プロセスを記述します。

3.1 空間と距離の定義#

$t_{SC}$ 個の強結合ねじれモードによって張られる配置空間を $\mathbb{T}^{t_{SC}}$ とします。ねじれ角は周期的であるため、これは $t_{SC}$ 次元トーラス（あるいは周期境界条件を持つユークリッド空間）と同相です。

この空間内の距離 $d(\mathbf{x}, \mathbf{y})$ は、ねじれ角の周期性を考慮したユークリッドノルム（Euclidean norm）として定義されます。

3.2 母点の集合#

区別可能な構造の集合に加え、対称性 $\sigma_\tau$ に由来する「区別不可能なレプリカ」も含めた母点集合 $S$ を定義します。

S = \{ \mathbf{x}_k \in \mathbb{R}^{t_{SC}} \mid k = 1, \dots, \tilde{J} \}

ここで $\tilde{J}$ は、区別可能な構造数 $J$ よりも大きくなります（例：3回対称のメチル回転を含む場合、物理的に等価な3つの点が空間内に配置されます）。

3.3 ボロノイセルの定義#

構造 $j$ に対応するボロノイセル $V_j$ は、以下の不等式を満たす点 $\mathbf{x}$ の集合として定義される凸多面体です。

V_j = \{ \mathbf{x} \in \mathbb{T}^{t_{SC}} \mid \|\mathbf{x} - \mathbf{x}_j\| \le \|\mathbf{x} - \mathbf{x}_k\|, \forall k \neq j \}

3.4 有効周期性パラメータの導出#

各セル $V_j$ の超体積（Hypervolume） $\Omega_j^{SC}$ を、Convex Hull法（qhull等）を用いて数値的に算出します。全空間の体積は $(2\pi)^{t_{SC}}$ であり、以下の規格化条件が満たされます。

\sum_{j=1}^{\tilde{J}} \Omega_j^{SC} = (2\pi)^{t_{SC}}

この体積 $\Omega_j^{SC}$ を用いて、構造 $j$ における強結合ねじれ群の有効周期性パラメータ $M_j^{SC}$ は次式で与えられます。

M_j^{SC} = \frac{2\pi}{(\Omega_j^{SC})^{1/t_{SC}}}

3.5 物理的意味#

上式は、体積 $\Omega_j^{SC}$ を持つ $t_{SC}$ 次元超立方体の一辺の長さ $\Delta \phi_{eff}$ を求め、その逆数（ $2\pi / \Delta \phi_{eff}$ ）を周期性（Wellsの密度）とみなす操作に相当します。これにより、複雑なモード間結合を含む配置空間の情報を、スカラー量 $M$ に縮約し、式(15b)の1次元的な非調和性補正関数 $f_{j,\tau}$ に適用することが可能となります。

補遺2：強結合系におけるボロノイ分割とモード間結合（直観と数理）#

本稿で解説した多重構造（MS）法において、ねじれモード間の結合（Coupling）が著しい場合に導入される「ボロノイ分割（Voronoi Tessellation）」について、その物理的直観と数理的定義を補足する。

1. 概念的解釈：ポテンシャル井戸の「勢力圏」#

ねじれモード間の結合とは、あるねじれ角の変化が他のねじれ角の許容範囲を制限する現象である。これを直観的に理解するために、配座空間を「地図」、安定配座（極小点）を地図上の「拠点」に見立てたモデルを導入する。

弱結合（Weak Coupling）の場合: 各ねじれが独立であるならば、拠点は格子状に規則正しく配置され、各拠点が支配する領域（ポテンシャル井戸の広さ）は均一な長方形（直積）となる。
強結合（Strong Coupling）の場合: 立体障害等により、特定の角度の組み合わせが禁止されると、拠点の配置は歪み、不規則になる。この結果、「ある拠点の支配領域は広いが、別の拠点の領域は極端に狭い」といった不均衡が生じる。

ボロノイ分割は、この歪んだ空間において「任意の点から最も近い距離にある拠点はどれか？」という基準で境界線を引く幾何学的手法である。これにより、結合によって形状が複雑化した各配座の「勢力圏（Basin of Attraction）」を、恣意性なく一意に定義することが可能となる。この「勢力圏の体積」こそが、その配座のエントロピー的な重み（分配関数への寄与）に対応する。