DFT code in one hour with Matlab - Part 1: The Hydrogen Atom with finite difference
Channel: Mikael Kuisma
Date: 03/10/2012

I was asked to make a tutorial about creating a density functional theory code with Matlab. We'll go very fast, since we have only one hour to code a complete periodic finite difference DFT code. I recommend stopping and using Google for important keywords such as 'finite difference stencil' or 'atomic units' etc. The code is going to be able to calculate real systems, but I have not yet decided will I go for the pseudo potentials or PAW approach, but we will see as the code progresses. Stay tuned for Part 2: The helium atom where we calculate also the exchange correlation energy and solve the poisson equation for Hartree potential.

Note: I seem to talk about row vector the full tutorial, but I mean a column vector.

http://en.wikipedia.org/wiki/Finite_difference_method
http://en.wikipedia.org/wiki/Kronecker_product
http://en.wikipedia.org/wiki/Laplace_operator

Excercises:
1. Calculate the convergence properties of the calculation with respect to grid spacing. Make a logarithmic plot of the error and estimate n of the error bound expession O(h^n) when using this 1 -2 1 stencil. Try to eliminate the error from finite size box.
2. Use a higher order finite difference stencil and repeat excercise 1.

% DFT code in one hour with Matlab - Part 1: The hydrogen atom
g = 60; g3 = g^3;
p=linspace(-7,7,g);
[X,Y,Z] = meshgrid(p,p,p);
h=p(2)-p(1);
X=X(:);Y=Y(:);Z=Z(:);
R=sqrt(X.^2+Y.^2+Z.^2);
Vext = -1./R;
e=ones(g,1);
L=spdiags([e -2*e e], -1:1, g, g)/h^2;
I = speye(g);
L3 = kron(kron(L,I),I) + kron(kron(I,L), I) + kron(kron(I,I),L);
E = eigs(-0.5*L3+spdiags(Vext, 0, g3, g3), 1, 'sa');
disp(['Total energy for Hydrogen atom ' num2str(E*27.21,5) ' eV.']);