Libxc导数接口：自旋PBE、r²SCAN与差分核对

本篇适合需要调用Libxc C接口、理解交换关联势和响应核的用户。先运行完整小程序，再查看数组和数值；普通VASP使用者可以从Si与水的完整输入入手。

下载C源码、真实CSV、数值检查与三张图的绘图程序。

1 运行并找到结果#

将下载文件保存为libxc-example.zip，在允许使用的CPU环境中运行：

unzip libxc-example.zip -d libxc-response
cd libxc-response
bash run.sh

程序使用oneAPI 2024.2的icx，分别链接libraries.json列出的Libxc 6.2.2、7.0.0、7.1.2。每版输出2000行导数核对记录；整包还包括名称查询、PBE0组合及解析LDA对照。完整运行保留每项差分误差，不只输出一个“通过”。

head results/7.1.2/derivatives.csv
cat results/7.1.2/metadata.csv

文件各列依次为泛函ID、族、自旋通道数、局域点、步长编号、导数行、导数列、解析接口值、差分值和实际扰动量。row=-1表示一阶导数；其他行表示Hessian，即二阶偏导。

2 对哪一个能量求导#

Libxc的zk是每电子交换关联能epsilon。应对能量密度求导：

非自旋：F = rho × epsilon
自旋：  F = (rho_up + rho_down) × epsilon

所有局域输入使用原子单位，各导数的量纲随求导变量改变。不能把不同导数的误差直接相加来判断精度；本例对每个分量单独用1e-7 + 3e-5 × abs(解析值)作为误差比较尺度。

改变rho求偏导时，必须固定sigma、lapl、tau。本例的四组基准点虽按给定比例选择，各次有限差分只改一个变量，不沿着整组缩放关系求导。

3 自旋GGA：输入有5个独立变量#

以一个局域点为例：

double rho[2] = {0.12, 0.08};
double sigma[3] = {0.012, 0.003, 0.008};
/* sigma顺序：up-up, up-down, down-down */
/* sigma_up-down = grad(rho_up) · grad(rho_down) */

这里的交叉项就是梯度点积，不额外乘2。调用xc_gga_exc_vxc_fxc，分别给出交换或关联的每电子能、一阶势和二阶块。完整可编译代码在probe.c。

示例把上述块展开为5×5 Hessian，再逐列扰动输入、对一阶输出做中心差分。PBE交换101和关联130各用4个局域点，每点扫描5个步长。它们是局域接口演示，不是铁磁矩或材料磁响应计算。

4 非自旋meta-GGA：再加入lapl和tau#

r²SCAN交换497、关联498使用：

double rho = 0.2;
double sigma = 0.015;
double lapl = 0.0;
double tau = 0.5;
/* 输入/输出均为原子单位；tau为独立动能密度变量 */

通过xc_mgga_exc_vxc_fxc取得4个一阶分量和10个独立二阶分量，展开为4×4矩阵。求导变量顺序为rho, sigma, lapl, tau。

本例的r²SCAN不依赖lapl，因此对应一阶、二阶分量应为0。程序在调用前把所有输出初始化为0，之后逐项检查这些位置恰为0，避免把未写入的位置当作有效结果。其他依赖拉普拉斯的meta-GGA不能照搬“相关项为0”的处理结论。

本包只演示非自旋r²SCAN的一、二阶接口；自旋meta-GGA的更多数组排列需要另行处理，不能直接使用本例的4×4尺寸。

5 为什么步长不能越小越好#

python3 - <<'PY'
import json
a = json.load(open('results/VALIDATION.json'))
for version, row in a['summary'].items():
    print(version, row['step_max_normalized_error'])
PY

相对步长为1e-3, 1e-4, 1e-5, 1e-6, 1e-7，实际扰动为相对步长 × max(abs(当前分量), 0.01)。这使0值的拉普拉斯输入也有明确扰动量。

纵轴是每个分量误差除以其比较尺度后取最大值。r²SCAN在粗步长下有明显截断误差，在过细步长下又受数值舍入影响；本例1e-5附近比两端更合适。不要把整张图每个点都解释为达到同一精度。本例要求每个独立分量至少有两个步长满足比较尺度，并保存各步长的within_tolerance。

6 高阶函数存在，是否就能调用#

先查运行时标志：

const xc_func_info_type *info = xc_func_get_info(&f);
int flags = xc_func_info_get_flags(info);
if (flags & XC_FLAGS_HAVE_KXC) {
    xc_lda_kxc(&f, 1, &rho, &kxc);
}

本包直接链接的具体库，LDA交换实际给出的能力为：

lda-higher.csv未提供的位置写nan，不写成0，也不强行调用。独立xc-info程序与重新链接的小程序可能给出不同能力信息，编程时以当前加载库的标志和实际调用为准；不要将此表泛化为Libxc同版本的全部构建。

三阶结果用Slater交换的解析式独立核对：

F(rho) = -C × rho^(4/3)
C = 3/4 × (3/pi)^(1/3)
F'''(rho) = (8C/27) × rho^(-5/3)

在rho=0.01时，本次7.0.0输出KXC = 471.4597413098021，对应原子单位。20个密度点的一阶、二阶及可用三阶逐项比较；没有生成四阶数值结果。

7 与宿主计算怎样衔接#

本程序展示Libxc局域函数及导数。VASP/ABACUS/CP2K还需实现密度、轨道、精确交换、力或响应方程，不能因为库可返回二阶/三阶导数，就直接推断宿主的声子、光学或高阶响应可用。

本专题的Si与水内置/外部泛函对照来自真实VASP结果；ABACUS的泛函使用见ABACUS专题。这里新增的是可复跑的库接口示例，不新增材料自洽数据。

接口说明参考Libxc官方手册。本包完整签名与各安装目录include/xc.h对应；宿主组合及版本条件见VASP LIBXC1说明。

继续操作：共线自旋meta-GGA#

新增自旋r²SCAN九变量与完整二阶矩阵实操，提供完整C程序、真实CSV与2张图，展示上下自旋交换和零极化退化对照。原有非自旋MGGA、自旋GGA及高阶能力说明继续保留。