メインコンテンツへスキップ
見出し画像

書記が物理やるだけ#418 配置間相互作用法の計算

    鈴華書記（Writer Rinka）

    ここからは電子相関を考慮した手法について紹介していく。


    問題

    画像


    説明

    HF法が単一の電子配置（スレーター行列式）で波動関数を表すのに対し、配置間相互作用法（CI法）は複数の電子配置を相互作用させることで、より現実に近い波動関数を構築する。

    画像


    解答

    CI法では、分子の全電子波動関数を、HF法で得られる参照配置（通常はHFの基底状態）と、そこから電子を仮想軌道に励起させた励起配置（スレーター行列式）の線形結合として定義する。このCI波動関数をシュレディンガー方程式に代入し変分原理を適用することで、行列の固有値問題であるCI方程式が導かれ、この方程式を解くことで、最も低い固有値として正確な基底状態エネルギーが、対応する固有ベクトルとして基底状態のCI係数が得られる。

    画像


    CI方程式に現れる行列Hの要素について、以下の性質が示される。非対角要素 Hij  がゼロでない（配置 Φi  と Φj  の間に相互作用がある）ことは、スレーター・コンドン則によって規定される。

    画像


    CI波動関数 Ψ を用いた1粒子密度行列 (1-RDM) は、分子軌道（MO）基底で定義され、電子の空間分布に関する情報を含む。

    画像


    画像


    CISD（一・二重励起配置間相互作用法）計算をPySCFで実行するコードの例を以下に示す。

    import numpy as np
from pyscf import gto, scf, ci

# --- 1. 分子構造と基底関数の定義 ---
mol = gto.Mole()
mol.atom = '''
O 0.00000 0.00000 0.1173
H 0.00000 0.75720 -0.4692
H 0.00000 -0.75720 -0.4692
'''
mol.basis = 'sto-3g'
mol.verbose = 3
mol.build()

# --- 2. Hartree-Fock (RHF) 計算 ---
print("\n--- Hartree-Fock Calculation ---")
mf = scf.RHF(mol)
mf.run()

# --- 3. CISD 計算 ---
print("\n--- CISD Calculation ---")
myci = ci.CISD(mf)
myci.run()

# ==========================================================
# 結果の解析
# ==========================================================

# --- A. エネルギー (Energy) ---
e_tot = myci.e_tot
e_corr = myci.e_corr
print(f"\n[A] Energy Results")
print(f"  Total Energy (CISD)     : {e_tot:.9f} Hartree")
print(f"  Correlation Energy      : {e_corr:.9f} Hartree")

# --- B. 波動関数 (Wave Function Analysis) ---
print(f"\n[B] Wave Function Analysis")

# 1. 生のCIベクトルと軌道情報を取得
raw_ci = myci.ci
nocc = myci.nocc
nmo = myci.nmo
nvir = nmo - nocc

# 2. 手動でベクトルを C0, C1, C2 に分解する
# PySCFのCISDベクトル構造: [c0, c1(nocc*nvir), c2(nocc*nocc*nvir*nvir)]
# c0: 参照状態 (1要素)
c0 = raw_ci[0]

# c1: 1電子励起 (インデックス 1 から 1 + nocc*nvir まで)
c1_end = 1 + nocc * nvir
c1 = raw_ci[1:c1_end].reshape(nocc, nvir)

# c2: 2電子励起 (残りの部分)
c2 = raw_ci[c1_end:].reshape(nocc, nocc, nvir, nvir)

# 3. 重み（Weight = 係数の二乗ノルム）を計算
w_c0 = c0**2
w_c1 = np.linalg.norm(c1)**2
w_c2 = np.linalg.norm(c2)**2
w_sum = w_c0 + w_c1 + w_c2

print(f"  Reference Weight (C0^2) : {w_c0:.5f}")
print(f"  Singles Weight (C1^2)   : {w_c1:.5f}")
print(f"  Doubles Weight (C2^2)   : {w_c2:.5f}")
print(f"  Normalization Check     : {w_sum:.5f}")

# --- C. 自然軌道 (Natural Orbitals) ---
print(f"\n[C] Natural Orbitals Analysis")
# CISDの1粒子密度行列(MO基底)を作成
dm1_mo = myci.make_rdm1()

# 対角化して自然軌道の占有数(occ)と係数(u)を得る
occ, u = np.linalg.eigh(dm1_mo)

# 占有数の大きい順に並べ替え
idx = np.argsort(occ)[::-1]
occ = occ[idx]
u = u[:, idx]

print("  Natural Orbital Occupations (Top 8):")
print(f"  {occ[:8]}")

# 自然軌道の係数をAO基底に変換 (C_nat = C_hf * U)
nat_orb_coeff = np.dot(mf.mo_coeff, u)


# --- D. 電荷 (Mulliken Population Analysis) ---
print(f"\n[D] Population Analysis (CISD Density)")

# MO基底の密度行列をAO基底に変換: P_ao = C * P_mo * C^T
dm1_ao = np.dot(mf.mo_coeff, np.dot(dm1_mo, mf.mo_coeff.T))

# Mulliken解析を実行 (密度行列dmにCISDのものを使用)
pop, chg = mf.mulliken_pop(dm=dm1_ao, verbose=False)

print(f"  Mulliken Charges per atom:")
for i, atom_chg in enumerate(chg):
    print(f"    Atom {i:>2} ({mol.atom_symbol(i)}): {atom_chg: .5f}")


    計算結果について：

    --- Hartree-Fock Calculation ---
converged SCF energy = -74.9630231384612

--- CISD Calculation ---
E(RCISD) = -75.01187316934728  E_corr = -0.04885003088605327

[A] Energy Results
  Total Energy (CISD)     : -75.011873169 Hartree
  Correlation Energy      : -0.048850031 Hartree

    CISD計算により、HF法では無視されていた電子相関が考慮され、系が約 0.049 Hartree (約 30.7 kcal/mol) 安定化している。


    [B] Wave Function Analysis
  Reference Weight (C0^2) : 0.97451
  Singles Weight (C1^2)   : 0.00016
  Doubles Weight (C2^2)   : 0.02464
  Normalization Check     : 0.99931

    電子相関による補正の大部分（約 2.5%）は、2電子励起 (C22 ) に由来する。これは、動的相関（電子の瞬間的な位置関係による相関）が主に2電子励起によって記述されるという量子化学の一般的な原則と一致している。


    [C] Natural Orbitals Analysis
  Natural Orbital Occupations (Top 8):
  [1.99999776 1.99849189 1.998032   1.97806234 1.97504092 0.02548022
 0.02489488]

    占有軌道 (Occupied): 最初の5つの軌道（水分子の価電子軌道）の占有数は、2.0 よりもわずかに小さい値（例: 1.975, 1.978）になっている。

    仮想軌道 (Virtual): それに対応して、HFでは占有数ゼロと見なされていた仮想軌道が、わずかに占有された値（例: 0.025, 0.024）を持っている。

    この「2.0からの減少」と「0.0からの増加」は、電子相関の結果である。電子が相互作用することで、HFでは完全に占有軌道に押し込められていた電子密度が、わずかに仮想軌道へと漏れ出し、波動関数がより柔軟な記述になっていることを示す。


    [D] Population Analysis (CISD Density)
  Mulliken Charges per atom:
    Atom  0 (O): -0.33583
    Atom  1 (H):  0.16791
    Atom  2 (H):  0.16791

    酸素原子は -0.336 の負電荷を、水素原子はそれぞれ +0.168 の正電荷を帯びている。これは、電気陰性度の高い酸素が、共有結合の電子を強く引きつけ、分極を起こしていることを示す。


    本記事のもくじはこちら：


    あなたへのおすすめ