均匀电子气（Homogeneous Electron Gas, HEG）的简约密度矩阵泛函理论
（Reduced Density Matrix Functional Theory, RDMFT）C++ 程序。本项目以模块化
的方式实现了若干个常见的 RDMFT 关联泛函（Hartree–Fock、Müller、Power 系列、
BBC1 …），并将得到的关联能 E_c(r_s) 与量子蒙特卡洛 (QMC) 的 PW92 参数化
结果作对比，便于后续设计新泛函并复现常见基准图。
上图为本程序在默认参数下生成的 figures/correlation_energy.png，与文献中
RDMFT 经典对比图的整体形貌一致：
| 泛函 | r_s = 1 处 E_c (Ha) |
r_s = 5 处 E_c (Ha) |
|---|---|---|
| QMC (PW92) | −0.0598 | −0.0282 |
| Power(α=0.55) | −0.0618 | −0.0346 |
| Power(α=0.58) | −0.0465 | −0.0240 |
| Müller (α=0.5) | −0.103 | −0.074 |
| HF | ~0 | ~0 |
对自旋无极化（顺磁）HEG，自然轨道为平面波，其每自旋占据数 n(k) ∈ [0, 1]
仅依赖 |k|。在原子单位下，每单位体积的能量为
T/V = (1 / (2π²)) ∫₀^∞ k⁴ n(k) dk
E_xc/V = -(1 / (2π³)) ∬₀^∞ k k' K(n(k), n(k'))
· ln|(k+k')/(k-k')| dk dk'
ρ = (1 / π²) ∫₀^∞ k² n(k) dk
其中两体核 K(n_i, n_j) 决定了具体的 RDMFT 泛函。Power 族给出
K = f(n_i) f(n_j)，f(n) = n^α：
- HF（精确交换） ⇔
α = 1 - Müller / BB ⇔
α = 1/2 - Power（Sharma 等） ⇔
α ∈ (0, 1)，常用0.55 ~ 0.58 - GU（Goedecker–Umrigar 1998） ⇔ Müller 的
α=1/2形式 外加去除单轨道自相互作用；在 HEG（平面波）极限下与 Müller 相同。
不可分（非乘积）核也已实现：
-
BBC1（Gritsenko–Pernal–Baerends 2005）：弱占据态对加负号。
-
CGA（Csányi–Goedecker–Arias，Phys. Rev. A 65, 032510 (2002)）：
K_CGA(n_i, n_j) = (1/2) * [ n_i n_j + sqrt(n_i (2 - n_i)) * sqrt(n_j (2 - n_j)) ]为 HF 交换加上空穴/统计相关贡献；因子
1/2与驱动里自洽方程的加权 一致（见CGAFunctional/solve_rdmft中的el_scale）。 -
CHF（Csányi–Arias 型 corrected-HF 核，Phys. Rev. B 61, 7348 (2000)）：
K_CHF(n_i, n_j) = n_i n_j + sqrt(n_i (1 - n_i)) * sqrt(n_j (1 - n_j))命令行中
CHF与CHF等价（历史别名），输出中的 functional 列为CHF。 -
GEO（本项目新增）：多次幂"几何均值"型核
K_GEO(n_p, n_q) = [ n_p n_q + (n_p n_q)^{1/2} + 2 (n_p n_q)^{3/4} ] / 4即 HF 核 (
alpha = 1)、Mueller 核 (alpha = 1/2) 与 Power(3/4) 核 以权重1 : 1 : 2等比例混合，并归一使K_GEO(1, 1) = 1（与 HF 的 饱和极限一致）。该核非可分，但 Euler–Lagrange 方程的左侧在n_i上 是单调递减函数（前提是U_alpha,i >= 0），因此可对每个网格点用 一维二分把n_i反解出来；外层照常用mu二分满足密度约束。 -
optGeo（optimized geometric-channel）：与 GEO 相同的三条通道
n_i n_j、(n_i n_j)^{1/2}、(n_i n_j)^{3/4}，但混合权重为w1=a^2、w2=b^2、w3=c^2，其中(a,b,c)在单位球面上（程序会对输入 做归一化，使a^2+b^2+c^2=1，从而K(1,1)=1）。CLI 使用分号分隔角度 （因为--funcs列表本身用逗号分隔），例如权重w1,w2,w3 = 0.00675,0.64213,0.35112时取OptGeo@-0.0821547206643049;0.8013311419635793;0.5925545538598386（即(-√w1,√w2,√w3)，和为 1 时已在单位球上）。三角度仅通过OptGeo@传入。 可用python3 scripts/optimize_optGeo.py在 PW92 参考下拟合(a,b,c)（脚本会输出OptGeo@...，并在build/optgeo_best/写入校验 TSV；默认先做少量粗网格 prescreen 再松收敛；--quiet减少输出，--verbose打印 C++ 标准输出；仅需 NumPy，可选 SciPy）。 -
HybOpt（图中简称 hybopt；HF 与 Power(α) 的凸组合）：两体核
K = (1−λ) n_i n_j + λ n_i^α n_j^α，即 (1−λ)·HF + λ·Power(α)（与单独选HF/Power@α时同一 JK 约定）。λ必须位于[0,1]，α必须位于(0,1)；越界或非有限参数会被拒绝，以保证 两通道 1D 二分 占据更新单调。不可因子化成单一f(n_i)f(n_j)。CLI：HybOpt@lambda;alpha（一个分号、两个浮点； bash 请加引号）。兼容旧键名OptGM@lambda;alpha（视为 HybOpt）。python3 scripts/optimize_hybopt.py（SciPy）在 PW92/QMC 参考Ec下拟合(λ, α)，输出推荐键值并写入build/hybopt_best/。 -
Beta：将
sqrt(n(1-n))空穴通道推广为可调指数K_beta(n_i, n_j) = n_i n_j + [ n_i (1 - n_i) * n_j (1 - n_j) ]^beta默认运行三个
beta值（0.45 / 0.55 / 0.65）以演示其在 HF 与 CGA 之间的插值行为：beta -> +infinity⇒ 纯 HF；beta = 1/2⇒ 与 CHF 相同的n(1-n)空穴括号；beta < 1/2⇒ 空穴贡献变强，相关性更强（Müller 极限方向）。
最大的数值难点是 ln|(k+k')/(k-k')| 在 k = k' 处的可积对数奇异性。
直接的梯形/Simpson 积分在该奇异点附近收敛极慢。本项目采用
乘积积分（product integration）：
- 在分段网格上把
k' f(n(k'))视为分段线性函数。 - 对每个区间上的
ln(k_i + k')与ln|k_i − k'|关于线性帽函数的积分用解析原函数t ln|t| − t、½ t² ln|t| − ¼ t²直接计算。
由此得到一个稠密矩阵 W[i, j]，使得
V_i = Σ_j W[i, j] k_j K(n_i, n_j),
E_xc/V = −(1 / (2π³)) Σ_i w_i k_i V_i.
对于 HF 阶跃占据数，这种积分方法在默认生产网格（801 点）下相对于解析结果
（e_x = −3 k_F / (4π)）的误差已优于 10⁻⁴。
基本算法是 Euler–Lagrange 方程
δ(T + E_xc)/δn_i = μ · δρ/δn_i ( 0 < n_i < 1 )
并辅以约束 n_i ∈ [0, 1]。
- 对 Power 族（核可分离），一步 Euler–Lagrange 解析地给出
n_i， 外层用二分法把μ调整到密度约束。 - 对一般的非可分核，使用投影梯度（projected gradient）作为兜底。
两种分支都通过混合（mixing）+ 二分 μ 的方式收敛到稳态。
按照惯例，本程序在输出文件中的 Ec_per_N 列由
E_c^RDMFT = E^RDMFT − E^HF（每电子，单位 Hartree）给出，HF 解析能量为
e_HF = (3/10) k_F² − (3/(4π)) k_F, k_F = (9π/4)^{1/3} / r_s.
QMC 参考由 PW92 拟合（include/QMC.hpp）给出。
.
├── include/
│ ├── HEG.hpp # 几何/解析公式（k_F、ρ、HF 能等）
│ ├── Grid.hpp # 一维 k 空间网格 + 积分权重
│ ├── ExchangeKernel.hpp # log-奇异核的乘积积分
│ ├── Functional.hpp # 泛函抽象基类 + HF / Müller / Power / BBC1
│ ├── Energy.hpp # T/V、E_xc/V、ρ、ε_i 的计算
│ ├── Solver.hpp # 自洽求解（含 μ 二分）
│ └── QMC.hpp # PW92 关联能参数化
├── src/main.cpp # 命令行驱动程序
├── src/MomentumDistributionGZ.cpp # 见 §5.1，Gori-Giorgi/Ziesche n(k, r_s)
├── tools/dump_gz_grid.cpp # 把 GZ n(k, r_s) 倒进 stdout 的小驱动（plot-gz 使用）
├── tools/reference/nk_GZ.f # Gori-Giorgi 发布的原始 FORTRAN 参考实现（仅存档）
├── tools/reference/README.md # 出处、SHA-256、FORTRAN ↔ C++ 例程对照表
├── tests/test_hf_exchange.cpp
├── tests/test_gz_momentum.cpp # GZ 端点 + 粒子数总和规则回归测试
├── scripts/plot_common.py # 关联能 / n(k) 图共用曲线列表与样式
├── scripts/plot_results.py # 读取 data/*.tsv 画关联能（固定几条泛函 + PW92）
├── scripts/plot_nk.py # 读取 data/nk/*.tsv 画 n(k)（同上泛函，多 r_s）
├── scripts/plot_nk_optgeo.py # 仅 optGeo：不同 r_s 的 n(k) 叠在同一张图
├── scripts/plot_nk_hybopt.py # 仅 hybopt：不同 r_s 的 n(k) 叠在同一张图
├── scripts/optimize_optGeo.py # 拟合 optGeo 三角度 vs PW92
├── scripts/optimize_hybopt.py # 拟合 HybOpt@lambda;alpha vs PW92（需 SciPy）
├── data/ # 每个泛函一个 .tsv 文件（HF.tsv / GEO.tsv / ...）
├── figures/ # 生成的对比图
├── Makefile # 简单 make 构建
└── CMakeLists.txt # CMake 构建（含单元测试）
任选其一：
# 选项 A：Makefile（最简单）
make # 编译 build/rdmft_heg 与 build/test_hf_exchange（默认 -fopenmp）
make USE_OPENMP=0 # 无 OpenMP（若环境不支持 libgomp）
make test # 运行单元测试
make run # 增量扫描：只处理 Makefile 里 FUNCS 列出的泛函（缺 TSV 才算）
make rerun # 对同一批 FUNCS 强制重算（--force）
make geo # 仅重算 GEO -> data/GEO.tsv
make optgeo # 仅重算 Makefile 中配置的 OptGeo@... -> data/*.tsv
make plot # correlation_energy.png、nk.png、nk_optgeo.png、nk_hybopt.png（需 data/ 与 make nk-data）
make clean-data # 删除所有 data/*.tsv；下一次 make run 会按 FUNCS 重新生成性能说明（相对旧实现）：预计算交换矩阵 W 的同时生成转置 W^T，使
∂E_xc/∂n_i 中按列访问 W 时内存连续；因子化泛函（HF / Müller / Power / GU）
的 V_inner 与 deps_xc 用一次矩阵–向量乘积代替 N² 次 kernel() 调用；所有
O(N²) 外层循环默认用 OpenMP 并行（make USE_OPENMP=0 可关）。GEO / optGeo /
HybOpt（hybopt）/ CGA / CHF / Beta 等非因子化核仍保持原物理公式，仅受益于连续访问与并行。
每个泛函的结果单独写到 data/<name>.tsv（例如 data/HF.tsv、
data/GEO.tsv、data/Power_0.55.tsv）。这样新增/修改一个泛函时只需
重跑该泛函（make geo 或 ./build/rdmft_heg --funcs <name> --force），
无需重复跑其他泛函的全部 r_s 扫描。scripts/plot_results.py 只叠 PW92 与
plot_common.py 中列出的几条泛函（与 plot_nk.py 一致）。make nk-data 导出
data/nk/ 后，make plot 里的 plot_nk.py 默认 --rs auto：只画目录里已有 nk 数据对应的 r_s（无长串缺文件告警）。
# 选项 B：CMake
cmake -S . -B build
cmake --build build -j
ctest --test-dir build --output-on-failure
./build/rdmft_heg --help
# 直接跑二进制时必须带 --funcs，例如：
# ./build/rdmft_heg --funcs Mueller,CGA --rs 2 --out-dir data编译需要 C++17 编译器（推荐
g++≥ 9）。scripts/plot_results.py依赖numpy和matplotlib，可通过pip install numpy matplotlib安装。
rdmft_heg [选项]
--rs <列表> 逗号分隔的 r_s 值，如 0.5,1,2,5
--funcs <列表> **必填**，逗号分隔泛函，如 HF,Mueller,Power@0.55,OptGeo@a;b;c、`HybOpt@lambda;alpha`（`OptGM@...` 为别名；含 `;` 时请整段加引号）
--N <整数> k 方向网格点数（奇数，默认 801）
--kmax <浮点> k_max 取 (factor × k_F(r_s))，默认 3
--out-dir <目录> 每个泛函一个 .tsv 文件的输出目录，默认 data
--force 覆盖已存在的 .tsv（默认跳过）
--verbose 打印自洽迭代日志
输出每个 data/<name>.tsv 的列依次为：
rs functional E_per_N Ec_per_N Ec_QMC T/N Exc/N mu rho_err converged iters
只需在 include/Functional.hpp 中继承 rdmft::Functional：
// f(n) = ln(1 + n) 仅作示意
class MyFunctional : public Functional {
public:
std::string name() const override { return "MyFunc"; }
double f (double n) const override { return std::log(1.0 + n); }
double df(double n) const override { return 1.0 / (1.0 + n); }
};如果你的核 K(n_i, n_j) 不是 f(n_i) f(n_j) 形式，再覆盖
kernel 与 kernel_grad 即可（参考 BBC1Functional 的实现）。
随后在 make_functional / main.cpp::make 注册一个标识符即可被命令行使用：
if (key == "MyFunc") return std::make_unique<MyFunctional>();新泛函会自动走通用路径：
EnergyEvaluator用kernel计算能量；Solver自动检测是否是已知的可分形式；不是则自动切到投影梯度；- QMC 与 HF 的对比无需改动。
如要使用 Power 族的解析占据更新，只要让你的类继承 PowerFunctional
（覆盖 name()），或修改 Solver::solve_rdmft 中的 dynamic_cast
分支判断。
加了新泛函之后只需要跑这一项即可，不必重算其他泛函：
# 把 MyFunc 加到 Makefile 的 FUNCS 列表（或直接传 --funcs MyFunc）
make run # 增量：只跑 data/MyFunc.tsv 还不存在的项
# 或者只重算单个：
./build/rdmft_heg --funcs MyFunc --force --out-dir data
make plot # 重画关联能与 n(k)；新泛函入图需改 scripts/plot_common.py 的 WANTED_SERIEStests/test_hf_exchange.cpp 用 401 点网格在 r_s = 2 处对比：
- 数值动能 vs
(3/10) k_F² - 数值交换能 vs
−(3/(4π)) k_F - HF 自洽能 vs 解析 HF 能
均要求相对误差 ≤ 5×10⁻³。该测试用 CTest 自动注册。
进一步的精度可以通过 --N 增加网格点数控制；驱动默认 801 用于生产数据，
单元测试仍用较粗的 401 点以加快 CI。
include/MomentumDistributionGZ.hpp / src/MomentumDistributionGZ.cpp
提供了 P. Gori-Giorgi 与 P. Ziesche 在 Phys. Rev. B 66, 235116 (2002)
（arXiv:cond-mat/0205342v3）中给出的均匀电子气参考动量分布
n(k, r_s) 的解析参数化。该分布并不来自 RDMFT 自洽，而是从 TY
有效势计算 + QMC 数据 + RPA/Wigner 极限拟合而来，可用作 RDMFT 计算的
"金标准" 对照（适用范围 r_s ≲ 12）。
-
复现论文 Fig. 6 上图（多个
r_s的n(k/k_F)）：make plot-gz # 写入 figures/nk_gz.png该目标会构建
build/dump_gz_grid（源码tools/dump_gz_grid.cpp）并调用scripts/plot_gz.py。 -
单元测试
tests/test_gz_momentum.cpp同时验证端点n(0)/n(1±)与论文 Eq. (3) 的粒子数总和∫ 3 x² n(x, r_s) dx = 1，因此构成make test/ctest的一部分。
已修复的一个论文级 bug：论文公式 Eq. (19) 的 a(r_s) 分母最后一项
印为 p_6 · r_s^6，作者本人发布的 FORTRAN 参考实现
（nk_GZ.f，Gori-Giorgi 在 csclub.uwaterloo.ca 上公开的版本；本仓库
亦在 tools/reference/nk_GZ.f 存档，并见 tools/reference/README.md
对照表）使用的是 r_s^4。先前版本的 C++ 代码忠实地拷贝了论文中的 r_s^6，导致
a(r_s=10) ≈ 0.012（应为 ≈ 0.48），论文 Eq. (3) 的归一化总和在
r_s = 5 处偏离 11%，在 r_s = 10 处偏离 35%，由此画出的 n(k, r_s)
曲线与论文 Fig. 6 在 r_s ≥ 5 时显著不符。改回 r_s^4 后归一化在整个
r_s ∈ [0.5, 12] 上偏差小于 1%，a(r_s) 也与 FORTRAN 一致。详情见
a_coeff 注释。
- J. P. Perdew & Y. Wang, Phys. Rev. B 45, 13244 (1992) — PW92 关联能。
- P. Gori-Giorgi & P. Ziesche, Phys. Rev. B 66, 235116 (2002) — 均匀电子气
动量分布
n(k, r_s)参数化（本项目include/MomentumDistributionGZ.hpp实现）；arXiv:cond-mat/0205342。 - P. Gori-Giorgi & J. P. Perdew, Phys. Rev. B 64, 155102 (2001) — 上式所
用的 on-top pair density
g(0, r_s)（GP 2001，本项目g0_on_top实现）。 - A. M. K. Müller, Phys. Lett. A 105, 446 (1984).
- S. Sharma, J. K. Dewhurst, N. N. Lathiotakis, E. K. U. Gross, Phys. Rev. B 78, 201103(R) (2008) — Power 泛函。
- O. V. Gritsenko, K. Pernal, E. J. Baerends, J. Chem. Phys. 122, 204102 (2005) — BBC1。
- M. Lubasch et al., Phys. Rev. A 88, 062512 (2013) — RDMFT 在 HEG 上的实现细节、与 QMC 的对比。