qsolver：一个用于精确求解束缚态和共振态的耦合径向薛定谔方程的程序

《Computer Physics Communications》：qsolver: a program for the accurate solution of coupled radial Schr?dinger equations for bound and resonant states

　　**摘要**
Fortran程序`qsolver`旨在求解具有束缚能级和共振（准束缚）能级的耦合径向薛定谔方程组。为了寻找共振态的能量和宽度，该程序利用了外部复数缩放方法。势能矩阵$V(r)$和可选的径向耦合矩阵$B(r)$由用户提供，它们可以是径向坐标$r$的实函数或复函数。对于共振态计算，`qsolver`还可以自动将实值的乘法函数$V_{mn}(r)$和$B_{mn}(r)$外推到复平面。一个精确的解析变换$\xi=\phi(\tau)$将修改后的耦合通道（CC）方程从半无限区间$r \in [a, \infty)$或无限区间$\tau \in (-\infty, +\infty)$映射到有限域$x \in [a, b]$上，从而可以使用具有指数收敛率的多项式配置方法对其进行数值求解。该程序可以通过检查光谱基中的本征函数来识别虚假（“幽灵”）共振态，并能够计算实数和复数本征值关于哈密顿参数（包括约化质量$\mu$）的偏导数。`qsolver`基于现代Fortran 2003接口构建，其中大多数参数都是可选的，因此使用起来非常直接，也可以轻松集成到用户的代码中。数值测试表明，该程序能够以高达13-14位有效数字的精度在双精度（real*8）运算中找到共振态的能量和宽度。

**程序概述**
**程序名称：** `qsolver`
**CPC库链接到程序文件：** https://doi.org/10.17632/x4dgx6k5g8.1
**开发者仓库链接：** https://github.com/vlavlam/qsolver
**许可规定：** MIT
**编程语言：** Fortran 2003
**外部例程/库：** BLAS, LAPACK, ARPACK

**问题性质：**
精确求解耦合径向薛定谔方程可能会遇到困难，特别是在处理接近解离和/或电离极限的密集能级时。当需要确定共振（准束缚）态的能量和宽度时，这些挑战变得更加严重。

**解决方法：**
采用了复数缩放技术，可以将束缚态和共振态统一作为具有零边界条件的复数缩放方程的本征态来处理。然后通过积分变量的解析变换，将耦合通道（CC）方程精确映射到有限域上，以便应用有效的数值方法。通过对变换后的问题进行多项式配置，该方法显式地处理了外部复数缩放引入的导数不连续性，从而保持了光谱收敛率。最终得到的实数或复数矩阵本征值问题使用标准线性代数库LAPACK和ARPACK来解决。

**其他说明（包括限制和不寻常的特性）：**
该方法的一个限制是它依赖于平滑的乘法径向函数；对于不够平滑的势能和耦合函数，会失去指数收敛率。此外，目前密集的哈密顿矩阵的显式构造和存储限制了可以处理的通道数量。该程序的一个独特之处在于它能够仅依赖原始的实值势能矩阵，通过复数缩放技术同时计算共振能量和宽度。矩阵会自动外推到复平面；不需要用户提供额外的解析延展。为了抑制通常出现在真实共振附近的虚假根，程序内部有一个过滤器来丢弃非物理本征值。程序还提供了计算出的能量和宽度关于任何哈密顿参数的一阶解析偏导数，这有利于基于梯度的优化或敏感性分析。`qsolver`作为一个单一的、自包含的Fortran 2003模块提供，具有一个子程序公共接口，可以通过直接包含源文件而不是作为单独构建的库来集成到用户的项目中。

**引言**
单通道和耦合径向薛定谔方程在原子和分子物理学的多个领域中起着关键作用，从对孤立状态的传统绝热处理到对非绝热耦合电子态流的全面微扰分析[1]，包括共振电子-原子散射和Rydberg态[2]、Feshbach共振和超冷分子形成[3]等专门领域。

**特别挑战：**
在解决逆光谱问题时会出现特定挑战，该问题需要从可用的实验数据中直接恢复势能曲线和非绝热耦合矩阵元素。这种方法需要重复求解大量的振动-转动态的径向薛定谔方程。当考虑密切相邻的电子态之间的非绝热相互作用时，计算分子的共振能量和预解离率尤其具有挑战性。

**已经开发了许多用于解决非耦合（单通道）径向薛定谔方程的强大程序**。例如，Bob LeRoy开发的广泛使用的`LEVEL`程序[4]可以找到单个径向方程的束缚能级和共振能级，并评估障碍势的共振能级宽度。然而，这些宽度是在半经典近似下获得的，这限制了它们的精度。其他程序如`Duo`[5]可以计算双原子分子的振动-电子能级，`BOUND`[6]专用于相互作用原子和分子的束缚态，以及`LILIX`包[7]。所有这些程序都仅涉及束缚态，尽管它们也可以用于稳定方法[8]框架内的共振态迭代搜索，但该方法只能提供共振能量和宽度的近似估计。共振参数也可以通过散射计算提取，例如使用`MOLSCAT`程序[9]，但只能通过计算截面的能量依赖性间接获得。`SpectrumSDT`程序[10]处理三原子系统的振动-转动计算中出现的大型耦合径向方程组，并使用复杂的吸收势计算共振能量和宽度。在外场中的分子碰撞计算中也会遇到耦合径向方程；对于这类问题，已经提出了一种对数导数方法来定位场可调的Feshbach共振[11]。

**工作介绍**
本文介绍了Fortran 2003程序`qsolver`，用于求解单通道和耦合通道径向薛定谔方程。它在单一框架内处理束缚态和共振（准束缚）态，并直接且高精度地计算共振能量和宽度。这使得`qsolver`非常适合上述具有挑战性的问题，特别是非绝热耦合电子态的逆光谱分析。

**部分代码示例**
**问题性质：**
耦合径向薛定谔方程组表示如下：
$$
-\frac{\hbar^2}{2m}\psi_{\parallel r}^2\mathbf{r}\cdot\mathbf{r} + \frac{\hbar^2}{2m}\mathbf{V}(\mathbf{r})\cdot\mathbf{r} + \frac{\hbar^2}{2m}[\mathbf{V}(\mathbf{r})\cdot\mathbf{r} - \lambda\mathbf{r}\cdot\mathbf{r}(\mathbf{V}^\dagger(\mathbf{r})]\mathbf{r}\cdot\mathbf{V}(\mathbf{r}) + \mathbf{E}_{\text{potential}}(\mathbf{r}) = \mathbf{E}_{\text{orbit}}(\mathbf{r})
$$
其中$\mu$是粒子的约化质量，$\psi_{\parallel r}(r)$是M分量通道向量，其分量为$\psi_m(r)$，$B(r)$和$V(r)$分别是M×M矩阵，其元素分别为$B_{mn}(r)$和$V_{mn}(r)$。矩阵$G(r)$是对角的，其元素为$G_m(r)$；其存在使方程(1)具有广义本征值问题的形式。如果未指定$G(r)$，则方程(1)保持原始形式。

**复数缩放方法概述**
一种众所周知的解决这个问题的技术是复数缩放方法[15]，它允许$r$取复数值，从而将系统(1)扩展到复平面。通过沿着复数轮廓延伸到无穷大可以避免实数轴，在实数轴上波函数倾向于发散。轮廓通过新的实数坐标$y$定义如下：
$$
\begin{cases}
\psi_{\parallel r}(y) = \psi_{\parallel r}(r) & \text{if } y \leq r_c \\
\psi_{\parallel r}(y) = \psi_{\parallel r}(r_c + y e^{i\theta} & \text{if } y > r_c
\end{cases}
$$
其中$rc$是轮廓离开实轴的任意选定点，$\theta$是进入复平面的旋转角度。

**解析映射到有限域**
如果方程(1)定义在具有无限端点的区间上，`qsolver`可以将其映射到有限域，从而便于使用标准数值方法。这是通过对变量$y$应用坐标映射实现的。在没有复数缩放的情况下，$y = r$。在外部复数缩放下，$y$参数化复轮廓$\xi = \phi(\tau)$。替换表示为：
$$
\begin{align}
\xi &= \psi_{\parallel r}(\xi) \\
\phi_{\tau}(x) &= \sqrt{\frac{\xi'}{\xi'} \cdot \psi_{\parallel r}(\xi)}
\end{align}
$$
其中$\xi' = e^{i\theta}\psi_{\parallel r}(\xi)$，$m = 1, \cdots, M$。

**通过配置方法进行数值求解**
对于修改后的薛定谔方程(14)的数值求解，程序使用光谱配置方法[20]，这是一类基于插值的算法，对于具有解析解的问题，随着网格的细化通常会呈指数收敛。

**共振态的虚假本征值过滤**
一旦通过光谱配置对微分方程进行离散化，得到的矩阵问题可能会包含非物理（虚假）本征值，这些本征值不对应于原始问题的真实解。对于不涉及共振态的问题，虚假本征值很容易与物理本征值区分开来，因为真正的束缚态本征值位于势能的解离限制以下，而非物理本征值位于该阈值以上。

**自动外推到复平面**
为了实施复数缩放方法，必须计算耦合通道系统(1)的系数函数在$r$超出$rc$时的复数值。这些是势能矩阵$V(r)$、非绝热耦合矩阵$B(r)$和对角矩阵$G(r)$的元素，以及在请求能量导数时的它们相对于模型参数的偏导数。程序支持两种模式。第一种模式接受用户提供的函数。

**计算本征值关于哈密顿参数的导数**
如果耦合通道系统(1)依赖于参数$p_k$，可能需要计算计算出的能量关于这些参数的导数。例如，在直接势能拟合方法中，这些导数对于执行非线性优化程序是必不可少的。`qsolver`提供了计算这些导数的功能。当应用复数缩放时，程序使用了Hellmann–Feynman定理的非厄米类似物[30]。

**基本用法**
与Fortran 77相比，现代Fortran大大简化了`qsolver`子程序的使用。它的大多数参数都是可选的，允许用户仅指定与当前任务特别相关的参数，而不必受不必要的细节的困扰。例如，可以使用以下调用计算量子谐振子的前10个本征值：
```fortran
qcvt calculated_eigenvalues(got_energy, NH, ns)
```
在这个例子中，空间边界没有明确给出，这意味着$\tau_{\min} = -\infty$和$\tau_{\max} = \infty$。

**数值示例**
程序`qtest.f90`包含30多个测试示例，探讨了`qsolver`实现的不同方面。这些示例涵盖了从简单的、可解析解决的问题（如盒子中的粒子、谐振子和库仑势）到更复杂的任务（包括涉及搜索共振能级的多通道情况）。这些示例应作为`qsolver`的实用指南。`qtest`的典型输出包含在程序分发中。

**关于ECS半径的实用选择**
上述数值示例为正确选择ECS半径值$rc$提供了一个有用的实用指南。轮廓应优先在共振波函数至少大致达到其外向波渐近形式后才进入复平面。如果$rc$放置得太靠近主要相互作用区域，波函数的非渐近部分的解析延展可能在复平面分支上产生较大振幅，从而恶化计算性能。

**计算成本和当前限制**
设$M$为通道数，$N$为每个通道的配置点数。离散化哈密顿矩阵的阶数为$\mathcal{O}(NH^2)$。当前的`qsolver`实现显式地构建和存储这个密集矩阵，因此其主要存储需求为$\mathcal{O}(NH^2H)$。假设每个双复数元素占16字节，单个$NH \times NH$矩阵大约占用3.6 GB内存；实数矩阵所需的内存量减半。

**结论性评论**
Fortran 2003程序`qsolver`可以计算给定一般形式的单通道和耦合径向薛定谔方程的束缚态和共振态。它已经在从孤立障碍势到文献中的耦合通道（CC）模型的各种基准测试中得到了验证，在这些测试中展示了指数收敛性和高达13-14位有效数字的绝对精度（使用real*8浮点运算）。

**VRediT作者贡献声明**
Vladimir Meshkov：撰写——原始草稿、可视化、验证、软件、方法论、调查、形式分析、概念化。Asen Pashov：撰写——审阅与编辑、验证。安德烈·斯托利亚罗夫：撰写、审稿与编辑、验证、监督、资金筹措、概念构思。
利益冲突声明：
作者声明不存在可能影响本文所述工作的已知财务利益冲突或个人关系。
弗拉基米尔·梅什科夫 | 阿森·帕绍夫 | 安德烈·斯托利亚罗夫
莫斯科国立罗蒙诺索夫大学化学系，列宁斯基耶戈雷街1/3号，莫斯科119991，俄罗斯