使用自洽场方法求解 Kohn-Sham 方程一. DFT简述(使用Atomic units)

密度泛函理论主要由 Kohn 和沈吕九在半个世纪之前创造。用于求解多电子体系基态性质，其主要思想是这样的：

首先：假设一组波函数 {ψi(x⃗)} 描述了电子们的状态

然后：由 {ψi(x⃗)} 导出体系能量的表达式(包含电子动能，库伦能等)：1. 电子动能：Tel=−12∑ni=1∫ψ∗i(x⃗)∇2ψi(x⃗)d3x2. 电子与原子之间相互作用能： Vext=∫n(x⃗)Vnuc(x⃗)d3x3. 电子受其他电子的库伦能：VH=12∫ϕ(x⃗)n(x⃗)d3x，其中ϕ(x⃗) 为电子电荷密度形成的库伦式，可以通过求解 Poisson 方程(∇2ϕ=−4πn)得到 4. 对库伦能的修正(交换能)：Ex=∫fx(n(x⃗))dV

注：其中：n(x⃗)=∑iψ∗iψi 为电子密度；对于交换能Ex，有各种近似方法，这里使用最简单的局域密度近似。具体介绍可以在这里找到。寻找更好的交换能函数，是凝聚态中一个十分重要的研究课题。

综上，可以得到总能量表达式为：

E[{ψi(x)}]=Tel+Vext+VH+Ex

最后：利用变分法，求解使体系能量最小时波函数需要满足的条件，就可以得到著名的Kohn-Sham 方程：

[−ℏ22m∇2+Vext(x⃗)+ϕ(x⃗)+Vx(x⃗)]ψi(x⃗)=ϵiψ(x⃗)

其中Vx 为交换势, 具体表达式见这里(局域密度近似)

令人吃惊的是，它与单电子 Schrodinger 方程是如此的相似。不同之处只是在于，由于电子之间相互作用，Hamiltonian 中的势能项包含了电子密度。这使得K-S方程成为一个非线性方程(哈密顿量与波函数有关)，与我们熟知的本征值问题不太一样。接下来，介绍如何编程解这个方程。二. 求解 Kohn-Sham 方程的自洽场(self-consistent field method SCF)算法

Initialization: 1. 确定电子个数(N)2. 用外势能近似总势能，即 Vtot=Vext，得到近似Hamiltonian；Iteration：1. 求Hamiltonian最小的 N/2 个本征值，及对应的本征函数 ψi(x⃗)(每个态上占据两个电子)2. 由得到的本征函数集 {ψi(x⃗)} 求交换势(ϕ(x⃗))和库伦势(Vx(x⃗))3. 更新Hamiltonian: H=−ℏ22m∇2+Vext(x⃗)+ϕ(x⃗)+Vx(x⃗)4. 判断当前结果是否满足要求，如果满足，就跳出循环三. 算法的 Matlab 实现

使用条件及近似方式：只考虑电子成对占据某一能态的原子；使用LDA近似；使用空间离散化的方法求解Hamiltonian的本征值；使用Dirichlet边界条件(边界处概率密度为0)；以4个电子为例。

主程序：

%For double occupationN=4;% num of enectronsg=50% num of latticesg3=g^3;p=linspace(-5,5,g);% one dimensiton space lattice[X,Y,Z]=meshgrid(p,p,p);% three dimension space latticeh=p(2)-p(1);% latice spacingX=X(:);Y=Y(:);Z=Z(:);% all elements of arraty as a single columnR=sqrt(X.^2+Y.^2+Z.^2);% distance from the centerVext=-N./R;% potential energy(2 protons)e=ones(g,1);L=spdiags([e-2*ee],-1:1,g,g)/h^2;% 1D finite difference Laplacian (with 0 boundary condition)I=speye(g);L3=kron(kron(L,I),I)+kron(kron(I,L),I)+kron(kron(I,I),L);% extend Laplacian to 3 DVtot=Vext;%initial guessncomp=exp(-R.^2/2);%compensation charge(for poisson equation)ncomp=-N*ncomp/sum(ncomp)/h^3;ncomppot=-N./R.*erf(R/sqrt(2));%solution of poisson eq. of compensation charge%%%%%%%%initial guess for N = 4%%%%%%%%%%%%%%%%%%E=[-4-1];PSI=[exp(-3.7*R)exp(-0.17*R)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%foriter=1:5%Do the main loop not for everH=-0.5*L3+spdiags(Vtot,0,g3,g3);% Hamiltonian of Helium[PSI,E]=lowestNEigen(H,PSI,E,N,5);fori=1:(N/2)PSI(:,i)=PSI(:,i)/norm(PSI(:,i));endPSI=PSI/h^(3/2);%normalize PSIn=0;fori=1:(N/2)%calculate density of electronn=n+2*PSI(:,i).^2;endVx=-(3/pi)^(1/3)*n.^(1/3);%exchange potantial (LDA)Vh=cgs(L3,-4*pi*(n+ncomp),1e-7,400)-ncomppot;%Hartree potantial(solution of poisson eq.: L3 Vh = -4*pi*n)Vtot=Vx+Vh+Vext;%total potantialT=0;fori=1:(N/2)%calculate Kinetic energyT=T+2*PSI(:,i)'*(-0.5*L3)*PSI(:,i)*h^3;%Kinetic enerty(expactation value of Kinetic energy oprator)endEext=sum(n.*Vext)*h^3;%external energyEh=0.5*sum(n.*Vh)*h^3;%hartree energyEx=sum(n.*Vx*(3/4))*h^3;%How come the 3/4????????????????Etot=T+Eext+Eh+Ex;moreoff;%see the disp  in the loop, but why?E%  disp(['Kinetic energy ' num2str(T,5) ]);%  disp(['Exchange energy ' num2str(Ex,5) ]);%  disp(['External energy ' num2str(Eext,5) ]);%  disp(['Potential energy ' num2str(Eh,5) ]);disp(['Total energy'num2str(Etot,5)]);end%scatter3(X(1:10:g3),Y(1:10:g3),Z(1:10:g3),n(1:10:g3)*1000);scatter3(X,Y,Z,n*4);//showelectrondensity

用Davidson method 求解本征值和本征向量：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  Usage: Get lowest eigenvalue and cooresponding Evetor by Davidson's method%  H: discrete Hamiltonian%  PSI, E： initial guess of eigenvectors and eigenvalue%  N: number of eigen pair needed%  iter: iteration number%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[PSI, E] =lowestNEigen(H, PSI, E, N, iter)fori=1:iter%"for loop": for faster convergentRR=H*PSI-PSI*diag(E);%All the Residual vectors at oncePSIe=[PSIRR];%SubsbaceHH=PSIe'*H*PSIe;%For transform subspace to a space inwhich H is diagnolizedSS=PSIe'*PSIe;%For transform subspace to a space inwhich H is diagnolizedHH=HH+HH';SS=SS+SS';%Ensure they are Hermition matrix, so that the eigen value will return in order[U,E]=eig(HH,SS);%For transform subspace to a space inwhich H is diagnolizedE=diag(E);%Diagnal Matix to vectorPSIe=PSIe*U;%%SIe' * H * PSIePSI=PSI(:,1:(N/2));E=E(1:(N/2));%Pick lowest N/2 eigen vectors and valuesend结果

Total energy: -10.793

电子密度分布图：

接下来

误差较大，如何升级？

如何在 DFT 中考虑空间角动量，自旋角动量，由此研究其磁性？参考资料：