tm - Troullier-Martins 伪化器
本模块实现 Troullier-Martins (TM) 方法，将全电子轨道伪化为光滑无节点的伪轨道，满足范数守恒条件。
Troullier-Martins (TM) 伪化器
实现 TM 方法生成赝轨道，满足范数守恒和连续性约束。
主要功能
tm_pseudize : 对给定 AE 轨道进行 TM 伪化 TMResult : 存储伪化结果的数据类
参考文献
Troullier & Martins, PRB 43, 1993 (1991)
算法文档
详见 docs/source/algorithm/tm_method.rst
- class atomppgen.tm.TMResult(u_ps, a_coeff, rc, eps, l, norm_error, continuity_orders, continuity_check, solver_info)[源代码]
基类：
objectTM 伪化结果
- 参数:
- u_ps
伪轨道 u(r)（径向波函数，已拼接内外区）
- Type:
np.ndarray
- a_coeff
TM 多项式系数 [a_0, a_2, a_4, ...] 内区形式：u = r^{l+1} exp(Σ a_{2i} r^{2i})
- Type:
np.ndarray
- norm_error
范数守恒相对误差 |Q_PS - Q_AE| / Q_AE
- Type:
- continuity_check
rc 处的连续性检查结果 {'u': {...}, 'du': {...}, 'd2u': {...}}
- Type:
Dict
- solver_info
求解器收敛信息 {'ier': int, 'nfev': int, 'mesg': str, 'fallback': bool}
- Type:
Dict
- atomppgen.tm.tm_pseudize(r, w, u_ae, eps, l, rc, continuity_orders=2)[源代码]
使用 Troullier-Martins 方法构造伪轨道
- 在 r ≤ rc 内区，伪轨道形式为：
u(r) = r^{l+1} · exp(Σ_{i=0}^{N} a_{2i} r^{2i})
约束条件： 1. 在 r = rc 处，u 及其导数（到 continuity_orders 阶）与 AE 轨道连续 2. 内区范数守恒：∫_0^{rc} |u|^2 dr = ∫_0^{rc} |u_ae|^2 dr
- 参数:
- 返回:
包含伪轨道、系数、范数误差等信息
- 返回类型:
- 抛出:
ValueError -- 如果 rc 不在网格范围内 如果 continuity_orders 不是 2 或 4
备注
TM 方法关键点： 1. 系数数量由连续性条件决定：
continuity_orders = 2: 需要 4 个系数 [a_0, a_2, a_4, a_6]
continuity_orders = 4: 需要 6 个系数 [a_0, ..., a_10]
数值稳定性：使用对数形式避免指数溢出
范数守恒通过非线性方程组保证（scipy.optimize.fsolve）
网格推荐：强烈建议使用 exp_transformed 或 log 网格类型。 linear 网格的导数精度不足，可能导致 TM 非线性方程组不收敛。
引用
Troullier & Martins, PRB 43, 1993 (1991), 方程 (14)-(18)
主要函数
- atomppgen.tm.tm_pseudize(r, w, u_ae, eps, l, rc, continuity_orders=2)[源代码]
使用 Troullier-Martins 方法构造伪轨道
- 在 r ≤ rc 内区，伪轨道形式为：
u(r) = r^{l+1} · exp(Σ_{i=0}^{N} a_{2i} r^{2i})
约束条件： 1. 在 r = rc 处，u 及其导数（到 continuity_orders 阶）与 AE 轨道连续 2. 内区范数守恒：∫_0^{rc} |u|^2 dr = ∫_0^{rc} |u_ae|^2 dr
- 参数:
- 返回:
包含伪轨道、系数、范数误差等信息
- 返回类型:
- 抛出:
ValueError -- 如果 rc 不在网格范围内 如果 continuity_orders 不是 2 或 4
备注
TM 方法关键点： 1. 系数数量由连续性条件决定：
continuity_orders = 2: 需要 4 个系数 [a_0, a_2, a_4, a_6]
continuity_orders = 4: 需要 6 个系数 [a_0, ..., a_10]
数值稳定性：使用对数形式避免指数溢出
范数守恒通过非线性方程组保证（scipy.optimize.fsolve）
网格推荐：强烈建议使用 exp_transformed 或 log 网格类型。 linear 网格的导数精度不足，可能导致 TM 非线性方程组不收敛。
引用
Troullier & Martins, PRB 43, 1993 (1991), 方程 (14)-(18)
数据类
- class atomppgen.tm.TMResult(u_ps, a_coeff, rc, eps, l, norm_error, continuity_orders, continuity_check, solver_info)[源代码]
TM 伪化结果
- 参数:
- u_ps
伪轨道 u(r)（径向波函数，已拼接内外区）
- Type:
np.ndarray
- a_coeff
TM 多项式系数 [a_0, a_2, a_4, ...] 内区形式：u = r^{l+1} exp(Σ a_{2i} r^{2i})
- Type:
np.ndarray
- norm_error
范数守恒相对误差 |Q_PS - Q_AE| / Q_AE
- Type:
- continuity_check
rc 处的连续性检查结果 {'u': {...}, 'du': {...}, 'd2u': {...}}
- Type:
Dict
- solver_info
求解器收敛信息 {'ier': int, 'nfev': int, 'mesg': str, 'fallback': bool}
- Type:
Dict
使用示例
基础用法
from atomppgen import solve_ae_atom, tm_pseudize
import numpy as np
# 1. 获取全电子解
ae_result = solve_ae_atom(
Z=13, # Al
spin_mode="LDA",
lmax=0,
grid_type="exp_transformed",
grid_params={"n": 800, "rmax": 100.0},
scf_params={"tol": 1e-6}
)
# 2. 提取 3s 轨道
r = ae_result.r
w = ae_result.w
u_ae_3s = ae_result.u_by_l[0][2]
eps_3s = ae_result.eps_by_l[0][2]
# 3. TM 伪化
result = tm_pseudize(
r=r,
w=w,
u_ae=u_ae_3s,
eps=eps_3s,
l=0,
rc=2.1, # Bohr
continuity_orders=2,
)
# 4. 查看结果
print(f"范数守恒误差: {result.norm_error:.2e}")
print(f"TM 系数: {result.a_coeff}")
# 5. 绘图对比
import matplotlib.pyplot as plt
plt.plot(r, u_ae_3s, label='AE')
plt.plot(r, result.u_ps, '--', label='PS')
plt.axvline(result.rc, color='k', ls=':', alpha=0.5, label=f'rc={result.rc}')
plt.legend()
plt.show()
不同连续性阶数
# 二阶连续性（4 个系数）
result_2nd = tm_pseudize(
r=r, w=w, u_ae=u_ae, eps=eps, l=0, rc=2.1,
continuity_orders=2
)
print(f"系数数量: {len(result_2nd.a_coeff)}") # 4
print(f"范数误差: {result_2nd.norm_error:.2e}")
# 四阶连续性（6 个系数，更严格）
result_4th = tm_pseudize(
r=r, w=w, u_ae=u_ae, eps=eps, l=0, rc=2.1,
continuity_orders=4
)
print(f"系数数量: {len(result_4th.a_coeff)}") # 6
print(f"范数误差: {result_4th.norm_error:.2e}")
多通道伪化
# Al 的 s, p, d 通道
ae_result = solve_ae_atom(Z=13, spin_mode="LDA", lmax=2)
rc_by_l = {0: 2.1, 1: 2.2, 2: 2.4} # 推荐截断半径
tm_results = {}
for l in range(3): # s, p, d
# 价层轨道索引（Al: 3s, 3p, 3d）
n_valence = {0: 2, 1: 1, 2: 0} # 对应 u_by_l 的索引
u_ae = ae_result.u_by_l[l][n_valence[l]]
eps = ae_result.eps_by_l[l][n_valence[l]]
result = tm_pseudize(
r=ae_result.r, w=ae_result.w,
u_ae=u_ae, eps=eps, l=l,
rc=rc_by_l[l],
continuity_orders=2
)
tm_results[l] = result
print(f"l={l}: norm_error={result.norm_error:.2e}")
连续性检查
result = tm_pseudize(r=r, w=w, u_ae=u_ae, eps=eps, l=0, rc=2.0)
# 检查各阶导数的连续性
for key in ['u', 'du', 'd2u']:
check = result.continuity_check[key]
print(f"{key}:")
print(f" PS = {check['ps']:.6e}")
print(f" AE = {check['ae']:.6e}")
print(f" 相对误差 = {check['rel_error']:.2e}")
技术细节
截断半径选择指南
推荐值（单位：Bohr）：
元素 |
通道 |
rc 推荐值 |
|---|---|---|
Al (Z=13) |
s |
2.0-2.2 |
p |
2.1-2.3 |
|
d |
2.3-2.5 |
选择原则：
最小值：rc 必须大于最外层节点半径（确保内区无节点）
最大值：rc 应小于共价半径的 1.5 倍（避免过软）
平衡点：较小的 rc 产生更硬的赝势（可转移性好），但计算成本高
连续性阶数对比
continuity_orders |
系数数量 |
约束条件 |
典型精度 |
|---|---|---|---|
2 |
4 |
u, u', u'', 范数 |
1e-6 |
4 |
6 |
u, u', ..., u'''', 范数 |
1e-7 |
说明：
**二阶**（推荐）：平衡精度与求解速度，适合大多数情况
四阶：更平滑的伪势，适合高精度计算或重元素
数值稳定性
TM 伪化内部采用以下保护措施：
指数溢出保护：在 exp(p(r)) 中，如果 p(r) > 100，自动调整初值
范数积分精度：使用 Simpson 法则，精度 O(h^4)
非线性求解器容差：xtol=1e-10，自动重试不收敛情况
如果遇到收敛问题：
尝试调整 rc（±0.1 Bohr）
检查 AE 解是否正确（converged=True）
降低 continuity_orders（从 4 降到 2）
范数守恒验收标准
教学标准：norm_error < 1e-5
生产标准：norm_error < 1e-6（推荐用于发表）
典型误差水平：
rc (Bohr) |
范数误差（Al 3s） |
|---|---|
1.8 |
~1e-7 |
2.0 |
~5e-7 |
2.2 |
~2e-6 |
2.4 |
~8e-6 |
说明：较大的 rc 会增加范数守恒误差，但通常仍在可接受范围内。
参考文献
原始论文：N. Troullier and J. L. Martins, Efficient pseudopotentials for plane-wave calculations, Phys. Rev. B 43, 1993-2006 (1991). DOI: 10.1103/PhysRevB.43.1993
QE 实现指南：P. Giannozzi, Notes on pseudopotential generation (2019), https://www.quantum-espresso.org/wp-content/uploads/2022/03/pseudo-gen.pdf