| 查看: 1141 | 回复: 2 | |||
| 【奖励】 本帖被评价2次，作者fangyongxinxi增加金币 1.6 个 | |||
[资源]
一个HF方法计算H2键长的matlab代码（一百多行）
|
|||
|
% Downloadable MATLAB code, copyright Amritanshu Palaria % H2 ground state from a non-relativistic Hartree Fock treatment using STO-4G basis set and Born Oppenheimer approximation % Program written by Amritanshu Palaria @ NCN, Apr 12, 2006 % For theoretical reference, please see "Computational Physics" by Thijssen % erf(x) = 2/sqrt(pi) * integral from 0 to x of exp(-t^2) dt clear all; % initialization Z=1,S=1,A=0,F=0,V=0,D=0,g=0,Vecp=0,Vec=0,Val=0,EVal=0; for y=1:41 clear T,S,A,F,V,D,g,Vecp,Vec,Val,EVal; % y = 11; disp('Step '); disp(y); r = 0.45+0.05*y; % distance between nucleii in a.u.(units of a0, ie. Bohr radius)... since the problem is essentially 1-D, this is all that is needed disp('H-H bondlength'); disp(r); NNr=1/r; % alphas for the 4 Gaussian basis functions for each H atom alpha1 = [13.00773 1.962079 0.444529 0.1219492]; % sto-4g, Phi = 0.03492*exp(-13.00773*r^2)+0.2347*exp(-1.962079*r^2)+0.8138*exp(-0.444529*r^2)+1*exp(-0.1219492*r^2) alpha = [alpha1 alpha1]; asize=size(alpha); % one H nucleus at origin, the other at r ra=[0 0 0 0 r r r r]; % the location of the atom with the ith basisfunction % first 4 basis fns come from the atom at origin % the next 4 come from the one at r % the S (overlap), T (K.E.) and A (Coulomb) matrices for i=1:8 % run over all basis functions for j=1:i % do not run over all basis functions to exploit symmetry % intermediate variables as = alpha(i)+alpha(j); ap = alpha(i)*alpha(j); rat=ap/as; % location of the gaussian resulting from product of % the 2 gaussian basis functions rp=(alpha(i)*ra(i)+alpha(j)*ra(j))/as; S(i,j)=(pi/as)^1.5*exp(-rat*(ra(i)-ra(j))^2); S(j,i)=S(i,j); % using symmetry T(i,j)=0.5*rat*(6-4*rat*(ra(i)-ra(j))^2)*S(i,j); T(j,i)=T(i,j); % using symmetry if (rp==0) F0=1+sqrt(pi)/2*erf(sqrt(as*(rp-r)^2))/sqrt(as*(rp-r)^2); % erf(x) = 2/sqrt(pi) * integral from 0 to x of exp(-t^2) dt elseif (rp==r) F0=sqrt(pi)/2*erf(sqrt(as*(rp-0)^2))/sqrt(as*(rp-0)^2)+1; else F0=sqrt(pi)/2*(erf(sqrt(as*(rp-0)^2))/sqrt(as*(rp-0)^2)+erf(sqrt(as*(rp-r)^2))/sqrt(as*(rp-r)^2)); end A(i,j)=-2*pi*Z/as*exp(-rat*(ra(i)-ra(j))^2)*F0; A(j,i)=A(i,j); % using symmetry end end S T A [V,D]=eig(S); % gives eigenvectors in columns of V, this diagnolizes S for i=1:asize(2) V(:,i)=V(:,i)/(D(i,i)^0.5); end % the two-electron integrals P taking care of all symmetries % it is noted that the gaussian functions evaluate to REAL values % and the g operator is real and symmetric w.r.t. 1 and 2 for i=1:8 for j=1:i % do not run over all basis functions to exploit symmetry as1=alpha(i)+alpha(j); ap1=alpha(i)*alpha(j); rp=(alpha(i)*ra(i)+alpha(j)*ra(j))/as1; for k=1:i-1 % do not run over all basis functions to exploit symmetry for l=1:k % do not run over all basis functions to exploit symmetry as2=alpha(k)+alpha(l); ap2=alpha(k)*alpha(l); rq=(alpha(k)*ra(k)+alpha(l)*ra(l))/as2; if ((rp-rq)==0) F0=1; else F0=sqrt(pi)/2*erf(sqrt(as1*as2*(rp-rq)^2/(as1+as2)))/sqrt(as1*as2*(rp-rq)^2/(as1+as2)); end g(i,j,k,l)=2*pi^2.5/as1/as2/(as1+as2)^0.5*exp(-ap1*(ra(i)-ra(j))^2/as1-ap2*(ra(k)-ra(l))^2/as2)*F0; g(k,l,i,j)=g(i,j,k,l); % using symmetry g(j,i,k,l)=g(i,j,k,l); % using symmetry g(i,j,l,k)=g(i,j,k,l); % using symmetry g(j,i,l,k)=g(i,j,k,l); % using symmetry g(k,l,j,i)=g(i,j,k,l); % using symmetry g(l,k,i,j)=g(i,j,k,l); % using symmetry g(l,k,j,i)=g(i,j,k,l); % using symmetry end k=i; for l=1:j as2=alpha(k)+alpha(l); ap2=alpha(k)*alpha(l); rq=(alpha(k)*ra(k)+alpha(l)*ra(l))/as2; if ((rp-rq)==0) F0=1; else F0=sqrt(pi)/2*erf(sqrt(as1*as2*(rp-rq)^2/(as1+as2)))/sqrt(as1*as2*(rp-rq)^2/(as1+as2)); end g(i,j,k,l)=2*pi^2.5/as1/as2/(as1+as2)^0.5*exp(-ap1*(ra(i)-ra(j))^2/as1-ap2*(ra(k)-ra(l))^2/as2)*F0; g(k,l,i,j)=g(i,j,k,l); % using symmetry g(j,i,k,l)=g(i,j,k,l); % using symmetry g(i,j,l,k)=g(i,j,k,l); % using symmetry g(j,i,l,k)=g(i,j,k,l); % using symmetry g(k,l,j,i)=g(i,j,k,l); % using symmetry g(l,k,i,j)=g(i,j,k,l); % using symmetry g(l,k,j,i)=g(i,j,k,l); % using symmetry end end end end % Self-consistent loop C=ones(1,asize(2)); %initial guess, 11 行几列的1矩阵 nor=C*S*C' % c' 转置 disp('Initial Guess'); C=C/nor Elast=1; EVal=0; P=C'*C; count=0; rplot=-2:0.01:3; % range in atomic units (bohr) for plotting the probability density % figure(); % Generating Fock matrix while((abs(EVal(1)-Elast))>0.000001) count=count+1; disp('Step #'); disp(count); Elast=EVal(1) for i=1:8 for j=1:8 J=0; for k=1:8 for l=1:8 J=J+P(k,l)*g(i,j,k,l); end end F(i,j)=T(i,j)+A(i,j)+J; end end % solve the generalized eigenvalue problem FC=ESC % inverse of V is its conjugate transpose because V is unitary Fp=V'*F*V; % modified Fock matrix [Vecp,Val]=eig(Fp); Vec=V*Vecp; EigVal=diag(Val,0) [EVal,index]=sort(EigVal); disp('Ground State eigenvalue from this step is (in a.u., i.e. Hartree):'); Erec(count)=EVal(1) % lowest eigenvalue GrCoeff=Vec(:,index(1)); % eigenvector corresponding to ground state disp('C matrix from this step is:'); C=GrCoeff.' %new C matrix, this is normalized w.r.t. S ie. C*S*C'=1 disp('Input density matrix for next step is:'); P=0.8*P+0.2*C'*C %new input density matrix GndWvFn=C(1:4)*exp(-(alpha(1:4)).'*(rplot.*rplot))+C(5:8)*exp(-(alpha(5:8)).'*((rplot-r).*(rplot-r))); % ground state wavefunction from STO-4G approximation % Uncomment the next 4 lines if you need a plot of probability with step of convergence % hold on; % plot(rplot,(abs(GndWvFn)).^2,'r--','LineWidth',1.5); % xlabel('r position'); % ylabel('Probability density'); % title('Plot of the probability density vs. r position with step of convergence'); end Qs=0; for i=1:8 for j=1:8 for k=1:8 for l=1:8 Qs=Qs+g(i,j,k,l)*C(i)*C(j)*C(k)*C(l); end end end end disp('The final total energy (including nuclear-nuclear repulsion) after convergence for this H-H bond is '); Eg(y)=2*C*(T+A)*C'+Qs+NNr; % the energy from h for each e, the e-e interaction and the nucleus-nucleus interaction added disp(Eg(y)); if Eg(y)==min(Eg) rmin = r; Erecmin = Erec; WvFnmin = GndWvFn; end end figure(); plot([0.5:0.05:2.5],Eg,'bo','LineWidth',2); title('Plot showing total energy from Hartree Fock treatment of H-H using Born Oppenheimer approximation'); xlabel('H-H bondlength'); ylabel('Total energy (in a.u.)'); disp('The least energy bond length (A) for H-H is '); disp(rmin*0.529177); disp('The energy (Hartree) of this bond is '); disp(min(Eg)); figure(); plot(Erecmin,'ro','LineWidth',2); axis([1 count min(Erecmin) max(Erecmin)+0.1]); title('Plot showing convergence of Fock level with every step for the least energy H-H bondlength '); xlabel('Step of convergence'); ylabel('Fock energy level (in a.u.)'); %plot of final probabiliy density figure(); plot(rplot,(abs(WvFnmin)).^2,'LineWidth',2.5); xlabel('r position'); ylabel('Probability density'); title('Plot of the probability density of electrons vs. r position for the least energy H-H bond'); axis tight; |
回复此楼» 收录本帖的淘帖专辑推荐
 量化 |
» 猜你喜欢
ELISA试验中不容忽视的细节盘点
已经有0人回复
-
ELISA试验中不容忽视的细节盘点(二)
已经有0人回复
-
物理化学论文润色/翻译怎么收费?
已经有163人回复
-
求助火焰封管的时候管子炸了
已经有1人回复
-
细胞培养，这22个细节一定要注意！（一）
已经有0人回复
-
CSC访学博后项目获批，外方学校暂停合作该怎么办？
已经有51人回复
-
请问四氢呋喃溶解的聚合物用甲醇沉淀时，如何使沉淀过程加速？
已经有2人回复
-
七嗪类物质合成求助
已经有0人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 求回归方程参数估计的最优化算法matlab代码
已经有9人回复
- 求助大侠帮忙，有关matlab模拟数据
已经有29人回复
- matlab 计算过程中工作空间变量保存及读取 求解决内存不足的方法
已经有5人回复
- 如何用MATLAB 实现化学反应方程式（写程序代码）?
已经有20人回复
- 计算晶体怎样减少键长键角的相对误差？怎样确定中心原子周围键的键组成？
已经有4人回复
- 【求助】我想计算一个分子中C原子的p轨道键长，用什么软件计算比较好。
已经有6人回复
- 【求助】【如何计算H-H的键能与键长之间的关系】
已经有6人回复
- 【讨论】大家讨论下计算弱键的方法-氢键，卤键
已经有9人回复
- 【求助】Weickert的各向异性扩散方程的滤波方法的matlab程序代码
已经有8人回复
2楼2011-12-05 15:48:44
★★★★★ 五星级,优秀推荐
|
计算结果（matlab 2011ra X64): The final total energy (including nuclear-nuclear repulsion) after convergence for this H-H bond is -1.0327 The least energy bond length (A) for H-H is 0.7408 The energy (Hartree) of this bond is -1.1265 [ Last edited by dhqdqk on 2011-12-6 at 12:54 ] |
3楼2011-12-06 11:53:05