Libxc自旋meta-GGA：九变量数组与r²SCAN实例
完整C程序演示自旋r²SCAN的一阶势、9×9二阶矩阵、自旋互换和零极化对照。
阅读本篇目录
上一份导数接口教程使用自旋PBE和非自旋r²SCAN。本篇继续展示共线自旋r²SCAN：输入从4个独立变量增加为9个，不能再按非自旋数组长度调用。只需要通过VASP使用r²SCAN的用户，无需编写这些代码，可先看Libxc的VASP实操。
1 下载、执行、查看#
将资料包保存为libxc-spin.zip，在允许使用的CPU环境中执行：
unzip libxc-spin.zip -d libxc-spin
cd libxc-spin
bash run.sh
入口加载oneAPI 2024.2，用icx分别链接libraries.json里的Libxc 6.2.2、7.0.0、7.1.2既有私有库。程序串行运行，不需要GPU，不修改宿主库。需在新目录运行，已有results时退出；绘图默认使用集群ASE Python中的Matplotlib，也可通过PLOTPY指定自己的绘图Python。
head results/7.1.2/derivatives.csv
cat results/7.1.2/spin-controls.csv
ls figures
每版得到3600行导数数据、8行自旋对照。reference保存实际结果；results/VALIDATION.json列出17,088项检查和每个步长的误差。这里计算的是局域泛函函数值，不是某个材料的自洽能量或磁化率。
2 九个变量如何排列#
本例使用以下顺序，所有量采用原子单位：
| 索引 | 变量 | 基准点示例 |
|---|---|---|
| 0、1 | rho_up, rho_down |
0.12、0.08 |
| 2、3、4 | sigma_uu, sigma_ud, sigma_dd |
0.012、0.003、0.008 |
| 5、6 | lapl_up, lapl_down |
0、0 |
| 7、8 | tau_up, tau_down |
0.3、0.2 |
sigma_uu和sigma_dd分别为两个自旋密度梯度的平方；sigma_ud为两个梯度的点积，不能额外乘2。tau为相应自旋通道的动能密度，不能用一个长度为1的数组代替。
double rho[2] = {0.12, 0.08};
double sigma[3] = {0.012, 0.003, 0.008};
double lapl[2] = {0.0, 0.0};
double tau[2] = {0.3, 0.2};
xc_func_type f;
xc_func_init(&f, 497, XC_POLARIZED);
/* 完整输出数组和调用见spin_mgga.c */
497是r²SCAN交换，498是r²SCAN关联，程序分别计算，不把单个分量称作完整XC。四个局域点使用缩放因子0.5、0.8、1、1.5构造；每次求偏导只扰动一个变量，其他变量固定。
3 二阶输出为什么有这么多块#
xc_mgga_exc_vxc_fxc分别返回9个一阶分量和下列二阶块：
| 二阶块 | 输出长度 |
|---|---|
| rho–rho | 3 |
| rho–sigma | 6 |
| rho–lapl、rho–tau | 各4 |
| sigma–sigma | 6 |
| sigma–lapl、sigma–tau | 各6 |
| lapl–lapl、tau–tau | 各3 |
| lapl–tau | 4 |
合计45个独立分量，展开后为9×9对称矩阵。同组变量的对称块按上三角逐行展开，混合块按行优先展开。例如rho–sigma的6个数，先给出rho_up对应三个sigma的导数，再给出rho_down对应的三个。
完整程序的rect和sym只负责按照接口排列展开，随后用一阶输出的独立有限差分核对9×9全部位置，不以展开后的对称性本身证明正确。
r²SCAN不使用拉普拉斯，因此本例索引5、6相关的一阶、二阶位置都应为0。调用前输出数组全部清零，调用后对这些位置再做严格零检查。其他依赖拉普拉斯的meta-GGA仍需保留对应导数。
4 数值差分核对什么#
一阶导数应来自能量密度：
F = (rho_up + rho_down) × epsilon_xc
库返回的exc是每电子能量，不是这个能量密度。二阶则由一阶势对各变量做中心差分得到。相对步长从1e-3扫描到1e-7；实际扰动为相对步长 × max(abs(变量值), 0.01)。
每个独立分量至少有两个步长满足1e-7 + 3e-5 × abs(解析值)的比较尺度。粗步长的截断误差、很小步长的舍入误差都会增大；图中保留所有步长的结果，不表示每一个点都同样精确。
python3 - <<'PY'
import csv
p = 'results/7.1.2/derivatives.csv'
for row in csv.DictReader(open(p)):
if row['id'] == '497' and row['point'] == '2':
if row['step'] == '2' and row['row'] == '0' and row['col'] == '1':
print(row)
PY
这条命令提取交换分量基准点上rho_up–rho_down的交叉二阶导数，以及相对步长1e-5的差分值。
5 两个自旋对照，帮助排除数组混用#
第一组把上下自旋完全互换：密度、拉普拉斯和动能密度对应交换，sigma_uu与sigma_dd交换，sigma_ud保留。没有外部自旋方向偏好的局域泛函，其每电子能量应不变。本例在三个版本的全部8个交换/关联局域点逐项比较，绝对差不超过1e-12 Ha。
第二组将同一个非自旋输入拆成零极化自旋输入：
rho_up = rho_down = rho / 2
sigma_uu = sigma_ud = sigma_dd = sigma / 4
lapl_up = lapl_down = lapl / 2
tau_up = tau_down = tau / 2
三个sigma项都是四分之一，因为梯度各为总密度梯度的一半。程序分别调用非自旋和自旋接口，本例能量差均在1e-12 Ha以内；真实结果在spin-controls.csv，不是预先填入的参考常数。
6 怎样读二阶结果#
图中只画rho–rho块，三条曲线具有相同量纲Ha·bohr³，可直接比较。横轴表示所选局域输入点的总密度；sigma和tau也随各基准点变化，因此曲线不是“固定其他局域量、连续对rho求导”的路径。
rho–sigma、sigma–tau等块的量纲不同，不与rho–rho混在同一个物理纵轴上。矩阵分量也不直接等于材料的磁化率、声子或介电常数；宿主软件还需完成相应响应方程。
本篇覆盖共线自旋r²SCAN的局域一阶/二阶接口，不延伸为非共线、自旋流或任意泛函的专有参数。高阶导数的运行时查询及宿主边界见Libxc导数接口，名称与编号见版本和组合泛函。接口定义参考Libxc官方手册及各安装前缀的include/xc.h。