Kohn-Sham方程及其解法

1Kohn-Sham方程及其解法1.Kohn-Sham方程如果原子核不动，材料可以看成是“外场下的非均匀电子气”，体系的基态性质是其电子密度的唯一泛函，而该电子密度)(rn满足Kohn-Sham方程：写在一起就是：Niirrn12)()(对所研究的体系解出该Kohn-Sham方程，就可得到其电子密度，而体系的性质由该电子密度决定：物理量F=F[n]2.Kohn-Sham方程中的各项：第1项：动能项（电子的动能，原子核不动）第3项：称为Hartree势（哈特利势），可以类比为库仑势。(')'21122()()()()()'([()]())()(()()())nrHreffiiiNiieffxcxcxcrrvrdrvrvEnVnrVrrVrrrrVrnrr)()()]}('')()([21{2rrnVrdrrrnrViiixcext2第2和第4项需要很多的说明。第2项：外势项。由原子核（或原子芯）的空间排列（即材料的结构）构成。原子由：{原子核＋全部的电子}构成all-electronscal.或{原子芯＋价电子}构成pseudo-potentialcal.由于全电子的计算工作量大（波函数在靠近原子核的地方振荡很厉害），非全电子的计算通常有优势。我们这里就将使用非全电子的计算（VASP程序包）。所以，需要有“赝势”的概念：赝势方程：如果不考虑原子的芯电子，则原子就成为“赝原子”（原子芯＋价电子）。这时，价电子运动受处的势场就相当于来自原子芯的“赝势”（原来是原子芯内所有电子提供的势场）。可以证明，将薛定鄂方程中的势能换成赝势，EVH)(pspspsEVH)(则存在相应的赝波函数，使体系的本征值不变。固体物理中的“正交化平面波”法，实际上对此做了证明:3从头赝势：赝势下，本征值是真的还是不够的（还不能研究与波函数或电荷密度相关的信息），我们希望波函数还要是真的。所谓从头赝势，就是能够在某个rc半径之外，使赝原子的能量本征值以及赝波函数都和“整个原子”时的解一致！在rc半径之内，赝势应该尽量的平滑（则使波函数的振荡很小），赝波函数没有节点。4如何构造从头赝势：参考文献也附上。对绝大多数的原子，赝势都已经有人构造完成。5PAWrepresentation:现在大家使用ProjectorAugmented-Wave（PAW），它结合了赝势和缀加平面波法。参考文献附上。（我们也只要直接调用即可，如果不管细节的话。）第4项：称为交换-关联势（exchange-correlation势）：通常的两种近似处理方式（1）LDA近似：LocalDensityApproximation那么，方程中的交换关联势近似为)])([)(()(][)]([rnrnrdndnnErnVxcxcxc实际的应用中，需要采用参数化的办法，例如：交换能..458.014943)(343)]([3/123/13/1uarrrnrnssx其中3/143nrs。6关联能，常用的是T.P.Perdew和A.Zunger根据D.M.Ceperley和B.L.Alder的用最精确的Monte-Carlo方法计算的均匀电子气的结果：)1(ln0040.00232.0ln0622.00960.0)1()3334.00529.11(2846.0ssssssssscrrrrrrrrr（2）GGA近似：GeneralizedGradientApproximation介绍VASP程序包中常用的两类GGA函数：1）Perdew-Wang’91（PW91）交换关联函数：45211210043211homsinh1)()(sinh12sasasaseaasasasxx其中19645.01a，7956.72a，2743.03a，1508.04a，004.05a。tsnnHLDACC,,210021042242)(121ln2,,scccetCnCCtAAtAtttsnH1/)(2122rnCeA7其中09.0，0667263212.0，7559.150cC，003521.01cC，nkrnts2)(而2/14Fskk，)()(homrnrnncc。2）Perdew-Burke-Ernerhof（PBE）交换关联函数：rdsnFrnrnrnEXCxxc,,)()()(hom其中)(r是局域密度，是相对自旋极化率，)(2)(rnkrnsF，则：211)(ssFx其中21951.032，而066725.0是与二级梯度展开有关的。对所有的r都有804.1)(sFx，则804.0，Perdew-Burke-Ernzerhof采用的是804.0。关联能可以写成与Perdew-Wang’91类似的形式，即：rdtrHrrnnnEsscGGAc,,),()(],[hom其中42222302111ln,,tAAtAttaetrHs这里)(2)(rnkrnts，2/104akkFs是Thomas-Fermi屏蔽波矢，211)(3/23/2是自旋放大系数，的值与交换项中的相同，即066725.0，22ln1，函数A的形式如下：1023hom1expaenAc。*****现在大家通常都使用GGA近似来计算。实际操作中，也只8要选择恰当的近似方法：什么LDA和什么GGA即可，如果不关心细节的话。方程的解法：3．Kohn-Sham方程是一个自洽方程：方程：（是一个自洽方程）或写成：)()()]([rrrnHii,其中.1)()(*)(occiiirrrn.即在哈密顿量H中含有需要求解的未知“波函数（这里是Kohn-Sham轨道）”（即：未知的需要求解的电荷密度或“波函数”被嵌套在必须已知的哈密顿量中），故方程是一个自洽方程，必须做自洽求解：自洽解法，常见的步骤：(a)从一个随意给定的0出发，构造电荷密度：.100)()(*)(occiorrrn,从而得知哈密顿量H[n0]（这样哈密顿量就确定了，但通常还不是系统真正的H）就可以解方程：)()()]([0rrrnH………eq.(1)得到1（这样得到的1一般说还不是体系的解，因为刚才的哈密顿量)()()]}()([21{2rrnrViiixc9还是猜测的）(b)但现在可以有了更好的出发点：1,可以再构造密度：.1111)()(*)(occirrrn,从而得知H[n1]再解方程：)()()]([1rrrnH………eq.(2)可以得到2。（2应该比1更加趋近于最后的解）【为了数值求解上的收敛，实际的做法是：1n是n与1n的恰当混合。】……重复以上过程，直到自洽为止（即1n与n相差很小）！可见，在以上整个的自洽求解过程中，实际上“自洽方程”的求解问题最后可以归结为求解：一个已知哈密顿量H的方程。4．已知H的“Kohn-Sham方程”的两种常见解法1.矩阵的对角化（标准的）2．迭代法（1）矩阵对角化：对于一个已知其哈密顿量H的Kohn-Sham或薛定谔方程：标准做法：(1).可用一组正交归一的完整集}{i来展开：1)(iiirC)()(ˆrrH10【实际中为了可以处理，必须做切断，以便数值解，即：NiiirC1)(,N取到足够大为止】(2).代入方程，则：NiiiNiiirCrCH11)()((3).两边同乘)(*rj,再对r空间积分，则：rdrrCrdrHrCNiijiNiiji11)()(*)(ˆ)(*在已知哈密顿量和你自己选择的基函数的情况下，以上积分都是确定值。记ijijHrdrHr)(ˆ)(*,而ijijrdrr)()(*（这里假设基函数是正交归一的）则：jiNiiijNiiCHC110][1jiijNiiHC……….(A)以上线性方程组有非零解的条件是其系数行列式为零：0jiijH也即：NNNNNNHHHHHHHHH.....................212222111211=0这样，K－S方程或薛定鄂方程的解转化成一个标准的“矩阵对角化”的数学问题。这至少可以使用标准的计算机程序来完成。11对上面矩阵进行对角化，可解出N个本征值i,每个本征值都可以代回方程（A）【方程（A）就成为一个已知系数的线性方程组】，就可以解出一组},1{NiCi，,即本征函数（波函数）。对于正交归一的完整集}{i的说明：目前有许多种基集的选择方式，也不一定要正交归一，或完整集。当展开的函数集是平面波时，则称：平面波法……..是APW是，就叫：APW法……LCAO，LAPW，LMTO方法练习1：一般地1)(iiirC，实际中必须做切断，以便数值解：NiiirC1)(，现假设N=2,请使用上述矩阵对角化方法求体系的本征值。练习2：以上使用正交归一的基函数ijijrdrr)()(*，如果基函数不是正交归一的，试推出其久期方程（即矩阵）。（2）迭代法Iterativemethods简介：[哈密顿量已知]还是可用一组正交归一的完整集}{i来展开：)()(ˆrrH121)(iiirC【一样做切断，以便数值解：NiiirC1)(】迭代法：从随意猜测的一组}{0出发。(进行bandbyband的计算，对每个k点，每个bandkn：)(1)随意猜测方程的解为1，则一般：11H就有“剩余”矢量residualvector111H[1是一个1×N矩阵，或说是一个N维矢量](2)下一次的猜测可用：112＋＝，则一般说，仍然有：22H，但是residualvector222H应该变小。(c)重复以上过程，只要方法合适，residualvector应该越来越小，…..，到residualvector近似零时，就是方程的解。bandbyband计算之后，再构筑.)()(*)(occnknknkrrrn，再重复进行。目前，迭代法对处理大的体系非常有用。因为可以不必存储N×N个矩阵元，当系统很大时，矩阵元数目是个海量。而且对目前的算法水平，迭代法的速度较矩阵对角化方法快很多。特别是并行计算技术的快速发展，更使迭代法的优势得到发挥。【矩阵对角化方法的并13行效率很低】***当然，若要不存储N×N个矩阵元，理论处理上必须能够使00'H，不能回归到使用H矩阵元。我们将使用的VASP程序包采用迭代法。