红外与拉曼光谱

1. 振动光谱的物理基础

分子的振动光谱——红外（IR）光谱与拉曼（Raman）光谱—— 是研究分子结构、官能团和化学键信息的核心工具。 两者都起源于分子的振动运动，但物理机制截然不同： IR 依赖于振动过程中偶极矩的变化， Raman 依赖于振动过程中极化率的变化。

1.1 偶极矩与红外活性

分子在外部光场（电磁波）的振荡电场 $\mathbf{E}(t)$ 中， 其固有偶极矩 $\vec{\mu}$ 与外场发生相互作用，引起能量吸收。 根据量子力学跃迁矩阵元理论，振动跃迁 $v_i \to v_f$ 的红外吸收概率正比于：

在谐振子近似下，将偶极矩对正则坐标 $Q_k$ 做 Taylor 展开：

利用谐振子选择定则 $\Delta v = \pm 1$，只有展开式中一次项对跃迁有贡献，因此：

1.2 极化率与拉曼活性

拉曼散射是光子与分子振动的非弹性散射。 入射光（频率 $\nu_0$）诱导分子产生瞬时偶极矩：

其中 $\boldsymbol{\alpha}$ 是 $3\times3$ 极化率张量。 将极化率对正则坐标展开：

类似地，拉曼活性条件为：

1.3 互斥定则（Mutual Exclusion Rule）

2. 红外吸收强度的计算

2.1 偶极矩对正则坐标的导数

设分子有 $3N-6$（非线性）或 $3N-5$（线性）个简正模式， 第 $k$ 个简正模式的正则坐标为 $Q_k$，约化质量为 $\mu_k$（单位：amu）。 以质量加权笛卡尔位移 $\Delta\mathbf{x}$ 表示分子构型变化：

其中 $\mathbf{L}_k$ 是第 $k$ 个简正模式的质量加权本征向量 （Gaussian 输出中的 "normal mode eigenvectors"，单位为 amu$^{-1/2}$）。

偶极矩导数的数值差分近似（中心差分，误差 $O(\Delta Q^2)$）：

2.2 强度公式与单位

红外积分强度（单位：km/mol）由偶极矩三个方向分量的导数平方之和决定：

其中换算常数 $C_{\mathrm{IR}}$（将 Debye/(amu$^{1/2}$·Å) 换算为 km/mol）：

2.3 数值差分方法详解

中心差分法相比前向/后向差分精度更高，误差阶为 $O(\Delta Q^2)$：

步长 $\Delta Q$ 的选择影响计算精度：

3. 拉曼散射活性的计算

3.1 极化率张量

极化率张量 $\boldsymbol{\alpha}$ 是描述分子在外电场作用下 诱导偶极矩响应的 $3\times3$ 实对称矩阵：

由于对称性，只有 6 个独立分量。 Gaussian 输出按以下顺序给出（原子单位，1 a.u. = 1.48185$\times$10⁻²⁵ cm³）：

Polarizability (field-independent basis, a.u.):
   XX=    39.294   YX=    -0.049   YY=   156.727
   ZX=     0.000   ZY=     0.000   ZZ=   112.948

程序内部按顺序存储为： $[\alpha_{xx},\,\alpha_{yx},\,\alpha_{yy},\,\alpha_{zx},\,\alpha_{zy},\,\alpha_{zz}]$。

3.2 各向同性与各向异性导数

对极化率张量的每个独立分量计算对正则坐标 $Q_k$ 的导数：

由此定义两个关键量：

各向同性极化率导数（Mean polarizability derivative）

各向异性极化率导数（Anisotropy derivative）

3.3 拉曼活性公式

拉曼活性（Raman scattering activity）由各向同性和各向异性导数共同贡献：

3.4 去偏振比（Depolarization Ratio）

对于 $90°$ 散射，平面偏振激发光的去偏振比 $\rho_p$ 定义为：

4. Gaussian 数值差分计算流程

4.1 输入文件设置

步骤 1–2：几何优化 + 频率计算

%chk=H2O.chk
%mem=2GB
%nprocshared=4
#p b3lyp/6-311+g(d,p) opt freq=raman

H2O 几何优化与频率计算（含 Raman 活性）

0 1
 O  0.0   0.00000000  -0.11081188
 H  0.0  -0.78397589   0.44324751
 H  0.0   0.78397589   0.44324751

步骤 4：对每个正则模式生成 ±ΔQ 位移的单点输入

%chk=H2O_mode1_plus.chk
#p b3lyp/6-311+g(d,p) iop(1/21=1) scf=tight prop=(dipole,polar) nosymm

H2O mode 1  +DeltaQ  (正则坐标位移 +0.01 amu^(1/2)*Angstrom)

0 1
 O  [displaced coords]
 H  [displaced coords]
 H  [displaced coords]

4.2 数据提取

从 .log 文件提取偶极矩

Dipole moment (field-independent basis, Debye):
   X=     0.0000    Y=     0.0000    Z=    -0.0052  Tot=     0.0052

从 .log 文件提取极化率

Polarizability (field-independent basis, a.u.):
   XX=    39.294   YX=    -0.049   YY=   156.727
   ZX=     0.000   ZY=     0.000   ZZ=   112.948

Python 正则表达式提取

import re

def extract_dipole(logfile):
    """提取偶极矩 (Debye)，返回 [X, Y, Z]"""
    pattern = r'X=\s+([\d\.\-]+)\s+Y=\s+([\d\.\-]+)\s+Z=\s+([\d\.\-]+)'
    with open(logfile) as f:
        text = f.read()
    matches = re.findall(pattern, text)
    if matches:
        last = matches[-1]  # 取最后一次出现（收敛后的结果）
        return [float(x) for x in last]
    return None

def extract_polarizability(logfile):
    """提取极化率张量 (a.u.)，返回 [XX, YX, YY, ZX, ZY, ZZ]"""
    pattern = (r'XX=\s+([\d\.\-]+)\s+YX=\s+([\d\.\-]+)\s+YY=\s+([\d\.\-]+)\s+'
               r'ZX=\s+([\d\.\-]+)\s+ZY=\s+([\d\.\-]+)\s+ZZ=\s+([\d\.\-]+)')
    with open(logfile) as f:
        text = f.read()
    matches = re.findall(pattern, text)
    if matches:
        return [float(x) for x in matches[-1]]
    return None

4.3 H₂O 示例

H₂O 属于 $C_{2v}$ 点群，有 $3N-6 = 3$ 个简正模式，均具有 IR 和 Raman 活性：

5. Python 程序：NumIR_Raman

NumIR_Raman.py 位于 program/ 子目录， 实现了从 Gaussian 输出数据出发，通过数值差分计算 IR 和 Raman 强度的完整流程。

5.1 程序结构与主要函数

5.2 使用示例

计算 IR 强度

from NumIR_Raman import calculate_ir_intensity

# 从 Gaussian 输出读取的偶极矩分量 (Debye)
# 以 H2O mode1 为例（仅 z 分量有值）
dipole_plus  = [0.0, 0.0, -0.0052]   # Q + ΔQ
dipole_minus = [0.0, 0.0,  0.0052]   # Q - ΔQ
reduced_mass = 1.0803   # amu

ir = calculate_ir_intensity(dipole_plus, dipole_minus, reduced_mass)
print(f"IR 强度: {ir:.4f} km/mol")
# → IR 强度: 3.8762 km/mol

计算拉曼活性

from NumIR_Raman import calculate_raman_derivatives, calculate_raman_intensity

# 极化率张量独立分量 [XX, YX, YY, ZX, ZY, ZZ]（原子单位）
alpha_plus  = [39.294, -0.049, 156.727, 0.000, 0.000, 112.948]
alpha_minus = [39.294,  0.049, 156.727, 0.000, 0.000, 112.948]
reduced_mass = 1.0803  # amu

# 计算极化率导数
dalpha_dQ = calculate_raman_derivatives(alpha_plus, alpha_minus, reduced_mass)

# 计算拉曼活性
S_k, ai2, ga2 = calculate_raman_intensity(dalpha_dQ)
print(f"拉曼活性  S_k = {S_k:.4f} Å⁴/amu")
print(f"各向同性  α'² = {ai2:.4f}")
print(f"各向异性  γ'² = {ga2:.4f}")
# 去偏振比
rho = 3*ga2 / (45*ai2 + 4*ga2)
print(f"去偏振比  ρ   = {rho:.4f}")

直接运行（使用内置示例数据）

python program/NumIR_Raman.py

输出示例：

======================================================================
数值计算红外(IR)和拉曼(Raman)强度
======================================================================

红外强度计算结果
----------------------------------------------------------------------
模式        dμ/dQ (au)     IR强度 (km/mol)   相对强度
----------------------------------------------------------------------
mode1       -0.000527      2.8454            1.0000
mode6       -0.002384      24.3745           8.5671
mode14      -0.005692      138.6566          48.7392

5.3 物理常数与单位换算

IR 强度换算因子

6. 常见问题与调试

Q1：计算的 IR 强度与 Gaussian 直接输出不一致？

Q2：拉曼活性数值偏差较大？

Q3：某振动模式 IR 强度为 0，但预期有活性？

Q4：如何快速验证 Gaussian 频率计算中的 IR 强度？