3. 力常数与动力学矩阵构建

在上一章中，我们把晶格振动的物理图像归结为“原子在平衡位置附近做微小位移，回复力由原子间势能曲率决定”。要把这一图像变成可计算的声子谱，核心任务只有两步：第一，把势能曲率写成可操作的力常数；第二，把力常数组装成动力学矩阵，再对角化得到频率与特征向量。本章围绕这条主线展开，重点讲清楚力常数的定义、对称性如何约束力常数、动力学矩阵如何组装，以及频率和特征向量如何求解与判读。

一、力常数的定义：从势能展开到回复力

设晶体中第 l 个原胞内第 κ 个原子沿 α 方向的位移为 u_α(lκ)。将体系总势能 Φ 在平衡位置附近做泰勒展开，一次项因平衡条件为零，二阶项为：

Φ = Φ₀ + ½ Σ Φ_{αβ}(lκ, l'κ') u_α(lκ) u_β(l'κ') + O(u³)

其中二阶导数 Φ_αβ(lκ, l'κ') 就是力常数，它描述“l'κ' 原子沿 β 方向位移一个单位时，lκ 原子沿 α 方向受到的力”。更常用的写法是力常数矩阵：

Φ_{αβ}(lκ, l'κ') = ∂²Φ / ∂u_α(lκ) ∂u_β(l'κ')

由定义可直接得到两条重要性质。其一是对称性：Φ_αβ(lκ, l'κ') = Φ_βα(l'κ', lκ)，即交换原子对与方向指标后矩阵转置相等。其二是平移不变性：晶体整体刚性平移时原子间相对位置不变，因此每个原子所受合力为零，这给出“声学求和规则”：

Σ_{l'κ'} Φ_{αβ}(lκ, l'κ') = 0

这条规则在数值计算中极其重要。如果从第一性原理或经验势得到的力常数不满足它，动力学矩阵在 Γ 点会出现非零的声学频率，通常表现为三个接近零但不为零的虚频或小实频，提示计算存在截断误差或对称性处理不当。

二、对称性约束：把独立力常数压到最少

晶体点群操作 R 作用在原子坐标和位移分量上，力常数必须满足：

Φ_{αβ}(lκ, l'κ') = Σ_{γδ} R_{αγ} R_{βδ} Φ_{γδ}(R⁻¹lκ, R⁻¹l'κ')

这条关系意味着：一旦知道某个“代表原子对”的力常数，就可以通过对称操作生成所有等价原子对的力常数。以立方晶系为例，沿 x、y、z 三个方向的最近邻力常数张量只有两个独立分量：纵向分量 Φ_∥ 和横向分量 Φ_⊥。若不做对称性约化，一个含 8 个原子的原胞、考虑 5 个近邻壳层，独立力常数可达数千个；利用点群与平移对称性后，往往只剩几十个甚至十几个独立参数。

实际操作中，对称性约束通常通过以下方式落地：先由空间群操作生成所有等价原子对，再对每个等价类只保留一个代表，最后用最小二乘或线性约束求解。这样既减少计算量，也保证动力学矩阵严格满足晶体对称性，避免出现“物理上不该有的分裂”。

三、动力学矩阵组装：从力常数到本征问题

对每个波矢 q，定义动力学矩阵：

D_{αβ}(κκ' | q) = (1/√(m_κ m_κ')) Σ_{l'} Φ_{αβ}(0κ, l'κ') exp[i q · (r(l'κ') − r(0κ))]

其中 m_κ 为原子质量，求和遍历所有原胞 l'。由于力常数满足平移不变性，D(q) 是厄米矩阵，其本征方程为：

D(q) e(qν) = ω²(qν) e(qν)

本征值 ω²(qν) 给出第 ν 支声子频率，本征向量 e(qν) 给出各原子在该模式下的位移方向与相对相位。组装时需要注意三点：一是质量加权必须与力常数定义一致；二是相位因子中的原子位置要用实际坐标，不能只取原胞内坐标；三是对于极性晶体，长程库仑相互作用需要在动力学矩阵中单独加入非解析项，否则 LO-TO 分裂会缺失。

四、频率与特征向量的求解与判读

对每个 q 点，动力学矩阵是 3N × 3N 的厄米矩阵（N 为原胞内原子数），对角化后得到 3N 个本征值。若 ω² 为正，频率为实数，对应稳定振动模式；若出现负值，则 ω 为虚数，表示该模式在动力学上不稳定，晶体可能发生结构相变。因此，检查全布里渊区内是否存在虚频，是判断结构稳定性的第一道关卡。

特征向量 e(qν) 是一个 3N 维复向量，其分量给出各原子在模式 ν 中的位移方向与相对相位。以双原子链为例，在 Γ 点附近，声学支的特征向量中两个原子同向运动，光学支中两个原子反向运动。通过观察特征向量，可以判断某个模式是“键伸缩型”还是“键弯曲型”，也可以识别出哪些原子参与振动最强烈。对于极性材料，纵向光学支的特征向量还会伴随宏观电场，需要在计算中单独处理。

五、一个可执行的最小示例

下面给出一段用 Python 组装一维双原子链动力学矩阵并求解的示例代码。该示例仅用于说明流程，参数为教学假设值，不代表任何真实材料。

import numpy as np

# 教学假设参数：一维双原子链，晶格常数 a=1，质量 m1=1.0, m2=2.0
# 力常数：最近邻弹簧 K=1.0，次近邻忽略
a = 1.0
m1, m2 = 1.0, 2.0
K = 1.0

def dynamical_matrix(q):
    # 2x2 动力学矩阵，基元含两个原子
    D = np.zeros((2, 2), dtype=complex)
    # 对角项：每个原子与左右最近邻的力常数之和
    D[0, 0] = 2 * K / m1
    D[1, 1] = 2 * K / m2
    # 非对角项：相位因子 exp(i q a)
    D[0, 1] = -K / np.sqrt(m1 * m2) * (1 + np.exp(-1j * q * a))
    D[1, 0] = np.conj(D[0, 1])
    return D

for q in [0.0, 0.5 * np.pi / a, np.pi / a]:
    D = dynamical_matrix(q)
    eigvals, eigvecs = np.linalg.eigh(D)
    freqs = np.sqrt(np.abs(eigvals)) * np.sign(eigvals)
    print(f"q={q:.3f}, omega={freqs}")

运行后可以观察到：在 q=0 处，两个本征值分别为 0 和 √(2K(1/m₁+1/m₂))，对应声学支和光学支；在 q=π/a 处，声学支频率升高，光学支频率降低，两支之间出现带隙。若把 K 改为负值或破坏声学求和规则，Γ 点会出现非零虚频，提示模型不稳定。

总结来说，力常数是晶格动力学的“输入参数”，对称性约束是“压缩参数空间”的关键，动力学矩阵是“组装平台”，而频率与特征向量是最终“输出结果”。四者环环相扣：力常数定义不清，后续矩阵就没有物理意义；对称性处理不当，矩阵会出现虚假分裂；组装时质量加权或相位因子写错，色散曲线就会整体偏移；求解后不检查虚频和声学模，就无法判断结构是否稳定。掌握这条链路，是后续声子谱计算与热力学性质分析的基础。