计算化学公社

标题: 量子化学学习用RHF方法的code，更新到RHF/6-311G解第二周期元素 [打印本页]

作者
Author: Shannon 时间: 2014-10-20 22:04
标题: 量子化学学习用RHF方法的code，更新到RHF/6-311G解第二周期元素
本帖最后由 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));

作者
Author: Shannon 时间: 2014-11-2 02:56
本帖最后由 Shannon 于 2014-11-26 00:05 编辑

次楼层的code是RHF法解H2分子。源文件来自网上，有很多修改，外加了一个可视化电子密度的小程序。
截图如下

作者
Author: Shannon 时间: 2014-11-2 03:00
画图的小程序在此楼：附件是code中用到的各个高斯函数积分，matlab绘图的底层文件等。
xvals = -3:0.1:3;
yvals = -3:0.1:6;
zvals = -3:0.1:2;
nel = numel(xvals);


orbit=2;
psi1=integrate_wave(xvals,yvals,zvals,d,c,spreads,centers,1,orbit);
psi2=integrate_wave(xvals,yvals,zvals,d,c,spreads,centers,2,orbit);
edensity=psi1.^2+psi2.^2;
%  surf(xvals,yvals',psi2+psi1)
maximus=max(max(max(edensity)));
ratio=0.5;
figure;
   %%
   %Volume visualization
   currentaxis=gca;
[X,Y,Z]=meshgrid(yvals,xvals,zvals);
colordata= isovalue(edensity);
   hold on

   for ii=1:5 %循环用来画圈圈
%             clear f;
         [v,f] = isosurf(edensity, [],  (0.5+(ii*0.1))*maximus, 0, 0); v=v'; f=f';%此函数来自底层mex程序，也就是c/c++编的，矩阵排列和matlab习惯不同 %改成matlab的排列方式
      v= axis_dim (X,Y,Z,v); %修改坐标到xyz的尺度
p(i) = patch('Vertices', v, 'Faces', f, 'EdgeColor', 'none','FaceVertexCData', colordata, ...
      'FaceColor', 'red','facelighting','flat','facealpha',ii*0.1);
   drawnow
   end
  %画个密度投影图，然后移动到底下
surface1 = pcolor(yvals,xvals,sum(edensity,3));
   set(surface1,'zdata',-ones(61,91));
  drawmolecule(currentaxis,Rnuc);
   view(3)
camlight

作者
Author: 516518 时间: 2014-11-3 08:44
学习，好贴

作者
Author: Shannon 时间: 2014-11-13 11:05
刚写完的HF解第二周期元素的code。 code的双电子积分器还在完善中，解CO的能量和高斯接近，势能矩阵，重叠矩阵和动能矩阵完全一致，双电子积分不完全对的上。解Ne的能量的结果以及各项双电子积分结果和高斯结果是完全一样的。希望数学高手大大们能帮助完善下双电子积分器代码……已经花了接近20小时Debug，双电子积分还是与高斯输出的结果有差别。
算法来自附件中的pdf文件。

作者
Author: Shannon 时间: 2014-11-18 02:28
更改了单双电子积分的解法，用了McMurchie-davidson的算法。结果已经和高斯完全一致了。算法来自helgaker的ppt。

作者
Author: helpme 时间: 2014-11-18 11:27
人才啊。

作者
Author: Shannon 时间: 2014-11-20 05:06
大幅优化了双电子积分，目前计算完CO的能量需4秒左右（CPU-core i7 3.4GHz）。
对重叠矩阵强-制使用对称性，已经可以得到稳定的势能面了。（S矩阵微小的不对称都会对得到的轨道产生巨大影响，必须得强-制使用对称性）

Co的势能面如图。

作者
Author: psfan 时间: 2014-11-25 08:04
厉害！下一步进行代码并行优化

作者
Author: Shannon 时间: 2014-11-25 23:40

psfan 发表于 2014-11-25 08:04
厉害！下一步进行代码并行优化

我的matlab没买并行工具包……学校没钱，哈哈。Matlab工具包里面有一堆好东西，可惜卖的实在是贵了点，一个100刀，有点黑

作者
Author: Shannon 时间: 2014-11-25 23:47
本帖最后由 Shannon 于 2014-11-28 12:25 编辑

增加了一个画波函数的小程序，另外增加了一个cauchy-schwartz筛选。稍微提高了点速度。算CO需要2.9s。 CO2需要14秒的样子。CO2的LUMO

HOMO

能量第三高的轨道

附件是此程序，另外的一个附件是6-311G计算CO的玩具代码，结果和高斯是完全一致的，速度灰常慢，约需4分钟。

作者
Author: psfan 时间: 2014-11-27 16:15

Shannon 发表于 2014-11-25 23:40
我的matlab没买并行工具包……学校没钱，哈哈。Matlab工具包里面有一堆好东西，可惜卖的实在是贵了点，一 ...

别的可以放一放，并行包应该买一套

作者
Author: sobereva 时间: 2014-11-27 17:52
并不并行现阶段来说无所谓。
真要并行的话就说明已经进入实用阶段了，那就肯定得写成C/Fortran了，到时候直接用openmp并行就行了，相当方便。

作者
Author: Shannon 时间: 2014-11-28 02:31

sobereva 发表于 2014-11-27 17:52
并不并行现阶段来说无所谓。
真要并行的话就说明已经进入实用阶段了，那就肯定得写成C/Fortran了，到时候 ...

确实是得用Fortran/C，刚才试用了下用matlab转C的工具包（免费试用版），速度加快了2.5 倍多。

作者
Author: trewqaz 时间: 2014-12-14 02:17
人才啊，可以搞中国的Gaussian了

作者
Author: Shannon 时间: 2014-12-14 03:47

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

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

作者
Author: trewqaz 时间: 2014-12-14 04:11
高斯的积分核心算法是“棱镜”

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

作者
Author: qwoop 时间: 2014-12-30 15:26
大赞楼主，学习量化编程很好的资料

作者
Author: ximi1986 时间: 2015-7-1 16:26
大赞一个！传说内联函数，比匿名函数会快些~~~有空也好想试试啊~~~~~~

作者
Author: 一声叹息010 时间: 2016-10-24 21:12

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

h2.m中忘记声明zetaH2了

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

作者
Author: sobereva 时间: 2017-7-4 01:51

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

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

作者
Author: zkw 时间: 2017-7-4 13:55

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

O(∩_∩)O谢谢，可是没学过Fortran

有MATLAB版的吗？理论方面有什么资料可以推荐一下么

作者
Author: zkw 时间: 2017-7-4 14:07

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

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

作者
Author: Warm_Cloud 时间: 2017-7-4 15:15

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

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

作者
Author: zkw 时间: 2017-7-4 17:16

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

谢谢~希望我能尽快看懂

作者
Author: sobereva 时间: 2017-7-5 06:22

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

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

欢迎光临计算化学公社 (http://bbs.keinsci.com/)