A common approximation for a ``quantum dot,'' a relatively small number of electrons trapped in a potential, is to model the trapping potential as a harmonic oscillator (which actually is the potential V(r) which you have currently in ``dft.m'') and to treat the electron-electron interactions within density-functional theory. To perform such a calculation, replace the code block in ``dft.m'',
%# Converge using steepest descents
W=lm(W,20); %# 20 iterations of simple lm() to get close to the minimum
[Wlm, Elm]=lm(W,50); %# 50 iterations of lm() from W, while recording results
.
.
.
semilogy(its,Elm-18,its,Epclm-18,its,Ecg1-18,its,Ecg2-18,its,Ecg3-18);
%# Set W to best result for plotting
W=Wcg3;
pause;
which compares convergence of various minimization methods
(now inappropriately assuming that the converged result is 18) with%# Be sure to show lots of digits! format long %# Converge a little using steepest descents W=sd(W,10); %# Converge more with preconditioned line minimization W=pclm(W,20); %# Finish convergence with full preconditioned conjugate gradients W=pccg(W,50,3);Be sure to remove any ``pause;'' statements so that the code will also run through the part which gets the final wave functions and plots them.
Run your code. You should find a converged value of the total energy for this quantum dot of 8 electrons (4 states, with two electrons each) near Etot=43.337. You should also find states with similar s and p characters as to what you had before, but now with eigenvalues of approximately 5.509, 6.949, 6.949, 6.949.