KSSOLV之一(KS方程求解器-matlab版)_奇趣物理实验

Today, I find a paper,titled "A MATLAB Toolbox for Solving the Kohn-Sham Equations.pdf" . As I am very familiar with matlab, not fortran, it is easy for me to read the code. It seems that this toolbox is suitable for education, instead of a calculation.

Besides, the developer of the software pocket is a chinese:)

KSSOLV's developer

In the file "A MATLAB Toolbox for Solving the Kohn-Sham Equations", there are some example to run KSSOLV in matlab. I test "sih4_setup.m" and find it is very easy to follow:)

1. run sih4_setup.m

2. typing mol on the command line

3. typing [Etotvec, X, vtot, rho] = scf(mol) on the command line

4. typing isosurface(rho) on the command line

补记：

KSSOLV is targeted primarily towards users who are interested in the numerical analysis aspects of electronic structure calculations. The developers focus is on numerical algorithms and how they can be easily prototyped within KSSOLV.

Exc(ρ) is known as the exchange-correlation energy, which is a correction term used to account for energy that the non-interacting reference system fails to capture. KSSOLV developers use LDA proposed in [Kohn and Sham 1965]. (2009-09-16)

—————————————————————————————————————————————

What I do not fully anderstand ：（

   1:When KSSOLVE try to calculate the integral "FT(f) = int_{0}^{rst} f(r) e^{i*g^Tr} d^3r" in pseudoinit.m file, you write the following content at 223 line.
          My question is: How do it try to calculate the integral?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            r1 = r(idx);

r2 = r(ip1);

% construct a 2-D grid (r,g1) so that the integration

% can be vectorized

X1 = r1*g1;

X2 = r2*g1;

X21 = X2-X1;

COSX12 = cos(X1)-cos(X2);%cos[(r1+delt_r)*g1]-cos(r1*g1)

SINX21 = sin(X2)-sin(X1);%sin[(r1+delt_r)*g1]-sin(r1*g1)

b1 = vloc(idx).*r1;

b2 = vloc(ip1).*r2 - b1;%[vloc(ip1).*r2-vloc(idx).*r1]为delt[vloc(r).*r]

B1 = repmat(b1,1,mnq);

B2 = repmat(b2,1,mnq)./X21;

C1 = COSX12.*B1 - X21.*B2.*cos(X2) + B2.*SINX21;

s = sum(C1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2:In pseudoinit.m file, KSSOLVE write the following content at 262 line.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% now add the analytical solution of the second (long range)

% integral

s=s-ch*4*pi*cos(g1*r(ifirst))./g1.^2;

s1=s1.*g124pi;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

          3:In pseudoinit.m file, you write the following content at 368 line.
             My question is: Why you use sin(X)./X instead of besselj function?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              A = sin(X)./X;
%             A = sinc(X/pi);
%             A = besselj(0,X);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

          4:In getvion.m file, KSSOLVE write the following content at 106 line.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xvec = (qvec-qi2(iqvec))./(qi2(iqvec+1)-qi2(iqvec));

% some kind of interpolation) ....

f1vec = 1-xvec-xvec.*(1-xvec)/2;

f2vec = xvec + xvec.*(1-xvec);

f3vec = -xvec.*(1-xvec)/2;

yvec = vq(iqvec,itype).*qi2(iqvec).^2.*f1vec ...

+ vq(iqvec+1,itype).*qi2(iqvec+1).^2.*f2vec ...

+ vq(iqvec+2,itype).*qi2(iqvec+2).^2.*f3vec;

y2vec = vqT(iqvec,itype).*qi2(iqvec).^2.*f1vec ...

+ vqT(iqvec+1,itype).*qi2(iqvec+1).^2.*f2vec ...

+ vqT(iqvec+2,itype).*qi2(iqvec+2).^2.*f3vec;

iql = find(qvec>=1e-6);

yvec(iql) = yvec(iql)./qvec(iql).^2;

y2vec(iql) = y2vec(iql)./qvec(iql).^2;

y1vec = rhoq(iqvec,itype).*f1vec+rhoq(iqvec+1,itype).*f2vec+rhoq(iqvec+2,itype).*f3vec;

vionf(inz) = vionf(inz) + yvec.*ccvec/vol;

vionTf(inz) = vionTf(inz)+ y2vec.*ccvec/vol;

rhof(inz) = rhof(inz) + y1vec.*ccvec/vol;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

5.I notice there is not getwq4solid.m file. I try to calculate sibulk_setup.m with Ham(mol).
Therefore, I get an error "Undefined command/function 'getwq4solid'."

阅读┊ 收藏 ┊ 喜欢 ▼ ┊打印┊举报/Report

前一篇：GW theory II

后一篇：UNIFIED PSEUDOPOTENTIAL FORMAT 1

新浪BLOG意见反馈留言板　欢迎批评指正