近日，来自北京大学应用物理与技术研究中心在读博士生刘千锐和陈默涵研究员在国产开源密度泛函理论软件ABACUS中实现了随机密度泛函理论（Stochastic Density Functional Theory，简称SDFT）和混合密度泛函理论（MDFT）计算方法，该方法使用了平面波基组与模守恒赝势，并且引入了布里渊区的多k点采样方法来处理对周期性体系的模拟。该方法已实现了大规模并行，并成功地应用于高效预测极端高温下的碳和硅体系的能量、受力、压强、态密度等物理性质。该工作题为“Plane-Wave-Based Stochastic-Deterministic Density Functional Theory for Extended Systems”，已被Physics Review B期刊发表（DOI：10.1103/PhysRevB.106.125132）

过去几十年中，Kohn-Sham密度泛函理论（KSDFT）由于较好的平衡了计算精度与效率，被广泛应用于化学、材料、生物等各个领域。然而传统KSDFT方法往往采用对角化体系哈密顿量的方法来得到电子能级和波函数，使得KSDFT的求解时间正比于体系电子数的3次方，即的标度。因此，使用传统的KSDFT方法模拟大尺寸的原子体系（例如数百到上千个原子）会十分困难。

过去十几年，基于有限温度的KSDFT方法也得到了大力发展，成为了研究极端高温（几个eV以上）和高压下的重要工具，例如稠密等离子体、温稠密物质的状态方程、X射线汤普森散射、光学性质和输运性质等等（有兴趣的可以查阅文中的参考文献）。然而这些计算在模拟极端高温时效率会变得十分低下。一方面，随着温度升高，更多的电子被电离，在赝势方法中需要描述更多的价电子来描述电子与离子的作用，这会极大提升需要的能量截断值；另一方面，随着温度升高，需要计算的Kohn-Sham电子波函数随体系温度增加，这使得计算时间随迅速增加。

近年来，来自以色列希伯来大学、美国加利福尼亚大学洛杉矶分校、以色列特拉维夫大学的Baer与Cytter、Neuhauser、Rabani等科学家提出了有限温度下的随机密度泛函理论（stochastic density functional theory，简称SDFT）[1,2]，该方法利用切比雪夫展开以及随机轨道求迹的方法得到基态电子密度，从而通过省去对角化矩阵得到波函数的步骤来直接获得体系电子密度，实现了线性标度的效率。此后，该方法被用于结合有限温度的DFT方法，成功模拟温稠密物质，并在极端高温情况下体现出了相较于传统有限温度KSDFT方法效率上的绝对优势。2020年，来自美国洛斯·阿拉莫斯国家实验室的White与Collins又提出了混合随机密度泛函理论（mixed stochastic-deterministic density functional theory，简称MDFT）方法[3]，该方法在低能级采用KS轨道，而高能级使用随机轨道，因此可以大大降低SDFT的随机误差。以上前人的工作为解决较大尺度和极端高温条件下传统KSDFT方法低效的问题带来了新的曙光。

传统KSDFT的流程是先初猜一个基态电子密度，使用该密度构建有效势和哈密顿算符，之后需要对角化哈密顿量，用对角化得到的本征能量与本征波函数得到新的基态电子密度，这样不断地进行所谓的电子迭代，直到电子密度收敛（见图1）。SDFT方法则是放弃了对角化的过程，取而代之的是用构建“准完备”的随机波函数通过求迹来得到基态电子密度，其余电子密度的迭代过程与传统KSDFT相同。

然而，实际程序运行时，随机波函数并不是严格完备的，这就导致其有随机的误差，这个误差会反比于随机轨道数量的平方根，因此可以通过增加随机轨道数量来减小误差。然而由于平方根的存在，这就导致当随机轨道数量较多时，再增加随机轨道的数量对减小误差的收益并不高，此时就可以采用混合了KS轨道与随机轨道的MDFT。

MDFT先解出部分KS轨道，然后将生成的随机轨道与解出的KS轨道正交，得到正交后的随机轨道，利用KS轨道与正交后的随机轨道求迹同样得到新的基态电子密度进而进行密度的迭代。事实上，当MDFT解出的KS轨道个数为零时，就是SDFT；反之，当使用的随机轨道个数为零时，就变成了传统KSDFT。

在MDFT与SDFT中，均需要计算密度算符的迹，也就是需要计算哈密顿量的费米-狄拉克函数形式作用在任意波函数上的结果:，这里就需要使用切比雪夫展开这个函数，通过计算哈密顿的幂次作用在波函数上结果完成计算，而随着温度的升高需要展开的阶数越小，计算越快。因此，MDFT很好的结合了KSDFT与SDFT的优势，与KSDFT相比，MDFT温度越高计算越快，轨道数不随温度增加而迅速增多；与SDFT相比，MDFT引入KS轨道可以迅速减小随机误差。

图1. (a)传统KSDFT对角化(DG)方法、(b)SDFT与(c)MDFT流程对比图

该工作采用平面波基组与模守恒赝势实现SDFT与MDFT的计算，并且首次对SDFT与MDFT算法引入对布里渊区的多k点的采样，这使得周期性边界条件下的计算结果更为准确。对于SDFT与MDFT的电子密度迭代、系统总能量、压强、原子受力和应力的公式推导到实现均已在国产开源密度泛函理论软件ABACUS完成，并且已在2.3版本中发布。接下来我们看看使用SDFT与MDFT计算的结果。

图2.电荷密度对比图

图2对比了温度为10 eV，无序结构的8个硅原子产生的电子密度。其中(a)是传统KSDFT使用1000条KS轨道计算的电子密度，这样求得每条KS波函数所带来的计算量是十分巨大的，而且所需的波函数个数还会随着温度上升而快速增加。（b-d）依次表示的是使用64条随机轨道的SDFT、256条随机轨道的SDFT、和256条随机轨道加32条KS轨道的MDFT这三种方法得到的电子密度，将这些电子密度和图a中得到的KSDFT电子密度做差所剩下的电子密度差值。可以看到，(b)中使用64条随机轨道计算得到的电子密度有明显的随机误差，而随着随机轨道数量的增加，误差会快速减小。最后，当采用混合KS轨道的MDFT方法时，误差就几乎不见了。

图3.平均原子自由能、压强、原子受力误差随KS轨道数量变化图

图3定量计算了温度为10 eV和密度为10.0 g/cm³硅的平均原子自由能、压强、原子受力的误差，图中横坐标表示MDFT中使用KS轨道的数量，不同的曲线使用的随机轨道数不同，红色空心圆圈使用2*2*2的布里渊区k点采样，其余均是Γ点。首先随着随机数量的增加，不同的线会向下平移，代表着误差的减小；此外，由于2*2*2的k点结果会对不同的k点的随机轨道取平均，这会进一步减小随机误差，导致算出的自由能、压强、原子受力误差也会减小；最后，在随机轨道数量保持不变时，增加KS轨道可以迅速减小随机误差，仅仅需要几十条KS轨道就能得到1000条轨道的KSDFT相近的结果了。

在SDFT与MDFT中求迹过程中，由于不同随机轨道的操作是独立的，因此可以很好地进行并行计算，图4(a)展示了对512碳原子使用16到256个CPU核时，单位时间内电子迭代的步数，可以发现对于SDFT，其几乎是一条直线，这代表着其良好的并行效率，而MDFT由于需要计算部分KS轨道，计算效率略慢于SDFT。SDFT与MDFT的效率也表征了其可以模拟高温大体系，而这对KSDFT来说是几乎不可能做到的。图4(b)对比了SDFT与传统KSDFT不同温度的计算时间，随着温度升高KSDFT用时会迅速升高，而SDFT却呈现时间减小的趋势。这是由于SDFT需要用切比雪夫展开费米-狄拉克函数，而随着温度升高，展开的阶数减小，这就导致了计算时间的下降，因此在高温计算时SDFT相比于传统KSDFT具有绝对优势。

图4. 不同核数，不同温度计算效率对比图

以上结果表明，在ABACUS中实现的SDFT与MDFT，可以得到与传统KSDFT相近的结果，并且在模拟大体系与高温体系以及多核并行上具有明显优势。

目前国产开源密度泛函理论软件ABACUS 从v2.3.0版本开始已支持SDFT与MDFT的计算。

https://github.com/deepmodeling/abacus-develop

在INPUT输入文件中，只需要将“calculation”参数设“sto-scf”或“sto-md”即可进行MDFT的电子自洽迭代（SCF）或分子动力学（MD）计算。在参数中“nbands”代表传统KS轨道数量，“nbands_sto”为随机轨道数量，“nche_sto”代表切比雪夫展开的阶数，一般模拟的温度越低需要的展开阶数越高。目前ABACUS的MDFT计算仅支持费米-狄拉克分布，即“smearing_method”要为“fd”，“smearing_sigma”为电子温度，单位为Ry。（1 Ry≈13.6 eV）其余的INPUT参数均与KSDFT的参数代表含义相同，可以参考网上文档。

（https://github.com/deepmodeling/abacus-develop/blob/develop/docs/input-main.md） 

欢迎有兴趣的老师和同学前来尝试~