量子化学学习用RHF方法的code，更新到RHF/6-311G解第二周期元素

Shannon

下面是 氦原子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));

sobereva

zkw 发表于 2017-7-4 13:55
O(∩_∩)O谢谢，可是没学过Fortran有MATLAB版的吗？理论方面有什么资料可以推荐一下么

就是个语法的事。搞懂了原理、算法，用哪个语言写都差不多。

zkw

Warm_Cloud 发表于 2017-7-4 15:15
因为写电子积分的难度更大，更复杂，写快也更难。如果想看自洽场的可以看看这个 http://bbs.keinsci.com/ ...

谢谢~希望我能尽快看懂

Warm_Cloud

zkw 发表于 2017-7-4 14:07
而且不太明白为什么都是与积分有关的，还有想请教一下怎么用自洽场的方法求解波函数啊

因为写电子积分的难度更大，更复杂，写快也更难。如果想看自洽场的可以看看这个 http://bbs.keinsci.com/forum.php?mod=viewthread&tid=6139

zkw

zkw 发表于 2017-7-4 13:55
O(∩_∩)O谢谢，可是没学过Fortran有MATLAB版的吗？理论方面有什么资料可以推荐一下么

而且不太明白为什么都是与积分有关的，还有想请教一下怎么用自洽场的方法求解波函数啊

zkw

sobereva 发表于 2017-7-4 01:51
仔细看
利用wfn文件计算电子密度的代码的编写方法
http://sobereva.com/182

O(∩_∩)O谢谢，可是没学过Fortran有MATLAB版的吗？理论方面有什么资料可以推荐一下么

sobereva

zkw 发表于 2017-7-3 15:40
楼主，我的期末MATLAB大作业准备画分子的电子云图，但是我还没有学量子力学，只是初步的了解了一下还有很多 ...

仔细看
利用wfn文件计算电子密度的代码的编写方法
http://sobereva.com/182

zkw

楼主，我的期末MATLAB大作业准备画分子的电子云图，但是我还没有学量子力学，只是初步的了解了一下还有很多不懂的地方，想具体了解一下楼主用的什么方法，能给一些参考文献或者资料让我学习一下吗

一声叹息010

Shannon 发表于 2014-11-2 02:56
次楼层的code是RHF法解H2分子。源文件来自网上，有很多修改，外加了一个可视化电子密度的小程序。
截图如 ...

h2.m中忘记声明zetaH2了

ximi1986

大赞一个！传说内联函数，比匿名函数会快些~~~有空也好想试试啊~~~~~~

qwoop

大赞楼主，学习量化编程很好的资料

Shannon

将双电子积分中用到的递归算法全部暴_力_拆开。 现在代码效率能搞定高角动量的极化函数了（虽然目前还没写进去）。附件是写程序的程序，用来解析 递归代码成 简单的加减乘除代码的代码。
迭代很多时候会计算 重复 东西，导致效率不高，尤其是高角动量时这种效率浪费特别大。 全部拆开之后速度又提高了很多。目前计算CO，STO2G需要1S。

trewqaz

高斯的积分核心算法是“棱镜”

Shannon

trewqaz 发表于 2014-12-14 02:17
人才啊，可以搞中国的Gaussian了

这个code真的只是用来学量化的玩具…高斯G6311（d）计算苯的 单点能量才只需要几秒钟（这个matlab代码大约需要以小时计算的时间），写过code之后方才发现高斯里面的数学优化量大的 简直丧-病。高角动量基底函数和 contracted Gaussian 的积分速度，高斯快的简直不可思议。
写成matlab一个 大的好处是用来学CI/NBO，了解轨道正交性之类 是怎么搞得非常有用。

trewqaz

人才啊，可以搞中国的Gaussian了