本帖最后由 Shannon 于 2014-11-28 12:26 编辑
下面是 氦原子HatreeFock的code。 用
guess=[0.7,1,2 ]
[E,xxx]=fminsearch(@Helium,guess,optimset('Display','on','MaxIter',10000,'PlotFcns',@optimplotfval ))
fminseach 方法:
(Lagarias, J.C., J. A. Reeds, M. H. Wright, and P. E. Wright, "Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions," SIAM Journal of Optimization, Vol. 9 Number 1, pp. 112-147, 1998.)
搜索。 能够得到E=-2.86167的HF最好解。
function E = Helium_final(input)
c1=input(1);c2=0.180687;
% c2=(1-c1^2)^0.5;
�=0.5;c2=0.5
Z=2; % 核电荷数。2,不需要改
% zeta(1)=1.45363;
% zeta(2)=2.91093;
phi = @(z,r) sqrt((z^3)/pi)*exp(-z*r); % 在此调教波函数
% c1=input(1);
% c2=input(2);
zeta(1)=input(2);
zeta(2)=input(3);
sfunc= @(a,b,r) phi(zeta(a),r).*phi(zeta(b),r).*(r.^2);
it=1
while (it<10) % 用内循环求解C2的值,确保函数 normalization。
%方法是newton 法
S= @(a,b) 4*pi*quad(@(r) sfunc(a,b,r),0,10);
% sreal=8*(zeta(1)^3*zeta(2)^3)^0.5/(zeta(1)+zeta(2))^3;
% hreal=(sreal/2)*(zeta(1)*zeta(2)-Z*(zeta(1)+zeta(2)));
% normalizatoin= quad2d (@(r1,r2) bigphi(c1,c2,r1,r2),0,10,0,10 )
% 这个是normalization 的函数: 必须等于1 (c1^2+2*c1*c2*S(1,2)+c2^2)
Deltac2 = (2*c1*S(1,2)+2*c2)\(1-(c1^2+2*c1*c2*S(1,2)+c2^2));
it=it+1;
c1^2*S(1,1)+2*c1*c2*S(1,2)+c2^2*S(2,2)
c2
c2=c2+Deltac2;
if (abs(Deltac2)<0.001)break;
end
end
hfunc=@(a,b,r) phi(zeta(a),r).*phi(zeta(b),r).*(4*pi*r);
h = @(a,b) -(zeta(b)^2*S(a,b))/2-(Z-zeta(b))*quad(@(r) hfunc(a,b,r),0,10);
% 下面是双电子积分,6-dimension,I think, at least 3 dimension need to be killed to lower
% computational cost
% 估计code效率 很低,但比较傻瓜式
gfunccomp = @(z1,z2,z3,z4,r1,r2) phi(zeta(z1),r1).*phi(zeta(z2),r2).*(1./abs(max(r1,r2))).*phi(zeta(z3),r1).*phi(zeta(z4),r2)*(16*pi^2).*r1.^2.*r2.^2;
gpiece = @(x1,x2,x3,x4) quad2d(@(r1,r2) gfunccomp(x1,x2,x3,x4,r1,r2), 0,10,0,10);
g = @(gi,gj) c1^2*gpiece(gi,1,gj,1)+c1*c2*gpiece(gi,1,gj,2)+c2*c1*gpiece(gi,2,gj,1)+c2^2*gpiece(gi,2,gj,2);
E=(2*((c1^2)*h(1,1)+2*c1*c2*h(1,2)+(c2^2)*h(2,2))+(c1^2)*g(1,1)+2*c1*c2*g(1,2)+c2^2*g(2,2)); |