Квантовая теория молекул (Часть 1). Шунина В.А - 27 стр.

                                       27
14.   (* одноцентровые параметры связывания *)
15.   Be={8.2,12.8,5.4,5.4};
16.   (* одноцентровые кулоновские интегралы *)
17.   Ga={10.017,13.663,12.848,12.848};
18.   (* координаты атомов в ангстремах *)
19.   R={{0.000,0.000,0.000},{0.000,1.146,0.000},{-0.933,-
      0.577,0.000},{+0.933,-0.577,0.000}};

20. DatIn={na,R,no,nocc,nbnlms,zt,Zcore,IP,Be,Ga};

21. CNDO2[DatIn_List]:=Module[{flag=0,na,R,no,nocc,nbnlms,zt,Z,U,Be,
22. Ga,β ,γ,fz,p,S,f,h,e,c,ct,Res,Etot,Etot0, δ =0.00001,Pat,summa,it,
23.  itmax=100,a,b,µ,ν,z,W,ce=27.212},
24.  {na,R,no,nocc,nbnlms,zt,Zcore,IP,Be,Ga}=DatIn;
25.  R=R/0.529167;IP=IP/ce;Be=Be/ce;Ga=Ga/ce;
26.  f=h=c=ct=p=Table[0.,{no},{no}];e=Table[0.,{no}];
27. γ=Table[0.,{na},{na}];Q=W=Pat=Table[0.,{na}];
28. β=Table[(Be[[a]]+Be[[b]])/2.,{a,na},{b,na}];
29.  Enuc=0.;
30.  For[a=1,a≤ na,a++,
31.   For[b=1,b<=a,b++,
32.    If[a b,γ[[a,a]]=Ga[[a]];Continue[]];
33.    Rab=R[[a]]-R[[b]];Rab=Sqrt[Rab.Rab];
34.    Enuc=Enuc+Zcore[[a]]*Zcore[[b]]/Rab;
35.    z=0.5*Rab*(Ga[[a]]+Ga[[b]]);fz=1./Sqrt[1.+z^2];
36.    γ[[a,b]]=0.5*(Ga[[a]]+Ga[[b]])*fz;γ[[b,a]]=γ[[a,b]]]];
37.  For[a=1,a≤ na,a++,
38.   W[[a]]=Sum[(1-KroneckerDelta[a,b])*Zcore[[b]]*γ[[a,b]],
                                                   {b,1,na}]];
39.   For[µ=1,µ≤ no,µ++,
40.    a=IntegerPart[nbnlms[[µ]]/10000];
41.    h[[µ,µ]]=IP[[µ]]-W[[a]]-(Zcore[[a]]-0.5)*γ[[a,a]];
42.    For[ν =1,ν ≤µ,ν++,
43.     b=IntegerPart[nbnlms[[ ν]]/10000];
44.     If[a b && µ≠ ν,h[[µ,ν]]=h[[ν,µ]]=0.;Continue[]];
45.     If[a!=b,
46.      S=Over[nbnlms[[µ]],zt[[µ]],R[[a]],
                 nbnlms[[ν ]],zt[[ν]],R[[b]]];
47.      If[S>10.,Return[{2,Etot,e,c,Q}]];
48.      h[[µ,ν]]=-β[[a,b]]*S;h[[ ν,µ]]=h[[µ,ν]]]]];
49.   For[it=0,it≤itmax,it++,
50.    If[it 0,
51.     Do[f[[µ,ν]]=If[µ ν,IP[[µ]],h[[µ,ν]]],{µ,no},{ν,no}],
52.     Etot0=Etot;
53.     For[b=1,b≤na,b++,
54.      summa=0.;
55.      For[µ=1,µ≤no,µ++,
56.       a=IntegerPart[nbnlms[[ µ]]/10000];
57.       If[a b,summa=summa+p[[ µ,µ]]]];
58.     Pat[[b]]=summa];
59.     For[µ=1,µ≤no,µ++,
60.      a=IntegerPart[nbnlms[[ µ]]/10000];