分子波函数与振动光谱模拟

1. 分子薛定谔方程

1.1 全分子哈密顿量

对含 $M$ 个原子核（质量 $M_I$，电荷 $Z_I$）和 $N$ 个电子的分子体系， 在非相对论、时间无关的近似下，全哈密顿量（原子单位制）为：

1.2 Born-Oppenheimer 近似

由于原子核质量远大于电子质量（$M_I \gg m_e$），核的运动远慢于电子， 可将全波函数分离为电子部分和核部分的乘积：

分离后得到两个独立方程：

1.3 电子薛定谔方程与电子哈密顿量

在固定核构型下，电子哈密顿量为：

求解电子薛定谔方程可得到：

• 电子能量 $E_e(\mathbf{R})$ 和波函数 $\psi_e$
• 电子密度分布、静电势、偶极矩
• 电离能（Ionization Potential）
• 分子光谱（UV-Vis、IR、Raman）
• 势能面（PES）与优化几何构型

2. 势能面扫描 — Gaussian 实践

通过系统地改变分子内坐标（键长、键角）并在每个点做单点能计算， 可以绘制绝热势能面。 Gaussian 的 scan 关键词可自动完成这一任务。

2.1 H₂O PES 扫描输入

%chk=test0076
#p rhf/3-21g scan

Water RHF scan

0,1
o
h,1,r
h,1,r,2,a

r .96 2 .05       ← O-H 键长从 0.96 Å 开始，扫描 2 步，步长 0.05 Å
a 104.5 3 1.0     ← H-O-H 键角从 104.5° 开始，扫描 3 步，步长 1.0°

2.2 断点续算（Restart）

对于耗时的扫描任务，若计算中断可使用 scan=restart 从断点继续：

%chk=test0076_break_4
%kjob l301 4          ← 在第 4 个 l301 作业后终止（模拟中断）
#p rhf/3-21g scan

Water RHF scan

0,1
o
h,1,r
h,1,r,2,a

r .96 2 .05
a 104.5 3 1.0

--Link1--
%chk=test0076_break_4
%nosave
#p rhf/3-21g scan=restart   ← 从检查点文件续算

3. 谐振子近似与振动频率

3.1 泰勒展开与力常数矩阵

在平衡构型 $\mathbf{R}_0$ 附近对势能面做泰勒展开（Born-Oppenheimer 绝热近似下）， 忽略三阶及以上项（谐振子近似）：

其中 $q_i = R_i - R_{i,0}$ 为原子位移，$k_{ij}$ 为力常数矩阵 元素（Hessian 矩阵，Gaussian 关键词 freq 计算此矩阵）。 一阶项为零（平衡点梯度为零）。

牛顿运动方程（矩阵形式）：

其中 $\mathbf{M} = \text{diag}(m_1,m_1,m_1,m_2,\ldots)$ 为质量矩阵， $\mathbf{K}$ 为力常数矩阵。 方程耦合（$\mathbf{K}$ 非对角），需通过对角化解耦。

3.2 质量加权坐标

引入质量加权位移坐标（Mass-weighted displacement）：

定义质量加权力常数矩阵（Mass-weighted Hessian）：

运动方程化为标准形式：

3.3 简正模式与对角化

对 $\mathbf{F}$ 进行正交对角化，求本征值与本征向量：

$\mathbf{L}$ 的各列为本征向量（简正模式方向）， $\boldsymbol{\Lambda}$ 的对角元为本征值。 定义简正坐标（Normal coordinates）：

每个简正坐标满足独立的谐振子方程，振动圆频率为：

3.3b 简正模式的物理图像

简正模式（Normal Mode）是分子集体振动的基本单元： 在一个简正模式中，所有原子以相同的频率和固定的相位关系 同时运动，各原子位移的相对大小和方向由本征向量 $\mathbf{L}_k$ 决定。

零频模式（Zero-frequency Modes）

对整个分子体系，$3N$ 个自由度中有若干模式对应 整体平动和转动，其势能不随运动改变， 故力常数为零，频率为零（或数值上因舍入误差出现极小虚频）：

剩余 $3N-6$（非线性）或 $3N-5$（线性）个模式为真正的振动简正模， 频率均为正实数。 Gaussian 在输出时会自动剔除这些平动/转动模式， 只报告振动频率。

H₂O 的三个简正模式

以 H₂O（$N=3$，$C_{2v}$ 对称性）为例，三个简正模式为：

3.3c Gaussian 日志中的简正模式数据

Gaussian freq 计算完成后，日志文件（.log） 中包含完整的简正模式信息。以下是典型输出块的逐行注释：

                      1                      2                      3
  ← 模式编号（振动模式按频率从低到高排列）

                   A1                     A1                     B2
  ← 对称性标签（点群不可约表示）；A1 = 全对称，B2 = 非全对称

 Frequencies --  1622.8153             3811.3037             3912.7373
  ← 振动频率（cm⁻¹）；负值表示虚频（过渡态或未收敛）

 Red. masses --     1.0824                1.0440                1.0812
  ← 约化质量 μk（amu）；= (Σ Lik²/mi)⁻¹，决定零点能和同位素效应

 Frc consts  --     1.6770                0.8938                0.9736
  ← 有效力常数 k̃ = (2πcν̃)² μk（mDyne/Å = 100 N/m）；反映键的强弱

 IR Inten    --    72.8370                5.3941               52.8490
  ← 红外强度（km/mol）；正比于 |∂μ/∂Qk|²；0 = 红外非活性

  Atom  AN    X      Y      Z       X      Y      Z       X      Y      Z
  ← 各原子在该简正模式下的笛卡尔位移向量（列归一化，无量纲）
     1   8  0.00   0.00   0.07    0.00   0.00  -0.07    0.00   0.07   0.00
  ← O（原子序号8）：模式1中沿 z 轴位移 0.07；模式3中沿 y 轴位移 0.07
     2   1  0.00   0.58  -0.56    0.00   0.58   0.44    0.00   0.41   0.58
  ← H₁（原子序号1）：位移向量较大，轻原子运动幅度大
     3   1  0.00  -0.58  -0.56    0.00  -0.58   0.44    0.00  -0.41   0.58
  ← H₂（原子序号1）：与 H₁ 对称

位移向量的含义与归一化约定

日志中输出的每列位移向量 $\mathbf{d}_k$（共 $3N$ 个分量，每原子 3 个） 采用列归一化（Column-normalized）约定：

它与质量加权本征向量 $\mathbf{L}_k$ 的关系为：

从位移向量到原子运动幅度

对于热振动，第 $k$ 个模式的经典振动幅度正比于 $1/\sqrt{\mu_k\omega_k^2}$。 原子 $i$ 的实际位移幅度（Å）为：

轻原子（H）的 $d_{k,i}/\sqrt{m_i}$ 较大， 所以在动画中氢原子的箭头总是比重原子长—— 这是质量效应的直接体现。

如何在日志文件中快速定位简正模式数据

# Linux/Mac
grep -A 50 "Frequencies --" water.log | head -80

# Windows PowerShell
Select-String -Path water.log -Pattern "Frequencies" -Context 0,15

# Python
import re
with open("water.log") as f:
    text = f.read()
blocks = re.findall(r"Frequencies --.*?(?=\nFrequencies|\Z)", text, re.DOTALL)
for b in blocks:
    print(b[:300])

3.4 归一化坐标与约化质量

质量加权坐标 $\tilde{Q}_k$ 有量纲（$\text{amu}^{1/2}\cdot\text{Å}$）， 非无量纲也非归一化。为此引入约化质量（Reduced mass）$\mu_k$ 和归一化简正坐标 $\xi_k$：

对于双原子分子，约化质量退化为熟悉的形式：

正交变换 vs 非正交变换

4. Gaussian 频率计算实践

4.1 输入文件示例 — H₂O 优化 + 频率

%chk=water
#p opt freq rhf/3-21g

Water RHF

0,1
o
h,1,r
h,1,r,2,a

r .96
a 104.5

4.2 其他分子计算任务

4.3 输出文件关键词

5. 光物理过程与 Jablonski 图

5.1 Jablonski 能级图

Jablonski 图描述分子吸收光子后在不同电子态之间的跃迁路径：

S2  ──────────────────────────────  (高激发单重态)
        │   IC（内转换，快速）
S1  ──────────────────────────────  (第一激发单重态)
        │\__  Fluorescence     │ ISC（系间窜越）
        │    (荧光, ~ns)       ↓
S0  ──────────────────   T1  ──────── (三重态)
  (基态)    ↑ Absorption           │ Phosphorescence
            hν                     │ (磷光, ~μs-ms)
                                   ↓
                              S0 ──────────────

量子产率（Quantum Yield）

$k_r$ 为辐射速率常数，$k_{nr}$ 为非辐射速率常数之和（包括 IC、ISC、淬灭等）。 $\Phi_f \in [0, 1]$，越大表示发光效率越高。

5.2 跃迁偶极矩与 Franck-Condon 展开

两个电子态 $i\to f$ 之间跃迁的跃迁偶极矩：

在 Born-Oppenheimer 近似下，$\Psi = \psi_e\chi_n$， 对偶极算符关于简正坐标做泰勒展开（Herzberg-Teller 展开）：

跃迁偶极矩分解为：

5.3 爱因斯坦系数与辐射速率

对于二能级体系（能量差 $\hbar\omega = E_2 - E_1$）， Einstein 于 1917 年引入三个系数描述辐射跃迁：

辐射速率常数 $k_r = A_{21}$，自发辐射的荧光寿命 $\tau_r = 1/A_{21}$。

6. 红外光谱（Infrared Spectrum）

6.1 偶极矩（Electric Dipole Moment）

偶极矩描述分子中正负电荷中心的分离程度， 是决定分子能否与红外光相互作用的核心物理量。 对含 $M$ 个原子核、$N$ 个电子的分子，偶极矩算符为：

量子力学期望值（即分子的永久偶极矩）：

其中 $\rho(\vec{r}) = N\int|\psi_e|^2\,d\vec{r}_2\cdots d\vec{r}_N$ 是电子密度。

Gaussian 如何计算偶极矩

在 HF 或 DFT 计算中，Gaussian 将波函数展开为基函数， 偶极矩通过下式计算：

其中 $P_{\mu\nu}$ 是密度矩阵，$\langle\phi_\mu|r_\alpha|\phi_\nu\rangle$ 是基函数的偶极矩积分（解析可求）。 Gaussian 的任何单点计算均自动输出偶极矩， 无需额外关键词：

#p b3lyp/6-311+g(d,p)        ← 输出中含 Dipole moment (Debye)

#p b3lyp/6-311+g(d,p) polar  ← 同时输出极化率张量（见第 7 节）

6.2 红外光谱基本原理与选择定则

红外光谱源于分子振动能级之间的跃迁。 当红外光子的频率与某振动模式的频率共振时（$h\nu_{\rm photon} = \hbar\omega_k$）， 分子吸收光子，振动量子数改变：$v = 0 \to 1$。

物理图像

分子沿简正坐标 $Q_k$ 振动时，若偶极矩随之变化， 振荡的偶极矩就像一根"振荡天线"， 能与同频率的电磁波（红外光）发生能量交换。 因此红外活性的选择定则为：

量子力学上，跃迁矩阵元为：

当 $(\partial\mu_\alpha/\partial Q_k)_0 = 0$ 时矩阵元为零，跃迁禁戒， 该模式红外非活性。

6.3 偶极矩导数：沿简正坐标的有限差分

Gaussian 在频率计算中对每个简正模式 $k$ 计算偶极矩导数 $(\partial\vec{\mu}/\partial Q_k)_0$。 GaussView 中可以动画演示的"沿简正坐标移动原子"， 正是这一计算的直观体现。

简正坐标位移与笛卡尔位移的关系

第 $k$ 个简正模式的本征向量 $\mathbf{L}_k$（质量加权，维数 $3N$） 给出了各原子的振动方向。 将简正坐标从平衡值 $Q_k = 0$ 移动到 $Q_k = \pm\delta$ 时， 第 $i$ 个原子的笛卡尔坐标变化为：

其中 $\mathbf{L}_k^{(i)}$ 是本征向量中对应原子 $i$ 的 $(x,y,z)$ 三分量， $m_i$ 为该原子质量，$\delta$ 为步长（通常约 $0.01\,\text{Å}\cdot\text{amu}^{1/2}$）。

有限差分计算偶极矩导数

对 $3N-6$（非线性）或 $3N-5$（线性）个简正模式各做一次 $\pm$ 位移，共需 $2(3N-6)$ 个单点能计算。 这正是 Gaussian 频率计算中 Numerical frequencies using 2 points per step 的含义。

6.4 红外强度公式

模式 $k$ 的红外吸收强度（积分摩尔吸收系数）为：

Gaussian 输出的 IR Intensities 单位为 km/mol， 正是上式乘以适当换算因子后的结果， 它直接对应实验中的摩尔消光系数（积分面积）。

6.5 Gaussian 计算与输出解读

%chk=water.chk
#p b3lyp/6-311+g(d,p) opt freq

Water B3LYP frequency

0,1
O   0.000000   0.000000   0.119748
H   0.000000   0.756950  -0.478993
H   0.000000  -0.756950  -0.478993

日志文件中频率输出块示例（H₂O 的三个模式）：

                     1                      2                      3
                  A1                     A1                     B2
 Frequencies --  1622.8153             3811.3037             3912.7373   (cm⁻¹)
 Red. masses --     1.0824                1.0440                1.0812   (amu)
 Frc consts  --     1.6770               0.8938                0.9736   (mDyne/Å)
 IR Inten    --    72.8370               5.3941               52.8490   (km/mol)

7. 拉曼光谱（Raman Spectrum）

7.1 极化率张量（Polarizability Tensor）

在外电场 $\mathbf{E}$ 作用下，分子的电子云发生形变，产生感应偶极矩：

$\boldsymbol{\alpha}$ 是 $3\times3$ 的极化率张量（对称，$\alpha_{ij}=\alpha_{ji}$）， 在分子主轴系中为对角形式。 其六个独立分量为：

两个常用的不变量：

Gaussian 如何计算极化率

极化率是能量对外电场的二阶导数：

Gaussian 提供两种计算途径：

#p b3lyp/6-311+g(d,p) polar    ← 静态极化率

输出示例：

 Dipole polarizability, Alpha (input orientation):
 Eigenvalues:   5.8312   8.6045  10.2334
 Isotropic =    8.2230   Anisotropy =    3.9751

 Alpha(1,1)=  10.2189  Alpha(2,1)=   0.0000  Alpha(3,1)=   0.0000
 Alpha(2,2)=   8.6045  Alpha(3,2)=   0.0000  Alpha(1,3)=   0.0000
 Alpha(3,3)=   5.8312                                             (Bohr³)

7.2 拉曼散射原理与选择定则

拉曼散射基于光的非弹性散射。 当频率为 $\nu_0$ 的激光照射分子时， 分子极化率随振动发生周期性变化， 导致感应偶极矩中出现新的频率分量。

经典图像

将极化率对简正坐标 $Q_k$ 展开，并令 $Q_k = Q_k^0\cos(\omega_k t)$：

感应偶极矩 $\mu_i^{\rm ind} = \alpha_{ij} E_j^0\cos(2\pi\nu_0 t)$ 中包含：

因此拉曼活性的选择定则为：

7.3 极化率导数：沿简正坐标的有限差分

Raman 活性要求计算每个简正模式的极化率导数张量 $(\partial\boldsymbol{\alpha}/\partial Q_k)_0$， 这是一个 $3\times3$ 对称矩阵（6 个独立元素）。 GaussView 中"沿简正模式移动原子"就是为了数值地估算这个导数。

有限差分步骤

每个模式需 2 次极化率计算；对 $3N-6$ 个模式，共需 $2(3N-6)$ 次 polar 计算。 这正是 freq=raman 比普通 freq 耗时更多的原因。

7.4 拉曼活性与强度公式

由极化率导数张量 $\boldsymbol{\alpha}'$ 构造两个旋转不变量：

Gaussian 输出的 Raman Activity（单位 Å⁴/amu）定义为：

实验可观测的差分散射截面为：

最后一项为玻尔兹曼热修正因子， 在室温下对高频模式接近 1，对低频模式 Anti-Stokes 峰有较大增强。

此外，退偏振比（Depolarization Ratio）定义为：

7.5 Gaussian 计算与输出解读

%chk=water_raman.chk
#p b3lyp/6-311+g(d,p) opt freq=raman

Water Raman calculation

0,1
O   0.000000   0.000000   0.119748
H   0.000000   0.756950  -0.478993
H   0.000000  -0.756950  -0.478993

日志文件中拉曼输出示例：

                     1                      2                      3
                  A1                     A1                     B2
 Frequencies --  1622.8153             3811.3037             3912.7373
 IR Inten    --    72.8370               5.3941               52.8490   (km/mol)
 Raman Activ --     4.7251              95.1832               33.6279   (Å⁴/amu)
 Depolar (P) --     0.1532               0.1137                0.7500
 Depolar (U) --     0.2654               0.2041                0.8571