分子转动光谱模拟

1. 环上粒子模型（Particle-on-a-ring）

环上粒子模型（Particle-on-a-ring）是箱中粒子（Particle-in-a-box）在 环形构型下的推广，是描述分子转动量子化的最简模型。 粒子（质量 $m$）被限制在半径为 $r$ 的一维环上运动， 以方位角 $\phi$ 描述其位置。

1.1 薛定谔方程

由线性薛定谔方程出发，令 $x = r\phi$（弧长），$\phi$ 为角坐标：

令 $b^2 = \dfrac{2mr^2 E}{\hbar^2}$，方程化为标准形式：

通解为：

1.2 边界条件与量子数

环形结构要求波函数满足周期性边界条件（Periodic BC）： 粒子绕环一周 $2\pi$ 后回到出发点，波函数值不变：

这意味着波函数在整个圆周上恰好容纳整数个完整振荡（de Broglie 驻波条件）， 因此 $b$ 必须为整数。令 $l = b$，得到转动量子数：

1.3 能级与波函数

将 $l$ 代入能量表达式，得到能级公式：

其特征：

• $l=0$：$E_0 = 0$，非简并
• $l \neq 0$：$+l$ 与 $-l$ 能量相同，二重简并（顺/逆时针转动）
• 能级间隔随 $|l|$ 增大而增大：$\Delta E \propto l$

归一化波函数（取实数形式）：

$\delta$ 为相位常数。取复数形式更常用： $\psi_l(\phi) = \dfrac{1}{\sqrt{2\pi}}e^{il\phi}$。

1.4 与箱中粒子的对比

2. 固态物理拓展——Born–von Kármán 边界条件

3. 分子转动的角动量理论

3.1 角动量与惯量张量

对由 $N$ 个原子组成的分子，第 $k$ 个原子的角动量贡献为：

对刚性分子，整体以角速度 $\boldsymbol{\omega}$ 转动时， 各原子速度为 $\vec{v}_k = \boldsymbol{\omega}\times\vec{r}_k$， 总角动量为：

其中 $\mathbf{I}$ 为惯量矩阵（Inertial Matrix），$3\times3$ 对称矩阵：

对角元素（转动惯量 Moments of Inertia）

非对角元素（惯量积 Products of Inertia）

3.1b 惯量积的物理意义

惯量矩阵的非对角元 $I_{\alpha\beta}$（$\alpha \neq \beta$）称为 惯量积（Products of Inertia），定义为：

以 $I_{xy}$ 为例：若分子中有质量偏向 $x > 0,\, y > 0$ 和 $x < 0,\, y < 0$ 两个象限，则 $\sum_k m_k x_k y_k > 0$，$I_{xy} < 0$，惯量积非零。

惯量积的物理含义

惯量积反映质量分布相对于坐标轴的非对称性。 其物理后果是：当分子绕某坐标轴转动时，如果惯量积不为零， 产生的角动量方向会偏离转动轴。

例如，若仅有 $\omega_z \neq 0$，但 $I_{xz} \neq 0$，则 $L_x = I_{xz}\omega_z \neq 0$—— 角动量在 $x$ 方向出现分量，导致转动轴"进动"，这正是陀螺偏转的根源。

何时惯量积为零？

若分子对某坐标平面具有对称性，则惯量积自动消失：

只有当坐标轴方向恰好对应惯量主轴时，惯量矩阵才变为对角矩阵， 此时转动轴与角动量轴重合，运动最简单。 这就是需要将分子坐标变换到惯性主轴坐标系的动机。

3.1c 变换到惯性主轴坐标系

对任意初始坐标下的分子，通过以下两步变换，可将坐标系转化为 惯性主轴坐标系（Principal Axis System, PAS）： 质心在原点，三个坐标轴分别指向惯量主轴， 且沿 $x$、$y$、$z$ 方向的转动惯量依次增大（$I_x \leq I_y \leq I_z$）。

第一步：平移——将质心移至原点

计算分子的质心坐标：

对所有原子做平移：

平移后，$\sum_k m_k \vec{r}_k' = \vec{0}$，质心在新坐标原点。 此时用 $\vec{r}_k'$ 代入惯量公式，计算得到质心系惯量矩阵 $\mathbf{I}'$。

第二步：旋转——将坐标轴转到惯量主轴方向

对质心系惯量矩阵 $\mathbf{I}'$（实对称矩阵）求解本征值问题：

得到三个本征值 $I_1 \leq I_2 \leq I_3$（即主转动惯量） 和对应的单位本征向量 $\mathbf{n}_1,\,\mathbf{n}_2,\,\mathbf{n}_3$（即主轴方向）。 由于 $\mathbf{I}'$ 是实对称矩阵，三个本征向量两两正交： $\mathbf{n}_i \cdot \mathbf{n}_j = \delta_{ij}$。

构造旋转矩阵 $\mathbf{R}$，以本征向量为行：

对每个原子施加旋转变换：

此时惯量矩阵变为对角形式：

变换结果：惯性主轴坐标系的性质

数值实现（Python 示例）

import numpy as np

def to_principal_axis_system(coords, masses):
    """
    coords : (N, 3) array，原始笛卡尔坐标（Å）
    masses : (N,)   array，原子质量（amu）
    返回：新坐标 coords_pas，主转动惯量 [Ix, Iy, Iz]（amu·Å²）
    """
    M = masses.sum()

    # Step 1: 平移到质心
    r_cm = (masses[:, None] * coords).sum(axis=0) / M
    coords_cm = coords - r_cm

    # Step 2: 构建惯量矩阵
    I = np.zeros((3, 3))
    for m, r in zip(masses, coords_cm):
        I += m * (np.dot(r, r) * np.eye(3) - np.outer(r, r))

    # Step 3: 对角化，取得本征值（主转动惯量）和本征向量（主轴）
    eigenvalues, eigenvectors = np.linalg.eigh(I)
    # eigenvalues 已按升序排列：Ix ≤ Iy ≤ Iz
    # eigenvectors 的列为本征向量

    # Step 4: 旋转矩阵（本征向量转置 = 行为主轴方向）
    R = eigenvectors.T          # shape (3, 3)

    # Step 5: 变换坐标
    coords_pas = coords_cm @ eigenvectors  # 等价于 (R @ coords_cm.T).T

    return coords_pas, eigenvalues  # eigenvalues = [Ix, Iy, Iz]

3.2 转动动能

对惯量矩阵进行对角化（选取主轴坐标系）， 可使所有非对角元为零，得到三个主转动惯量 $I_A \leq I_B \leq I_C$：

对角化等价于将分子（或坐标轴）旋转至惯量主轴方向， 这一变换通过求惯量矩阵的本征向量来实现。

3.3 转动常数（Rotational Constants）

转动能级

对线性分子（或一般分子沿某一主轴 $B$）， 将转动动能用角动量量子数 $J$ 表达，得到量子化的转动能级：

相邻能级之间的跃迁能量（$J \to J+1$）为：

由能量到频率——转动常数的来历

光谱实验中观测到的是电磁波的频率而非能量。 根据 de Broglie（普朗克-爱因斯坦）关系 $E = h\nu$， 将转动能量直接除以普朗克常数 $h$，即得到对应的频率：

将其中与 $J$ 无关的常数部分单独提出，定义为转动常数 $B$：

于是转动能级对应的频率简洁地写为：

类似地，对三个主轴各定义一个转动常数：

也可转换为波数（cm⁻¹）表示，即用光速 $c$ 除以频率：

4. 用 Gaussian 计算转动常数

Gaussian 是量子化学领域最广泛使用的计算软件。 通过几何优化或单点计算，Gaussian 会自动输出分子的转动常数。

4.1 Gaussian 输入文件格式

Gaussian 输入文件（.gjf 或 .com）的基本结构：

%chk=filename.chk          ← 检查点文件（可选，用于后续读取）
%mem=2GB                   ← 内存（可选）
%nprocshared=4             ← 处理器数（可选）
#p <方法>/<基组> <关键词>   ← 路由段（Route section）

分子标题（任意文字）        ← 空行分隔

电荷  自旋多重度            ← 如 0 1 = 中性单重态
元素符号  X  Y  Z          ← 原子坐标（Å）
...
                           ← 末尾空行（必须！）

常用计算方法与基组

获取转动常数的关键词

单点能计算（#p hf/6-31g）默认输出转动常数。 若需几何优化后给出，使用：

#p b3lyp/6-311+g(d,p) opt freq

opt：优化几何构型； freq：计算振动频率（同时输出零点能和更精确的转动常数）。

4.2 示例一：H₂ 分子转动常数

输入文件

%chk=H2.chk
#p hf/6-31g

H2

0 1
 H  0.0  0.0  0.3
 H  0.0  0.0 -0.3

Gaussian 输出结果

Rotational constants (GHZ):
   0.0000000        2785.8616915        2785.8616915
AtmWgt=  1.0078250  1.0078250

结果解读

4.3 示例二：H₂O 分子转动常数

输入文件

%chk=H2O.chk
#p hf/6-31g

H2O

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

Gaussian 输出结果

Rotational constants (GHZ):
   919.6759406      407.9403341      282.5913764
AtmWgt=  15.9949146   1.0078250   1.0078250

结果解读

4.4 Gaussian 输出文件关键信息

在 Gaussian 日志文件（.log）中，搜索以下关键词：

5. 编程作业：任意分子转动常数计算

从 Gaussian 日志文件中读取原子坐标，手动计算惯量矩阵并求转动常数， 与 Gaussian 输出验证。

1. 从 Gaussian 日志读取原子坐标
   在 .log 文件中找到 Standard orientation: 段落， 提取原子序号、原子类型、$x,y,z$ 坐标（Å）及原子质量（AtmWgt）。

   import re
   
   def read_gaussian_log(filename):
       atoms = []
       with open(filename) as f:
           lines = f.readlines()
       # 找 Standard orientation 块
       for i, line in enumerate(lines):
           if "Standard orientation" in line:
               start = i + 5
               for j in range(start, len(lines)):
                   if "---" in lines[j]:
                       break
                   parts = lines[j].split()
                   if len(parts) == 6:
                       atoms.append({
                           'Z': int(parts[1]),
                           'x': float(parts[3]),
                           'y': float(parts[4]),
                           'z': float(parts[5])
                       })
       return atoms
   

2. 平移坐标原点到质心
   计算质心坐标，将所有原子坐标减去质心： $$ \vec{r}_{\rm cm} = \frac{\sum_k m_k\,\vec{r}_k}{\sum_k m_k}, \qquad \vec{r}_k' = \vec{r}_k - \vec{r}_{\rm cm} $$

   import numpy as np
   
   def to_center_of_mass(atoms, masses):
       coords = np.array([[a['x'], a['y'], a['z']] for a in atoms])
       m = np.array(masses)
       r_cm = np.sum(m[:, None] * coords, axis=0) / np.sum(m)
       return coords - r_cm
   

3. 构建惯量矩阵并对角化（转到主轴系）
   $$ I_{\alpha\beta} = \sum_k m_k\bigl(|\vec{r}_k'|^2\delta_{\alpha\beta} - r_{k\alpha}' r_{k\beta}'\bigr) $$

   def inertia_tensor(coords, masses):
       I = np.zeros((3, 3))
       for m, r in zip(masses, coords):
           r2 = np.dot(r, r)
           I += m * (r2 * np.eye(3) - np.outer(r, r))
       return I
   
   # 对角化得主转动惯量
   eigenvalues, eigenvectors = np.linalg.eigh(I)
   # eigenvalues 即 I_A, I_B, I_C（amu·Å²，需换算 SI 单位）
   

4. 计算转动常数
   $$ B = \frac{h}{8\pi^2 I_B} \quad\text{(Hz)} \qquad B[\text{GHz}] = \frac{h[\text{J·s}]}{8\pi^2\,I_B[\text{kg·m}^2]} \times 10^{-9} $$

   h  = 6.62607015e-34   # J·s
   amu_to_kg = 1.66053906660e-27
   A_to_m    = 1e-10
   
   def rotational_constants_GHz(I_eigenvalues_amu_A2):
       consts = []
       for I in I_eigenvalues_amu_A2:
           if I < 1e-10:
               consts.append(0.0)
           else:
               I_SI = I * amu_to_kg * A_to_m**2
               B = h / (8 * np.pi**2 * I_SI) * 1e-9  # GHz
               consts.append(B)
       return sorted(consts, reverse=True)  # A >= B >= C