SymPy 氢原子波函数与能级：从 `sympy.physics.hydrogen` 到解析量子力学计算

SymPy 氢原子波函数与能级：从sympy.physics.hydrogen到解析量子力学计算

【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy

本文围绕 SymPy 的物理模块sympy.physics.hydrogen（对应文档 doc/src/modules/physics/hydrogen.rst）展开，系统讲解氢原子（以及类氢离子）径向波函数、角向球谐函数、完整波函数与非相对论/相对论能级的符号计算 API。读者将掌握R_nl、Y_lm、Z_lm、Psi_nlm、E_nl、E_nl_dirac六个函数的参数语义、调用方式、归一化验证与底层实现原理，可直接用于量子化学教学、轨道可视化数据准备与能级理论计算的符号推导。

模块概览：一个纯 Python 实现的氢原子解析工具箱

sympy.physics.hydrogen是 SymPy 物理子包中专门处理氢原子问题的模块，源码位于 sympy/physics/hydrogen.py，全部内容以Hartree 原子单位制（Hartree atomic units）表达：约化普朗克常数、电子质量、电荷量均为 1，能量单位为 Hartree，长度单位为 Bohr 半径。该模块共导出 6 个公开函数：

函数	功能	主要参数
R_nl(n, l, r, Z=1)	径向波函数 R_{nl}	主量子数 n、角量子数 l、径向坐标 r、原子序数 Z
Y_lm(l, m, phi, theta)	复球谐函数 Y_l^m	角量子数 l、磁量子数 m、方位角 phi、极角 theta
Z_lm(l, m, phi, theta)	实球谐函数 Z_l^m	同上
Psi_nlm(n, l, m, r, phi, theta, Z=1)	完整波函数 ψ_{nlm}	n、l、m、r、phi、theta、Z
E_nl(n, Z=1)	非相对论能级	n、Z
E_nl_dirac(n, l, spin_up=True, Z=1, c=...)	Dirac 相对论能级	n、l、自旋方向、Z、光速 c

在文档体系中，hydrogen.rst通过 Sphinx 的automodule指令自动抓取该模块 docstring 生成 API 文档，因此模块内的每一个 docstring 都同时充当文档正文与可执行示例（doctest），文章后续所有示例均可在 Python 会话中直接运行验证。

径向波函数R_nl：主量子数与角量子数的解析表达

径向波函数描述电子波函数随到核距离 r 的变化，是氢原子波函数的核心因子。

参数语义与取值范围

R_nl(n, l, r, Z=1)接受四个参数：

• n（主量子数）：正整数，可取 1, 2, 3, 4, ...；
• l（角动量量子数）：取值范围 0 到 n-1（对应光谱符号 s、p、d、f...）；
• r（径向坐标）：可以是符号或数值；
• Z（原子序数）：默认 1（氢），2 为氦离子 He⁺，依此类推，用于描述类氢离子。

典型调用示例

>>> from sympy.physics.hydrogen import R_nl >>> from sympy.abc import r, Z >>> R_nl(1, 0, r, Z) 2*sqrt(Z**3)*exp(-Z*r) >>> R_nl(2, 0, r, Z) sqrt(2)*(-Z*r + 2)*sqrt(Z**3)*exp(-Z*r/2)/4 >>> R_nl(2, 1, r, Z) sqrt(6)*Z*r*sqrt(Z**3)*exp(-Z*r/2)/12

对于氢原子可直接使用默认的 Z=1：

>>> R_nl(1, 0, r) 2*exp(-r) >>> R_nl(2, 0, r) sqrt(2)*(2 - r)*exp(-r/2)/4 >>> R_nl(3, 0, r) 2*sqrt(3)*(2*r**2/9 - 2*r + 3)*exp(-r/3)/27

需要说明的是，不同教材对径向函数中 r 的量纲约定可能不同（是否已除以玻尔半径），SymPy 此处约定 r 本身以 Bohr 半径为单位。

任意原子序数的类氢离子

通过修改 Z 即可推广到类氢离子。以银离子（Z=47）为例：

>>> R_nl(1, 0, r, Z=47) 94*sqrt(47)*exp(-47*r) >>> R_nl(2, 0, r, Z=47) 47*sqrt(94)*(2 - 47*r)*exp(-47*r/2)/4 >>> R_nl(3, 0, r, Z=47) 94*sqrt(141)*(4418*r**2/9 - 94*r + 3)*exp(-47*r/3)/27

归一化验证：径向积分恒为 1

径向波函数满足归一化条件 ∫₀^∞ |R_{nl}|² r² dr = 1。用 SymPy 的符号积分可以直接验证（r² 来自三维球坐标体积元的径向部分）：

>>> from sympy import integrate, oo >>> integrate(R_nl(1, 0, r)**2 * r**2, (r, 0, oo)) 1 >>> integrate(R_nl(2, 0, r)**2 * r**2, (r, 0, oo)) 1 >>> integrate(R_nl(2, 1, r)**2 * r**2, (r, 0, oo)) 1

该性质对任意原子序数均成立：

>>> integrate(R_nl(1, 0, r, Z=2)**2 * r**2, (r, 0, oo)) 1 >>> integrate(R_nl(2, 0, r, Z=3)**2 * r**2, (r, 0, oo)) 1 >>> integrate(R_nl(2, 1, r, Z=4)**2 * r**2, (r, 0, oo)) 1

源码实现要点

从 hydrogen.py 的实现可以看到 R_nl 的构造逻辑：

n, l, r, Z = map(S, [n, l, r, Z]) # 统一符号化 n_r = n - l - 1 # 径向量子数 a = 1/Z # Bohr 半径（原子单位下） r0 = 2 * r / (n * a) # 无量纲化径向变量 C = sqrt((S(2)/(n*a))**3 * factorial(n_r) / (2*n*factorial(n + l))) return C * r0**l * assoc_laguerre(n_r, 2*l + 1, r0).expand() * exp(-r0/2)

关键点有三处：

1. 关联拉盖尔多项式：R_nl 的节点结构由assoc_laguerre(n_r, 2*l + 1, r0)决定（assoc_laguerre定义于 sympy/functions/special/polynomials.py），其阶数正是径向量子数 n_r = n - l - 1，因此轨道径向节点数恰好为 n - l - 1；
2. 归一化系数：系数 C 由阶乘项精确给出，源码注释中还保留了一个等价的教科书形式C = S(2)/n**2 * sqrt(1/a**3 * factorial(n_r) / (factorial(n+l)))，两种写法数值等价；
3. 渐进指数因子：exp(-r0/2)与幂因子r0**l共同保证波函数在原点的正则行为与无穷远处的指数衰减。

测试文件 sympy/physics/tests/test_hydrogen.py 中的test_wavefunction将R_nl对 (n,l) = (1,0)、(2,0)、(2,1)、(3,0)、(3,1)、(3,2)、(4,0)...(4,3) 的解析结果与手工推导的标准表达式逐一作simplify(...) == 0断言，覆盖了低量子数的全部组合；test_norm则对 n ≤ 2 的所有 (n, l) 组合验证归一化积分恒为 1。

角向部分：复球谐函数Y_lm与实球谐函数Z_lm

角向因子决定波函数在角度方向的分布。模块提供两套约定：复球谐Y_lm与实球谐Z_lm。

Y_lm：复球谐函数

Y_lm(l, m, phi, theta)的参数中phi 为方位角（azimuthal angle），theta 为极角（polar angle），l 取 0 到 n-1，m 取 -l 到 l。函数遵循Condon-Shortley 相位约定，根据 l、m 的取值可能返回实值或复值（涉及exp(I*phi)项）。

>>> import sympy >>> from sympy.physics.hydrogen import Y_lm >>> phi, theta = sympy.symbols('phi theta', real=True) >>> Y_lm(1, 0, phi, theta) sqrt(3)*cos(theta)/(2*sqrt(pi)) >>> Y_lm(2, -1, phi, theta).simplify() sqrt(30)*exp(-I*phi)*sin(2*theta)/(8*sqrt(pi)) >>> Y_lm(1, -1, phi, theta).simplify() sqrt(6)*exp(-I*phi)*sin(theta)/(4*sqrt(pi)) >>> Y_lm(1, 1, 0, 1).evalf() -0.290723302201011

一个实用的衍生能力是求角向节点：对轨道角向概率密度求零值点，可得节点所在极角。例如对 d_z² 轨道（l=2, m=0），其角向分布为Y_lm(2, 0, phi, theta)**2，节点角度（弧度）为：

>>> dz2 = Y_lm(2, 0, phi, theta)**2 >>> [sol.evalf() for sol in sympy.solve(dz2)] [4.0969092717143, 5.32786868905508, 2.18627603546528, 0.955316618124509]

Z_lm：实球谐函数

Z_lm返回实球谐函数 Z_l^m，同样遵循 Condon-Shortley 相位约定。如果结果中虚部没有按预期消去，将 phi、theta 声明为实数符号即可（如上面示例所示）。

>>> from sympy.physics.hydrogen import Z_lm >>> Z_lm(1, 0, phi, theta) sqrt(3)*cos(theta)/(2*sqrt(pi)) >>> Z_lm(2, -1, phi, theta).simplify() -sqrt(15)*sin(phi)*sin(2*theta)/(4*sqrt(pi)) >>> Z_lm(1, -1, phi, theta).simplify() -sqrt(3)*sin(phi)*sin(theta)/(2*sqrt(pi)) >>> Z_lm(1, 1, 0, 1).evalf() -0.411144836870562

Z_lm同样支持节点求解，对 d_z² 轨道（l=2, m=0）给出的节点极角与Y_lm版本完全一致：

>>> dz2 = Z_lm(2, 0, phi, theta)**2 >>> [sol.evalf() for sol in sympy.solve(dz2)] [4.0969092717143, 5.32786868905508, 2.18627603546528, 0.955316618124509]

底层实现与约定

Y_lm与Z_lm都是对 sympy/functions/special/spherical_harmonics.py 中Ynm/Znm的封装，实现仅一行：

return Ynm(l, m, theta, phi).expand(func=True) # Y_lm return Znm(l, m, theta, phi).expand(func=True) # Z_lm

注意参数顺序的差异：Y_lm的签名为(l, m, phi, theta)，而底层Ynm的签名为(n, m, theta, phi)，封装时完成了两者互换。Ynm的定义式为 Y_n^m(θ,φ) = sqrt((2n+1)(n-m)!/(4π(n+m)!)) · exp(imφ) · P_n^m(cosθ)；Znm则按 m>0、m=0、m<0 三种情形由 Ynm 的实部/虚部线性组合构造实函数（见 spherical_harmonics.py）。Ynm还内置了对称性化简规则，例如对阶数 Y_n^{-m} 与角度翻转的简化，这保证了expand(func=True)后能得到干净的正交函数展开式。测试用例test_y_lm/test_z_lm（test_hydrogen.py）对照球谐函数表验证了 l ≤ 2 各 m 取值的显式结果。

完整波函数Psi_nlm：径向与角向的乘积

氢原子定态波函数 ψ_{nlm} 是径向波函数 R_{nl} 与球谐函数 Y_l^m 的乘积。Psi_nlm(n, l, m, r, phi, theta, Z=1)一次调用即给出完整表达式，参数语义与前述函数一致：n 为正整数，l 取值 0 到 n-1，m 取值 -l 到 l。

>>> from sympy.physics.hydrogen import Psi_nlm >>> from sympy import Symbol >>> r = Symbol("r", positive=True) >>> phi = Symbol("phi", real=True) >>> theta = Symbol("theta", real=True) >>> Z = Symbol("Z", positive=True, integer=True, nonzero=True) >>> Psi_nlm(1,0,0,r,phi,theta,Z) Z**(3/2)*exp(-Z*r)/sqrt(pi) >>> Psi_nlm(2,1,1,r,phi,theta,Z) -Z**(5/2)*r*exp(I*phi)*exp(-Z*r/2)*sin(theta)/(8*sqrt(pi))

物理合法性校验

与纯函数式的R_nl不同，Psi_nlm在计算前会对量子数做物理合法性检查（hydrogen.py）：

• n 必须是正整数，否则抛出ValueError("'n' must be positive integer")；
• 必须满足 n > l，否则抛出ValueError("'n' must be greater than 'l'")；
• 必须满足 |m| ≤ l，否则抛出ValueError("|'m'| must be less or equal 'l'")。

由于这些检查依赖is_integer/ 比较运算，只有传入具体整数或带有整数假设的符号时才会触发校验。

全空间归一化验证

ψ_{nlm} 的全空间归一化需要乘上三维球坐标的 Jacobian 因子 r² sinθ，再对 r、φ、θ 三重积分。模块 docstring 给出了完整推导式：

>>> from sympy import integrate, conjugate, pi, oo, sin >>> wf = Psi_nlm(2,1,1,r,phi,theta,Z) >>> abs_sqrd = wf*conjugate(wf) >>> jacobi = r**2*sin(theta) >>> integrate(abs_sqrd*jacobi, (r,0,oo), (phi,0,2*pi), (theta,0,pi)) 1

这一三重积分同时验证了径向归一化与球谐函数的正交归一性，是理解"波函数模方为概率密度"这一物理诠释的符号计算示例。测试test_psi_nlm（test_hydrogen.py）还对 (1,0,0)、(2,1,-1)、(3,2,1) 三个态（含 Z=2 情形）的展开式做了逐项断言。

非相对论能级E_nl

E_nl(n, Z=1)返回态 (n, l) 在 Hartree 原子单位下的能量，公式为经典的玻尔能级：

E = -Z²/(2n²)

值得注意的是，能量不依赖角量子数 l——这正是库仑势场中"偶然简并"（l 简并）的体现。

>>> from sympy.physics.hydrogen import E_nl >>> from sympy.abc import n, Z >>> E_nl(n, Z) -Z**2/(2*n**2) >>> E_nl(1) -1/2 >>> E_nl(2) -1/8 >>> E_nl(3) -1/18 >>> E_nl(3, 47) -2209/18

实现上，E_nl同样做了符号化与合法性检查：若 n 为整数且 n < 1 则抛出ValueError("'n' must be positive integer")（hydrogen.py）。测试test_hydrogen_energies（test_hydrogen.py）覆盖了 n = 1..100 及 Z = 47 的情形，并验证 n=0 时抛错。

相对论能级E_nl_dirac：Dirac 方程的精细结构

E_nl_dirac(n, l, spin_up=True, Z=1, c=...)基于 Dirac 方程计算态 (n, l, spin) 的相对论能量（Hartree 单位，不含静质量能），能自然产生自旋-轨道耦合导致的能级精细结构分裂。

参数说明

• n、l：主量子数与角量子数，约束为 l ≥ 0 且 n > l，否则抛出ValueError；
• spin_up：电子自旋向上（默认 True）或向下；
• Z：原子序数；
• c：原子单位下的光速，默认值137.035999037，取自 2010 年 arXiv 论文（arXiv:1012.3627）的 CODATA 推荐值。

精细结构示例

>>> from sympy.physics.hydrogen import E_nl_dirac >>> E_nl_dirac(1, 0) -0.500006656595360 >>> E_nl_dirac(2, 0) -0.125002080189006 >>> E_nl_dirac(2, 1) -0.125000416028342 >>> E_nl_dirac(2, 1, False) -0.125002080189006 >>> E_nl_dirac(3, 0) -0.0555562951740285 >>> E_nl_dirac(3, 1) -0.0555558020932949 >>> E_nl_dirac(3, 1, False) -0.0555562951740285 >>> E_nl_dirac(3, 2) -0.0555556377366884 >>> E_nl_dirac(3, 2, False) -0.0555558020932949

对比同 (n, l) 下的非相对论值（如 E_nl(2) = -1/8 = -0.125、E_nl(3) = -1/18 ≈ -0.0555556）可以看到，相对论修正量级约为 α² 量级，且l 越大简并分裂越小、自旋方向的影响也越小——这正是精细结构实验观测到的规律。

实现要点与限制

实现（hydrogen.py）先做合法性校验（l ≥ 0、n > l、且l = 0 时 spin_up 必须为 True，否则抛ValueError("Spin must be up for l==0.")——s 轨道无轨道角动量，不存在自旋-轨道分裂），随后构造：

skappa = -l - 1 if spin_up else -l # 有效量子数 κ beta = sqrt(skappa**2 - Z**2/c**2) return c**2/sqrt(1 + Z**2/(n + skappa + beta)**2/c**2) - c**2

测试test_hydrogen_energies_relat（test_hydrogen.py）从三个层面验证该函数：对小的 c 值（c=1, 2, 3）与手工推导的精确闭式解比对；对 c=137 与含 sqrt 的精确表达式比对；以及在真实光速下与E_nl的相对误差检查（n ≤ 4，Z=1, 2, 3 时误差阈值逐步放宽），并覆盖三个异常分支（n=0、l 为负、l=0 自旋向下）。

综合实战：从波函数到可量化数据

将上述 API 组合起来，可以完成一条完整的"量子态 → 解析表达式 → 数值数据"工作流：

1. 构造轨道波函数：用Psi_nlm(n, l, m, r, phi, theta, Z)直接得到指定量子态 (n,l,m) 的解析形式；
2. 分离径向与角向：需要单独考察径向分布时用R_nl，考察角度分布（如轨道形状、节面）时用Y_lm或实函数Z_lm；
3. 求解节点：对概率密度表达式用sympy.solve求零值点，得到径向节点位置或角向节点极角（弧度），再evalf()转为数值；
4. 能级对照：用E_nl得非相对论基准值，用E_nl_dirac考察精细结构分裂，二者差值即相对论修正；
5. 归一化核验：用integrate验证径向归一化与全空间归一化，作为教学或数值代码的解析参照。

所有这些函数都接受符号参数，因此既能对具体量子数求数值结果，也能保留 Z、r 等符号做一般性推导；模块与测试文件 test_hydrogen.py 中的断言本身即是可复用的回归用例。若需要进一步了解球谐函数的对称性与化简规则，可直接查阅其底层实现 spherical_harmonics.py 的Ynm/Znm/Ynm_c类。

小结

sympy.physics.hydrogen以极简的 API 覆盖了氢原子量子力学中最常用的解析对象：径向波函数R_nl（关联拉盖尔多项式 + 归一化系数）、角向球谐Y_lm/Z_lm（Condon-Shortley 约定）、完整波函数Psi_nlm（内置量子数合法性校验）、非相对论能级E_nl与 Dirac 相对论能级E_nl_dirac。全模块遵循 Hartree 原子单位制并支持符号参数，配合 SymPy 的integrate、simplify、solve即可完成归一化验证、节点求解与精细结构分析等典型任务，是符号计算进入量子力学教学的实用入口。

【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy