一. DFT简述（使用Atomic units）密度泛函理论主要由 Kohn 和沈吕九在半个世纪之前创造。用于求解多电子体系基态性质，其主要思想是这样的：
首先：假设一组波函数 {ψi(x⃗ )} 描述了电子们的状态
然后：由 {ψi(x⃗ )} 导出体系能量的表达式(包含电子动能，库伦能等)：
1. 电子动能：Tel=−12∑ni=1∫ψ∗i(x⃗ )∇2ψi(x⃗ )d3x
2. 电子与原子之间相互作用能： Vext=∫n(x⃗ )Vnuc(x⃗ )d3x
3. 电子受其他电子的库伦能：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 中考虑空间角动量，自旋角动量，由此研究其磁性？
参考资料：
1h DFT code in matlab
Laplace 算子的离散化方法及离散空间中的分离变量法
快速求解矩阵最小（或最大）本征值本征向量的Davidson method
Gaussian compensation charge
