last_modified: 2026-01-07
生成AIによる自動生成記事に関する免責事項: 本記事は、学術論文（J. M. Bowman, 1978 および R. B. Gerber & M. A. Ratner, 1979）の内容に基づき、大規模言語モデル（生成AI）によって作成された解説記事です。数式の導出やアルゴリズムの記述においては、学術的な正確性を期しておりますが、実用的な実装や研究への応用にあたっては、必ず原著論文および標準的な教科書を参照してください。なお、参照している論文本文や図表の再配布を許諾するものではありません。
1. 序論：振動問題における平均場近似の必要性
分子の振動スペクトル解析において、もっとも基礎的かつ広範に用いられている手法は、ポテンシャルエネルギー曲面（PES）の平衡点周りでの二次展開に基づく 調和近似（Harmonic Approximation） である。正規モード解析（Normal Mode Analysis）は、結合された振動子系を独立した調和振動子の集合へと分離することを可能にし、低エネルギー領域においては極めて良好な記述を与える。
しかしながら、分子の高振動励起状態、あるいは結合解離や大振幅振動（Large Amplitude Motion）を伴う系においては、ポテンシャルエネルギーの 非調和性（Anharmonicity） およびモード間結合（Mode Coupling）が無視できない寄与を持つ 。これらの効果を取り入れるために、電子状態計算におけるHartree-Fock（HF）法のアナロジーとして発展したのが、振動自己無撞着場（Vibrational Self-Consistent Field: VSCF）法である 。
本稿では、Bowman およびGerber, Ratnerら による先駆的な研究に基づき、VSCF法の理論的枠組み、変分原理に基づく導出、そして計算機実装に資する詳細なアルゴリズムについて詳述する。
2. 歴史的背景と位置づけ
2.1 摂動論から変分法へ
非調和振動を取り扱う古典的な手法としては、二次以上のポテンシャル項を摂動として扱う**振動摂動論（Vibrational Perturbation Theory: VPT）**が存在する。VPTは計算コストが低く、多くの分光学的応用に供されてきたが、モード間結合が強い系や縮退に近い系においては、摂動展開の収束性が悪化するという本質的な問題を有している。
これに対し、1970年代後半、Bowmanらは、電子構造理論におけるSCFの概念を核の振動運動に適用することを提唱した 。Bowmanは、結合された非調和振動子系に対し、単純積波動関数（Hartree積）を用いた変分計算を行い、正確な量子力学解（Exact Quantum: EQ）と比較して極めて良好な結果が得られることを示した 。
2.2 半古典論との融合
同時期、GerberとRatnerは、VSCFの枠組みに半古典（Semiclassical）近似を導入したSC SCF法を提案した 。これは、各モードの波動方程式を厳密に解く代わりに、WKB近似に基づく量子化条件を用いるものであり、計算の簡便さと精度を両立させた手法である 。彼らの手法では、波動関数そのものではなく、エネルギーと古典的なターニングポイント（Turning Points）のみを自己無撞着に決定する 。これらの研究は、後のVibrational Configuration Interaction (VCI) や Vibrational Coupled Cluster (VCC) といった、より高精度な「振動相関」理論の礎となった。
3. 数学的定式化（Bowman 1978 に基づく）
本節では、VSCF方程式の実装可能な形式への導出を行う。
3.1 ハミルトニアンの定義
個の振動自由度を持つ系を考える。一般化座標を とし、ハミルトニアン を以下のように、各モード独立の項 と、モード間結合ポテンシャル の和として記述する 。
ここで、単一モードのハミルトニアン は通常、運動エネルギー項と対角ポテンシャル項（例：調和ポテンシャルやMorseポテンシャル）を含む 。
3.2 変分原理とHartree積
VSCF法では、全振動波動関数 を、各モードごとの一粒子波動関数（Modal） の積（Hartree積）で近似する 。ここで は各モードの量子数を表す。
この波動関数に対し、変分原理 を適用する 。ただし、各モード波動関数の規格直交性 を拘束条件として課す 。
3.2.1 コラム：なぜSlater行列式ではなくHartree積なのか？——「区別できる粒子」とコード上の表現
VSCF法を学ぶ際、電子状態計算（Hartree-Fock法）に馴染みがある読者は、「なぜ波動関数をSlater行列式（反対称化積）にしないのか？」という疑問を抱くかもしれない。この違いは、対象とする「粒子（モード）」の物理的性質と、それがプログラムコード上でどう扱われているかを対比することで明確に理解できる。
1. 物理的背景：空間の共有 vs 専用の座標軸
電子は全宇宙で区別がつかないフェルミ粒子であり、すべての電子が同一の3次元空間（）を共有している。そのため、パウリの排他律を満たすためにSlater行列式が必要となる。一方、振動モード（Vibrons）は「区別できる（Distinguishable）」粒子として扱われる。
これは、各モードが互いに独立した専用の座標軸を持っているためである。
- モード1（）: 例として「C-H結合の伸縮」。 軸上に住んでいる。
- モード2（）: 例として「H-O-H結合の変角」。 軸上に住んでいる。
これらは物理的に混ざり合うことがなく、ラベル（1番、2番）によって明確に区別される。したがって、波動関数の符号反転（反対称化）を考慮する必要はなく、単純な積（Hartree積） で記述することが正当化される。
2. コードとの対応関係(コードの中身の全体はあとで示す。)
この「交換項がない（Exchange-free）」という特徴は、実装コードの有効ポテンシャル計算部分に色濃く反映されている。電子のHartree-Fock法では、ハミルトニアン行列の構築に際して**クーロン項（）と交換項（）**の両方が現れるが、VSCFのコードには「交換項」に相当する処理が存在しない。
# VSCFのコード抜粋：有効ポテンシャルの構築（交換項がない！）
# モード2の波動関数を使って、モード2の位置の期待値（平均場）を計算
# これは電子計算における「クーロン積分」のようなもの（相手がどこにいるかの確率）
exp_q2_2 = self.get_expectation_value(self.C2, n2, self.Q2_2)
# モード1のハミルトニアンに「平均場」を足す
# H_eff = H_bare + (結合定数 * 相手の平均位置) * 自分の演算子
# ここに「自分と相手を入れ替える（Exchange）」ような項は一切入っていない
H1_eff = self.H1_0 + (self.lam * exp_q2_2) * self.Q1_2
上記のコードにおいて、exp_q2_2 は単なるスカラー値（数）であり、モード1のハミルトニアン H1_0 に対して外場（有効ポテンシャル）として加算されているだけである。もしSlater行列式が必要であれば、ここに行列のインデックスを入れ替える複雑な演算（交換積分）が必要になるが、VSCFではその必要がない。つまり、「異なる座標軸 を持つため、単純に相手の確率分布を『背景の壁（ポテンシャル）』として扱えばよい」 というHartree積の物理的描像が、「交換項を含まない単純な行列加算」 というアルゴリズムに直結しているのである。
3.3 VSCF方程式の導出
汎関数 を特定のモード の波動関数 で変分をとることにより、以下の方程式が得られる 。
ここで、括弧内の第二項は、モード 以外のすべてのモード がそれぞれの波動関数 で分布していると感じる平均場ポテンシャル（Mean-Field Potential）、あるいは有効ポテンシャル（Effective Potential）である 。これを と定義する。
最終的に、解くべきVSCF方程式（一粒子Schrödinger方程式）は以下の形となる 。
この方程式は、各モードの有効ポテンシャルが他モードの波動関数の形状に依存し、その波動関数はまた有効ポテンシャルに依存するという、相互依存関係（Self-Consistency）を持っている。したがって、反復法（Iterative Method）による数値解法が必要となる 。
3.4 全エネルギー（BowmanやGerberらが扱った2モード系の場合）
VSCF近似における全エネルギー は、固有値 の単純和とはならず、平均場ポテンシャルによる相互作用エネルギーの二重カウントを補正する必要がある。 定義通りに記述すれば、全エネルギーは全ハミルトニアンの期待値である 。
BowmanやGerberらが扱った （2モード系）の場合、全エネルギーは以下のように表される。
ここで、固有値 には既に平均場 の寄与が含まれているため、単純和 では相互作用エネルギー が2回カウントされてしまう。そのため、 を1回差し引くことで正しい全エネルギーが得られる。
3.4.1 全エネルギー（一般化）
VSCF法において、各モードの方程式から得られる固有値 は、そのモードが感じる「運動エネルギー ＋ 対角ポテンシャル ＋ 平均場ポテンシャル」のエネルギー和である。したがって、全系の総エネルギー を求める際に、単にこれらの固有値を足し合わせるだけでは、モード間の相互作用エネルギーが重複してカウントされてしまう。任意のモード数 および任意の相互作用形式 に対して常に成立する全エネルギーの定義は、得られたSCF波動関数 を用いた全ハミルトニアンの期待値である。
実装上は、各モードの「一粒子エネルギー期待値（運動＋対角ポテンシャル）」の総和に、全相互作用エネルギーの期待値を一度だけ加算するのが最も確実な手順となる。なお、相互作用 がペアポテンシャルの和 で記述される一般的なケース（多くの力場や分子モデル）においては、固有値の総和 と全エネルギー の関係は以下のようになる。
これは、各ペア相互作用 が、モード の有効ポテンシャルとモード の有効ポテンシャルの両方で計算され、固有値の和の中では「2重カウント」されているためである。Bowmanらの論文 1 で扱われている2モード系（）もこのケースに該当し、 で全エネルギーが与えられる。
※なお、多体相互作用項（3体以上）を含む場合や、モードごとに異なる分解を用いる定式化では、固有値と全エネルギーの関係式はこの限りではなく、常に期待値定義に立ち返る必要がある。
4. 数値的実装の詳細
理論を実際のプログラムコードとして実装する際の具体的な手順と、行列要素の計算手法について解説する。ここでは、Bowman論文で扱われている結合調和振動子モデルを例にとる。
4.1 具体的なハミルトニアンモデル
2つの結合した振動子（モード1, 2）を考える。ポテンシャルには3次の非調和項が含まれるとする 。
ここで単一モード項は調和振動子である 。
（※本実装では原子単位系 を仮定する ）
4.2 平均場ポテンシャルの構築
モード1に対する有効ポテンシャル を導出する。結合項 のうち、 に依存する部分について、モード2の波動関数 で期待値をとる 。
は に依存しないためそのまま残り、 は となる。
同様に、モード2に対する有効ポテンシャル は以下のようになる 。
（注： の項は に依存しない定数項となり、エネルギー固有値 をシフトさせるが、波動関数の形状には影響しない。ただし全エネルギー計算には寄与する。）
4.3 アルゴリズムの実装例
以下に、Bowman (1978) の Table I の結果を再現するためのPythonコードを示す。ここでは、VSCF法によって導かれる有効1次元Schrödinger方程式の解法として、調和振動子基底関数（Harmonic Oscillator Basis）による展開（変分法）を用いている。
import numpy as np
from scipy.linalg import eigh
from scipy.special import eval_hermite
import math
import matplotlib.pyplot as plt
class VSCF_RigorousSolver:
"""
結合調和振動子系に対するVSCF法の厳密な実装
Hamiltonian: H = h1(q1) + h2(q2) + lambda * (q1 * q2^2 + eta * q1^3)
Reference: J. M. Bowman, J. Chem. Phys. 68, 608 (1978).
"""
def __init__(self, omega1, omega2, lam, eta, n_basis=20):
self.w1 = omega1
self.w2 = omega2
self.lam = lam
self.eta = eta
self.n_basis = n_basis
# --- 1. 演算子行列の構築 (厳密解) ---
# 調和振動子基底における位置演算子 q の行列要素 <n|q|m>
# q = sqrt(1/2w) * (a + a_dag) を利用 (hbar=1, m=1)
self.Q1 = self._build_position_matrix(self.w1, self.n_basis)
self.Q2 = self._build_position_matrix(self.w2, self.n_basis)
# べき乗演算子の行列 Q^2, Q^3
self.Q1_2 = self.Q1 @ self.Q1
self.Q1_3 = self.Q1_2 @ self.Q1
self.Q2_2 = self.Q2 @ self.Q2
# 運動エネルギー + 調和ポテンシャル項 (対角行列)
# H0 = hbar * w * (n + 1/2)
self.H1_0 = np.diag([(i + 0.5) * self.w1 for i in range(self.n_basis)])
self.H2_0 = np.diag([(i + 0.5) * self.w2 for i in range(self.n_basis)])
# 初期密度（波動関数係数）：初期状態は非摂動の基底状態とする
self.C1 = np.zeros((self.n_basis, self.n_basis))
self.C2 = np.zeros((self.n_basis, self.n_basis))
np.fill_diagonal(self.C1, 1.0) # 単位行列
np.fill_diagonal(self.C2, 1.0)
def _build_position_matrix(self, omega, n_dim):
"""生成・消滅演算子による位置行列の解析的構築"""
mat = np.zeros((n_dim, n_dim))
coeff = np.sqrt(1.0 / (2.0 * omega))
# <n|q|n+1> = sqrt(n+1), <n+1|q|n> = sqrt(n+1)
for i in range(n_dim - 1):
elem = np.sqrt(i + 1)
mat[i, i+1] = elem
mat[i+1, i] = elem
return mat * coeff
def get_expectation_value(self, C_matrix, state_idx, Operator_matrix):
"""ある状態における演算子の期待値 <psi|O|psi>"""
vec = C_matrix[:, state_idx] # 列ベクトルが固有ベクトル
return vec.T @ Operator_matrix @ vec
def solve(self, state_idx=(0, 0), max_iter=100, tol=1e-9):
"""
VSCF方程式を反復法で解く
state_idx: ターゲットとする量子数 (n1, n2)
"""
n1, n2 = state_idx
E_total_prev = 0.0
print(f"--- VSCF Calculation for State (n1={n1}, n2={n2}) ---")
for it in range(max_iter):
# ---------------------------------------------------
# Step 1: モード1の有効ハミルトニアンを解く
# ---------------------------------------------------
# モード2の期待値計算: <q2^2>
exp_q2_2 = self.get_expectation_value(self.C2, n2, self.Q2_2)
# 有効ポテンシャル項 (行列形式)
# V_eff(q1) = lambda * <q2^2> * q1 + lambda * eta * q1^3
V1_eff_matrix = (self.lam * exp_q2_2) * self.Q1 + \
(self.lam * self.eta) * self.Q1_3
# 対角化
H1_eff = self.H1_0 + V1_eff_matrix
eps1_all, C1_new = eigh(H1_eff)
eps1 = eps1_all[n1] # 固有エネルギー epsilon_1
# 更新 (混合なしで直接更新)
self.C1 = C1_new
# ---------------------------------------------------
# Step 2: モード2の有効ハミルトニアンを解く
# ---------------------------------------------------
# モード1の期待値計算: <q1>, <q1^3>
exp_q1 = self.get_expectation_value(self.C1, n1, self.Q1)
exp_q1_3 = self.get_expectation_value(self.C1, n1, self.Q1_3)
# 有効ポテンシャル項
# V_eff(q2) = lambda * <q1> * q2^2 ( + 定数項: lambda * eta * <q1^3>)
V2_eff_matrix = (self.lam * exp_q1) * self.Q2_2
# 定数シフト項 (Energy shift term): lambda * eta * <q1^3>
# モード2の波動関数形状には影響しないが、エネルギー固有値には寄与する
const_shift_2 = self.lam * self.eta * exp_q1_3
# 対角化
H2_eff = self.H2_0 + V2_eff_matrix
eps2_all, C2_new = eigh(H2_eff)
# 固有値に定数シフトを加算して本来の epsilon_2 にする
eps2 = eps2_all[n2] + const_shift_2
# 更新
self.C2 = C2_new
# ---------------------------------------------------
# Step 3: 全エネルギーの計算と収束判定
# ---------------------------------------------------
# E_total = eps1 + eps2 - <W>
# <W> = < psi1 psi2 | lambda q1 q2^2 + lambda eta q1^3 | psi1 psi2 >
# = lambda <q1><q2^2> + lambda eta <q1^3>
# この補正項は相互作用エネルギーの二重カウントを防ぐために必要
W_mean = self.lam * exp_q1 * exp_q2_2 + \
self.lam * self.eta * exp_q1_3
E_total = eps1 + eps2 - W_mean
diff = abs(E_total - E_total_prev)
if diff < tol:
print(f"Converged at Iter {it+1}. Total Energy = {E_total:.8f}")
return E_total
E_total_prev = E_total
print("Warning: Max iteration reached.")
return E_total
def get_wavefunction_1d(self, mode_index, state_idx, grid_points):
"""実空間グリッド上の1次元波動関数 psi(x) を計算する"""
if mode_index == 1:
omega = self.w1
coeffs = self.C1[:, state_idx]
else:
omega = self.w2
coeffs = self.C2[:, state_idx]
psi_x = np.zeros_like(grid_points)
alpha = np.sqrt(omega) # スケール因子
for n, c_n in enumerate(coeffs):
if abs(c_n) < 1e-10: continue
xi = alpha * grid_points
# 正規化係数 (math.factorialを使用)
normalization = np.sqrt(1.0 / (2**n * math.factorial(n))) * \
(omega / np.pi)**0.25
phi_n = normalization * np.exp(-0.5 * xi**2) * eval_hermite(n, xi)
psi_x += c_n * phi_n
return psi_x
# --- メイン実行部 ---
if __name__ == "__main__":
# Bowman (1978) のパラメータ (Table I footnote a)
# w1^2 = 0.29375 -> w1 = 0.541987
# w2^2 = 2.12581 -> w2 = 1.458016
w1 = np.sqrt(0.29375)
w2 = np.sqrt(2.12581)
lam = -0.1116
eta = 0.08414
# ソルバーの初期化
solver = VSCF_RigorousSolver(w1, w2, lam, eta, n_basis=30)
# 基底状態 (0,0) の計算
E00 = solver.solve(state_idx=(0, 0))
print(f"Reference Exact Energy (0,0) from paper: 0.9916 ")
print(f"Calculated VSCF Energy (0,0): {E00:.5f}")
# --- 波動関数の可視化 ---
x_grid = np.linspace(-6, 6, 200)
psi1 = solver.get_wavefunction_1d(mode_index=1, state_idx=0, grid_points=x_grid)
# ゼロ次（調和振動子）の波動関数との比較用
alpha1 = np.sqrt(w1)
psi1_ho = (w1/np.pi)**0.25 * np.exp(-0.5 * (alpha1*x_grid)**2)
plt.figure(figsize=(6, 5))
plt.title("Mode 1 Wavefunction (Ground State)")
plt.plot(x_grid, psi1_ho, '--', label='Harmonic (0th order)', color='gray')
plt.plot(x_grid, psi1, '-', label='VSCF (Optimized)', color='blue', linewidth=2)
plt.xlabel('Coord q1')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)
plt.show()
5. 成果と考察：Bowman論文の結果に基づく検証（原著結果の再現）
Bowman（1978）の研究では、上述のモデルハミルトニアンを用いて、VSCF法の結果を厳密な量子力学計算（Exact Quantum: EQ）、および調和近似（Harmonic Oscillator: HO）と比較している。表1のデータ 1 に基づき、VSCF法の性能を評価する。
5.1 基底状態エネルギー
- 厳密解 (EQ): 0.9916
- VSCF: 0.9925
- 調和近似 (HO): 1.0000
調和近似がエネルギーを過大評価（+0.0084）しているのに対し、VSCF法は厳密解に対して +0.0009 程度の誤差に収まっている。これは、平均場近似が非調和相互作用の主要部分（特に静的な変形効果）を効果的に取り込んでいることを示唆している。
5.2 励起状態エネルギー および
- EQ (1, 0): 1.5159, VSCF: 1.5190
- EQ (0, 1): 2.4188, VSCF: 2.4214
励起状態においても、VSCFは調和近似（それぞれ1.5420, 2.4580）と比較して圧倒的に高い精度を示している 3。特に注目すべきは、モード間結合定数 を変化させた際の挙動である。結合が強くなる（が増加する）につれて、調和近似の破綻は顕著になるが、VSCFは依然として厳密解に近い値を追随する。
5.3 波動関数の形状変化
VSCF法の大きな利点は、エネルギーだけでなく、最適化された波動関数（Modal）が得られる点にある。図1に示されるように、SCF収束後の波動関数は、ゼロ次の調和振動子波動関数と比較して、ピーク位置のシフトや広がりの変化が見られる。具体的には、3次項のポテンシャル（）の影響により、波動関数の対称性が崩れ、ポテンシャルの軟らかい方向へ核の存在確率密度が染み出している様子が記述される。これは、「最適化された座標系」を自動的に見つけ出していると解釈できる。
6. VSCF法の限界と発展
VSCF法は強力な手法であるが、定義上、 平均場近似（相関の無視） に基づいているため、以下の限界が存在する。
- 動的相関の欠如: 瞬時的なモード間の相関運動（あるモードが伸びた瞬間に別のモードが縮むといった、時間の関数としての相関）は記述できない。電子状態理論における「電子相関」と同様の問題である。
- トンネル効果や共鳴の記述: 縮退または近縮退したモード間では、平均場描像が破綻しやすい（Fermi共鳴など）。これらの問題を克服するため、VSCFを参照関数とした事後補正法が開発されている。
- VCI (Vibrational Configuration Interaction): VSCF波動関数（配置）の線形結合をとる手法。Bowman自身も、SCF波動関数の線形結合（CI）が急速に収束することを指摘している。
また、GerberとRatnerが示したように、半古典論（Semiclassical Theory）との融合は、高密度な状態密度を持つ巨大分子系において、完全な量子計算の代替として計算コストを大幅に削減する道を開いた。彼らのSC SCF法は、波動関数そのものではなく、エネルギーと古典的なターニングポイント（Turning Points）のみを自己無撞着に決定する手法であり、量子SCFの結果を驚くべき精度で再現した。
7. 結論
VSCF法は、多原子分子の非調和振動解析において、調和近似と厳密な量子計算の中間に位置する実用的かつ強力なアプローチである。数学的には、Hartree積近似と変分原理に基づく連立非線形方程式の解法に帰着される。BowmanやGerberらの研究は、この手法が単なる数値計算テクニックにとどまらず、結合振動子系の物理的描像（有効ポテンシャル中の独立粒子運動）を理解するための堅固な理論的基盤を提供していることを示している 。現代の計算化学において、VSCFおよびその拡張手法は、気相クラスター、生体分子、表面吸着種などの赤外スペクトル予測において不可欠なツールとなっている。
参考文献
- Bowman, J. M. “Self-consistent field energies and wavefunctions for coupled oscillators,” J. Chem. Phys. 1978, 68, 608-610.
- Gerber, R. B.; Ratner, M. A. “A Semiclassical Self-Consistent Field (SC SCF) Approximation for Eigenvalues of Coupled-Vibration Systems,” Chem. Phys. Lett. 1979, 68, 195-198.
- Carney, G. D.; Sprandel, L. L.; Kern, C. W. “Variational Approaches to Vibration-Rotation Spectroscopy for Polyatomic Molecules,” Adv. Chem. Phys. 1978, 37, 305.
補足: VSCF法を簡易的に実装したもう一つの例
本文では非調和性が弱く、あまりVSCF法の強みがわかりにくいため、以下に非調和性とモード間結合の強い系に対する使用例を追加した。
import numpy as np
from scipy.linalg import eigh
import matplotlib.pyplot as plt
import math
class StrongAnharmonicVSCF_2D:
"""
2次元強結合非調和振動子系のVSCFおよび厳密解ソルバー
H = h1(q1) + h2(q2) + gamma * q1^2 * q2^2
h_i = p^2/2 + 0.5*w^2*q^2 + alpha * q^4
"""
def __init__(self, w1, w2, alpha, gamma, n_basis=25):
self.w1 = w1
self.w2 = w2
self.alpha = alpha
self.gamma = gamma
self.nb = n_basis
# --- 行列構築 (調和振動子基底) ---
self.Q1 = self._build_position_matrix(self.w1, self.nb)
self.Q2 = self._build_position_matrix(self.w2, self.nb)
# べき乗演算子
self.Q1_2 = self.Q1 @ self.Q1
self.Q1_4 = self.Q1_2 @ self.Q1_2
self.Q2_2 = self.Q2 @ self.Q2
self.Q2_4 = self.Q2_2 @ self.Q2_2
# 非摂動ハミルトニアン (対角成分 + 4次項)
# H0 = 0.5*w*Q^2 + alpha*Q^4 (運動エネルギー項は対角成分に吸収済みとみなす)
# 正確には H_ho = T + V_ho = diag(E_n). ここでは H_model = H_ho + alpha*Q^4 とする。
self.H1_0 = np.diag([(i + 0.5) * self.w1 for i in range(self.nb)]) + self.alpha * self.Q1_4
self.H2_0 = np.diag([(i + 0.5) * self.w2 for i in range(self.nb)]) + self.alpha * self.Q2_4
# 初期密度行列
self.C1 = np.eye(self.nb)
self.C2 = np.eye(self.nb)
def _build_position_matrix(self, w, n):
mat = np.zeros((n, n))
coeff = np.sqrt(1.0 / (2.0 * w))
for i in range(n - 1):
val = np.sqrt(i + 1)
mat[i, i+1] = val
mat[i+1, i] = val
return mat * coeff
def get_expectation(self, C, state_idx, Op):
vec = C[:, state_idx]
return vec.T @ Op @ vec
def solve_vscf(self, target_state=(0,0), max_iter=100, tol=1e-9):
n1, n2 = target_state
E_prev = 0.0
for it in range(max_iter):
# --- Mode 1 Optimization ---
# 平均場: <q2^2> を計算
exp_q2_2 = self.get_expectation(self.C2, n2, self.Q2_2)
# 有効ハミルトニアン: H1_eff = H1_0 + gamma * <q2^2> * q1^2
H1_eff = self.H1_0 + (self.gamma * exp_q2_2) * self.Q1_2
e1, self.C1 = eigh(H1_eff)
eps1 = e1[n1]
# --- Mode 2 Optimization ---
# 平均場: <q1^2> を計算
exp_q1_2 = self.get_expectation(self.C1, n1, self.Q1_2)
# 有効ハミルトニアン: H2_eff = H2_0 + gamma * <q1^2> * q2^2
H2_eff = self.H2_0 + (self.gamma * exp_q1_2) * self.Q2_2
e2, self.C2 = eigh(H2_eff)
eps2 = e2[n2]
# --- Total Energy ---
# 二重カウントの補正: E = eps1 + eps2 - <W>
W_avg = self.gamma * exp_q1_2 * exp_q2_2
E_total = eps1 + eps2 - W_avg
if abs(E_total - E_prev) < tol:
break
E_prev = E_total
return E_total, self.C1[:, n1], self.C2[:, n2]
def solve_exact(self, target_state_idx=0):
# Full CI (厳密対角化)
I = np.eye(self.nb)
# H_total = H1_0(x)I + I(x)H2_0 + gamma * Q1^2(x)Q2^2
H_total = np.kron(self.H1_0, I) + np.kron(I, self.H2_0) + \
self.gamma * np.kron(self.Q1_2, self.Q2_2)
eigvals, eigvecs = eigh(H_total)
return eigvals[target_state_idx], eigvecs[:, target_state_idx]
def _get_ho_func(self, n, w, grid):
"""調和振動子基底関数の値 (正規化定数のオーバーフロー対策済み)"""
# math.sqrtを使ってPythonの多倍長整数演算を行う
denom = math.sqrt(2**n * math.factorial(n))
norm = (w/np.pi)**0.25 / denom
herm = np.polynomial.hermite.Hermite.basis(n)(np.sqrt(w)*grid)
return norm * np.exp(-0.5 * w * grid**2) * herm
def get_vscf_wavefunction(self, vec, w, grid):
"""VSCFの1次元波動関数を構築"""
psi = np.zeros_like(grid)
for n, c in enumerate(vec):
if abs(c) > 1e-8:
psi += c * self._get_ho_func(n, w, grid)
return psi
def get_exact_marginal_wavefunctions(self, full_vec, grid):
"""
厳密解(2D)から、モード1およびモード2それぞれの周辺確率振幅を計算
psi_marginal_1(x1) = sqrt( integral |Psi(x1, x2)|^2 dx2 )
psi_marginal_2(x2) = sqrt( integral |Psi(x1, x2)|^2 dx1 )
"""
# 係数ベクトルを行列形式に戻す C[n1, n2]
C_matrix = full_vec.reshape((self.nb, self.nb))
# グリッド上の基底関数値
chi1 = np.array([self._get_ho_func(n, self.w1, grid) for n in range(self.nb)])
chi2 = np.array([self._get_ho_func(n, self.w2, grid) for n in range(self.nb)])
# 2次元波動関数 Psi[i, j] の再構成
# Psi(x1_i, x2_j) = sum_{n,m} chi1_n(x1_i) * C_nm * chi2_m(x2_j)
Psi_2d = np.einsum('ni, nm, mj -> ij', chi1, C_matrix, chi2)
dx = grid[1] - grid[0]
# モード1の周辺密度 (x2方向=axis1 で積分)
dens1 = np.sum(Psi_2d**2, axis=1) * dx
# モード2の周辺密度 (x1方向=axis0 で積分)
dens2 = np.sum(Psi_2d**2, axis=0) * dx
return np.sqrt(dens1), np.sqrt(dens2)
# --- 計算設定 ---
# 2つのモードの違いを見るため、振動数を非対称にする
w1 = 1.0
w2 = 1.3 # モード2の方を少し「硬い」調和ポテンシャルにする
alpha = 0.5
gamma = 0.5
solver = StrongAnharmonicVSCF_2D(w1, w2, alpha, gamma, n_basis=25)
# 計算実行
E_harm = 0.5*w1 + 0.5*w2
E_vscf, vec1, vec2 = solver.solve_vscf()
E_exact, vec_exact_full = solver.solve_exact()
print(f"Total Energy Results:")
print(f" Harmonic: {E_harm:.6f}")
print(f" VSCF : {E_vscf:.6f}")
print(f" Exact : {E_exact:.6f}")
# --- プロット用データ作成 ---
x = np.linspace(-3, 3, 200)
# 1. 波動関数 (VSCF vs Exact)
psi1_vscf = solver.get_vscf_wavefunction(vec1, w1, x)
psi2_vscf = solver.get_vscf_wavefunction(vec2, w2, x)
psi1_exact, psi2_exact = solver.get_exact_marginal_wavefunctions(vec_exact_full, x)
# 符号合わせ (位相補正)
def align_sign(target, ref):
# ピーク位置付近での符号を比較して合わせる
idx = np.argmax(np.abs(ref))
if np.sign(target[idx]) != np.sign(ref[idx]):
return -target
return target
psi1_exact = align_sign(psi1_exact, psi1_vscf)
psi2_exact = align_sign(psi2_exact, psi2_vscf)
# 調和近似の波動関数
psi1_harm = (w1/np.pi)**0.25 * np.exp(-0.5 * w1 * x**2)
psi2_harm = (w2/np.pi)**0.25 * np.exp(-0.5 * w2 * x**2)
# 2. 有効ポテンシャル
# モード1用: V_eff_1(q1) = 0.5*w1^2*q1^2 + alpha*q1^4 + gamma * <q2^2> * q1^2
exp_q2_2 = vec2.T @ solver.Q2_2 @ vec2
V1_eff = 0.5 * w1**2 * x**2 + alpha * x**4 + gamma * exp_q2_2 * x**2
# モード2用: V_eff_2(q2) = 0.5*w2^2*q2^2 + alpha*q2^4 + gamma * <q1^2> * q2^2
exp_q1_2 = vec1.T @ solver.Q1_2 @ vec1
V2_eff = 0.5 * w2**2 * x**2 + alpha * x**4 + gamma * exp_q1_2 * x**2
V1_harm = 0.5 * w1**2 * x**2
V2_harm = 0.5 * w2**2 * x**2
# --- 描画 (2パネル) ---
# 下部に凡例用のスペースを空けるため、figsizeを少し調整し、subplots_adjustを入れる
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 7))
plt.subplots_adjust(bottom=0.2) # 下側20%を凡例用に空ける
# モード1のプロット
# 修正: f'...' -> rf'...' にして \omega の警告を回避
ax1.set_title(rf'Mode 1 ($q_1$): $\omega_1={w1}$')
ax1.set_xlabel('$q_1$')
ax1.set_ylabel('Potential / Amplitude')
ax1.set_ylim(0, 5)
# ポテンシャル
ax1.plot(x, V1_harm, 'k:', alpha=0.4, label='Harmonic Pot.')
ax1.plot(x, V1_eff, 'b-', alpha=0.3, linewidth=2, label='VSCF Eff. Pot.')
# 波動関数
ax1_w = ax1.twinx() # 右軸
ax1_w.set_yticks([]) # 右軸の目盛りは消す
ax1_w.plot(x, psi1_harm, 'k--', alpha=0.4, label='Harmonic Approx')
ax1_w.plot(x, psi1_vscf, 'r-', linewidth=2, alpha=0.8, label='VSCF')
ax1_w.plot(x, psi1_exact, 'k:', linewidth=2.5, label='Exact (Marginal)')
# モード2のプロット
# 修正: rf'...' を使用
ax2.set_title(rf'Mode 2 ($q_2$): $\omega_2={w2}$')
ax2.set_xlabel('$q_2$')
ax2.set_ylim(0, 5)
# ポテンシャル
ax2.plot(x, V2_harm, 'k:', alpha=0.4)
ax2.plot(x, V2_eff, 'b-', alpha=0.3, linewidth=2)
# 波動関数
ax2_w = ax2.twinx()
ax2_w.set_yticks([])
ax2_w.plot(x, psi2_harm, 'k--', alpha=0.4)
ax2_w.plot(x, psi2_vscf, 'r-', linewidth=2, alpha=0.8)
ax2_w.plot(x, psi2_exact, 'k:', linewidth=2.5)
# 凡例 (図の下部に配置)
# モード1の凡例ハンドルを取得して代表させる（モード2も同じ色使いなので）
lines, labels = ax1.get_legend_handles_labels()
lines_w, labels_w = ax1_w.get_legend_handles_labels()
# まとめて1つの凡例にする
fig.legend(lines + lines_w, labels + labels_w,
loc='lower center',
bbox_to_anchor=(0.5, 0.05), # 図の底辺中央より少し上
ncol=5, # 横並びの列数
frameon=True, # 枠線をつける
fancybox=True, shadow=True)
plt.show()
補足2: 本文中のコードをIRスペクトルの計算に応用した例
VSCF法の真価は、非調和性による波動関数の「歪み」を記述できる点にある。これにより、調和近似では禁止されている倍音や結合音の遷移強度（IRスペクトル）を計算することが可能になる。 以下に、Bowman (1978) のモデル系に対し、VSCF計算を行ってエネルギー準位と遷移双極子モーメント（IRスペクトル）を算出するPythonコードを示す。
import numpy as np
from scipy.linalg import eigh
from scipy.special import eval_hermite
import math
import matplotlib.pyplot as plt
class VSCF_SpectrumSolver:
"""
VSCF法を用いて振動エネルギーと遷移強度（IRスペクトル）を計算するクラス
[Disclaimer / 免責事項]
This code is intended for educational purposes.
It demonstrates state-specific VSCF calculations and qualitative IR spectra
using a simplified model Hamiltonian and dipole moment surface.
It is not meant for quantitative spectroscopy or production-level calculations.
Hamiltonian: H = h1 + h2 + lambda*(q1*q2^2 + eta*q1^3)
Dipole Model: mu = q1 + q2 + 0.1*q1*q2 (Qualitative model)
"""
def __init__(self, w1, w2, lam, eta, n_basis=20):
self.w1 = w1
self.w2 = w2
self.lam = lam
self.eta = eta
self.nb = n_basis
# 演算子行列の構築
self.Q1 = self._build_pos_mat(self.w1, self.nb)
self.Q2 = self._build_pos_mat(self.w2, self.nb)
self.Q1_2 = self.Q1 @ self.Q1
self.Q1_3 = self.Q1_2 @ self.Q1
self.Q2_2 = self.Q2 @ self.Q2
# 非摂動ハミルトニアン
self.H1_0 = np.diag([(i + 0.5) * self.w1 for i in range(self.nb)])
self.H2_0 = np.diag([(i + 0.5) * self.w2 for i in range(self.nb)])
def _build_pos_mat(self, w, n):
"""調和振動子基底における位置演算子行列"""
mat = np.zeros((n, n))
c = np.sqrt(1.0 / (2.0 * w))
for i in range(n - 1):
v = np.sqrt(i + 1)
mat[i, i+1] = v
mat[i+1, i] = v
return mat * c
def get_exp(self, C, idx, Op):
"""期待値計算 <psi|Op|psi>"""
return C[:, idx].T @ Op @ C[:, idx]
def solve_vscf(self, target_state=(0,0), max_iter=50):
"""State-Specific VSCF calculation"""
n1, n2 = target_state
C1 = np.eye(self.nb)
C2 = np.eye(self.nb)
E_total = 0.0
for it in range(max_iter):
# Mode 1 SCF
exp_q2_2 = self.get_exp(C2, n2, self.Q2_2)
H1_eff = self.H1_0 + (self.lam * exp_q2_2) * self.Q1 + \
(self.lam * self.eta) * self.Q1_3
e1, C1 = eigh(H1_eff)
eps1 = e1[n1]
# Mode 2 SCF
exp_q1 = self.get_exp(C1, n1, self.Q1)
exp_q1_3 = self.get_exp(C1, n1, self.Q1_3)
H2_eff = self.H2_0 + (self.lam * exp_q1) * self.Q2_2
const_shift = self.lam * self.eta * exp_q1_3
e2, C2 = eigh(H2_eff)
eps2 = e2[n2] + const_shift
# Total Energy
W_mean = self.lam * exp_q1 * exp_q2_2 + self.lam * self.eta * exp_q1_3
E_curr = eps1 + eps2 - W_mean
if abs(E_curr - E_total) < 1e-9:
break
E_total = E_curr
return E_total, C1[:, n1], C2[:, n2]
def compute_dipole_matrix_element(self, vec1_i, vec2_i, vec1_f, vec2_f):
"""遷移双極子モーメント計算"""
S1 = vec1_f.T @ vec1_i
S2 = vec2_f.T @ vec2_i
q1_fi = vec1_f.T @ self.Q1 @ vec1_i
q2_fi = vec2_f.T @ self.Q2 @ vec2_i
# mu = q1 + q2 (線形近似) + 0.1*q1*q2 (電気的非調和性)
# ※調和近似との比較を公平にするため、ここでは電気的非調和項はあえて小さく/無視するか、
# あるいは「調和近似側」では厳密に選択則を適用する。
term_q1 = q1_fi * S2
term_q2 = S1 * q2_fi
return term_q1 + term_q2
# --- パラメータ設定 (Bowmanモデル) ---
w1 = np.sqrt(0.29375) # ~0.542
w2 = np.sqrt(2.12581) # ~1.458
lam = -0.1116
eta = 0.08414
solver = VSCF_SpectrumSolver(w1, w2, lam, eta, n_basis=20)
# --- 計算対象の状態 ---
states_to_calc = [(0,0), (1,0), (0,1), (2,0), (1,1), (0,2)]
results = {}
print("Calculating VSCF states...")
for s in states_to_calc:
E, v1, v2 = solver.solve_vscf(target_state=s)
results[s] = {'E': E, 'v1': v1, 'v2': v2}
# --- VSCFスペクトルデータの生成 ---
transitions_vscf = []
ground = results[(0,0)]
for s in states_to_calc:
if s == (0,0): continue
final = results[s]
freq = final['E'] - ground['E']
mu_fi = solver.compute_dipole_matrix_element(
ground['v1'], ground['v2'], final['v1'], final['v2']
)
intensity = mu_fi**2
transitions_vscf.append({'freq': freq, 'int': intensity, 'label': str(s)})
# --- 調和近似 (Harmonic) スペクトルデータの生成 ---
# Double Harmonic Approximation:
# 1. Potential is Harmonic (freq = w)
# 2. Dipole is Linear (mu = q), implies strictly selection rule delta_n = +/- 1
transitions_harm = []
# (1,0) Fundamental
freq_h1 = w1
int_h1 = 1.0 / (2.0 * w1) # <0|q|1>^2 = 1/2w
transitions_harm.append({'freq': freq_h1, 'int': int_h1, 'label': '(1,0)'})
# (0,1) Fundamental
freq_h2 = w2
int_h2 = 1.0 / (2.0 * w2)
transitions_harm.append({'freq': freq_h2, 'int': int_h2, 'label': '(0,1)'})
# Overtones/Combinations are strictly zero intensity in Double Harmonic approx.
# We add them with zero intensity just for frequency comparison
transitions_harm.append({'freq': 2*w1, 'int': 0.0, 'label': '(2,0)'})
transitions_harm.append({'freq': w1+w2, 'int': 0.0, 'label': '(1,1)'})
transitions_harm.append({'freq': 2*w2, 'int': 0.0, 'label': '(0,2)'})
# --- プロット (比較) ---
plt.figure(figsize=(10, 6))
# 1. Harmonic (Gray Dashed)
freqs_h = [t['freq'] for t in transitions_harm if t['int'] > 0]
ints_h = [t['int'] for t in transitions_harm if t['int'] > 0]
plt.stem(freqs_h, ints_h, linefmt='gray', markerfmt='D', basefmt=" ", label='Harmonic (Double Approx)')
# ゼロ強度の位置にも点線を引いて「位置のズレ」を示す
for t in transitions_harm:
plt.axvline(t['freq'], color='gray', linestyle=':', alpha=0.5)
# 2. VSCF (Blue/Red)
freqs_v = [t['freq'] for t in transitions_vscf]
ints_v = [t['int'] for t in transitions_vscf]
plt.stem(freqs_v, ints_v, linefmt='b-', markerfmt='bo', basefmt=" ", label='VSCF (Anharmonic)')
# ラベル付け
for t in transitions_vscf:
plt.text(t['freq'], t['int'], f" {t['label']}", rotation=90, verticalalignment='bottom', color='blue')
plt.title("IR Spectrum Comparison: Harmonic vs VSCF")
plt.xlabel("Frequency (Energy Difference)")
plt.ylabel("Intensity (Dipole Strength)")
plt.xlim(0, max(freqs_v)*1.1)
plt.ylim(0, max(ints_v)*1.3)
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
VSCF法による計算結果は調和近似と比較して低波数シフトを示しており、これは非調和性による ポテンシャルの軟化（softening） が適切に考慮された結果であるといえる。