電荷分割解析の理論

電荷分割解析の理論#

収束済みの波動関数には原子がまったく入っていません — 分子密度と非局在化した分子軌道の集合だけが あります。化学者が報告するあらゆる「原子電荷」は、したがって分割スキームの結果です:本質的に 分子的な量を原子核に戻して分ける方法の選択であり、異なるスキームは真に異なる選択をします、単に 1 つの「真の」原子電荷への異なる数値的近似というわけではありません(厳密に原子的な系を離れれば そのような単一の真理は存在しません)。本章は qc-rs の qc.prop.chrg/qc.prop.bond 名前空間が 実装する電荷分割解析法のすべて — Mulliken、Löwdin、Hirshfeld ストックホルダー族(素の版・反復版・ CM5・ADCH)、MBIS、NPA/NAO、VDD、そして Mayer/Wiberg 結合指数 — を導出します。qc-prop クレートの ソースに直接基づきます(それぞれが引用された論文からのクリーンルーム実装であり、JANPA のような 参照コードに対して交差検証されています)。

AO 分割電荷:Mulliken と Löwdin#

最も単純な分割スキームは完全に AO インデックス空間で動作し、実空間にはまったく触れません。どちらも 同じ対象 — AO 基底の占有行列 — から出発しますが、和を取る前に密度をどの基底で表現するかが 異なります。

Mullikenは密度行列を重なりと縮約した和を取り、各非対角の AO 対占有を、関与する 2 つの AO の 原子に完全に帰属させます(異なる場合は均等に分割します):

\[ P_\mu = (DS)_{\mu\mu}, \qquad N_A = \sum_{\mu\in A}P_\mu, \qquad q_A = Z_A - N_A. \]

Mulliken のよく知られた弱点は基底依存性です:結合対のどちらかの原子により拡散した関数を加えると、 実際の物理的電子分布をまったく変えずに両者間の重なり占有をシフトさせるので、\(q_A\) は同じ分子・ 幾何でも基底によって符号すら含めて大きく変わりえます。

Löwdin電荷はこれを、まず AO 基底を対称的に直交化する — 対角を取る前に密度の両側を \(S^{1/2}\) で 変換する — ことで修正し、恣意的に分割すべき非対角の重なり寄与を残しません:

\[ P_\mu = \bigl(S^{1/2}DS^{1/2}\bigr)_{\mu\mu}, \qquad N_A = \sum_{\mu\in A}P_\mu, \qquad q_A = Z_A - N_A. \]

Löwdin 電荷はまさにこの理由で Mulliken より著しく基底に対して安定です — そもそも恣意的に分割 すべき重なりが残っていません — しかしどちらのスキームも根本的には AO インデックス分割です: 実際に電子密度が空間のどこにあるかについては何も述べておらず、ある基底あるいは別の基底で表現 された密度行列をどう分けるかだけです。

Hirshfeld ストックホルダー族#

Hirshfeld分割は代わりに実空間で動作し、前分子(promolecule)— 球面平均化された孤立中性 自由原子密度を実際の分子幾何に配置した重ね合わせ — を、各点で分子の実際の密度をどれだけ「気前 よく」各原子に帰属させるべきかの参照として使います:

\[ w_A(\mathbf r) = \frac{\rho_A^0(\mathbf r)}{\sum_B\rho_B^0(\mathbf r)}, \qquad N_A = \int w_A(\mathbf r)\,\rho(\mathbf r)\,d\mathbf r, \qquad q_A = Z_A - N_A, \]

\(\rho_A^0\) は \(A\) を中心とする自由中性原子の密度です。構成上 \(\sum_Aw_A(\mathbf r)=1\) が至る ところで成り立つので、これは真の単位分割です — 空間のすべての点が完全な密度重みを何らかの原子の 組み合わせに寄与します — しかし参照前分子は分子の実際の結合状態にかかわらず常に中性の自由原子 から構築されます、これが素の Hirshfeld 電荷が化学的直感や双極子フィット電荷が示唆するより絶対値が 小さくなりがちであるという文書化された理由です:実際の分子中で真に帯電した原子でも、依然として 中性参照と比較されます。

Hirshfeld-Iは前分子を自己無撞着にすることでまさにこの恣意性を除去します:常に中性原子と 比較する代わりに、原子自身の現在の分数占有 \(N_A\) を挟む 2 つの整数電子数の間で孤立原子参照密度を 補間します:

\[ \rho_A^0[N_A](\mathbf r) = \bigl(\lceil N_A\rceil-N_A\bigr)\,\rho_A[\lfloor N_A\rfloor](\mathbf r) + \bigl(N_A-\lfloor N_A\rfloor\bigr)\,\rho_A[\lceil N_A\rceil](\mathbf r), \]

そして \(N_A=\int w_A\rho\,d\mathbf r\) を自己無撞着まで反復します — 各サイクルの前分子は前サイクル の分数占有から構築され、変化しなくなるまで続けます。これは素の Hirshfeld の「常に中性参照」バイアス を直接除去し、系統的により大きな絶対値の、よりよく収束した電荷を与えます(Bultinck ら, 2007)、 代償として(中性原子のものだけでなく)帯電原子参照密度の小さなライブラリ(挟む整数電子数)が必要に なります。

CM5は Hirshfeld-I とは逆の戦略を取ります:背後の分割スキームの自己無撞着性を固定する代わりに、 実験/高精度理論の分子双極子モーメントを再現するようフィットされた、純粋に経験的で幾何ベースの 事後補正を素の Hirshfeld 電荷に直接適用します:

\[ q_A^{\text{CM5}} = q_A^{\text{HPA}} + \sum_{B\ne A}T_{AB}\,B_{AB}(r_{AB}), \qquad B_{AB} = \exp\bigl[-\alpha(r_{AB}-R_A-R_B)\bigr], \]

\(\alpha=2.474\ \text{Å}^{-1}\)、\(R_A,R_B\) は表になった共有結合半径、\(T_{AB}\) は Hirshfeld から「真の」 双極子への大規模な参照集合の対応に対してフィットされた原子対固有の(ほとんどの元素対については 単純な原子ごとの差の)補正係数です。CM5 の補正は距離とともに指数関数的に減衰し(結合しているか それに近い対だけが有意に互いを補正します)、まったく自己無撞着解を必要としません — 単一の素の Hirshfeld 計算の後に適用される一発の閉形式補正です。

ADCHは関連するが構造的に異なる問題を解きます:電荷自体を経験的に補正する代わりに、Hirshfeld 電荷だけでは分子双極子モーメントを再現できないという事実を補正します、Hirshfeld が各原子自身の 原子双極子— その原子周りで実際の密度が球対称からどれだけずれているかの 1 次モーメント、 \(\boldsymbol\mu_A=-\int w_A(\mathbf r)(\mathbf r-\mathbf R_A)\,\rho(\mathbf r)\,d\mathbf r\) — を捨てているからです。ADCH はこの捨てられた原子双極子を、原子の隣接原子上に置かれた点電荷補正の 集合(Thole-Duijnen 式再分配)に展開し、2 つの厳密な制約を同時に満たす補正電荷 \(q_{AB}\) を解きます — 再分配の電荷中性と、1 次モーメントとしての原子双極子の厳密な再現です:

\[ \sum_Bq_{AB}=0, \qquad \sum_Bq_{AB}\,\mathbf r_B=\boldsymbol\mu_A, \qquad Q_A = q_A^{\text{HPA}}+\sum_Bq_{BA}. \]

これらは(外部参照データへのフィットではなく、局所重み付き共分散行列 \(\boldsymbol\Lambda=\sum_Bv_{AB} (\mathbf r_B-\langle\mathbf r\rangle)(\mathbf r_B-\langle\mathbf r\rangle)^{\mathsf T}\) を通じて解かれる) 厳密な線形制約なので、ADCH の最終点電荷から計算される分子双極子は構成上厳密に QM 双極子モーメント を再現します— 較正された近似ではなく確固たる数学的保証であり、qc-rs 実装が使う自然な正しさの 検定です。CM5 と異なり、ADCH はまったく元素固有のフィットパラメータを必要としません。

MBIS:自己無撞着性への異なる経路#

MBIS(最小基底反復ストックホルダー)は 2 つ目の自己無撞着ストックホルダースキームですが、 Hirshfeld-I とは異なる機構で自己無撞着性に到達します:表になった孤立原子参照密度の間を補間する 代わりに、各原子の前原子密度を原子中心 Slater(指数)殻の最小展開として直接表現します:

\[ \rho_A^0(\mathbf r) = \sum_iN_{Ai}\,f_{Ai}(\mathbf r), \qquad f_{Ai}(\mathbf r) = \frac{1}{8\pi\sigma_{Ai}^3}\exp\Bigl(-\frac{|\mathbf r-\mathbf R_A|}{\sigma_{Ai}}\Bigr), \]

各殻は個別に 1 電子に正規化されます(\(\int f_{Ai}=1\))。ストックホルダー重みはこれらの殻から Hirshfeld とまったく同様に構築されます、\(w_A=\rho_A^0/\sum_B\rho_B^0\)、しかし今度は殻占有 \(N_{Ai}\) と幅 \(\sigma_{Ai}\) の両方が、外部の帯電原子参照ライブラリなしに、分子自身の密度に 対して直接自己無撞着まで精密化されます:

\[ N_{Ai} = \int\rho(\mathbf r)\,\frac{\rho_{Ai}^0(\mathbf r)}{\rho_A^0(\mathbf r)}\,d\mathbf r, \qquad \sigma_{Ai} = \frac{1}{3N_{Ai}}\int\rho(\mathbf r)\,\frac{\rho_{Ai}^0(\mathbf r)}{\rho_A^0(\mathbf r)}\,|\mathbf r-\mathbf R_A|\,d\mathbf r, \qquad q_A = Z_A-\sum_iN_{Ai}. \]

これが Hirshfeld-I との本質的な構造的違いです:Hirshfeld-I の自己無撞着性はどの事前計算済み参照 密度と比較するか(小さな外部テーブルにわたる補間問題)に宿るのに対し、MBIS の自己無撞着性は 前原子自身の形状を手元の分子に直接フィットすること(反復の種となる初期中性原子殻占有以外に 外部参照密度を必要としない最適化問題)に宿ります。

NPA/NAO:自然軌道を介した基底頑健な占有#

Mulliken と Löwdin 電荷は AO インデックス分割であり、どちらも(Löwdin は Mulliken ほど深刻では ないが測定可能な程度に)使われた特定の AO 基底に敏感なままです、どちらの構成も、化学的に意味のある AO 寄与と、隣接原子の関数と大きく重なるだけの数値的に冗長で拡散した寄与を区別しないからです。 **自然電荷分割解析(NPA)**は、まず密度自体から直交規格な、原子ごとの、対称適応された基底 — 自然原子軌道(NAO)— を構築し、その後になって初めて占有を読み取ることで、任意に外部から選ばれた AO 基底を直接分割するのではなくこれを根本から解決します。

中心対象は占有行列 \(\mathbf P=\mathbf S\mathbf D\mathbf S\) です(裸の密度 \(\mathbf D\) では ありません)— 密度単独ではなく完全な \(SDS\) サンドイッチを使うことが、まさに結果の占有を基底頑健に するものです、AO 基底自身の重なり構造が素朴な密度行列の読みをどう膨らませたり縮めたりするかを 説明するからです。NAO の構成は一連の異なる代数的段階(Reed-Weinstock-Weinhold の元の手順の JANPA 再導出)であり、それぞれがきちんと定義された局所問題を解きます:

アルゴリズム 7 (自然原子軌道の構成(概略))

入力: AO 密度 \(D\)、AO 重なり \(S\)、原子ごと/角運動量ごとの AO ラベル。 出力: 自然原子軌道とその占有数。

  1. 原子内自然化:各 \((atom,l)\) ブロック内で \(m\) にわたり角度平均を取り、一般化固有値問題 \(\bar P\mathbf c=n\bar S\mathbf c\) を解く — 固有ベクトルが前 NAO、固有値 \(n\) がその占有重み。

  2. 分割:各 \((atom,l)\) について最高占有の前 NAO を、自由原子最小基底数まで、自然最小基底 (NMB)とする;残りは自然 Rydberg 基底(NRB)を形成する。

  3. NMB の占有重み付き対称直交化— 高占有軌道を優先的に保存する、通常の Löwdin 直交化の重み 付き類似物。

  4. NRB の各関数を(今や直交化された)NMB に対してSchmidt 直交化する。

  5. NRB 集合を再自然化する(post-Schmidt NRB 内で再び原子内自然化)、その重みを精密化する。

  6. 精密化された NRB 集合の占有重み付き直交化。

  7. 完全な NMB+NRB 集合内での最終原子内自然化、最終的な NAO を与える。

NAO 占有数は、最終 NAO 基底 \(\mathbf C\) における \(\mathbf C^{\mathsf T}(\mathbf S\mathbf D\mathbf S)\mathbf C\) の対角です;原子ごとに合計すると、ここでの他のすべてのスキームとまったく同様に \(N_A\) と \(q_A=Z_A-N_A\) が得られます。NPA を基底頑健にするものは、ステップ 1〜2 が、それらが どんな直交化によっても混ざる前に「真に占有され、化学的に意味のある」軌道(NMB)と「基底には 形式的に存在するが本質的に占有されていない、拡散/分極関数の」軌道(NRB)を明示的に分離することです — より大きく、より拡散した基底は、生の Mulliken/Löwdin 分割を摂動させるのとは異なり、ほぼゼロの 占有を持つ NRB 性を主に追加するだけです。

VDD:硬分割の実空間電荷#

Voronoi 変形密度(VDD)はさらに別の実空間アプローチを取り、Mulliken/Löwdin の AO 基底分割も Hirshfeld 族のなめらかなストックホルダー重みも避けます。代わりに、空間は硬いVoronoi セルに 分割されます — すべての点は完全に最も近い原子核に属し、なめらかな重み関数はまったくありません — そして電荷は各セルにわたる変形密度(実際の分子密度が Hirshfeld が使うのと同じ中性原子前分子と どれだけ異なるか)の積分です:

\[ q_A = -\int_{\text{Voronoi}(A)}\bigl[\rho(\mathbf r)-\rho^0(\mathbf r)\bigr]\,d\mathbf r, \]

\(\rho^0=\sum_B\rho_B^0\) は Hirshfeld が構築するのと同じ前分子重ね合わせです。物理的解釈は直接的 です:VDD は「全密度のどれだけが原子 \(A\) に属するか」(Hirshfeld の問い)ではなく「原子 \(A\) 自身の 空間領域に、その同じ硬分割された領域での中性前分子ベースラインと比べて、どれだけ多い(あるいは 少ない)密度があるか」を問います — 真に異なる問いであり、真に異なる(なめらかではなく硬い)空間 分割を使っており、単なる Hirshfeld の数値的変種ではありません。

Mayer と Wiberg 結合次数#

電荷は「原子にどれだけの」密度が属するかに答えます;結合次数は関連するが異なる問いに答えます — 原子対の間にどれだけの共有された、共有結合的性格が存在するかです。Mayer 結合次数は AO 密度と重なりから直接構築され、教科書的な結合次数の発想を非直交 AO 基底へ一般化します:

\[ B_{AB} = 2\sum_{\mu\in A,\nu\in B}\Bigl[(P^\alpha S)_{\mu\nu}(P^\alpha S)_{\nu\mu} + (P^\beta S)_{\mu\nu}(P^\beta S)_{\nu\mu}\Bigr], \]

閉殻(\(P^\alpha=P^\beta=D/2\))については \(B_{AB}=\sum_{\mu\in A,\nu\in B}(DS)_{\mu\nu}(DS)_{\nu\mu}\) に帰着します。各原子の原子価 \(V_A=\sum_{B\ne A}B_{AB}\) は、他のすべての原子への結合次数を 合計します(対角 \(B_{AA}\)、原子内項は、この和からは除外されますが参照のため完全な行列には保持 されます)。Wiberg 結合指数は、Löwdin 対称直交化基底で評価されたまさに Mayer 構成です — まず密度の両側を \(S^{1/2}\) で変換します(そうすると Mayer 式に現れるはずの重なりが恒等元になり 消えます):

\[ \tilde P^\sigma = S^{1/2}P^\sigma S^{1/2}, \qquad W_{AB} = 2\sum_{\mu\in A,\nu\in B} \Bigl[(\tilde P^\alpha_{\mu\nu})^2+(\tilde P^\beta_{\mu\nu})^2\Bigr]. \]

これはまさに Löwdin 電荷が Mulliken 電荷に対して持つのと同じ関係です — 同じ背後の式を適用する前に 対称直交化された AO 集合への基底変換 — そして理由も同じです:素の(Mayer/Mulliken)非直交基底式が さもなくば行わなければならない、基底依存の恣意的な重なり分割の決定を除去します。

検証済み例:すべてのスキーム、1 つの分子#

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
m = qc.chk.new(atom=water, ao="def2-svp", unit="angstrom").scf(ref="r").run()

m.prop.chrg.mulliken()["charges"]      # [-0.3471, 0.1735, 0.1735]
m.prop.chrg.lowdin()["charges"]        # [-0.1555, 0.0778, 0.0778]
m.prop.chrg.hirshfeld()                # [-0.3277, 0.1639, 0.1639]
m.prop.chrg.hirshfeld_i()              # [-0.9181, 0.4590, 0.4590]
m.prop.chrg.cm5()                      # [-0.6642, 0.3321, 0.3321]
m.prop.chrg.adch()                     # [-0.7576, 0.3788, 0.3788]
m.prop.chrg.mbis()                     # [-0.8702, 0.4351, 0.4351]
m.prop.chrg.npa()["charges"]           # [-0.9115, 0.4557, 0.4557]
m.prop.chrg.vdd()                      # [-0.3023, 0.1518, 0.1504]

m.prop.bond.mayer()["valence"]         # [1.9948, 1.0040, 1.0040]
m.prop.bond.wiberg()["valence"]        # [2.3054, 1.2081, 1.2081]

すべてのスキームが、水の酸素は負に、両水素は等しく(逆符号で)正に帯電していることに一致します、 化学的に期待される通りです — しかし大きさは Löwdin の控えめな \(-0.16\) から Hirshfeld-I/NPA/ MBIS のはるかに大きな \(\sim-0.87\)〜\(-0.92\) まで著しく変わります。この広がりは 1 つの真の数値に ついての数値的な不一致ではありません;まさに本章の冒頭で述べた点です — これらは「原子にどれだけの 電荷が属するか」の真に異なる、自己無撞着な定義であり、手法をまたいで電荷を比較することは、 1 つの手法内で関連する分子の系列にわたって一貫して比較するよりはるかに意味が薄いのです。 Mayer と Wiberg の原子価はどちらも正しく 2 つの O–H 単結合を識別します(各水素で原子価 \(\approx1\)、酸素で \(\approx2\))、Wiberg の Löwdin 基底構成は Mayer の生の非直交構成よりいくらか 大きな絶対値を与えます — Löwdin/Mulliken 電荷が示すのと同じ一般的な基底依存パターンと整合します。

練習 16

  1. Mulliken 電荷は Löwdin 電荷より基底依存性が大きいことが知られています。上記の AO 分割導出を 使って、Löwdin の構成が除去し Mulliken のものが除去しない具体的にどのステップか、そしてなぜ そのステップが基底依存性が入り込む場所なのかを一文で説明しなさい。

  2. Hirshfeld-I と MBIS はどちらも自己無撞着性に到達しますが、異なる機構によってです。各スキームで 反復ごとに何が変化するか(収束させている「未知数」は何か)を一文で説明しなさい。

  3. ADCH は構成上分子双極子モーメントの厳密な再現を保証しますが、CM5 は(フィットされた経験的補正 を通じて)近似的にしか再現しません。両方とも同じ素の Hirshfeld 電荷から出発します。この 保証対近似という違いを説明する、それぞれがこれらの電荷をどう補正するかの構造的違いは何ですか。

電荷と結合次数は密度の原子間分布を記述します;位相解析の章 は 代わりに密度自身の臨界点構造(QTAIM、ELF/LOL)を直接研究します — 同じ収束済み密度から化学的意味を 抽出する、真に異なる、分割不要な方法です。