计算化学公社

 找回密码 Forget password
 注册 Register
Views: 27162|回复 Reply: 27
打印 Print 上一主题 Last thread 下一主题 Next thread

[算法与编程] 量子化学学习用RHF方法的code,更新到RHF/6-311G解第二周期元素

[复制链接 Copy URL]

127

帖子

1

威望

1231

eV
积分
1378

Level 4 (黑子)

跳转到指定楼层 Go to specific reply
楼主
本帖最后由 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));

评分 Rate

参与人数
Participants 5
威望 +1 eV +19 收起 理由
Reason
snljty + 5 好物!
sobereva + 1 + 5
卡开发发 + 2 好物!
pwang123 + 2
xmseet + 5

查看全部评分 View all ratings

127

帖子

1

威望

1231

eV
积分
1378

Level 4 (黑子)

2#
 楼主 Author| 发表于 Post on 2014-11-2 02:56:32 | 只看该作者 Only view this author
本帖最后由 Shannon 于 2014-11-26 00:05 编辑

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





Baidu IME_2014-11-1_14-55-58.jpg (246.23 KB, 下载次数 Times of downloads: 132)

Baidu IME_2014-11-1_14-55-58.jpg

h2.m

6.3 KB, 下载次数 Times of downloads: 26

127

帖子

1

威望

1231

eV
积分
1378

Level 4 (黑子)

3#
 楼主 Author| 发表于 Post on 2014-11-2 03:00:49 | 只看该作者 Only view this author
画图的小程序在此楼:附件是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

HF_code_111.zip

17.65 KB, 下载次数 Times of downloads: 75

评分 Rate

参与人数
Participants 1
eV +4 收起 理由
Reason
sobereva + 4

查看全部评分 View all ratings

22

帖子

0

威望

342

eV
积分
364

Level 3 能力者

4#
发表于 Post on 2014-11-3 08:44:35 | 只看该作者 Only view this author
学习,好贴

127

帖子

1

威望

1231

eV
积分
1378

Level 4 (黑子)

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

HF_FUN_solver_v1.zip

13.8 KB, 下载次数 Times of downloads: 67

如何推倒高斯函数.pdf

155.23 KB, 下载次数 Times of downloads: 104

评分 Rate

参与人数
Participants 2
eV +9 收起 理由
Reason
卡开发发 + 4 好物!
sobereva + 5

查看全部评分 View all ratings

127

帖子

1

威望

1231

eV
积分
1378

Level 4 (黑子)

6#
 楼主 Author| 发表于 Post on 2014-11-18 02:28:01 | 只看该作者 Only view this author
更改了单双电子积分的解法,用了McMurchie-davidson的算法。结果已经和高斯完全一致了。算法来自helgaker的ppt。

HF_fun_slv_v1.1.zip

25.01 KB, 下载次数 Times of downloads: 67

Helgaker_molecular_integral.pdf

266.54 KB, 下载次数 Times of downloads: 98

评分 Rate

参与人数
Participants 2
eV +7 收起 理由
Reason
卡开发发 + 3 赞!
sobereva + 4

查看全部评分 View all ratings

296

帖子

1

威望

2588

eV
积分
2904

科音成员

7#
发表于 Post on 2014-11-18 11:27:26 | 只看该作者 Only view this author
人才啊。
华北电力大学数理学院,理论与计算化学,团簇、表面的结构与反应机理。(招第一性原理计算,量子化学计算方向的教师、硕士/博士研究生)

127

帖子

1

威望

1231

eV
积分
1378

Level 4 (黑子)

8#
 楼主 Author| 发表于 Post on 2014-11-20 05:06:20 | 只看该作者 Only view this author
大幅优化了双电子积分,目前计算完CO的能量需4秒左右(CPU-core i7 3.4GHz)。
对重叠矩阵 强-制使用对称性,已经可以得到稳定的 势能面了。(S矩阵微小的不对称都会对得到的轨道产生巨大影响,必须得强-制使用对称性)

Co的势能面如图。




CO_P_surface.jpg (46.98 KB, 下载次数 Times of downloads: 115)

CO_P_surface.jpg

HF_fun_slv_v1.3.zip

40.8 KB, 下载次数 Times of downloads: 51

评分 Rate

参与人数
Participants 1
eV +3 收起 理由
Reason
sobereva + 3 хорошо!

查看全部评分 View all ratings

39

帖子

0

威望

132

eV
积分
172

Level 3 能力者

9#
发表于 Post on 2014-11-25 08:04:20 | 只看该作者 Only view this author
厉害!下一步进行代码并行优化

127

帖子

1

威望

1231

eV
积分
1378

Level 4 (黑子)

10#
 楼主 Author| 发表于 Post on 2014-11-25 23:40:19 | 只看该作者 Only view this author
psfan 发表于 2014-11-25 08:04
厉害!下一步进行代码并行优化

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

127

帖子

1

威望

1231

eV
积分
1378

Level 4 (黑子)

11#
 楼主 Author| 发表于 Post on 2014-11-25 23:47:17 | 只看该作者 Only view this author
本帖最后由 Shannon 于 2014-11-28 12:25 编辑

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


HOMO

能量第三高的轨道

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

co21.png (202.3 KB, 下载次数 Times of downloads: 107)

co21.png

co22.png (33.59 KB, 下载次数 Times of downloads: 123)

co22.png

co23.png (26.37 KB, 下载次数 Times of downloads: 122)

co23.png

McMerchie_devidson_helgaker_v1.35.zip

56.33 KB, 下载次数 Times of downloads: 58

basis6311_info.m

2.55 KB, 下载次数 Times of downloads: 35

myCO6311.m

13.49 KB, 下载次数 Times of downloads: 35

评分 Rate

参与人数
Participants 1
eV +6 收起 理由
Reason
sobereva + 6

查看全部评分 View all ratings

39

帖子

0

威望

132

eV
积分
172

Level 3 能力者

12#
发表于 Post on 2014-11-27 16:15:46 | 只看该作者 Only view this author
Shannon 发表于 2014-11-25 23:40
我的matlab没买并行工具包……学校没钱,哈哈。Matlab工具包里面有一堆好东西,可惜卖的实在是贵了点,一 ...

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

6万

帖子

99

威望

5万

eV
积分
120179

管理员

公社社长

13#
发表于 Post on 2014-11-27 17:52:40 | 只看该作者 Only view this author
并不并行现阶段来说无所谓。
真要并行的话就说明已经进入实用阶段了,那就肯定得写成C/Fortran了,到时候直接用openmp并行就行了,相当方便。
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办极高质量的各种计算化学类培训:初级量子化学培训班中级量子化学培训班高级量子化学培训班量子化学波函数分析与Multiwfn程序培训班分子动力学与GROMACS培训班CP2K第一性原理计算培训班,内容介绍以及往届资料购买请点击相应链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的高速路!培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入北京科音微信公众号获取北京科音培训的最新消息,并避免错过网上有价值的计算化学文章!
欢迎加入人气极高、专业性特别强的理论与计算化学综合交流群思想家公社QQ群(群号见此链接),合计达一万多人。北京科音培训班的学员在群中可申请VIP头衔,提问将得到群主Sobereva的最优先解答。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(十分强大、极为流行的量子化学波函数分析程序)
Google Scholar:https://scholar.google.com/citations?user=tiKE0qkAAAAJ
ResearchGate:https://www.researchgate.net/profile/Tian_Lu

127

帖子

1

威望

1231

eV
积分
1378

Level 4 (黑子)

14#
 楼主 Author| 发表于 Post on 2014-11-28 02:31:48 | 只看该作者 Only view this author
sobereva 发表于 2014-11-27 17:52
并不并行现阶段来说无所谓。
真要并行的话就说明已经进入实用阶段了,那就肯定得写成C/Fortran了,到时候 ...

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

21

帖子

0

威望

4

eV
积分
25

Level 2 能力者

15#
发表于 Post on 2014-12-14 02:17:54 | 只看该作者 Only view this author
人才啊,可以搞中国的Gaussian了

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2025-8-17 01:45 , Processed in 0.194914 second(s), 25 queries , Gzip On.

快速回复 返回顶部 返回列表 Return to list