绪论：计算物理的世界 & 一维 Schrödinger 方程的数值解

本讲导读

这是课程的第一讲，要同时完成两件看起来不太相干的事：

• 建立课上要用的世界观——计算物理到底能算什么、算到什么精度、为什么有时候算不动；
• 上手第一个真正的数值方法——用差分法把一维 Schrödinger 方程变成矩阵问题来求解。

这两件事其实是连着的：第二部分提到的所有宏大应用（OLED 发光效率、载流子迁移率、光谱预测）， 最终都要落回「解一个方程」这个动作上。而本讲给出的差分法，是最简单、最能亲手写出来的那一种解法。

一、课程概览与定位对应讲义 PPT 第 1–3 页

1.1 课程信息

1.2 课程目标

课程描述里给出的三条目标，可以概括成「讲方法、动手做、自己探索」：

1.3 这门课会学什么

二、计算物理能做什么对应讲义 PPT 第 4–23 页

2.1 算力与精度的矛盾

回顾计算工具的演化线：穿孔卡片机 → I386 个人计算机 → 计算机集群 → AI for Science。 这条线之所以重要，是因为它一直在推动同一件事：把「算不动」变成「算得动」。

在电子计算机出现之前（1950 年代以前），能用高精度方法处理的体系数量极其有限。 六七十年代，电子计算机虽然出现了，但那时是「几台非常昂贵、很难使用」的机器。 直到集群普及，上千原子的第一性原理计算才成为可能。

1998 年诺贝尔化学奖得主 Walter Kohn 与 John Pople 曾有过一段著名的论断： 「在真实世界里，这最终可能意味着大多数化学实验将在芯片的硅片中进行，而不是在实验室的玻璃器皿中。 关掉本生灯吧，十年之内它就不再被需要了。」

这两位获奖者的贡献正好对应了计算化学的两条主线：

2.2 有机光电材料：OLED 与光物理

OLED（有机发光二极管）的基本结构：玻璃衬底上做出阳极，阴极在另一侧， 中间夹着发光层与导电层的有机分子或聚合物薄膜。器件的电致发光效率可以分解为：

上式中只有 $\eta_{\mathrm{out}}$ 靠光学设计，其余三项全部由分子本身的电子结构与振动耦合决定—— 这正是本课程后续要算的量。

光物理过程与 Jablonski 图

分子吸收光子后，能量通过几条相互竞争的通道耗散：

把每条箭头换成一个速率常数，光物理过程就变成可定量求解的速率方程组：

2.3 光谱、无辐射跃迁与 AIE

并苯系列的光谱

把蒽（Anthracene）→ 并四苯（Tetracene）→ 并五苯（Pentacene）放到一起算，可以看出清晰的规律：

为什么算这个：共轭长度是调控有机半导体带隙最直接的手段。算准了这一系列， 才能有信心去计算还没有合成出来的分子。

薁（Azulene）的反常发光

薁是一个著名的反例——它违背了 Kasha 规则，发射不是从最低激发态 S₁ 而是从 S₂ 出来的：

文献：2010, Yingli Niu, Qian Peng, Chunmei Deng, Xing Gao and Zhigang Shuai, J. Phys. Chem. A, 114, 7817–7831。

无辐射跃迁速率

发光效率通常不是被「发光」限制，而是被无辐射通道拖垮的：

其中包含两个因子：电子耦合项（非绝热耦合或自旋—轨道耦合） 与 Franck-Condon 加权密度（FCWD）。 要给学生的直觉是：分子越「软」、激发时结构变化越大，无辐射通道就越畅通。

聚集诱导发光（AIE）

反常现象：多数发光分子在溶液里很亮、聚集后变暗（ACQ）；但 AIE 分子正好相反—— 在溶液里几乎不发光，一聚集反而大亮。

• 溶液里不亮：分子可以自由旋转/振动，激发态能量通过分子内运动迅速无辐射耗散。
• 聚集后变亮：聚集把分子锁住，分子内转动被限制，无辐射通道被掐断。

计算设置：QM 部分用 B3LYP/6-31G*，MM 部分用 GAFF， 通过 Chemshell 3.5 把 Turbomole 与 DL_POLY 耦合起来。 文献：J. Phys. Chem. A, 2014, 118, 9094。这个例子展示了 QM/MM 混合方法的思路： 只有发光核心需要量子力学，周围的堆积环境用分子力学就够了。

重组能的内坐标分解

把总重组能按内坐标逐项分解，就能看出哪些键长/键角变化真正「吃掉」了能量。程序输出的表格大致长这样：

     n  type   |         int-def |       zd |   Reorganization energy 1 (cm-1)
                |                 |          |       diag   non-diag        sum
  ----------------------------------------------------------------------------
     1  BOND   |     C1 - C2     |   0.031  |       12.4       3.1         15.5
     2  BOND   |     C2 - C3     |   0.018  |        8.7       1.9         10.6
     3  ANGLE  |  C1 - C2 - C3   |   0.007  |        4.2       0.8          5.0
  ----------------------------------------------------------------------------
  Total (cm-1)  |                 |          |      176.3      41.2        217.5

diag 是该项自身弛豫贡献的能量，non-diag 是该项与其他坐标耦合产生的交叉项。 把贡献最大的几根键挑出来，就等于拿到了「该改哪里」的分子设计指南。

2.4 TADF 与电荷传输

热活化延迟荧光（TADF）

TADF 是第三代 OLED 发光机制的核心概念，它的逻辑只有三步：

1. 电致激发形成激子，单线态与三线态按 1 : 3 生成——按传统荧光机制，75% 的三线态被白白浪费。
2. 如果 S₁ 与 T₁ 的能隙 $\Delta E_{\mathrm{ST}}$ 足够小（< 0.2 eV）， 三线态可以反向系间窜越（RISC）回到单线态。
3. 再从 S₁ 发光，于是那 75% 的三线态也通过荧光通道发光——理论内量子效率可达 100%。

文献：J. Phys. Chem. C 2017, 121, 13448–13456。 计算的核心任务：算出 $\Delta E_{\mathrm{ST}}$ 与 RISC 速率 $k_{\mathrm{RISC}}$。

分子晶体中的电荷传输

电荷怎么走，取决于两个耦合强度的竞争：

其中 $V$ 是电子耦合（相邻分子间的电荷转移积分，$m \ne n$），$g$ 是电子—声子耦合。 理论框架是 Holstein-Peierls 哈密顿量。有机半导体室温下通常落在 Hopping 区—— 这就是为什么有机材料的迁移率远低于硅。

载流子迁移率的定义

在 Hopping 图像下，电荷传输被描述为一个扩散过程，与迁移率之间由 Einstein 关系连接： $\displaystyle \mu=\frac{eD}{k_BT}$。

随机行走模拟

真实材料里每个分子的环境都不完全一样（能量无序），于是速率也不一样。这时用随机行走模拟扩散：

1. 生成随机数 $r$（$0<r<1$）；
2. 计算各相邻跳跃的电荷转移速率 $k_{\mathrm{CT}}$；
3. 按速率归一化后抽样，决定跳向哪个格点；
4. 重复直到扩散距离超过晶格常数 2–3 个数量级；
5. 重复上千次后取平均，得到均方位移 $\langle r^2\rangle$ 与时间的线性关系。

三、计算资源与教材对应讲义 PPT 第 24–27 页

计算在校内计算集群上完成。集群提供多节点服务器，节点间由高速网络互联； 软件包括 Gaussian 16 等量子化学程序以及 Python / C 编译环境。 小体系也可以在本地 Windows 机器上用 Gaussian 16W 计算。

推荐教材

四、一维定态 Schrödinger 方程的数值解对应讲义 PPT 第 28–32 页

4.1 问题与原子单位

一维定态 Schrödinger 方程（能量本征值问题）：

为什么要数值解：解析解只对极少数特殊势能（无限深方势阱、谐振子、库仑势）成立。 一旦 $V(x)$ 的形状稍微复杂（双势阱、周期势、或来自拟合的势能曲线）， 就必须把方程离散化后用矩阵求解。

第一步：换成原子单位（atomic unit），令 $\hbar=m_e=e=1$：

4.2 边界条件

对于束缚态（bound state），波函数在无穷远处必须趋于零，等价于在有限区间两端取零：

4.3 中心差分与矩阵形式

把区间 $[0,L]$ 等分成 $N+1$ 段，步长 $h=L/(N+1)$，格点 $x_i=ih$， 记 $\psi_i=\psi(x_i)$。用中心差分近似二阶导数：

代入方程、两边同乘 $2h^2$：

写成矩阵形式：

4.4 收敛性验证

建议先做的三个例子：

验收方法：把格点数 $N$ 加倍重算，看前几个本征值是否稳定。 这一步不能省——很多「结果看起来合理」的计算，只是碰巧在未收敛的格点密度下给出的。

五、一维含时 Schrödinger 方程的数值解对应讲义 PPT 第 33–44 页

5.1 三种时间离散格式

含时 Schrödinger 方程（原子单位）：$i\dfrac{\partial\psi}{\partial t}=H\psi$。 对时间做差分有几种写法，它们的稳定性天差地别：

5.2 Cayley 形式与幺正性

把 Cayley 格式整理成传播子形式：

实际怎么算：Cayley 形式两边都含 $\psi$，不能直接显式推进。整理后得到线性方程组

右边可显式算出，左边是一个三对角方程组——于是每推进一步，就要解一次三对角系统。

5.3 追赶法解三对角方程组

无论是定态问题的矩阵对角化，还是含时问题的隐式推进，最终都要面对同一个结构：

一般的高斯消元要 $O(n^3)$，而利用三对角结构只需 $O(n)$：

1. 前向消元（追）：从第一行开始，把 $x_i$ 表示成 $x_i=\alpha_i x_{i+1}+\beta_i$，其中 $\alpha_i=\dfrac{-c_i}{a_i\alpha_{i-1}+b_i}$， $\beta_i=\dfrac{d_i-a_i\beta_{i-1}}{a_i\alpha_{i-1}+b_i}$。
2. 回代（赶）：从最后一行出发 $x_n=\beta_n$，再逐个往前代 $x_i=\alpha_i x_{i+1}+\beta_i$，直到 $x_1$。
3. 复杂度：两次扫描，总代价 $O(n)$。

这个算法又叫 Thomas 算法，是数值求解偏微分方程最常用的基础工具之一。 参考：追赶法解三对角矩阵。

5.4 初始波函数：高斯波包

含时问题必须先给一个初始态。最常用的构造是高斯波包：

六、本讲小结与练习

小结

练习

附：PPT 页码对照表

本讲为讲义体系的第 1 讲，数字页码与源课件 PPT 页码一一对应（1–44）， 没有插页，故放映文件共 44 页。 深链有两种写法：#p=N 按页码跳转（N = 1…44）， #s=M 按页序跳转（M = 1…44，两者本讲相同）。