c
c Index
c*****************************************************************
c********************** gbind.f *************************
c*****************************************************************
c*** ***
c*
c* this is a program for evaluation of g-functions (1983,1995)
c* for 1s, 2p, 2s initial and final 1s 2p 2s
c*
c*
c*****************************************************************
external f1s,ff1s,f2s,ff2s,f2p,ff2p
dimension vec(820),ved(820),bres(410)
dimension rgf(500)
integer chosen,yes,no
yes=1
no=0
write(*,*) 'Calculates Radial Matrix Elements'
write(*,*) 'of Coulomb interaction'
write(*,*) '.'
write(*,*) '.'
write(*,*) 'Mesh data: number, increment, odd number'
write(*,*)'m,h:'
read(*,*) m,h
write(*,*) 'ni li nf lf L'
read(*,*) ni,li,nf,lf,L
if(ni.GE.nf) then
n1=ni
n2=nf
l1=li
l2=lf
else
n2=ni
n1=nf
l2=li
l1=lf
endif
if(ni.EQ.nf) then
if(li.GT.lf) then
l1=li
l2=lf
else
l2=li
l1=lf
endif
endif
c
c writing out all possibilities
c
chosen = no
c
c including selection rules
c
c 1s 1s
if(n1.eq.1.and.l1.eq.0.and.n2.eq.1.and.l2.eq.0) then
if(L.eq.0) then
call subtt(l1,l2,L,ff1s,f1s,vec,ved,m,h)
chosen=yes
open(20,file='b1s1s')
endif
endif
c 2s 1s
if(n1.eq.2.and.l1.eq.0.and.n2.eq.1.and.l2.eq.0) then
if(L.eq.0) then
call subtt(l1,l2,L,ff2s,f1s,vec,ved,m,h)
chosen=yes
open(20,file='b2s1s')
endif
endif
c 2p 1s
if(n1.eq.2.and.l1.eq.1.and.n2.eq.1.and.l2.eq.0) then
if(L.eq.1) then
call subtt(l1,l2,L,ff2p,f1s,vec,ved,m,h)
chosen=yes
open(20,file='b2p1s')
endif
endif
c 2p 2s
if(n1.eq.2.and.l1.eq.1.and.n2.eq.2.and.l2.eq.0) then
if(L.eq.1) then
call subtt(l1,l2,L,ff2p,f2s,vec,ved,m,h)
chosen=yes
open(20,file='b2p2s')
endif
endif
c 2p 2p
if(n1.eq.2.and.l1.eq.1.and.n2.eq.2.and.l2.eq.1) then
if(L.eq.0.or.L.eq.2) then
call subtt(l1,l2,L,ff2p,f2p,vec,ved,m,h)
chosen=yes
if(L.eq.0) then
open(20,file='b2p2p0')
else
open(20,file='b2p2p2')
endif
endif
endif
c 2s 2s
if(n1.eq.2.and.l1.eq.0.and.n2.eq.2.and.l2.eq.0) then
if(L.eq.0) then
call subtt(l1,l2,L,ff2s,f2s,vec,ved,m,h)
chosen=yes
open(20,file='b2s2s')
endif
endif
c
if(chosen.eq.yes) then
write(20,303)
write(20,304) n1,l1,n2,l2,L
hst=2. *h
mst=m/2+1
n=m/2
ai1=.0
ai2=sumd(ved,1,m,h)
ai1=.0
bres(1) = ai2
k=0
m1=0
m2=0
do 10 k=1,n
m1=1+2*(k-1)
m2=2*k+1
ai1=ai1+sumd(vec,m1,m2,h)
ai2=ai2-sumd(ved,m1,m2,h)
r=(m2-1)*h
rgf(k)=r
bres(k+1)=ai1/(r**(l+1)) + ai2*r**l
if(l.gt.0) bres(1)=.0
write(20,*)rgf(k),bres(k+1)
10 continue
else
write(*,333) ni,li,nf,lf,L
333 format('UNSUPPORTED COMBINATION:',5I5)
endif
303 format('# g-functions bound states')
304 format('# ni=',i2,' li=',i2,' nf=',i2,' lf=',i2,' L=',I2)
stop
end
c**************************************************************
c
c constructing the overlaps for all points
c
c**************************************************************
c
subroutine subtt(li,lf,l,func,func2,vec,ved,m,h)
external func
external func2
dimension vec(820),ved(820)
data pi/3.141592653 /
c
c constructing g(l,r) calculates overlap function
c
open(21,file='test')
open(22,file='test2')
do 1 k=1,m
x=h*(k-1)
if(k.eq.1) x=1.d-002*h
f2=x*func2(x)
f=x*func(x)
write(21,*) x,f
write(22,*) x,f2
z=f2*f
vec(k)=x**l*z
ved(k)=z/x**(l+1)
1 continue
return
end
c**********************************************************
c
c simpson rule, from m1-th to m2-th point
c
c**********************************************************
c
function sumd(vec,m1,m2,xh)
c
dimension vec(820)
sumd=vec(m1)+vec(m2)
m3=m1+1
do 200 i=m3,m2,2
200 sumd=sumd+4.*vec(i)
if(m2-m1.le.2) go to 400
m3=m3+1
m4=m2-2
do 300 i=m3,m4,2
300 sumd=sumd+2.*vec(i)
400 sumd=sumd*xh/3.
return
end
c
c ***********************************************************
c * *
c * computes hydrogenic wavefunctions for the *
c * 1s, 2s, 2p case *
c * *
c * doubled, because of possible call to same function *
c * *
c ***********************************************************
c
function f1s(x)
f1s=2.*exp(-x)
return
end
c
c
function ff1s(x)
ff1s=2.*exp(-x)
return
end
c
function f2s(x)
f2s=(1. /sqrt(2. ))*(1. -0.5 *x)*exp(-0.5 *x)
return
end
c
function ff2s(x)
ff2s=(1. /sqrt(2. ))*(1. -0.5 *x)*exp(-0.5 *x)
return
end
c
function f2p(x)
f2p=(0.5 /sqrt(6. ))*x*exp(-0.5 *x)
return
end
c
function ff2p(x)
ff2p=(0.5 /sqrt(6.0 ))*x*exp(-0.50 *x)
return
end