implicit real*8(a-h,o-z) dimension rr(3) c=137.0359895d0 c2=c*c c4=c2*c2 write(6,12) write(6,11) 12 format(1x,'Slater-type',5x,'kinetically-balanced equation',15x, & 'Schroedinger equation') 11 format(4x,'z',2(9x,'exp',13x,'enr',9x,'S',4x), & 2x,'(s2/s1-1)*10^6') do 1 iz=1,137 esom=dsqrt(c4-iz**2*c2)-c2 f=0.0d0 q=1.0d0 a=a1opt(f,q,iz,eop) call overlap(s1,s2,iz,eb) exp=iz/(1.0d0-0.25d0*iz**2/c2) write(6,10)iz,a,eop,s2,exp,eb,s1,(s1/s2-1.0d0)*1.0d6 1 continue 13 format(1x,'Gauss-type, one function') 10 format(1x,i4,2(2f16.5,f10.6),f16.6) write(6,13) do 2 iz=1,137 call goverlap(s1,s2,en1,en2,exp1,exp2,iz) write(6,10)iz,exp2,en2,s2,exp1,en1,s1 2 continue stop end double precision function gauss(r,g,a,b) implicit real*8(a-h,o-z) e=-a*r-b*r**2 gauss=r**g*dexp(e) return end double precision function gint(r0,g,a,b) implicit real*8(a-h,o-z) r=r0 2 c=gauss(r,g,a,b) if(c.lt.1.0d-12)go to 1 r=r+1.0d-1 go to 2 1 h=r*5.0d-4 s1=0.0d0 s2=0.0d0 do 3 i=1,2000,2 r1=i*h r2=r1+h s1=s1+gauss(r1,g,a,b) s2=s2+gauss(r2,g,a,b) 3 continue gint=h*(4.0d0*s1+2.0d0*s2)/3.0d0 return end double precision function gn(n,q,f) implicit real*8(a-h,o-z) gn=1.0d0+2.0d0*f/(0.5d0+0.5d0*q)**n+f**2/q**n return end double precision function a1opt(f,q,iz,eop) implicit real*8(a-h,o-z) pi=4.0d0*datan(1.0d0) c=137.0359895d0 c2=c*c c f - expansion coefficent associated with G(a2). c It is assumed that the coefficient of G(a1) is 1.0d0. c q = a2/a1, where a2, a1 are the exponents in G(a2), G(a1) c a1opt - the exponent in G(a1) corresponding to the c extremum of the ground state energy c iz - the nuclear charge c eop - the ground-state energy ze=iz/c rp=0.5d0*(gn(2,q,f)/gn(3,q,f)+gn(0,q,f)/gn(1,q,f)) rm=0.5d0*(gn(2,q,f)/gn(3,q,f)-gn(0,q,f)/gn(1,q,f)) r0=gn(1,q,f)/gn(3,q,f) aa=ze**4*rm**2*(rp**2-rm**2)+ & ze**2*r0*(rp**2-2*rm**2)-r0**2 bb=ze*rm*(r0-ze**2*(rp**2-rm**2)) cc=ze**2*(rp**2-rm**2) delta=bb**2-aa*cc aop=(-bb-dsqrt(delta))/aa re=(1.0d0-ze*aop*rm)**2+aop**2*r0 eop=-ze*aop*rp-1.0d0+dsqrt(re) eop=c2*eop a1opt=c*aop return end subroutine overlap(s1,s2,i,eb) implicit real*8(a-h,o-z) dimension xnorm(3),xl(3),xs(3),xg(3),xb(3),ss(3,3) c=137.0359895d0 z=i*1.0d0 zz=z*z zc=z/c zc2=zc**2 c2=c*c esom=-c2*(1.0d0-dsqrt(1.0d0-zc2)) eb=-0.50d0*zz/(1.0d0-0.25d0*zc2) esh=-0.5d0*zz em=1.0d0/(1.0d0-0.25d0*zc2) gam=dsqrt(1.0d0-zc2) xnorm(1)=(z+z)**(gam+0.5d0)/ & dsqrt(2*s14aaf(2*gam+1.0d0,ifa)) xnorm(2)=2*(z*em)**1.5d0/dsqrt(1.0d0+0.25d0*zc2) xnorm(3)=(z/gam)**1.5d0*dsqrt(2*(1.0d0+gam)) xl(1)=dsqrt(1.0d0+gam) xl(2)=1.0d0 xl(3)=1.0d0 xs(1)=-dsqrt(1.0d0-gam) xs(2)=-0.5d0*zc xs(3)=-zc/(1.0d0+gam) xg(1)=gam xg(2)=1.0d0 xg(3)=1.0d0 xb(1)=z xb(2)=z*em xb(3)=z/gam do 2 k=1,3 do 2 l=1,3 ss(k,l)=xnorm(k)*xnorm(l)*(xl(k)*xl(l)+xs(k)*xs(l)) gama=s14aaf(xg(k)+xg(l)+1.0d0,ifa) beba=(xb(k)+xb(l))**(xg(k)+xg(l)+1.0d0) ss(k,l)=ss(k,l)*gama/beba 2 continue s1=ss(1,2) s2=ss(1,3) return end subroutine goverlap(s1,s2,en1,en2,exp1,exp2,iz) implicit real*8(a-h,o-z) dimension xnorm(3),xl(3),xs(3),xg1(3),xg2(3),xa(3), & xb(3),ss(3,3) c=137.0359895d0 z=iz*1.0d0 pi=4.0d0*datan(1.0d0) zz=z*z zc=z/c zc2=zc**2 c2=c*c esom=-c2*(1.0d0-dsqrt(1.0d0-zc2)) c em=1.0d0/(1.0d0-2.0d0*zc2/(9.0d0*pi)) gam=dsqrt(1.0d0-zc2) en1=-4.0d0*em*z**2/(3.0d0*pi) bb1=1.0d0/(2.0d0*em*c) xnorm(1)=(z+z)**(gam+0.5d0)/ & dsqrt(2*s14aaf(2*gam+1.0d0,ifa)) x=pi*(3.0d0*pi+2.0d0*zc2) xnorm(2)=16*(z*em)**1.5d0/(3.0d0*dsqrt(x)) exp1=8*em**2*z**2/(9*pi) tt=dsqrt(8.0d0/(9.0d0*pi))*zc c x=5.0d0/dsqrt(1.0d0-2*tt**2) aa=2.0d0*tt*(1.0d0+x)/(12.0d0+tt**2) exp2=c2*aa**2 x=(0.5d0*tt*aa-1.0d0)**2+3*aa**2 en2=(dsqrt(x)-1.0d0-2.5d0*tt*aa) en2=c2*en2 c bb2=en2/(3*exp2*c)+tt/dsqrt(exp2) x=8*exp2**1.50d0*dsqrt(2/pi) xnorm(3)=dsqrt(x/(1.0d0+3.0d0*bb2**2*exp2)) xl(1)=dsqrt(1.0d0+gam) xl(2)=1.0d0 xl(3)=1.0d0 xs(1)=-dsqrt(1.0d0-gam) xs(2)=-2.0d0*exp1*bb1 xs(3)=-2.0d0*exp2*bb2 xg1(1)=gam xg1(2)=1.0d0 xg1(3)=1.0d0 xg2(1)=gam xg2(2)=2.0d0 xg2(3)=2.0d0 xa(1)=z xa(2)=0.0d0 xa(3)=0.0d0 xb(1)=0.0d0 xb(2)=exp1 xb(3)=exp2 do 2 k=1,3 do 2 l=1,3 sn=xnorm(k)*xnorm(l) sg1=xg1(k)+xg1(l) sg2=xg2(k)+xg2(l) sa=xa(k)+xa(l) sb=xb(k)+xb(l) if(k.eq.1.and.l.eq.1)go to 5 r01=(-sa+dsqrt(sa**2+8*sb*sg1))/(4*sb) r02=(-sa+dsqrt(sa**2+8*sb*sg2))/(4*sb) go to 6 5 r01=sg1/sa r02=sg2/sa 6 ss(k,l)=xnorm(k)*xnorm(l)*(xl(k)*xl(l)*gint(r01,sg1,sa,sb) & + xs(k)*xs(l)*gint(r02,sg2,sa,sb)) 2 continue s1=ss(1,2) s2=ss(1,3) return end DOUBLE PRECISION FUNCTION S14AAF(X, IFAIL) C GAMMA FUNCTION DOUBLE PRECISION X INTEGER IFAIL DOUBLE PRECISION SRNAME DOUBLE PRECISION G, GBIG, T, XBIG, XMINV, XSMALL, Y INTEGER I, M INTEGER P01AAF DATA SRNAME /8H S14AAF / DATA XSMALL/1.0D-17/ DATA XBIG,GBIG,XMINV /57.0D0,7.1D+74,1.4D-76/ T = DABS(X) IF (T.GT.XBIG) GO TO 160 IF (T.LE.XSMALL) GO TO 140 M = X IF (X.LT.0.0D0) GO TO 80 T = X - DBLE(M) M = M - 1 G = 1.0D0 IF (M) 20, 120, 40 20 G = G/X GO TO 120 40 DO 60 I=1,M G = G*(X-DBLE(I)) 60 CONTINUE GO TO 120 80 T = X - DBLE(M-1) IF (T.EQ.1.0D0) GO TO 220 M = 1 - M G = X DO 100 I=1,M G = G*(X+DBLE(I)) 100 CONTINUE G = 1.0D0/G 120 T = 2.0D0*T - 1.0D0 Y= +8.86226925452758013D-1+T*( +1.61691987244425092D-2+ AT*( +1.03703363422075456D-1+T*( -1.34118505705967765D-2+ BT*( +9.04033494028101968D-3+T*( -2.42259538436268176D-3+ CT*( +9.15785997288933120D-4+T*( -2.96890121633200000D-4+ DT*( +1.00928148823365120D-4+T*( -3.36375833240268800D-5+ ET*( +1.12524642975590400D-5+T*( -3.75499034136576000D-6+ FT*( +1.25281466396672000D-6+T*( -4.17808776355840000D-7+ GT*( +1.39383522590720000D-7+T*( -4.64774927155200000D-8+ HT*( +1.53835215257600000D-8+T*( -5.11961333760000000D-9+ IT*( +1.82243164160000000D-9+T*( -6.13513953280000000D-10+ JT*( +1.27679856640000000D-10+T*( -4.01499750400000000D-11+ KT*( +4.26560716800000000D-11+ LT*( -1.46381209600000000D-11))))))))))))))))))))))) S14AAF = Y*G IFAIL = 0 GO TO 240 140 IF (T.LT.XMINV) GO TO 200 S14AAF = 1.0D0/X IFAIL = 0 GO TO 240 160 IF (X.LT.0.0D0) GO TO 180 IFAIL = P01AAF(IFAIL,1,SRNAME) S14AAF = GBIG GO TO 240 180 IFAIL = P01AAF(IFAIL,2,SRNAME) S14AAF = 0.0D0 GO TO 240 200 IFAIL = P01AAF(IFAIL,3,SRNAME) T = X IF (X.EQ.0.0D0) T = 1.0D0 S14AAF = DSIGN(1.0D0/XMINV,T) GO TO 240 220 IFAIL = P01AAF(IFAIL,4,SRNAME) S14AAF = GBIG 240 RETURN END SUBROUTINE X04AAF(I,NERR) INTEGER I, NERR INTEGER NERR1 DATA NERR1 /6/ IF (I.EQ.0) NERR = NERR1 IF (I.EQ.1) NERR1 = NERR RETURN END INTEGER FUNCTION P01AAF(IFAIL, ERROR, SRNAME) INTEGER ERROR, IFAIL, NOUT DOUBLE PRECISION SRNAME IF (ERROR.EQ.0) GO TO 20 CALL X04AAF (0,NOUT) IF (MOD(IFAIL,10).EQ.1) GO TO 10 WRITE (NOUT,99999) SRNAME, ERROR STOP 10 IF (MOD(IFAIL/10,10).EQ.0) GO TO 20 WRITE (NOUT,99999) SRNAME, ERROR 20 P01AAF = ERROR RETURN 99999 FORMAT (1H0, 38HERROR DETECTED BY NAG LIBRARY ROUTINE , A8, * 11H - IFAIL = , I5//) END c double precision function ga1opt(f,q,iz,eop) implicit real*8(a-h,o-z) dimension r1(4),r2(4),x(4) pi=4.0d0*datan(1.0d0) c=137.0359895d0 c2=c*c c f - expansion coefficent associated with G(a2). c It is assumed that the coefficient of G(a1) is 1.0d0. c q = a2/a1, where a2, a1 are the exponents in G(a2), G(a1) c a1opt - the exponent in G(a1) corresponding to the c extremum of the ground state energy c iz - the nuclear charge c eop - the ground-state energy r1(1)=8.0d0*q/(1.0d0+q)**2 r1(2)=4.0d0/(1.0d0+q) r1(3)=dsqrt(32.0d0)/(1.0d0+q)**1.5d0 r1(4)=dsqrt(128.0d0)*q/(1.0d0+q)**2.5d0 r2(1)=1.0d0 r2(2)=1.0d0/q r2(3)=1.0d0/q**1.5d0 r2(4)=1.0d0/dsqrt(q) do 1 i=1,4 1 x(i)=1.0d0+f*r1(i)+f*f*r2(i) w1=x(1)/x(4) w2=x(2)/x(3) w3=x(4)/x(3) t=dsqrt(8.0d0/(9.0d0*pi))*iz/c xm=w3 xl=t*(w1+1.5d0*w2) xk=t*(w1-1.5d0*w2) rt=1.0d0+(xk**2-xl**2)/(3.0d0*xm) ds=dsqrt(rt) a=(xl/ds-xk)/(3.0d0*xm+xk**2) a1opt=c2*a**2 e=(xk*xl+3.0d0*xm*ds)/(3.0d0*xm+xk**2) eop=(e-1.0d0)*c2 par1=xk*xl/(3.0d0*xm) par2=(xk*xk-xl*xl)/(3.0d0*xm) par3=xk*xk/(3.0d0*xm) c par4=1/t**2 write(6,100)q,w1,w2,w3 write(6,100)q,par1,par2,par3,eop stt=12.0d0/t**2 write(6,100)q,par1*stt,par2*stt,par3*stt c write(6,100)q,rt,a,e,dsqrt(1.0d0-iz*iz/c2) 100 format(4x,f7.2,4f14.8) return end