100 ! **********************************************************************
110 ! Program          ONELECT2.800     Utg}va  1.0      1984-04-23
120 ! av Torbj|rn Alm
130 ! Ins{nd av Torbj|rn Alm
140 ! F|r ABC800M ABC800C ABC802 ABC806
150 ! Testad p} ABC806
160 ! 
170 ! Screened by --"-TCP-"--
180 ! **********************************************************************
190 ! ------------------------------------------------------------------------
200 ! O N E L E C T R O N 2 - Compute one electron integrals for STO
210 ! After QCPE 117, by Peter Offenhartz
220 ! ABC-800 Basic II version 830102 Torbj|rn Alm, ABC-116
230 ! Input data:
240 ! NA     Primary Quantum number for Atom A
250 ! NB     Primary Quantum number for Atom B
260 ! Refers to Main Electron Shells, K=1, L=2, M=3 etc.
270 ! LA     Secondary Quantum number for Atom A
280 ! LB     Secondary Quantum number for Atom B
290 ! M      Magnetic Quantum number for Atoms A and B
300 ! M=0    Sigma, M=1  Pi,  M=2  Delta.
310 ! The program only accepts M=0-2.
320 ! The program works up to NA+NB=16, higher than any known atom
330 ! in the periodic system so far.
340 ! LA<NA, LB<NB, according to quantum chemical rules
350 ! X1     Exponent on orbital of Atom A
360 ! X1     Exponent on orbital of Atom B
370 ! R      Distace in Atomic Units between orbitals A and B
380 ! If R  <0, abs(R) is in ]ngstr|m
390 ! ------------------------------------------------------------------------
400 INTEGER : EXTEND 
410 Sqt.=SQR(3.) : Sq5.=SQR(5.)
420 DIM Raa.(20),Rbb.(20),Fack.(20)
430 Fack.(0)=1. : FOR I=1 TO 20 : Fack.(I)=Fack.(I-1)*I : NEXT I 
440 ; "One Electron 2" : ; " na,nb,la,lb,m,x1,x2,r"
450 INPUT Na,Nb,La,Lb,M,Ex1.,Ex2.,R. : IF R.<0. R.=-R./.519728
460 Alph.=.5*R.*(Ex1.+Ex2.) : Beta.=.5*R.*(Ex1.-Ex2.)
470 An.=FNFab.(Ex1.,Na) : Bn.=FNFab.(Ex2.,Nb)
480 Xami.=EXP(-Alph.) : Xbmi.=EXP(-Beta.)
490 Raa.(0)=Xami./Alph. : Xm.=Na+Nb+4
500 FOR I=1 TO Na+Nb+4 : Raa.(I)=(Xami.+I*Raa.(I-1))/Alph. : NEXT I 
510 IF ABS(Beta.)<=.368*Xm.+.184*LOG(Xm.)+.821 THEN 570
520 ! upward stable recurrence
530 Rbb.(0)=FNBfx.(0,Beta.,Xbmi.) : Xm.=1./Xbmi. : FOR I=1 TO Na+Nb+4
540   Xm.=-Xm. : Rbb.(I)=(Xm.-Xbmi.+I*Rbb.(I-1))/Beta.
550 NEXT I 
560 GOTO 620
570 Rbb.(Na+Nb+4)=FNBfx.(Na+Nb+4,Beta.,Xbmi.) : I=Na+Nb
580 IF I AND 1 THEN I=1 ELSE I=-1
590 Xm.=I/Xbmi. : FOR I=Na+Nb+4 TO 1 STEP -1
600   Xm.=-Xm. : Rbb.(I-1)=(Beta.*Rbb.(I)+Xbmi.-Xm.)/I
610 NEXT I 
620 ; '        Overlap=';FNOnel.(1)
630 ; ' Kinetic Energy=';FNOnel.(2)
640 ; '          1/r1 =';FNOnel.(4)
650 ; '          1/r2 =';FNOnel.(3)
660 ; '        Dipole =';FNOnel.(5)
670 ; '     Q-pole 22 =';FNOnel.(6)
680 ; '     Q-pole 02 =';FNOnel.(7)
690 ; '     Q-pole 20 =';FNOnel.(8)
700 GOTO 440
710 END 
720 DEF FNFab.(X.,N)=SQR((2.*X.)^(N*2+1)/Fack.(N*2)) 
730 DEF FNRmx.(N,M) LOCAL R.,S.,T.,I
740   IF M>N RETURN (0.)
750   IF M=N RETURN (1.)
760   R.=1. : I=1 : WHILE I<=N-M
770   R.=-(R.*(N-I+1))/I : I=I+1 : WEND 
780   RETURN (R.)
790 FNEND 
800 DEF FNRmp.(K) LOCAL R.,I
810   R.=0. : I=1 : WHILE I<=Jpx+1
820     R.=R.+FNRmx.(Jpx,Jpx-I+1)*ABS(FNRmx.(Ipx,K+I-1-Jpx))
830   I=I+1 : WEND 
840   RETURN (R.)
850 FNEND 
860 DEF FNRi.(M,N) LOCAL S.,K,L
870   S.=0. : K=M : L=M+Ipx+Jpx : WHILE K<=L
880     S.=S.+Raa.(K+Iaf)*Rbb.(L+N-K+Ibf)*FNRmp.(K-M)
890   K=K+1 : WEND : RETURN (S.)
900 FNEND 
910 DEF FNBfx.(N,Beta.,Xbmi.) LOCAL Ab.,Rn.,Qfn.,B.,F.,T1.,T2.,S1.,S2.,S3.,S4.
920   Ab.=ABS(Beta.) : Rn.=N : IF Ab.<.001 THEN IF N AND 1 THEN RETURN (0.) ELSE RETURN (1./(Rn.+.5))
930   IF N=0 THEN RETURN ((1./Xbmi.-Xbmi.)/Beta.)
940   B.=-Ab. : T1.=1. : T2.=1. : S1.=1. : S2.=1. : I=1 : WHILE I<=N
950     T1.=T1.*Ab./I : T2.=T2.*B./I : S1.=S1.+T1. : S2.=S2.+T2.
960   I=I+1 : WEND 
970   T1.=T1.*Ab./(N+1) : Qfn.=T1. : S3.=T1. : T2.=T2.*B./(N+1) : S4.=T2. : I=N+2
980   WHILE T1./S3.>1.E-08 : T1.=T1.*Ab./I : S3.=S3.+T1.
990   I=I+1 : WEND 
1000   I=N+2 : Rn.=N+1 : WHILE ABS(T2./S4.)>1.E-09
1010   T2.=T2.*B./I : S4.=S4.+T2. : I=I+1 : WEND 
1020   T2.=(S3.*S2.-S1.*S4.)/(Qfn.*Rn.) : IF Beta.<0. THEN RETURN (ABS(T2.)) ELSE RETURN (T2.)
1030 FNEND 
1040 DEF FNGx.(Jna,Jnb,Ncih) LOCAL K,G.
1050   K=3*La+Lb+2*M+1 : Ipx=Jna-La : Jpx=Jnb-Lb
1060   G.=.5*An.*Bn.*EXP(LOG(.5*R.)*(Jna+Jnb+Ncih+1))
1070   IF K=1 RETURN (G.*FNRi.(0,0))
1080   IF K=2 RETURN (G.*Sqt.*(FNRi.(0,0)-FNRi.(1,1)))
1090   IF K=3 RETURN (G.*(3.*(FNRi.(0,0)+FNRi.(2,2))-4.*FNRi.(1,1)-FNRi.(0,2)-FNRi.(2,0))*.5*Sq5.)
1100   IF K=4 OR K=10 OR K=12 RETURN (-1000.)
1110   IF K=5 RETURN (G.*3.*(FNRi.(0,0)-FNRi.(2,2)))
1120   IF K=6 RETURN (G.*.5*Sqt.*Sq5.*(3.*(FNRi.(0,0)+FNRi.(3,3))-FNRi.(1,1)-FNRi.(2,2)-FNRi.(2,0)-FNRi.(0,2)-FNRi.(3,1)-FNRi.(1,3)))
1130   IF K=7 RETURN (G.*1.5*(FNRi.(2,0)+FNRi.(0,2)-FNRi.(0,0)-FNRi.(2,2)))
1140   IF K=8 RETURN (G.*1.5*Sq5.*(FNRi.(3,3)+FNRi.(1,1)-FNRi.(3,1)-FNRi.(1,3)-FNRi.(2,2)+FNRi.(2,0)+FNRi.(0,2)-FNRi.(0,0)))
1150   IF K=9 RETURN (G.*1.25*(9.*(FNRi.(0,0)+FNRi.(4,4))-6.*(FNRi.(2,4)+FNRi.(4,2)+FNRi.(0,2)+FNRi.(2,0))+FNRi.(4,0)+FNRi.(0,4)+4.*FNRi.(2,2)))
1160   IF K=11 RETURN (G.*7.5*(FNRi.(4,4)-FNRi.(4,2)-FNRi.(2,4)+FNRi.(2,0)+FNRi.(0,2)-FNRi.(0,0)))
1170   IF K=13 RETURN (G.*1.875*(FNRi.(0,0)+FNRi.(4,4)+FNRi.(0,4)+FNRi.(4,0)+4.*FNRi.(2,2)-2.*(FNRi.(0,2)+FNRi.(2,0)+FNRi.(2,4)+FNRi.(4,2))))
1180 FNEND 
1190 DEF FNOnel.(Itax) LOCAL Kk,Nsw,Jna,Jnb,Cof.,Coff.,G.,Eex.
1200   Iaf=0 : Ibf=0
1210   IF Itax=1 RETURN (FNGx.(Na,Nb,0))
1220   IF Itax=3 RETURN (FNGx.(Na,Nb-1,0))
1230   IF Itax=4 RETURN (FNGx.(Na-1,Nb,0))
1240   IF Itax=5 Iaf=1 : Ibf=1 : RETURN (-FNGx.(Na,Nb,1))
1250   IF Itax=6 Iaf=2 : Ibf=2 : RETURN (FNGx.(Na,Nb,2))
1260   IF Itax=7 Iaf=2 : RETURN (FNGx.(Na,Nb,2))
1270   IF Itax=8 Ibf=2 : RETURN (FNGx.(Na,Nb,2))
1280   Kk=1-(Na>La+1)-2*(Nb>Lb+1) : Cof.=0. : Nsw.=0
1290   IF Kk=2 Nsw=1
1300   IF Kk=4 Cof.=-.5*(Na+La)*(Na-La-1)*FNGx.(Na-2,Nb,0)
1310   IF Nsw THEN G.=2.*Nb : Eex.=Ex2. : Jna=Na : Jnb=Nb-1 ELSE G.=2.*Na : Eex.=Ex1. : Jna=Na-1 : Jnb=Nb
1320   Coff.=.5*G.*FNGx.(Jna,Jnb,0)*Eex.
1330   RETURN (Coff.+Cof.-.5*Eex.*Eex.*FNGx.(Na,Nb,0))
1340 FNEND