来源论文: https://arxiv.org/abs/2606.13283v1 生成时间: Jun 13, 2026 12:53

ARCHÊ:基于自洽场无轨道密度泛函理论(OFDFT)的高效温稠密物质状态方程模拟平台

0. 执行摘要

在等离子体物理、天体物理以及惯性约束聚变(ICF)等前沿领域,准确构建宽流体范围(跨越数个数量级的温度与密度)的物质状态方程(EOS)具有至关重要的科学与工程意义。传统的 Kohn-Sham 密度泛函理论(KSDFT)由于受到波函数正交化带来的 $\mathcal{O}(N^3)$ 复杂度限制,在高温($T > 10$ eV)和大规模体系(大量原子)的模拟中面临严重的计算瓶颈——不仅是因为轨道数量随原子数增加而激增,更是由于费米-狄拉克分布导致的高能态热激发占据,使得所需计算的轨道数随温度升高而呈指数级或立方级增长。

为了攻克这一瓶颈,ARCHÊ 应运而生。作为由法国原子能和替代能源委员会(CEA)等机构开发的、用以替代传统 Fortran 版 ofmd 代码的新一代开源 C++/CUDA 分子动力学软件,ARCHÊ 基于无轨道密度泛函理论(OFDFT),直接以电子密度 $n(\vec{r})$ 作为基本变量,从而将计算复杂度成功降低至 $\mathcal{O}(N)$ 线性级别。

ARCHÊ 的核心创新和技术优势主要体现在:

  1. 独特的自洽场(SCF)求解路径:与绝大多数现存 OFDFT 代码采用共轭梯度(CG)直接极小化自由能的方法不同,ARCHÊ 独树一帜地采用自洽场(SCF)迭代配合牛顿法求解化学势,为引入先进的收敛加速算法奠定了数学基础。
  2. 双重收敛加速算法
    • 基于原子位置投影的电子密度初始化(Density Initialization):利用上一步分子动力学(MD)平衡时的自洽密度提取单中心原子剖面,并在新的原子坐标下重新组装,从而提供极度接近真实的初始密度猜测。
    • 基于微扰自由能极小化的自适应混合算法(More-style Density Mixing):将 More 等人最初用于平均原子模型的自适应混合方法推广至分子动力学,通过二阶泰勒展开极小化近似自由能,动态确定最完美的电荷密度混合比例。
  3. 非凡的计算效率与 GPU 加速:得益于 HeFFTe(百亿亿级傅里叶变换库)及 CUDA 架构的底层支持,单张 NVIDIA A100 GPU 相比于 256 核 CPU 展现出了 8 至 15 倍的加速性能。更具物理启示性的是,与 KSDFT 运行时间随温度升高而飙升相反,ARCHÊ 的运行时间随温度升高而显著减少,这使其成为研究极高温度等离子体和温稠密物质(WDM)的终极利器。

1. 核心科学问题、理论基础、技术难点与方法细节

1.1 核心科学问题:温稠密物质(WDM)与 KSDFT 的失效

温稠密物质(Warm Dense Matter, WDM)是指温度在几万度到数百万度(约 $1 \sim 100$ eV),密度接近固体密度的物质状态。它是凝聚态物理学向等离子体物理学过渡的桥梁。在 ICF 靶丸内爆、巨行星内部核幔边界等物理场景中,物质处于强耦合与半简并状态。为了模拟该体系的动力学演化并提取高精度的输运系数、不透明度及状态方程(EOS),必须采用量子力学级别的 ab initio 分子动力学。

然而,当使用传统的 KSDFT 进行模拟时,系统总自由能通过求解一电子自洽方程(Kohn-Sham 方程)获得:

$$\left( -\frac{1}{2}\nabla^2 + V_{\text{eff}}(\vec{r}) \right) \psi_i(\vec{r}) = \varepsilon_i \psi_i(\vec{r})$$

其中 $\psi_i$ 为单粒子轨道。由于必须保证这些轨道在空间上相互正交,计算复杂度包含对重叠矩阵的对角化或正交化步骤,其复杂度正比于轨道数 $N_{\text{orb}}$ 的三次方 $\mathcal{O}(N_{\text{orb}}^3)$。在有限温度下,电子遵循费米-狄拉克分布:

$$f_i = \frac{1}{e^{(\varepsilon_i - \mu)/k_B T_e} + 1}$$

随着电子温度 $T_e$ 增加,热激发使得高能态轨道被部分占据,为保证电荷中性及能量收敛,必须包含极其庞大的高能虚轨道,这直接导致传统一性原理分子动力学(FPMD)在高温下变得难以为继。

OFDFT 则通过构建仅依赖于电子密度 $n(\vec{r})$ 的非相互作用动能泛函 $F_0[n]$,彻底抛弃了波函数(轨道)的概念,其电子密度的电中性约束条件为:

$$N_e = \int_{\Omega} n(\vec{r}) d\vec{r}$$

这使得 OFDFT 的核心计算步骤转化为求解关于电子密度的欧拉方程或泊松方程,计算复杂度天然降至 $\mathcal{O}(N)$,为模拟更大规模、更高温度的 WDM 铺平了道路。

1.2 理论基础:自由能密度泛函与托马斯-费米模型

在有限温度下,根据 Mermin 定理,电子体系的平衡态由大热力学势泛函的极小值决定。体系的自由能泛函 $F[n, T]$ 写为:

$$F[n, T] = F_0[n, T] + E_H[n] + F_{\text{xc}}[n, T] + E_{\text{ext}}[n]$$

其中 $F_0[n, T]$ 为非相互作用电子气的自由能泛函,$E_H[n]$ 为电子自相互作用的哈特里(Hartree)能:

$$E_H[n] = \frac{1}{2} \int_{\Omega} \int_{\Omega} \frac{n(\vec{r})n(\vec{r}')}{|\vec{r} - \vec{r}'|} d\vec{r} d\vec{r}'$$

$F_{\text{xc}}[n, T]$ 为有限温度下的交换关联自由能。对大热力学势关于密度 $n(\vec{r})$ 求变分并施加粒子数守恒约束,得到欧拉-拉格朗日方程:

$$\frac{\delta F_0}{\delta n} + V_{\text{eff}}(\vec{r}) = \mu$$

其中有效势 $V_{\text{eff}}(\vec{r}) = V_H(\vec{r}) + V_{\text{xc}}(\vec{r}) + V_{\text{ext}}(\vec{r})$,$\mu$ 为化学势。在最经典的托马斯-费米(Thomas-Fermi, TF)近似下,局域非相互作用自由能泛函具有精确的解析形式,对应的电子密度 $n(\vec{r})$ 与有效势之间满足托马斯-费米解析关系:

$$n(\vec{r}, \mu) = \frac{(2k_B T_e)^{3/2}}{2\pi^2} I_{1/2}\left( \frac{\mu - V_{\text{eff}}(\vec{r})}{k_B T_e} \right)$$

其中 $I_{1/2}(x)$ 是 1/2 阶费米-狄拉克积分:

$$I_{j}(x) = \int_{0}^{\infty} \frac{t^j}{e^{t-x}+1} dt$$

通过对这一对偶关系进行自洽迭代,即可在不计算任何波函数的前提下确定空间电子密度分布。

1.3 技术难点一:网格计算中的库仑奇点与模守恒赝势

由于 ARCHÊ 采用三维周期性边界条件并在实空间/倒空间网格上对电子密度和电势进行离散化,原子核产生的纯库仑势 $V(r) = -Z/r$ 在核中心($r=0$)存在不可导的奇点。在网格尺寸有限的情况下,这种尖锐的奇点会导致哈特里能和电子-离子相互作用能的严重离散化误差,从而破坏自洽迭代的稳定性,导致力不守恒或非物理的能量漂移。

为了克服这一困难,ARCHÊ 引入了源自平均原子模型(Average-Atom Model)的模守恒赝势(Norm-Conserving Pseudopotential)。其核心思想是在截止半径 $r_{\text{cut}}$ 以内,用平滑的伪电荷密度 $\tilde{n}(r)$ 代替剧烈振荡的真实电荷密度 $n(r)$,而在 $r \ge r_{\text{cut}}$ 时二者严格相等。伪密度函数定义为:

$$\tilde{n}(r) = \begin{cases} \exp(a + br^2 + cr^4), & r < r_{\text{cut}} \\ n(r), & r \ge r_{\text{cut}} \end{cases}$$

参数 $a, b, c$ 通过匹配 $r_{\text{cut}}$ 处的密度值、一阶导数连续性以及截止半径内总电荷守恒来确定:

$$\tilde{n}(r_{\text{cut}}) = n(r_{\text{cut}})$$$$\left. \frac{\partial \tilde{n}}{\partial r} \right|_{r_{\text{cut}}} = \left. \frac{\partial n}{\partial r} \right|_{r_{\text{cut}}}$$$$\int_{0}^{r_{\text{cut}}} 4\pi r^2 \tilde{n}(r) dr = \int_{0}^{r_{\text{cut}}} 4\pi r^2 n(r) dr$$

在确定了平滑的伪密度 $\tilde{n}(r)$ 后,通过逆向求解欧拉方程得到平滑的单体赝势 $V_s(r)$:

$$V_s(r) = \mu - k_B T_e I_{1/2}^{-1} \left( \frac{2\pi^2 \tilde{n}(r)}{(2k_B T_e)^{3/2}} \right) - \frac{1}{r} \int_{0}^{r} 4\pi x^2 \tilde{n}(x) dx - \int_{r}^{r_{\text{cut}}} 4\pi x \tilde{n}(x) dx - V_{\text{xc}}[\tilde{n}]$$

该赝势不仅消除了核中心的数学奇点,而且通过模守恒性质完美保留了芯部电子的屏蔽效应。此外,由于赝势的使用改变了总能量的基准,代码在后处理阶段会自动施加一个基于平均原子模型计算的系统性能量修正项 $\Delta E_{\text{AA}}$(例如,在 $T=100$ eV 铝体系中,该项使总能量从赝势下的 $-2971$ eV 修正回真实库仑势下的 $-6550$ eV),以确保热力学内部能量的绝对准确。

1.4 技术难点二:自洽场(SCF)收敛机制与传统能量极小化的抉择

目前国际上几乎所有的 OFDFT 分子动力学代码(如 PROFESS, DFTpy)在求解电子密度时,均将自由能 $F[n]$ 作为优化目标,在保持电中性约束的流形上采用共轭梯度(CG)或准牛顿法(L-BFGS)直接寻找泛函的极小值点。这种直接极小化方法无需计算化学势,但在等离子体和极端高温下容易遇到 Hessian 矩阵病态、搜索步长难以自适应调整以及高频噪声振荡等收敛性缺陷。

ARCHÊ 另辟蹊径,将这一极值问题转化为求解自洽欧拉方程。其数学框架为: 对于给定的输入电势 $V_i(\vec{r})$(初始由输入密度 $n_i$ 通过泊松方程求解哈特里势得到),使用牛顿迭代法快速寻优一个全局统一的化学势 $\mu$,使得在该化学势下由托马斯-费米解析式计算得到的活性密度 $n_a(\vec{r}, \mu)$ 其积分值严格等于总电子数 $N_e$:

$$\int_{\Omega} n_a(\vec{r}, \mu) d\vec{r} - N_e = 0$$

由于 $\partial n_a / \partial \mu > 0$ 严格单调递增,牛顿法通常在 3 至 5 次迭代内即可精确收敛。随后,将生成的活性密度 $n_a$ 与输入密度 $n_i$ 按照特定物理规则混合,生成下一次迭代的输入密度,直至达到收敛标准:

$$\|n_i - n_a\|_{L^2} < \varepsilon$$

该 SCF 路径避免了高维 CG 搜索中的非物理负密度问题,且天生与各种先进的电荷混合加速算法相兼容。

1.5 方法细节:两种原创的 SCF 收敛加速算法

为了最大化 SCF 迭代的效率,ARCHÊ 首次在轨道自由分子动力学领域集成了以下两项突破性的加速技术:

1.5.1 基于原子位置投影的电子密度初始化

在分子动力学模拟中,相邻时间步长($t$ 到 $t+\Delta t$)之间的原子位移极小,这意味着电子密度分布具有极强的空间关联。若每个 MD 步长都使用均匀电荷分布(Uniform Density, $n_0 = N_e / |\Omega|$)作为 SCF 迭代的起点,将造成极大的计算资源浪费。

ARCHÊ 的做法是,将电子密度 $n(\vec{r})$ 近似看作各个单原子中心对称密度 profile 的线性叠加。在时间步 $n$ 收敛后,我们已知收敛密度 $n^n(\vec{r})$ 和对应的原子结构因子 $S_Z^n(\vec{k})$。假设单中心平均电荷密度 profile 的傅里叶变换为 $g(\vec{k})$,则满足:

$$n^n(\vec{k}) = g(\vec{k}) S_Z^n(\vec{k}) \implies g(\vec{k}) = \frac{n^n(\vec{k})}{S_Z^n(\vec{k})}$$

在下一个 MD 时间步 $n+1$,原子的位置发生移动,结构因子更新为 $S_Z^{n+1}(\vec{k})$。此时,我们直接利用保存的 $g(\vec{k})$ 组装出极具物理真实性的初始电荷密度预测值:

$$n^{n+1}_i(\vec{k}) = g(\vec{k}) S_Z^{n+1}(\vec{k}) = n^n(\vec{k}) \frac{S_Z^{n+1}(\vec{k})}{S_Z^n(\vec{k})}$$

通过这一简单而优雅的傅里叶变换投影,初始密度的平均误差大幅度降低,极大地缩短了后续 SCF 达到精度要求所需的步数。

1.5.2 基于自由能二阶微扰的自适应电荷混合(More 混合方案)

在经典的自洽迭代中,简单的线性电荷混合(Linear Mixing, $n_{\text{new}} = (1-w)n_i + w n_a$)需要人工微调一个极小的保守混合参数 $w \sim 0.05$,这不仅导致收敛极为缓慢,且在强非线性体系中极易发生密度电荷发散。而 ARCHÊ 创新性地引入了 More 混合方案,其通过分析体系总自由能泛函对密度波动的二阶变分,直接解析计算出最优的混合系数 $w$。

定义密度波动差值为 $\Delta n = n_i - n_a$。我们期望混合密度 $n_o = n_a + w \Delta n$ 能够使大热力学自由能泛函 $F[n]$ 达到极小值,即 $\frac{dF[n_o]}{dw} = 0$。将自由能泛函在活性密度 $n_a$ 处进行二阶泰勒展开:

$$F[n_o] \simeq F[n_a] + w \int_{\Omega} \frac{\delta F}{\delta n} \Delta n d\vec{r} + \frac{1}{2} w^2 \int_{\Omega} \int_{\Omega} \frac{\delta^2 F}{\delta n(\vec{r}) \delta n(\vec{r}')} \Delta n(\vec{r}) \Delta n(\vec{r}') d\vec{r} d\vec{r}'$$

利用变分极小化条件 $\frac{dF}{dw} = 0$,可以直接求得:

$$w = \frac{U}{K + U}$$

其中 $U$ 项表征了由电荷不匹配引起的一阶响应,而 $K+U$ 则为体系的二阶自洽硬度。利用拉格朗日导数关系:

$$\frac{\delta^2 F}{\delta n(\vec{r}) \delta n(\vec{r}')} = \left( \frac{\partial n_a}{\partial \mu_a} \right)^{-1} \delta(\vec{r} - \vec{r}') + \frac{1}{|\vec{r} - \vec{r}'|}$$

最终可以直接得到 $U$ 与 $K$ 的实空间积分闭合解析形式:

$$U = \int_{\Omega} \int_{\Omega} \frac{\Delta n(\vec{r}) \Delta n(\vec{r}')}{|\vec{r} - \vec{r}'|} d\vec{r} d\vec{r}'$$$$K = \int_{\Omega} \frac{[\Delta n(\vec{r})]^2}{\partial n_a(\vec{r}) / \partial \mu_a} d\vec{r}$$

其中 $\partial n_a / \partial \mu_a$ 在牛顿迭代中已作为副产物被计算。此算法巧妙地结合了静电排斥的二阶贡献(即哈特里势的贡献 $U$)与非相互作用动能的非局域贡献($K$),动态指引混合比例 $w$ 的演化。实测表明,More 自适应混合不仅彻底消除了收敛振荡,还与上述初始化方法形成强力协同。


2. 关键 Benchmark 体系、计算数据与性能表现

2.1 物理体系验证:铝(Al)的状态方程(EOS)

为了严谨评估 ARCHÊ 求解 WDM 体系状态方程的精确度,论文选择铝(Al)作为验证载体。对比的基准数据来自于目前量子化学界公认的、基于 Kohn-Sham DFT 框架且开启了 Extended DFT 模式(通过平面波对高能占据态进行近似以规避正交化)的高精度一性原理计算软件 ABINIT

模拟设置如下:

  • 网格离散:使用大小为 $256^3$ 的三维实空间网格。
  • 动力学系综:微正则系综(NVE),采用 Verlet 积分算法,时间步长满足等离子体振荡周期约束:$\Delta t < 0.1 \min(\tau_p, L/v_{\text{th}})$。
  • 相空间范围:压缩度范围跨越 $0.1\rho_0 \sim 10\rho_0$(其中 $\rho_0 = 2.7$ g/cm³),温度范围跨越 $10^4$ K $\sim 10^8$ K(相当于约 $0.86$ eV $\sim 8600$ eV)。

2.1.1 压强等温线(Pressure Isochores)对比分析

在极宽的等温线区间内,ARCHÊ 算得的压强与 ABINIT Extended 的高度一致(见图 6):

  • 在高温区($T \ge 10^5$ K,约 $10$ eV 以上),ARCHÊ 结果与 ABINIT 的相对偏差严格小于 1%,表现出令人惊叹的契合度。这是因为在极高温度下,量子壳层效应逐渐被热涨落抹平,理想费米气体的托马斯-费米近似(即一阶局域动能近似)趋于严格物理精确。
  • 在低温区($T \sim 10^4$ K),ARCHÊ 的压强数据出现轻微系统性偏低。这一偏差在物理上是完全合理的,因为在低能区,非局域动能效应(如 von Weizsäcker 梯度修正)以及强电子交换关联不可忽略,仅使用 Thomas-Fermi 泛函会在一定程度上低估壳层重组带来的排斥压。

2.1.2 内部能量等温线(Internal Energy Isotherms)对比分析

在施加了由平均原子模型得出的系统性能量漂移修正 $\Delta E_{\text{AA}}$ 之后,两款代码计算得到的总内部能量表现出高度契合(见图 7):

  • 在 $T > 10$ eV 范围内,能量曲线斜率与绝对值完全重合,误差在 1% 左右。
  • 实验证明,通过引入精细的局域密度近似(LDA)或有限温度 KSDT 交换关联泛函,这一符合度在低密度区能够得到进一步收敛。

2.2 两种收敛加速算法的协同加速表现

论文在 $T=300$ eV、$\rho = 4\rho_0$ 的超高压温条件下(包含 128 个传播原子,空间网格 $256^3$,进行 101 个 MD 时间步长模拟),定量剖析了上述两项创新算法的消减步数作用(见图 3):

加速方案配置SCF 迭代总次数相对收敛加速倍率收敛稳定性与鲁棒性
A. 朴素配置(均匀密度初始 + 固定系数线性混合 $w=0.05$)38381.00x极差,在某些 MD 步极易出现电荷震荡
B. 仅引入 More 混合(均匀密度初始 + 自由能自适应混合)13742.79x优异,无任何发散倾向
C. 仅引入原子位置投影初始化(不加 More 混合,仅普通混合)24751.55x一般,受限于小混合系数
D. 双剑合璧(More 混合 + 原子位置投影初始化)6495.91x完美,收敛轨迹极其平滑

结论:双重加速算法实现了物理上的深度协同。初始化投影提供了一个极其优秀的起点,而 More 自适应混合算法则通过宏观动力学分析快速纠正该起点的微小偏差,两者结合带来了高达 6倍 的绝对计算加速,使得自洽场 OFDFT 的效率达到了前所未有的高度。

2.3 算法复杂度与强弱可扩展性

2.3.1 实空间网格大小 $N_g$ 的可扩展性

理论上,基于 FFT 求解泊松方程和有效势的复杂度为 $\mathcal{O}(N_g \log N_g)$。但在实际测试中(见图 10,当网格从 $2^{17}$ 递增至 $2^{31}$),由于底层采用了百亿亿级高度优化的 HeFFTe 并行库,ARCHÊ 的实测总运行时间随网格增长呈现出几近完美的线性趋势。在 $2^{31}$(约 20 亿网格点)的极大规模下,代码仍能保持高效吞吐。

2.3.2 原子数 $N_i$ 的可扩展性

在固定网格大小下,改变系统中的原子数 $N_i$(由 32 个铝原子增至 500 个铝原子),总模拟时间展现出**绝对线性($\mathcal{O}(N_i)$)**特征(见图 11)。传统 KSDFT 所拥有的非线性二次/三次项在 ARCHÊ 中被完美消解,这极大有利于对宏观流体、界面扩散等需要千原子以上超大体系的直接物理建模。

2.3.3 独特的反温度运行时间趋势(Inverse Temperature Trend)

在分子动力学研究中,这是一个极其有趣的物理现象(见图 12):系统的温度越高,ARCHÊ 的总运行时间(或者说 SCF 收敛所需的平均迭代步数)反而越短!

  • 物理机制:在低温下,电子紧密凝聚在原子核周围,形成极高对比度的空间密度起伏,此时自洽场需要反复迭代才能精确刻画这种陡峭的局域边界;而随着温度急剧升高,热激发导致电子剧烈电离并趋向于弥散、平滑,空间密度分布逐渐退化为均匀电子气,这使得自洽欧拉方程的非线性特征大幅减弱,SCF 算法仅需极少次数即可轻松收敛。
  • 与之形成强烈对比的是,在 KSDFT 模拟中,温度每升高一倍,占据轨道数就随之呈暴涨之势,计算时间会呈指数式恶化。这一“反向运行趋势”确立了 ARCHÊ 在极端高温等离子体领域不可动摇的技术统治地位。

2.4 GPU 加速性能与多卡强扩展性

2.4.1 单卡 A100 与多核 CPU 强力对决

为了将硬件压榨至极限,开发团队将 ARCHÊ 最耗时的电子密度求解、Hartree 势泊松方程求解以及费米-狄拉克积分计算全部移植到了 CUDA 架构中。在单张 NVIDIA A100 GPU 相比 256 核 AMD Milan EPYC 7763 CPU 的强扩展性对比中(见图 14):

  • 在所有温度和压缩比测试中,单张 GPU 的运行速度均以压倒性优势战胜 256 核 CPU 节点。
  • 加速比(Total Time Ratio CPU/GPU)在低温区达到最大,接近 15 倍;即使在高温区(由于迭代步数变少导致非 GPU 部分占比相对增加),其加速比也稳定维持在 8 倍 以上。这意味着单张高端显卡即可提供数千核高性能计算集群的等效算力。

2.4.2 多卡通信瓶颈分析

在多卡强可扩展性测试(图 16,基于 216 个氧原子,网格大小跨越 $256^3 \sim 1024^3$)中,研究团队发现,当跨越节点进行多 GPU 强扩展时,通信延迟成为了最主要的性能瓶颈。由于 3D FFT 算法需要进行高频次、大带宽的全局全对全(All-to-All)跨卡数据交换,导致在小网格下多卡效率下降明显。因此,研究人员强烈建议:在单卡显存空间足够容纳网格阵列的前提下,优先采用单 GPU 进行 EOS 状态方程的生产计算,这能提供最佳的能效比和墙钟时间。


3. 代码实现细节与复现指南

3.1 软件架构设计与模块化哲学

ARCHÊ 放弃了前作 ofmd 繁重复杂的 Fortran 77/90 遗留代码,采用现代化 C++17 构建核心计算引擎,并遵循**自由函数(Free-function-based)**的设计范式。所有的能量项(动力学能、电学能、哈特里能、交换关联能)和相互作用力都被实现为互相独立、无状态、低耦合的模块。研究人员可以像搭积木一样,随意替换或增加新型的非局域动能泛函或交换关联泛函,而无需担心破坏整体的代码稳定性。

模块名称底层语言核心物理/计算职责
Core SCF LoopC++/CUDA负责通过牛顿法求解化学势 $\mu$,并实现 More 自适应电荷混合
HeFFTe WrapperC++统一调度底层傅里叶变换,在 CPU 架构上调用 FFTW,在 GPU 上无缝切换为 CuFFT
Force EngineCUDA高效并行的 Ewald Summation,解析计算离子-离子力和离子-电子力
XC InterfaceC++挂接开源 Libxc 库提供 0-K 关联能,并集成自研 KSDT 有限温度泛函
Front-end ToolPython快速生成初始晶格、伪电荷参数文件、动力学后处理轨迹分析

3.2 依赖环境部署

要成功编译并运行 ARCHÊ,高性能计算节点必须配备以下工具链:

  1. 编译器支持:支持 C++17 标准的 GCC (推荐 $\ge 9.3$) 或 Clang,以及 CUDA Toolkit (推荐 $\ge 11.0$)。
  2. 通信库:兼容 MPI 3.0+ 标准的实现(如 OpenMPI 或 MPICH)。
  3. 并行傅里叶变换库 HeFFTeHeFFTe (Highly Efficient FFT for Exascale) 是该代码的核心瓶颈解决器,编译时需激活 CUDA 选项以调用 cuFFT 后端。
  4. 交换关联泛函库 LibxcLibxc 用于支持零温下的交换关联能计算。

3.3 详细编译复现指南

以下是在配备 NVIDIA GPU 及 CUDA 环境的 Linux 超级计算机节点上,从零编译并运行 ARCHÊ 经典铝(Al)等离子体体系的保姆级步骤:

# 1. 克隆并进入构建目录(假设 ARCHÊ 已托管在开源平台,如 github)
git clone --recursive https://github.com/WilliamWeens/Arche.git arche_source
cd arche_source
mkdir build && cd build

# 2. 环境模块加载(以常见 HPC Lmod 系统为例)
module load gcc/10.2.0
module load openmpi/4.0.5
module load cuda/11.2

# 3. 运行 CMake 配置,激活 CUDA 选项并指定 HeFFTe 与 Libxc 路径
cmake .. \
  -DCMAKE_BUILD_TYPE=Release \
  -DARCHE_USE_CUDA=ON \
  -DARCHE_USE_MPI=ON \
  -DHeFFTe_DIR=/path/to/installed/heffte/lib/cmake/HeFFTe \
  -DLibxc_DIR=/path/to/installed/libxc/lib/cmake/Libxc

# 4. 编译并生成可执行文件
make -j16 arche_md

# 5. 运行内建的单元与集成测试以验证编译正确性
ctest --output-on-failure

输入文件(Parameters.in)典型配置示例

ARCHÊ 的分子动力学运行由一个简单直观的文本控制:

# --- 晶胞与原子几何配置 ---
BoxSizeX        10.45   # 晶胞 X 方向长度(a.u.)
BoxSizeY        10.45
BoxSizeZ        10.45
GridPointsX     256     # 实空间网格离散度
GridPointsY     256
GridPointsZ     256
TotalAtoms      128     # 体系总原子数
AtomType        Al      # 原子种类

# --- 热力学条件 ---
Temperature_eV  300.0   # 电子与离子温度 (300 eV)
ElectronDensity 4.0     # 压缩比 (4*rho_0)

# --- 分子动力学参数 ---
MD_Integrator   Verlet  # 选用经典 Verlet 积分器
TimeStep_au     0.05    # MD 步长(原子单位)
TotalSteps      1000    # 总演化步数

# --- SCF 算法开关 ---
UseDensityInit  true    # 开启原创的“原子位置投影”密度初始化!
UseMoreMixing   true    # 开启原创的“自适应二阶自由能”More 混合!
SCFTolerance    1e-6    # 自洽收敛绝对误差限

运行命令:

# 在单张 NVIDIA A100 GPU 上高效执行 128 原子的高温 WDM 动力学模拟
mpirun -np 1 ./arche_md Parameters.in

4. 关键引用文献与局限性评论

4.1 关键引用文献分析

ARCHÊ 的成功构建,深植于前人一系列经典的理论与算法突破之上:

  • Hohenberg-Kohn [50] & Mermin [51]:奠定了密度泛函理论从基态到有限温度热力学平衡态的完备物理合法性。
  • Lambert [40] (法国巴黎十一大博士论文):作为 ARCHÊ 的直接前身,提出了将模守恒赝势逆向反解并与 OFDFT 分子动力学成功耦合的完整数值方案,为克服网格奇点指明了方向。
  • More, Warren, Young & Zimmerman [41]:最初设计出这一基于二阶自由能泰勒展开的密度混合方法(QEOS 模型),旨在解决大尺度平均原子模型中的非线性不稳定性。ARCHÊ 团队极具远见地将其引入并推广到了多体三维分子动力学领域,解决了长期困扰 OFMD 收敛效率差的致命痛点。
  • Abinit Extended [13]:提供了高精度一性原理对比数据,确立了 ARCHÊ 作为快速、可靠 EOS 生成工具的学术公信力。
  • HeFFTe [42]:提供了极佳的跨节点并行 FFT 框架,使得 ARCHÊ 得以在万级处理器核上爆发强大算力。

4.2 局限性与改进方向批判

尽管 ARCHÊ 在高温和高效率方面展现出无可比拟的优势,但作为一门正处于快速迭代期的前沿代码,它在物理精度和并行架构上仍有以下局限性:

1. 动能泛函精度对低温区模拟的制约

目前,ARCHÊ 在主线验证中默认采用托马斯-费米(TF)局部非相互作用动能泛函。这种泛函完全忽略了电荷密度梯度的贡献。因此,在电子温度较低、非简并效应显著(如 $T < 10$ eV)的区域,其电子密度的壳层重构和原子间共价键结合能的描述存在先天缺陷,极易高估流体压缩性,导致 EOS 偏离实验曲线。未来的版本必须深度整合更具通用性的半局域(如带有 von Weizsäcker 修正的 GGA 动能泛函)甚至是基于核函数响应的非局域(Non-local)动能泛函。

2. 自洽场框架下引入非局域动能泛函的数学挑战

与现存 CG 直接能量极小化代码能较为方便地计算动能变分不同,在自洽场(SCF)欧拉方程框架下,若要引入极其复杂的非局域动能泛函,意味着无法再简单地通过托马斯-费米解析式和单体有效势来直接迭代求解电子密度。如何在大范围非局域算符作用下重构快速收敛的欧拉关系式,在数学和算法设计上是一个悬而未决的难题。

3. 强多卡强扩展性通信墙

正如 Benchmark 性能图 16 所示,HeFFTe 并行三维傅里叶变换的跨节点 All-to-All 组通信延迟在节点数量较多时会急剧上升。GPU 极高的本地计算速度使得这部分网络通信时间占总运行时间的比重不降反升(Amdahl 定律的又一佐证)。如何在底层算法中引入电荷密度的局域裁剪、实空间积分近似或者重叠通信与计算,以进一步降低多卡互连下的网络带宽依赖,是提升 ARCHÊ 工业界超级计算机大规模强扩展能力的重中之重。


5. 补充:无轨道 DFT(OFDFT)与 Kohn-Sham DFT(KSDFT)在温稠密物质领域的深层技术对比

为了更直观地展示为什么 OFDFT 才是处理温稠密物质(WDM)的首选,我们在此从物理机理、算法设计、并行扩展和计算复杂度等四个维度对两者进行全方位拆解对比:

特性维度Kohn-Sham DFT (KSDFT)Orbital-Free DFT (OFDFT) [ARCHÊ]物理/技术原因分析
基本计算变量单粒子波函数 $\{\psi_i(\vec{r})\}$空间电子密度 $n(\vec{r})$KSDFT 需跟踪数万个轨道的空间信息,而 OFDFT 仅需跟踪一个标量密度场。
理论计算复杂度$\mathcal{O}(N_{\text{orb}}^3)$ 三次方复杂度$\mathcal{O}(N_{\text{grid}} \log N_{\text{grid}})$ 近乎线性KSDFT 的瓶颈在于轨道的空间正交化过程;OFDFT 的核心瓶颈仅在于泊松方程的三维快速傅里叶变换。
温度可扩展性极其不友好。温度升高,高能态占据数激增,所需计算的波函数数量暴涨,计算耗时呈指数/立方恶化。极度友好。温度越高,计算速度反而越快(运行时间显著缩短)。高温下电子密度的空间分布趋于平滑、均匀,导致自洽欧拉方程的非线性变弱,SCF 迭代步数大幅度减少。
空间尺寸可扩展性只能局限于模拟数十到数百个原子的体系,计算极其昂贵。能够轻松支撑数千至数万个原子的超大物理胞模拟,线性级消耗。无轨道设计使得内存和处理器核心开销随原子数呈严格线性增长,消除二次方以上项。
动能描述精度绝对精确。通过单粒子波函数动能算符和占据数精确求和。近似描述。依赖于 Thomas-Fermi 或附加梯度修正的半局域/非局域动能泛函。失去了轨道的支持,无法直接求解动能算符,必须诉诸关于电荷密度梯度的宏观唯象泛函展开。
主要误差来源交换关联泛函(XC)的不确定性。动能泛函(KEDF)的不精确性(尤其在低 T 区)。OFDFT 的核心科学任务是在没有波函数的情况下寻找更加精确的密度动能描述。
并行化实现难点大量轨道的正交化过程需要繁重的高性能矩阵对角化库(如 ScaLAPACK, ELPA)。高维网格的域分解(Domain Decomposition)以及多卡跨节点大吞吐傅里叶变换通信。OFDFT 的并行设计极大简化了线性代数的复杂度,而完全聚焦于高效的数据搬运和一维/三维网格分解。

总结

ARCHÊ 的面世,不仅为科研界提供了一个极速、可靠、易编译的等离子体物质状态方程生成平台,更为基于自洽场(SCF)路径推进 orbital-free density functional theory 动能泛函的理论演化提供了一个坚实的现代化软件基石。通过双重自适应加速算法,ARCHÊ 成功证明了:在极端高温等离子体和温稠密物质这一传统 FPMD 的“计算禁区”中,无轨道自洽方法不仅在物理上是极其优秀的,而且在计算效率上拥有传统 Kohn-Sham 方法无法比拟的、革命性的领先优势。