Troullier-Martins 伪化方法

本节详细介绍 Troullier-Martins (TM) 方法的数学推导和数值实现。

数学基础

伪化目标

给定全电子（AE）轨道 \(\phi_{\text{AE}}(r)\) 和本征能量 \(\varepsilon_l\)， 构造伪轨道 \(\phi_{\text{PS}}(r)\)，满足：

  1. 外区一致性 （\(r > r_c\)）：

    \[\phi_{\text{PS}}(r) = \phi_{\text{AE}}(r), \quad r > r_c\]

  2. 内区光滑性 （\(r \leq r_c\)）：

    • 无节点（避免不必要的振荡）

    • 足够光滑（高阶导数连续）

  3. 范数守恒：

    \[\int_0^{r_c} |\phi_{\text{PS}}(r)|^2 dr = \int_0^{r_c} |\phi_{\text{AE}}(r)|^2 dr\]

径向波函数形式

使用径向波函数 \(u_{nl}(r) = r R_{nl}(r)\)，在 \(r \leq r_c\) 内采用 TM 参数化：

\[u_l(r) = r^{l+1} \exp\left( \sum_{i=0}^{N} a_{2i} r^{2i} \right)\]

说明： - 因子 \(r^{l+1}\) 保证原点处行为正确（\(u_l(0) = 0\)） - 指数形式 \(\exp(p(r))\) 保证正定性 - 偶次幂 \(r^{2i}\) 保证光滑性（所有导数在原点连续）

为什么选择指数形式？

TM 方法选择 \(u(r) = r^{l+1} \exp(\sum a_{2i} r^{2i})\) 而不是简单多项式，有三个原因：

  1. 正定性：\(\exp(\cdot) > 0\)，保证内区无节点。 多项式可能有零点，导致非物理的节点。

  2. 解析性：偶次幂 \(r^{2i}\) 保证所有导数在 \(r=0\) 连续。 如果包含奇次幂，奇数阶导数在原点会有奇异性。

  3. 柔性：指数的参数空间比多项式更“平坦”，非线性方程组更容易收敛。

直观理解：伪波函数在内区的形状像一个“隆起”——在原点是零（因为 \(r^{l+1}\)）， 向外先上升后下降，最终在 \(r_c\) 处平滑接到真实波函数。

TM 方法

连续性约束

在截断半径 \(r = r_c\) 处，要求伪轨道与 AE 轨道的函数值及前 \(2N-1\) 阶导数连续：

\[u_{\text{PS}}^{(k)}(r_c) = u_{\text{AE}}^{(k)}(r_c), \quad k = 0, 1, 2, \ldots, 2N-1\]

原始 TM 方法使用 \(N=3\) （6 个约束），我们简化为 \(N=2\) （4 个约束）：

  • 函数值匹配：\(u(r_c)\)

  • 一阶导数匹配：\(u'(r_c)\)

  • 二阶导数匹配：\(u''(r_c)\)

  • 范数守恒：\(\int_0^{r_c} |u|^2 dr\)

这给出 4 个方程，求解 4 个系数 \(a_0, a_2, a_4, a_6\)。

导数计算

设 \(p(r) = \sum_{i=0}^N a_{2i} r^{2i}\)，则：

\[\begin{split}u(r) &= r^{l+1} e^{p(r)} \\ u'(r) &= \left[(l+1) r^l + r^{l+1} p'(r)\right] e^{p(r)} \\ u''(r) &= \left[l(l+1) r^{l-1} + 2(l+1) r^l p'(r) + r^{l+1} (p'^2(r) + p''(r))\right] e^{p(r)}\end{split}\]

其中：

\[\begin{split}p'(r) &= \sum_{i=1}^N 2i \, a_{2i} r^{2i-1} \\ p''(r) &= \sum_{i=1}^N 2i(2i-1) a_{2i} r^{2i-2}\end{split}\]

范数积分

内区范数：

\[Q = \int_0^{r_c} u^2(r) dr = \int_0^{r_c} r^{2(l+1)} e^{2p(r)} dr\]

这是一个超越积分，需要数值求解（使用 Simpson 法则或 Gauss 积分）。

非线性方程组

定义残差函数 \(\mathbf{F}(\mathbf{a})\) 如下（\(\mathbf{a} = [a_0, a_2, a_4, a_6]^T\)）：

\[\begin{split}F_1 &= u_{\text{PS}}(r_c) - u_{\text{AE}}(r_c) \\ F_2 &= u'_{\text{PS}}(r_c) - u'_{\text{AE}}(r_c) \\ F_3 &= u''_{\text{PS}}(r_c) - u''_{\text{AE}}(r_c) \\ F_4 &= Q_{\text{PS}} - Q_{\text{AE}}\end{split}\]

求解 \(\mathbf{F}(\mathbf{a}) = \mathbf{0}\)，使用 scipy.optimize.fsolve （混合 Powell 算法）。

数值实现

初值选择

系数初值对收敛性影响很大。推荐初值：

\[\begin{split}a_0 &\approx \ln\left(\frac{u_{\text{AE}}(r_c)}{r_c^{l+1}}\right) \\ a_{2i} &= 0, \quad i = 1, 2, 3\end{split}\]

内外区拼接

完整伪轨道：

\[\begin{split}u_{\text{PS}}(r) = \begin{cases} r^{l+1} \exp\left( \sum_{i=0}^N a_{2i} r^{2i} \right), & r \leq r_c \\\\ u_{\text{AE}}(r), & r > r_c \end{cases}\end{split}\]

数值稳定性

  1. 指数溢出保护：

    在计算 \(e^{p(r)}\) 时，如果 \(p(r)\) 过大（\(> 100\)），可能溢出。 使用对数形式处理：

    \[\ln u(r) = (l+1) \ln r + p(r)\]

  2. 范数积分精度：

    使用自适应积分（如 scipy.integrate.quad）或高阶 Simpson 法则。

  3. 求解器容差：

    fsolve 的容差设为 xtol=1e-8, maxfev=1000。

数值实现指南

伪代码框架

def tm_pseudize(r, w, u_ae, eps, l, rc, N=2):
    # 1. 找到 rc 对应的网格索引
    i_rc = find_index(r, rc)

    # 2. 计算 AE 在 rc 处的导数
    u_rc, du_rc, d2u_rc = compute_derivatives(r, u_ae, i_rc)

    # 3. 计算 AE 内区范数
    Q_ae = simpson_integral(u_ae[:i_rc]**2, r[:i_rc])

    # 4. 定义残差函数
    def residuals(a):
        u_ps, du_ps, d2u_ps = eval_tm_at_rc(rc, l, a)
        Q_ps = compute_tm_norm(r[:i_rc], l, a)
        return [
            u_ps - u_rc,
            du_ps - du_rc,
            d2u_ps - d2u_rc,
            Q_ps - Q_ae
        ]

    # 5. 求解
    a_init = [log(u_rc / rc**(l+1)), 0, 0, 0]
    a_solution = fsolve(residuals, a_init)

    # 6. 拼接
    u_ps = splice(u_ae, a_solution, l, i_rc)

    return u_ps, a_solution

精度检验

范数守恒误差

定义相对误差：

\[\delta_{\text{norm}} = \frac{|Q_{\text{PS}} - Q_{\text{AE}}|}{Q_{\text{AE}}}\]

验收标准：\(\delta_{\text{norm}} < 10^{-5}\)

连续性误差

检查在 \(r_c\) 处的相对误差：

\[\delta_k = \frac{|u_{\text{PS}}^{(k)}(r_c) - u_{\text{AE}}^{(k)}(r_c)|}{|u_{\text{AE}}^{(k)}(r_c)|}\]

验收标准：\(\delta_k < 10^{-4}\) for \(k = 0, 1, 2\)

参考文献

  • 原始论文：N. Troullier and J. L. Martins, "Efficient pseudopotentials for plane-wave calculations," Phys. Rev. B 43, 1993-2006 (1991). 方程 (14)-(18) 定义了 TM 方法的核心公式。

  • QE 实现指南：P. Giannozzi, Notes on pseudopotential generation (2019), 第 2.2-2.3 节详细说明了 TM 伪化、截断半径选择与实现要点。 https://www.quantum-espresso.org/wp-content/uploads/2022/03/pseudo-gen.pdf