DFT+U Calculation

  • 为什么要进行DFT+U校正

    • DFT方法中包含了电子自相关作用项（这部分没有物理意义）造成对电子相关作用的高估
    • 特别是L(S)DA经常不能描述具有局域(强相关)d和f电子的系统(这主要表现为不现实的单电子能量)。在某些情况下，这可以通过以类似Hartree-Fock(屏蔽)的方式引入强原子内相互作用来补救。这种方法通常被称为L(S)DA+U方法。设置LDAU = .TRUE。
    • DFT（LDA和GGA）对于一般体系的计算结果是令人满意的，尤其是能带结构的计算，这些一般体系主要是金属体系或者是只包含前三周期元素的体系。但是，对于包含d电子或者f电子的体系，特别是过渡金属氧化物或者氮化物，DFT直接计算的结果往往是错误的，所以在金属/绝缘体的判定上常常出错。LDA和GGA往往会低估一些绝缘体或者半导体的带隙，甚至最高占据轨道（VBM或者HOMO）在Fermi面之上，变成金属。对于包含d或者f电子的体系，VBM或者HOMO往往是来自这些金属原子的d电子或者f电子，而DFT无法直接处理d轨道或者f轨道的强关联相互作用，目前广泛采用LDA+U的方法来处理d电子或者f电子的这种强关联相互作用。
    • 何时用何时不用
      • 对金属单质不用
      • 对含d电子的金属氧化物必须要用
      • 注意在不同的DFT+U条件下计算出的能量，没有可比性！
    • LDA＋U核心思想是：首先将研究体系的轨道分隔成两个子体系（subsystem），其中一部分是一般的DFT算法（如LSDA，GGA）等可以比 较准确描述的体系，另外是定域在原子周围的轨道如d或者f轨道，这些轨道在标准的DFT计算下不能获得正确的能量与占据数之间的关系（如DFT总是认为分 数占据是能量最小的，而不是整数占据）；对于d或者f轨道，能带模型采用Hubbard模型，而其他轨道仍然是按照Kohn－Sham方程求解；d以及f 轨道电子之间的关联能采用一个和轨道占据以及自旋相关的有效U表示；整体计算的时候需要将原来DFT计算过程中已经包含的部分关联能扣除，这部分一般叫 Double Counting part，并且用一个新的U来表示，最终的结果是在DFT计算的基础上新增加一个和d或者f轨道直接相关的分裂势的微扰项，这部分能量可以采用一般微扰理论计算

  • 有关参数：

    DFT+U Calculation
LDAU   = .TRUE.        (Activate DFT+U)
LDATYPE=  2            (Dudarev; only U-J matters)
LDAUL  =  2 -1         (Orbitals for each species)
LDAUU  =  2  0         (U for each species)
LDAUJ  =  0  0         (J for each species)
LMAXMIX=  4            (Mixing cut-off; 4-d, 6-f)

  • 说明：

    • LDAU flag：判断是否进行DFT+U计算

    • LDATYPE flag：进行+U修正的类型

      • 1：The rotationally invariant LSDA+U introduced by Liechtenstein et al.[3] $$ E{HF}=\frac{1}{2} \sum{\gamma}(U{\gamma_1 \gamma_3 \gamma_2 \gamma_4}-U{\gamma1 \gamma_3 \gamma_4 \gamma_2})\hat{n}{\gamma1 \gamma_2}\hat{n}{\gamma3 \gamma_4} $$ 其中，$\hat{n}{\gamma1 \gamma_2}=\lang \Phi^{s_2} |m_2\rang \lang m_1 | \Phi{s_1}\rang$ (PAW占据数)

        $U_{\gamma_1 \gamma_3 \gamma_2 \gamma_4}$为电子相关作用项

        最终总能形式 $$ E{total}(n,\hat{n})=E{DFT}(n)+E{HF}(\hat{n})-E{dc}(\hat{n}) $$ 其中，$E{dc}(\hat{n})=\frac{U}{2}\hat{n}{tot}(\hat{n}{tot}-1)-\frac{J}{2}\sum{\sigma}\hat{n}{tot}^{\sigma}(\hat{n}{tot}^{\sigma}-1)$

        U和J由参数LDAUU与LDAUJ指定

      • 2：The simplified (rotationally invariant) approach to the LSDA+U, introduced by Dudarev et al.[4]

        最终能量形式为： $$ E{LSDA+U}=E{LSDA}+\frac{(U-J)}{2} \sum{\sigma}[(\sum{m1}n{m1,m_2}^\sigma)-(\sum{m1,m_2}\hat{n}{m1,m_2}^\sigma \hat{n}{m_2,m_1}^\sigma)] $$ U和J由参数LDAUU与LDAUJ指定(注意这里只有U-J的值是有意义 的)

      • 4：和1一样不过是LDA+U而非LSDA+U

        最终能量形式为： $$ E{dc}(\hat{n})=\frac{U}{2}\hat{n}{tot}(\hat{n}{tot}-1)-\frac{J}{2}\sum{\sigma}\hat{n}{tot}^{\sigma}(\hat{n}{tot}^{\sigma}-1) $$

    • LDAUL flag:

      • 其值为对应轨道的主量子数 （d: 2, f:3）
      • 取-1时代表对此轨道不加U修正
      • 格式为LDAUL=num1 num2 ...(依次对应POSCAR中的原子顺序)

    • LDAUU以及LDAUJ flag：

      • U和J的值取得见仁见智，一般搜索文献取前人的结果即可

      • 这里列表如下，以方便后续查询（更新中）：

        计算中的DFT+U设置

    • LDAUPRINT = 0 | 1 | 2控制LDA+U计算的输出

      • LDAUPRINT=0: silent.
      • LDAUPRINT=1: Write occupancy matrix to the OUTCAR file.
      • LDAUPRINT=2: same as LDAUPRINT=1, plus potential matrix dumped to stdout

    • LMAXMIX = [integer](default:2) mixing cutoff

      若量子数大于此则不用混合，一般无需增加（设置）此值，除非：

      • L(S)DA+U calculations require in many cases an increase of LMAXMIX to 4 for d-electrons (or 6 for f-elements) in order to obtain fast convergence to the groundstate
      • The CHGCAR file will contain the one-center PAW occupancy matrices up to LMAXMIX. When the CHGCAR file is read and kept fixed in the course of the calculations (ICHARG=11), the results will be necessarily not identical to a self-consistent run. The deviations will be large for L(S)DA+U calculations. For the calculation of band structures within the L(S)DA+U approach it is strictly required to increase LMAXMIX to 4 for d-elements and to 6 for f-elements. （计算能带结构或读入CHGCAR并且不改变时）

  • PWmat存在可计算U值的方法：

    U值计算目前仍然是一个研究的难点和热点，最近PWmat实现了一种用线性响应理论计算U值的方法，下面我们将对这种方法进行详细介绍：

    首先U值表示定域轨道电子与电子之间的相互作用，其操作算符如下：

    img

    其中，ak, ai+表示湮灭和产生算符，i,j,k,l是原子轨道基组

img

其中V仅仅表示数字，需要注意的是，上式中索引指数i,j,k,l和变量r,r`，需要考虑存在自旋指数σ。

进一步，我们假设img是相同角动量L内的2L+1原子轨道，对于球谐函数YL,i（降低自旋指数），然后考虑对称性，V的形式可以简化为：

img

其中ah是Klein-Gordon系数，形式如下：

img

只有偶次的Fh不为零，例如F0，F2，F4和F6（f轨道）。人们习惯上假设J=（F2+F4）/14=0（d轨道），定义U=F0。

现在，如果我们插入单个Slater行列式波函数来计算Hint的期望值，并减去重复项，然后整理之后得到：

img

I遍历不同的原子，σ是自旋指数，img是一个矩阵，形式如下：

img

其中，img是Kohn-Sham波函数，img是它的占据数。

到目前为止，我们应该如何计算U值呢？我们可以通过下式计算：

img

但是这不包括周围环境对该原子的屏蔽左右，因此计算出来的U值偏大。另外人们也可以通过经验方法估算U值，得到任何我们想要的值。这里我们介绍一种用Koopman`s条件来得到U值【3】。作为原子轨道占据数的函数，在电子数N到N+1的变化时，总能应该是直线。所以我们可以调整U的值，直到最后的总能变成直线。

img

实际上，很难occupy所有img。事实证明，进行线性响应计算更容易。另外，决定哪个轨道占据也很困难，因此在线性响应计算中，我们将从所有轨道上获取电荷：

img

我们具有自旋σ的总和，所有最后，我们只有一个U值，而不是两个U，每个自旋一个。我们将通过线性响应计算来确定响应电荷来自哪个轨道（或者多个轨道）。

至此，我们得到一个新的Hamiltonian（线性能量项下的线性响应）：

img

其中，img仅仅是原子I上的一个数值。

在计算中，可以设置U=0，（从LDA开始），即：

img

对原子I设置一个最终的img，其他原子的img设置为零，然后进行SCF计算，现在就有E(img)，移除最后一项，还可以得到ELDA(img)和ELDA(nI)，最后得到U的形式为：

img

该公式可以理解为，当我们添加额外LDA+U项EU时，该项的二阶导数将抵消LDA项的二阶导数，从而使得最终结果成为nI的线性函数（Koopmans`s 条件），即

img

最终U的形式为：

img

事实证明，由于上述过程中nI的变化，总电荷是守恒的（实际上并未在原子I上增加一个额外的电荷），因此存在一个人为引入的项img。该项对应于SCF计算前img下nI的变化（无库伦相互作用），必须从上述U的计算公式中减掉，所有我们可以得到：

img

这里，偏导数应该理解为只有该原子具有nI的变化，其他原子没有变化，所以有：

img

其中，当我们将一个αI放在一个原子上，而将其他原子的αI设为零，ΔαI和ΔnI是从上述计算中得到的数值。但是其他原子上的ΔnI不为零。

为了计算img，我们首先计算：

img

然后得到：

img

注意，对于Koopmans条件，这实际上是在假设，当存在一个ΔnI时，所有其他ΔnI应该为零。实际上，也可以避免这种情况，也就是说，当存在一个ΔnI（有一个外部ΔαI和变化引起）时，在SCF响应中，其他原子的ΔnI时，所有其他可以是任意值。此时，U可以被定义为：

img

results matching ""

    No results matching ""