ABACUS使用DFT+U计算教程 |基础版

本教程是介绍ABACUS v3.6.0发布的全新DFT+U功能模块和理论背景和实例展示。ABACUS（Atomic-orbtial Based Ab-initio Computation at UStc，中文名原子算筹）是国产开源密度泛函理论软件，以下是一些常用地址：

ABACUS 在 DeepModeling 社区中的 GitHub 仓库地址为： https://github.com/deepmodeling/abacus-develop

ABACUS 的 Gitee 镜像仓库地址为： https://gitee.com/deepmodeling/abacus-develop

ABACUS 网站访问： http://abacus.ustc.edu.cn/

文档（包括安装方法、输入输出参数介绍、功能介绍、算例介绍、开发者须知等）： https://abacus.deepmodeling.com/en/latest/。

一、DFT+U算法介绍

在计算含过渡金属元素体系的性质时，使用密度泛函理论（Density Functional Theory，简称DFT）进行计算是一种常用方法。然而，DFT在描述过渡金属元素的电子关联和局域化效应方面存在一定的局限性，这可能导致计算结果偏离实际观测值。为了解决这个问题，人们引入了DFT+U（DFT with Hubbard-U）方法，通过引入Hubbard U项来修正过渡金属元素的电子关联和局域化效应，从而提高计算的准确性。

DFT+U方法的理论背景可以概括为以下几点：

1. $\mathbf{电子关联与局域化效应}$：过渡金属元素的特点是具有不完全填充的d轨道或f轨道，这些轨道上的电子之间存在较强的库伦排斥作用，导致电子关联和局域化效应。而传统的DFT方法通常采用单电子近似，难以描述这种强关联效应。
2. $\mathbf{Hubbard}\mathbf{模型}$：Hubbard模型是描述强关联电子体系的一个简单模型，它主要包括两个项：跃迁项（描述电子在格点之间跃迁的能量）和关联项（描述同一格点上两个电子的库伦排斥能）。在DFT+U方法中，Hubbard U项被引入到Kohn-Sham方程中，以修正电子关联效应。
3. $\mathbf{Double}\mathbf{Counting}\mathbf{问题}$：由于DFT方法已经包含了一定程度的电子关联信息，因此在引入Hubbard U项时需要避免双计数问题。DFT+U方法通过减去一个与U相关的能量项来消除双计数，从而使得修正后的能量更接近真实值。
4. $\mathbf{有效}\mathbf{U}\mathbf{值的选择}$：Hubbard U的值通常需要通过实验或者从头算的方法来确定。不同的过渡金属元素和体系会有不同的U值，因此在进行DFT+U计算时，选择合适数值的U非常重要。一般来说，通过比较DFT+U计算结果与实验数据或者高级理论方法的计算结果，可以得到比较合适的U值。 总之，DFT+U方法通过引入Hubbard U项，可以较好地描述过渡金属元素体系中的电子关联和局域化效应，从而提高计算的准确性。然而，这种方法需要选择合适的U值，并且在处理非局域化问题时可能存在一定局限性。因此，根据具体体系的特点，选择合适的计算方法和参数至关重要。

ABACUS中的DFT+U功能采用了实空间局域轨道投影的方法计算DFT+U方法中的投影密度矩阵 $${\mathrm{n}}_{m{m}^{\prime }}^{I\sigma }=\sum _{R,R^{\prime }}\sum _{\mu \nu }{\rho }_{\mu \nu }^{R^{\prime }-R,\sigma }⟨{\varphi }_{\mu }^R|{\alpha }_m^I⟩⟨{\alpha }_{m^{\prime }}^I|{\varphi }_{\nu }^{R^{\prime }}⟩$$

其中投影密度算符为 $P_{m{m}^{\prime }}^I=|{\alpha }_m^I⟩⟨{\alpha }_{m^{\prime }}^I|$, 投影轨道$|\alpha ⟩$来自ABACUS使用的数值原子轨道（NAOs:Numerical-Atomic-Orbitals）$|{\varphi }_{\zeta lm}^I⟩$的局域化调制，其中 ζ 为描述相同角动量 $l$ 轨道的序号，$m$ 为磁量子数。投影轨道选取$\zeta =0,l=$orbital_corr（章节2.3） 对应的NAOs 轨道调制到截断半径onsite_radius的相似局域轨道获得。$|{\varphi }_{\mu }^R⟩$和$|{\varphi }_{\nu }^{R^{\prime }}⟩$为NAOs基组，$R、R^{\prime }$为实空间对不同周期性晶格的遍历，$\mu ,\nu$为原子轨道基组的轨道指标，$I$为目标加U的原子，$\sigma$为电子自旋指标, ${\rho }_{\mu \nu }^{R^{\prime }-R,\sigma }$为实空间密度矩阵。

能量修正采用了Dudarev 1998[1]提出的形式

$$E_U\left[\left\{n_{m{m}^{\prime }}^{I\sigma }\right\}\right]=\frac{U}{2}\sum _I\sum _{m,\sigma }\left\{n_{mm}^{I\sigma }-\sum _{m^{\prime }}{n}_{m{m}^{\prime }}^{I\sigma }{n}_{m^{\prime }m}^{I\sigma }\right\}$$

二、重点参数介绍

2.1 onsite_radius

默认值： 5.0 ； 单位：Bohr

我们的DFT+U算法具有一个独特的特性，即支持对投影轨道半径和形状的调制。这可以通过调整参数onsite_radius来实现。投影轨道是局域轨道的一种近似描述，它在计算过程中起着重要作用，尤其是在描述局域化电子的关联效应时。通过参数onsite_radius基于ABACUS的数值原子轨道$|{\varphi }_{\zeta =0}^I⟩$调整$|{\alpha }^I⟩$投影轨道的截断半径和形状，可以有效改进我们的算法的数值稳定性，并确保在一定范围内都能得到合理的结果。

这种灵活性使得我们的DFT+U算法能够更好地适应各种不同的过渡金属元素体系，从而为用户提供更准确、更可靠的计算结果。我们建议用户在进行DFT+U计算时，根据具体体系的特点和需求，适当调整onsite_radius参数，以获得最佳的计算性能。同时，我们的技术团队也会为用户提供全方位的支持，帮助用户更好地利用这一特性，优化计算结果。

这里展示的图片为不同体系的onsite_radius参数对bandgap影响的测试结果：

图1：LiCoO2和NaCoO2体系取不同的onsite_radius（单位：Bohr）计算得到的带隙

如图1所示，不同的onsite_radius均能计算出大于0的gap，但是同一个U值下对应的gap大小会有不同，具体取值需要用户结合自身场景进行测试进行选择。

2.2 dft_plus_u

该参数用于选择是否开启DFT+U功能，

• 1 ： 使用新版DFT+U功能 （建议使用）
• 2 ： 使用旧版DFT+U功能 （仅建议测试使用，不支持参数onsite_radius、mixing_dmr）

2.3 orbital_corr

设置每种元素需要+U的具体轨道。-1表示不加U，0/1/2/3分别表示在s/p/d/f轨道上加U。

示例： 2 -1 第一个元素的d轨道加U，第二个元素不加U。

2.4 hubbard_u

设置每种元素需要+U的大小。单位 eV。

示例： 5.0 0 第一个元素的U值为 5.0 eV，第二个元素的U值为 0 eV。

2.5 yukawa_potential

设置是否使用yukawa算法[2]自动计算相应的U值，如果不知道目标的U值（hubbard_u未设置）可以使用该方法。

需要说明的是，该方法仅用于对目标U值进行估计，实际使用时请优先考虑手动设置hubbard_u

默认值：false

三、辅助收敛算法

我们的DFT+U算法支持独特的自洽场（SCF）收敛性算法，可以有效优化SCF收敛过程。这种优化主要依赖于两个参数：mixing_restart和mixing_dmr。以下是对这两个参数的详细描述。

3.1 mixing_restart

默认值：0.0

这个参数用于设定一个阈值，当电荷密度误差小于该阈值时，算法会自动清空历史电荷密度信息。这有助于消除由于历史信息积累而导致的收敛性问题。设定mixing_restart的值为$\delta$, 如果当前电荷密度与上一步的电荷密度之间的差异$\Delta \rho$满足以下条件：

$$\Delta \rho <\delta ,$$

则清空历史电荷密度信息。

3.2 mixing_dmr

默认值：false

这个参数用于开启密度矩阵辅助混合功能，可以进一步优化SCF收敛性。当设置mixing_dmr=1时，算法会在每一次迭代中计算密度矩阵的变化，并利用这个信息来辅助更新电荷密度。

这两个参数需要同时使用，以实现更好的SCF收敛性。当参数mixing_restart设定了一个合适的阈值，并且mixing_dmr=1时，我们的DFT+U算法能够有效地提高SCF收敛速度和稳定性，从而为用户提供更加高效、准确的计算结果。在进行DFT+U计算时，我们建议用户根据具体体系的特点和需求，适当调整这两个参数的值，以获得最佳的计算性能。

四、实际算例：FeO算例

请连接ABACUS 3.6.0镜像以及c16_m64_cpu的机器再进行此算例的运行。

此算例我们以FeO为例，介绍如何在ABACUS中使用DFT +U的功能。

我们会一起进行+U和不+U的计算，并对其计算得到的带隙，原子磁矩和PDOS进行对比。

4.1 计算前的准备

4.1.1 下载FeO的赝势和轨道

4.1.2 安装abacus-plot

ABACUS的代码仓库中提供了一个可以进行DOS/PDOS画图的python包abacus-plot，我们要先安装此软件。
更多详细说明请参考:https://github.com/deepmodeling/abacus-develop/blob/develop/tools/plot-tools/README.md

4.1.3 建立晶体结构

首先我们进入/root目录，并创建两个路径FeO_dft和FeO_dft_u分别用于进行不打开和打开+U功能的计算

接下来，我们需要准备ABACUS计算需要的构型文件STRU。

4.1.4 准备输入参数文件INPUT

当前ABACUS v3.6.0版本的DFT +U功能只支持LCAO基组，我们需要设置basis_type为lcao。
同时我们设置out_bandgap为1用于输出计算过程中的带隙；设置out_mul用于获取原子磁矩的相关信息。

为了进行+U的计算，我们需要额外以下的设置：

• dft_plus_u 1: 是否打开+U的功能。0：关闭DFT+U功能（默认值）；1：使用新版DFT+U；2：使用旧版DFT+U
• orbital_corr 2 -1: 设置每个原子类型需要在哪个轨道上+U。-1表示不加U，1/2/3分别表示在p/d/f轨道上加U。需要注意的是，这里的顺序需要与STRU文件中ATOMIC_SPECIES中设置的原子类型一一对应
• hubbard_u 4.0 0: 设置对应原子类型添加的hubbard U的大小（单位：eV）。需要注意的是，当前的ABACUS版本（v3.6.0）中还不支持设置hubbard J值的设置，此处实际设置的是有效U值
• onsite_radius 4: 设置投影轨道半径（单位：bohr）

在打开了+U的功能时，有些情况下SCF计算会更难收敛，我们推荐额外加上以下的参数用于加速SCF的收敛：

• mixing_restart 5e-3: 设置当两步之间的密度差（drho）小于此值时，SCF使用上一步输出的电荷密度作为初猜重新进行SCF的计算。通常SCF最开始的初猜会比较差，会导致SCF在很长的步数下drho都难以下降到理想的值。此时，使用SCF过程中一个相对好的电荷密度作为初猜，重新进行SCF的计算，往往可以更容易的收敛
• mixing_dmr 1: 设置SCF过程第n步满足mixing_restart条件时，采用电荷密度的混合系数对实空间的电荷密度进行混合

4.2 不开+U的计算

4.2.1 带隙

• 带隙：在OUT.ABACUS/running_scf.log中，通过搜索关键词E_bandgap查看带隙。需要注意的是E_bandgap只有在INPUT中设置了cal_band 1时，才会输出。ABACUS会在每步SCF过程中都输出E_bandgap的值，通常我们只需要最后一个SCF收敛后的结果。

----------------------------------------------------------
    Energy           Rydberg                 eV
----------------------------------------------------------
E_KohnSham          -537.2059357810     -7309.0617260437
E_KS(sigma->0)      -537.2058532385     -7309.0606029954
E_Harris            -537.2059397812     -7309.0617804696
E_band               -39.4906550285      -537.2979261402
E_one_elec          -234.1988547208     -3186.4388892774
E_Hartree            139.2708955046      1894.8777444254
E_xc                 -64.3507532326      -875.5369145554
E_Ewald             -377.9270582472     -5141.9614205398
E_entropy(-TS)        -0.0001650850        -0.0022460965
E_descf                0.0000000000         0.0000000000
E_exx                  0.0000000000         0.0000000000
E_Fermi                1.0669422120        14.5164935197
E_bandgap              0.0008010944         0.0108994484
----------------------------------------------------------

4.2.2 磁矩

• 磁矩：在OUT.ABACUS/running_scf.log中，通过搜索关键词total magnetism查看系统的总磁矩；搜索关键词absolute magnetism查看系统的磁矩绝对值的总和

          total magnetism (Bohr mag/cell) = -4.04024e-09
       absolute magnetism (Bohr mag/cell) = 7.07319
                        nelec for spin up = 22
                      nelec for spin down = 22

• 原子磁矩：
  • 基于Mulliken计算方法：在OUT.ABACUS/mulliken.txt中，会输出每个原子上的投影电子数以及磁矩。需要注意的是，只有在INPUT文件中设置了out_mul 1才会输出mulliken.txt文件。

>>> grep "Total" OUT.ABACUS/mulliken.txt
 Total charge of spin 1:        22
 Total charge of spin 2:        22
 Total charge:  44
Total Charge on atom:  Fe               15.16
Total Magnetism on atom:  Fe               3.544
Total Charge on atom:  Fe               15.16
Total Magnetism on atom:  Fe              -3.544
Total Charge on atom:  O               6.844
Total Magnetism on atom:  O                   0
Total Charge on atom:  O               6.844
Total Magnetism on atom:  O                   0

4.2.1 投影态密度（PDOS）

完成一次SCF的计算后，再进行一次NSCF的计算，可以在OUT.ABACUS文件夹中发现PDOS文件（在INPUT中设置了out_dos 1），此文件记录了电荷密度在各个轨道上投影的结果。使用abacus仓库中的abacus-plot工具可以对PDOS的结果进行作图（关于abacus-plot的安装以及使用说明请参考：https://github.com/deepmodeling/abacus-develop/blob/develop/tools/plot-tools/README.md）。

4.3 打开+U的计算

4.3.1 查看+U投影密度矩阵以及电子占据数和原子磁矩

如果进行计算时设置了dft_plus_u 1，则在OUT.ABACUS/running_scf.log文件中会同步输出每个+U的原子上相应+U轨道的电子占据数/投影密度矩阵/原子的磁矩：

//=========================L(S)DA+U===========================//
atom_type=0  L=2  chi=0    U=4.000000eV
Local occupation matrices
atoms  0
L  2
zeta  0
eigenvalues  0
0.99354375  0.98862906  0.99340911  0.99354375  0.99340911  4.96253479
spin  0
0.99354222  0.00000583  0.00000583  0.00000000  0.00001167
0.00000583  0.99181678 -0.00159386  0.00001011  0.00159386
0.00000583 -0.00159386  0.99181678 -0.00001011  0.00159386
0.00000000  0.00001011 -0.00001011  0.99354222  0.00000000
0.00001167  0.00159386  0.00159386  0.00000000  0.99181678
eigenvalues  1
0.17433776  0.98019229  0.04445079  0.17433776  0.04445079  1.41776938
spin  1
0.17256726  0.00614858  0.00614858  0.00000000  0.01229717
0.00614858  0.35754496  0.31132367  0.01064966 -0.31132367
0.00614858  0.31132367  0.35754496 -0.01064966 -0.31132367
0.00000000  0.01064966 -0.01064966  0.17256726 -0.00000000
0.01229717 -0.31132367 -0.31132367 -0.00000000  0.35754496
atomic mag: 0 3.54476541

上述的输出，表明:

• 第0号原子（Fe）的spin up的5个d轨道的电子占据数分别是：0.99354375 0.98862906 0.99340911 0.99354375 0.99340911；
• 第0号原子（Fe）的spin down的5个d轨道的电子占据数是：0.17433776 0.98019229 0.04445079 0.17433776 0.04445079。
• 最后的atomic mag输出了此原子上的总磁矩为3.54476541 ${\mu }_B$。

4.4 结果分析

对FeO（反铁磁）体系，我们分别进行关闭DFT+U和打开DFT+U的的计算。在DFT+U的计算中，Fe原子的d轨道上设置了4.0 eV的U值。相比普通的PBE计算，打开DFT+U功能，计算的FeO带隙从0.01 eV增大至2.1 eV。并且从PDOS的图上，我们可以很明显的看到：不开+U时，投影在Fe的d轨道上的态密度在费米能附近是连续的；打开DFT+U功能后，投影在Fe的d轨道上的态密度可以在费米能附近形成一个较大的带隙。

五、参考文献

1. Dudarev S L, Botton G A, Savrasov S Y, et al. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+ U study[J]. Physical Review B, 1998, 57(3): 1505.
2. Yue Chao Wang and Hong Jiang. Local screened Coulomb correction approach to strongly correlated d - electron systems. J. Chem. Phys., 2019, 150(15):16–20.