ccc PROGRAM GIANT ccc Hartree-Fock and RPA with Skyrme interaction, ready for isoscalar ccc modes of general multipole; written by Nguyen Van Giai (1986) ccc and slightly modified by G. Colo` (1992). Units used are fm for ccc length and MeV for energy. implicit real*8 (a-h,o-z) parameter nnn=200,nmx=16,nmy=nmx+1,nsp=250,ncf=380,ncf2=2*ncf ccc nnn=dimension of vectors in radial space; greater or equal than nri ccc nmx=greater or equal than nosc ccc nsp,ncf=greater or equal than the number of single particle states ccc and than the number of particle-hole configurations, respectively common/betalf/calf(7),cbet(7) common/bwu2/ih(ncf),ip(ncf),ecf(ncf),nig(ncf),nqqx(10,2), 1nqqy(10,2),ngr,irpa,jtot,itot,ipi common/bnri/nri common/bwu1/nconf,ifail,ntot,sqbel(ncf),s0,icoulb,pcent0 common/bdiam/aa(ncf,ncf),bb(ncf,ncf),evr(ncf),vecr(ncf2,ncf) common/bwf/wf(nnn,nsp),dwf(nnn,nsp) common/bwg/wg(nnn,nmx),ewg(nmx),xor(nmx),dwg(nnn,nmx) common/bqwf/nn(nsp),ll(nsp),lj(nsp),lt(nsp),ehf(nsp),xnor(nsp) common/bfg/f0(nnn),g0(nnn) common/bpas/h,cof common/bpot/xmn(nnn),xmp(nnn),vn(nnn),vp(nnn),yn(nnn),yp(nnn), 1vsn(nnn),vsp(nnn),vc(nnn),wnd(nnn),wpd(nnn) common/blecp1/t0,t1,t2,t3,t13,alfe,x0,x1,x2,x3,anucl,znucl common/blecp2/dt(nnn),taut(nnn),dt1(nnn),dt2(nnn),dp(nnn),dn(nnn) common/pot/u(nnn),rmass,eta,xka,xmsm(nnn),sreg(nnn),defa,sig, 1sirg(nnn),sgne,dreg(nnn),dirg(nnn),c2d,s2d,sigma common/bjost/nr,refp,aifp,refd,aifd common/brint/ri0,ris,rx1,rx6(2),rx7(2),rx23(2,2),rx45(2,2) common/bhf/nos,npo,ni dimension nord(ncf2),vec(ncf2),tjl1(15),opr(nnn),rmat(ncf), 1bel(ncf2) namelist/param1/anucl,znucl,jtot,itot,ipi,nos,nossp,nossn,irpa,cof,xkf namelist/param2/t0,t1,t2,t3,x0,x1,x2,x3,alfe,nosc,nri,iliu namelist/param3/nnna,nmxa,nspa,ncfa,afa0 100 format(15i3) 101 format(6e12.6) 102 format(4i3,e12.6) 103 format(2x,'number of s.p. states i=',i3,'is larger than', 1' reserved dimension nsp=',i3) 104 format(2x,'x**2-y**2 is negative for state number',i3) 105 format(2x,'single particle states',//3x,'i',4x,'n',4x,'l',4x,'j', 14x,'lt',6x,'e',13x,'norm') 106 format(1x,i3,3i5,'/2',i3,3x,e12.5,2x,e12.5) 107 format(2x,'p-h configurations',//4x,'i',5x,'ih',4x,'ip',6x, 1'eph',8x,'b(el)') 108 format(2x,i4,2x,i4,2x,i4,3x,4e12.5) 109 format(2x,'moments m-1=',e12.5,3x,'m0=',e12.5,3x,'m1=', 1e12.5,3x,'m3=',e12.5) 110 format(//2x,'state',5x,'energy',8x,'b(el)',4x,'fract(newsr)', 12x,'fract(ewsr)') 111 format(3x,i3,3x,4e13.5) 112 format(10e12.5) 113 format(2x,'solution') 114 format(120('*')) 116 format(2x,'cof,afa0 =', e12.6,2x,e16.5) 117 format(2x,'nri,ipunch,iliu,icoulb,h =',4i5,2x,f5.2) 118 format(2x,'jtot,itot,ipi,irpa =',4i5) 119 format(2x,'nos,npo,nossp,nossn,nosc,ni =',6i5) 120 format(4i3,f5.2) ccc CONSTANTS pi=3.14159265D0 nnna=nnn nmxa=nmx nspa=nsp ncfa=ncf hbc=197.3D0 dmshb=0.04823D0 bmas=0.5D0*dmshb*hbc*hbc ccc INPUTS read *,cof,afa0 ccc cof=0. for woods-saxon,1. for hf potential ccc afa0 must be set as: hbar*omega*mc2/(hbc**2); a parametrization ccc for hbar*omega is 41*A**(-0.333333), mc2 is the nucleon ccc mass in MeV, hbc is hbar*c (197.3 MeV*fm) write(1,116),cof,afa0 read *,nri,ipunch,iliu,icoulb,h ccc nri is the number of points used for radial integrals and h=step ccc if ipunch=0 RPA data are not punched out, they are if ipunch=1 ccc iliu=0 keeps sigma.sigma dependence in t0,t3, and iliu=1 ccc drops sigma.sigma dependence in t0,t3 ccc icoulb=0 if no coulomb, =1 with coulomb write(1,117),nri,ipunch,iliu,icoulb,h if(nri.gt.nnn)then write(1,*) 'Parameter nnn must be increased' stop end if nr=nri/2 read *,nos,npo,nossp,nossn,nosc,ni ccc nos=number of occupied states, npo=number of proton states (proton ccc must be put before neutrons in the com file) ccc nossp,nossn=first holes considered when building the particle-hole ccc configurations, of proton and neutron respectively ccc nosc=number of harmonic oscillator shells (10-15 shells are sufficient ccc even for heavy nuclei, cfr. Blaizot-Gogny, Nucl.Phys.A284, p.429) ccc ni=number of iterations for the Hartree-Fock equations write(1,119),nos,npo,nossp,nossn,nosc,ni if(nosc.gt.nmx)then write(1,*) 'Parameter nmx must be increased' stop end if ccc SINGLE PARTICLE STATES if(cof.lt.0.001D0) call wsax if(cof.gt.0.001D0) then call hf call lecpot end if dr=h if(cof.lt.0.001D0) rmass=(anucl+1.D0)/anucl if(cof.gt.0.001D0) rmass=(anucl-1.D0)/anucl read *,jtot,itot,ipi,irpa ccc jtot,itot,ipi=spin,isospin,parity of rpa states ccc irpa=1 if tda,=2 if rpa,=3 if hf write(1,118),jtot,itot,ipi,irpa bosc=1.d0/afa0 hbom=hbc*hbc/(bmas*bosc) bosc=dsqrt(bosc) lxmy=0 do 1 i=1,nos read 100,nn(i),ll(i),lj(i),lt(i) ccc quantum numbers n,l,j,t(z) of the single particle states (n,l have ccc their true values, lj=2*j, lt=1 if proton and =0 if neutron) noc=nn(i) lp=ll(i) jp=lj(i) lt1=lt(i)+1 call qwg(bosc,nosc,lp,jp,lt1) ehf(i)=ewg(noc) do 300 j=1,nri dwf(j,i)=dwg(j,noc) 300 wf(j,i)=wg(j,noc) xnor(i)=xor(noc) lxmy=max0(lxmy,ll(i)) 1 continue nos1=nos+1 i=nos 200 read 100,lp,jp,ltp,n1,n2 ccc n1,n2=number of nodes of first and last unoccupied state ccc considered, for each given set (lp,jp,ltp) if(lp.lt.0) go to 201 ltp1=ltp+1 lxmy=max0(lxmy,lp) ibb=(jp-2*lp+3)/2 nqqy(lp+1,ibb)=n1 nqqx(lp+1,ibb)=n2 call qwg(bosc,nosc,lp,jp,ltp1) do 202 ii= n1,n2 i=i+1 if(i.le.nsp) go to 203 write(1,103),i,nsp stop 203 nn(i)=ii ll(i)=lp lj(i)=jp lt(i)=ltp ehf(i)=ewg(ii) do 204 j=1,nri dwf(j,i)=dwg(j,ii) 204 wf(j,i)=wg(j,ii) 202 continue go to 200 201 ntot=i lxmy=lxmy+1 28 it1=itot+1 go to (4,5),it1 4 f1=0.75D0*t0 g1=-0.25D0*t0*(1.D0-2.D0*x0) vv1=3.d0*t1/8 vv2=t2*(5.d0+4.d0*x2)/4.d0 do 6 j=1,nri f0(j)=f1+(t3/48.D0)*(dt(j)**alfe)*(3.D0*(alfe+1.D0)*(alfe+2. 1D0)+alfe*(1.D0-alfe)*(1.D0+2.D0*x3)*(((dn(j)-dp(j))/dt(j))**2)) g0(j)=g1-(t3/24.D0)*(1.D0-2.D0*x3)*(dt(j)**alfe) 6 continue go to 7 5 f1=-0.25D0*t0*(1.D0+2.D0*x0) g1=-0.25D0*t0 vv1=-t1*(1.d0+2.d0*x1)/8.d0 vv2=t2*(1.d0+2.d0*x2)/4.d0 do 8 j=1,nri f0(j)=f1-(t3/24.D0)*(1.D0+2.D0*x3)*(dt(j)**alfe) g0(j)=g1-(t3/24.D0)* (dt(j)**alfe) 8 continue 40 continue 7 continue avv1=-3.d0*t1/32.d0 avv2=-2.d0*avv1 avv4=-0.25d0*t2*(x2+5.d0/4.d0) bvv1=-0.125d0*t1*(x1/2.d0-0.25d0) bvv2=-2.d0*bvv1 bvv4=-0.25d0*t2*(x2/2.d0+0.25d0) calf(1)=avv1 calf(2)=avv2 calf(3)=avv2 cbet(1)=bvv1 cbet(2)=bvv2 cbet(3)=bvv2 do 500 iy=4,7 iyy=1-2*(iy/6) calf(iy)=iyy*avv4 cbet(iy)=iyy*bvv4 500 continue if(iliu.eq.0) go to 305 do 306 j=1,nri 306 g0(j)=0.d0 305 continue l1g=iabs(jtot-1) l2g=jtot+1 l1=jtot if(jtot.eq.0) l1=2 do 234 j=1,nri xj=dr*j opr(j)=xj**l1 234 continue 5556 format(5x,4e15.5) 5555 continue ccc BUILDING CONFIGURATIONS i=0 do 205 ihh=1,nos if(ihh.lt.nossp)go to 205 if(ihh.gt.npo.and.ihh.lt.nossn)go to 205 l1=ll(ihh) j1=lj(ihh) lt1=lt(ihh) do 206 ipp=nos1,ntot l2=ll(ipp) j2=lj(ipp) lt2=lt(ipp) if(lt1.ne.lt2) go to 206 l12=l1+l2 l12=1-2*mod(l12,2) if(l12.ne.ipi) go to 206 j12m=iabs(j1-j2)/2 j12p=(j1+j2)/2 if((jtot-j12m)*(jtot-j12p)) 207,207,206 207 i=i+1 ih(i)=ihh ip(i)=ipp ecf(i)=ehf(ipp)-ehf(ihh) 206 continue 205 continue nconf=i nconf2=2*nconf write(1,114) write(1,105) do 222 i=1,ntot 222 write(1,106),i,nn(i),ll(i),lj(i),lt(i),ehf(i),xnor(i) write(1,114) write(1,1234),nconf 1234 format(2x,'nconf=',i4//) write(1,107) if(nconf.gt.ncf) then write(1,*) 'Parameter ncf must be increased' stop end if ccc BUILDING MATRICES AA AND BB sm1=0.D0 s0=0.D0 s1=0.D0 s3=0.D0 do 208 i1=1,nconf ip1=ip(i1) ih1=ih(i1) lp1=ll(ip1) lh1=ll(ih1) jp1=lj(ip1) jh1=lj(ih1) xp1=lp1*(lp1+1.d0) xh1=lh1*(lh1+1.d0) zlp1=lp1 zlh1=lh1 zjp1=0.5d0*jp1 zjh1=0.5d0*jh1 yp1=1.d0 yh1=1.d0 if(lp1.eq.0) yp1=0.d0 if(lh1.eq.0) yh1=0.d0 yph1=1.d0 k=(jtot-iabs(lp1-lh1))*(jtot-lp1-lh1) if(k.gt.0) yph1=0.d0 yl1=yl(lp1,jp1,lh1,jh1,jtot) do 233 lg=l1g,l2g 233 tjl1(lg)=tjl(lp1,jp1,lh1,jh1,jtot,lg) x=0.D0 do 235 j=1,nri 235 x=x+opr(j)*wf(j,ip1)*wf(j,ih1) rmat(i1)=dr*x*yl1 x=rmat(i1)**2 sm1=sm1+x/ecf(i1) s0=s0+x s1=s1+x*ecf(i1) s3=s3+x*(ecf(i1)**3) do 209 i2=i1,nconf ip2=ip(i2) ih2=ih(i2) jp2=lj(ip2) jh2=lj(ih2) aaa=0.d0 if(irpa.eq.3) go to 5559 nrip=nri call phme(ip1,ih2,ip2,ih1,jtot,nrip,dr,aaa) 5559 continue if (i1.eq.i2) aaapr=aaa if(i1.eq.i2) aaa=aaa+ecf(i1) aa(i1,i2)=aaa aa(i2,i1)=aaa bbb=0.d0 if(irpa.eq.1) go to 209 if(irpa.eq.3) go to 209 jbb=(jp2+jh2)/2 sgn=1.d0-2.d0*mod(jbb,2) call phme(ip1,ip2,ih2,ih1,jtot,nrip,dr,bbb) bbb=sgn*bbb bb(i1,i2)=bbb bb(i2,i1)=bbb if(i1.eq.i2) bbbpr=bbb 209 continue write(1,108),i1,ih1,ip1,ecf(i1),x,aaapr,bbbpr 208 continue write(1,109),sm1,s0,s1,s3 ccc DIAGONALIZATION AND NORMALIZATION nc=nconf irpa1=1 if(irpa.eq.2) irpa1=2 if(irpa.eq.2) nc=nconf2 nc1=nc+1 call diagma(nconf,irpa1) write(1,114) write(1,112),(evr(i),i=1,nconf) do 218 i=1,nconf if(evr(i).le.0.D0) go to 218 xno=0.d0 do 219 j=1,nc sg=1.D0 if(j.gt.nconf) sg=-1.D0 219 xno=xno+sg*(vecr(j,i)**2) if(xno.gt.0.D0) go to 220 write(1,104),i stop 220 xno=dsqrt(xno) do 221 j=1,nc 221 vecr(j,i)=vecr(j,i)/xno 218 continue ccc ORDERING OF STATES do 224 i=1,nconf 224 nord(i)=i do 225 i=1,nconf k=i x=evr(i) do 226 j=1,nc 226 vec(j)=vecr(j,i) do 227 ii=i,nconf if(x-evr(ii)) 227,227,228 228 k=ii x=evr(ii) do 229 j=1,nc 229 vec(j)=vecr(j,ii) 227 continue kk=nord(k) nord(k)=nord(i) nord(i)=kk evr(k)=evr(i) do 230 j=1,nc 230 vecr(j,k)=vecr(j,i) evr(i)=x do 231 j=1,nc 231 vecr(j,i)=vec(j) 225 continue ccc TRANSITION PROBABILITIES sm1=0.D0 s0=0.D0 s1=0.D0 s3=0.D0 do 236 i=1,nconf if(evr(i).le.0.d0) go to 236 x=0.d0 do 237 ii=1,nconf y=vecr(ii,i) if(irpa.eq.2) y=y-vecr(ii+nconf,i) 237 x=x+y*rmat(ii) sqbel(i)=x x=x*x bel(i)=x sm1=sm1+x/evr(i) s0=s0+x s1=s1+x*evr(i) s3=s3+x*(evr(i)**3) 236 continue write(1,114) write(1,109),sm1,s0,s1,s3 write(1,110) ii=0 do 238 i=1,nconf if(evr(i).lt.0.D0) go to 238 x=bel(i)/s0 ii=ii+1 nig(ii)=i 239 continue y=evr(i)*bel(i)/s1 write(1,111),i,evr(i),bel(i),x,y 238 continue ngr=ii if(ipunch.eq.0) go to 240 write(12,100) jtot,itot,ipi,nconf,nos,ntot,nosc,ngr,lxmy do 241 i=1,lxmy 241 write(12,100) (nqqx(i,ii),ii=1,2),(nqqy(i,ii),ii=1,2) do 242 i=1,ntot 242 write(12,100) nn(i),ll(i),lj(i),lt(i) write(12,100) (ih(i),i=1,nconf) write(12,100) (ip(i),i=1,nconf) do 243 ii=1,ngr i=nig(ii) write(12,101) evr(i) write(12,101) (vecr(j,i),j=1,nconf) 243 continue 240 continue stop end SUBROUTINE WSAX implicit real*8 (a-h,o-z) parameter nnn=200 common/bpot/xmn(nnn),xmp(nnn),vn(nnn),vp(nnn),yn(nnn),yp(nnn), 1vsn(nnn),vsp(nnn),vc(nnn),wnd(nnn),wpd(nnn) common/bpas/h,cof common/blecp1/t0,t1,t2,t3,t13,alfe,x0,x1,x2,x3,anucl,znucl common/blecp2/dt(nnn),taut(nnn),dt1(nnn),dt2(nnn),dp(nnn),dn(nnn) dimension dal(3),mash(3) 208 format(6e12.6) 209 format(15i3) read *,t0,t1,t2,t3,t13,x0 read *,alfe,x3,x1,x2 read *,anucl,znucl,h,(dal(i),i=1,3) read *,v0n,r0n,a0n,vson,rson,ason read *,v0p,r0p,a0p,vsop,rsop,asop read *,rcb write(1,102),v0n,r0n,a0n,vson,rson,ason write(1,103),v0p,r0p,a0p,vsop,rsop,asop write(1,104),rcb 100 format(2x,'woods-saxon parameters v0=',e12.5,5x,'r0=', 1e12.5,5x,'a0=',e12.5) 102 format(2x,'v0n,r0n,a0n,vson,rson,ason',3x,6e14.5) 103 format(2x,'v0p,r0p,a0p,vsop,rsop,asop',3x,6e14.5) 104 format(2x,'coulomb radius',3x,e14.5) rcb3=rcb**3 rcb2=rcb**2 e2=(znucl-1.d0)*197.3d0/137.d0 do 101 j=1,nnn x=h*j xmn(j)=(anucl+1.)/(anucl*0.04819) xmp(j)=xmn(j) yn(j)=0. yp(j)=0. x1=h*(j-1) x2=h*(j+1) wnd(j)=0.d0 wpd(j)=0.d0 z=(x-r0n)/a0n vn(j)=v0n/(1.d0+dexp(z)) z=(x-r0p)/a0p vp(j)=v0p/(1.d0+dexp(z)) z1=(x1-rson)/ason z2=(x2-rson)/ason z=(1.d0/(1.d0+dexp(z2))-1.d0/(1.d0+dexp(z1)))/(2.d0*h) vsn(j)=-vson*z/x z1=(x1-rsop)/asop z2=(x2-rsop)/asop z=(1.d0/(1.d0+dexp(z2))-1.d0/(1.d0+dexp(z1)))/(2.d0*h) vsp(j)=-vsop*z/x if(x.gt.rcb) go to 105 vc(j)=(3.d0*rcb2-x*x)*e2/(2.d0*rcb3) go to 101 105 vc(j)=e2/x 101 continue return end SUBROUTINE HF parameter nnn=200,neoc=38 ccc neoc=greater or equal than the number of occupied states common/bhf/nos,npo,ni common/bpas/h,cof common/der/fct(nnn),df(nnn) common rm(nnn),du(nnn),vc(nnn),u(nnn),vs(nnn),taun(nnn), +taup(nnn),taut(nnn),djn1(nnn),djp1(nnn),djt1(nnn),vn(nnn),vp(nnn), +vson(nnn),vsop(nnn),ecoulb(nnn),vqre(nnn),rm2(nnn),rm3(nnn), +hmen(nnn),hmep(nnn),dhmen(nnn),dhmep(nnn),d2hmen(nnn),d2hmep(nnn), +dn(nnn),dp(nnn),dt(nnn),dn1(nnn),dp1(nnn),dt1(nnn),dn2(nnn), +dp2(nnn),dt2(nnn),djn(nnn),djp(nnn),djt(nnn),dc(nnn),avp(nnn, +neoc),dunl(nnn,neoc),unl(nnn,neoc) dimension wf(10,150) dimension rmsm(150) dimension nn(neoc),ll(neoc),lj(neoc),lt(neoc),deg(neoc),ehf(neoc), +evso(neoc) dimension ov(neoc,neoc) dimension rec(8),unlosc(nnn,neoc) dimension dal(3),mash(3) character*6 sx(31),st(neoc) dimension nx(31),lx(31),jx(31) dimension vq(150),uq(150) real ma,mesmp(nnn),mesmn(nnn) data cp/.79577471e-01/,qp/12.566371/ data sx/' 1s1/2',' 1p3/2',' 1p1/2',' 1d5/2',' 2s1/2',' 1d3/2', 1' 1f7/2',' 2p3/2',' 1f5/2',' 2p1/2',' 1g9/2',' 2d5/2',' 1g7/2', 2' 3s1/2',' 2d3/2','1h11/2',' 2f7/2',' 1h9/2','1i13/2',' 3p3/2', 3' 2f5/2',' 3p1/2','1i11/2',' 2g9/2','1j15/2',' 3d5/2',' 4s1/2', 4' 2g7/2',' 3d3/2',' 3f5/2',' 3f7/2'/ data nx/1,1,1,1,2,1,1,2,1,2,1,2,1,3,2,1,2,1,1,3,2,3,1,2,1,3,4,2,3, 13,3/ data lx/0,1,1,2,0,2,3,1,3,1,4,2,4,0,2,5,3,5,6,1,3,1,6,4,7,2,0,4,2, 13,3/ data jx/1,3,1,5,1,3,7,3,5,1,9,5,7,1,3,11,7,9,13,3,5,1,11,9, 115,5,1,7,3,5,7/ 200 format( //13h iterations i2,10h xmu=f5.3//) 201 format(2x,a6,i2,5x,'E=',e11.4,5x,'deg=',e11.4,5x,' d=',e8.2) 203 format(/) 204 format(5x,a6,i2,5x,'not bound') 205 format(//' Nucleus with A=',f5.1,3x,'and Z=',i2/) 206 format(1h1) 209 format(8e10.3) 211 format(10(2x,e12.5)) 2110 format(2x,'t0,t1,t2,t13 =',5(2x,e12.6)) 2111 format(2x,'x0,x1,x2,x3,w1,w2 =',6(2x,e12.6)) 2112 format(2x,'alfe =',e12.6) 212 format(10e12.5) 213 format(//17h neutron density/) 214 format(//16h proton density/) 215 format(5e15.8) 216 format(//18h charge density /) 217 format(/6h ec/a=e15.8,7h ehf/a=e15.8,8h erea/a=e15.8,8h etot/a=e 1 15.8//) 226 format(/10h fichier i4//) 221 format(//4h rn=e13.5,4h rp=e13.5,4h rc=e13.5//) 400 format(20i3) 4001 format(2x,'ioc1,ideg1,ioc2,ideg2,ioc3,ideg3,ioc4,ideg4,ioc5,ideg5 1,ioc6,ideg6 =',/3x,12(2x,i2)) 4002 format(2x,'kop1,ipchpo,ipchwf,idcc,ipunch =',5(2x,i2)) 401 format(6e12.6) 403 format(8(f12.2,3x)) 1200 format(//2x,'mass density'/) 1201 format(120('*')) read *,t0,t1,t2,t3 read *,x0,x1,x2,x3 read *,w1,w2 read *,t13 read *,alfe ccc parameters of the Skyrme interaction (w2=w, w1=2*w, where w is the ccc spin-orbit parameter) write(1,2110),t0,t1,t2,t3,t13 write(1,2111),x0,x1,x2,x3,w1,w2 read *,ibox ccc if ibox=1 the system can be put in a finite box if(abs(alfe).lt.1.e-04) alfe=1. write(1,2112),alfe 1 read*,ioc1,ideg1,ioc2,ideg2,ioc3,ideg,ioc4,ideg4,ioc5,ideg5,ioc6,ideg6 ccc if ioc1=n, then the occupation factor of state n, instead of being ccc 2*j(n)+1, it is set to ideg1, and the same for the others ioc's, ideg's if(ioc1.lt.0) go to 1000 write(1,4001),ioc1,ideg1,ioc2,ideg2,ioc3,ideg3 1,ioc4,ideg4,ioc5,ideg5,ioc6,ideg6 read *,kop1,ipchpo,ipchwf,idcc,ipunch write(1,4002),kop1,ipchpo,ipchwf,idcc,ipunch write(1,1201) do 2 j=1,nnn 2 u(j)=0. ipre = 2 do 95 j=1,nnn rm(j) = h*float(j) rm2(j) = rm(j) * rm(j) rm3(j) = rm(j)*rm2(j) 95 continue uns4pi = cp nn1 = nnn- 1 nn2 = nnn- 2 nn3 = nnn- 3 nn4 = nnn- 4 w3 = w2 w0 = w1 cro = .25*(t1*(1.+.5*x1)+t2*(1.+.5*x2)) croq = .125*(t2*(1.+2.*x2) -t1*(1.+2.*x1)) txro = t0*(1. + .5*x0) txroq = -t0*(x0 +.5) txdro = .5*(t2*(1.+.5*x2)-t1*(1.+.5*x1)) txd2ro = .125*(t2*(1.+.5*x2) - 3.*t1*(1.+.5*x1)) txdroq = .25*(t1*(1.+2.*x1) + t2*(1.+2.*x2)) txd2rq = .0625*(3.*t1*(1. +2.*x1) + t2*(1.+2.*x2)) txtau = cro txtauq = croq c13=t13/24. do 344 i = 1,nos if (i-npo) 376,376,377 376 lt(i) = 1 j = i go to 378 377 lt(i) = 0 j = i - npo 378 nn(i) = nx(j) ll(i) = lx(j) lj(i) = jx(j) 344 st(i) = sx(j) ccc DEFINITION OF OCCUPATION FACTORS do 402 i = 1,nos 402 deg(i) = lj(i) + 1 if(ioc1.ne.0) deg(ioc1)=float(ideg1) if(ioc2.ne.0) deg(ioc2)=float(ideg2) if(ioc3.ne.0) deg(ioc3)=float(ideg3) if(ioc4.ne.0) deg(ioc4)=float(ideg4) if(ioc5.ne.0) deg(ioc5)=float(ideg5) if(ioc6.ne.0) deg(ioc6)=float(ideg6) jz = 0. ma = 0. do 100 i=1,nos jz = jz + lt(i)*deg(i) 100 ma = ma + deg(i) write(1,205),ma,jz dmshb = .04823*ma/(ma-1) hb = 1./dmshb ccc WOODS-SAXON r0 = 1.21 a1 = 55.5 a2 = 33.2 a3 = .36 a5 = .68 cwpi = 1.41 cw2 = cwpi*cwpi al = a5 r = r0*ma**(1./3.) ccc START iter = 0 go to(137,341,601), kop1 137 read 209, ((avp(j,i),j=1,nnn),i=1,nos),(ehf(i),i=1,nos) go to 107 341 continue do 94 j=1,nnn hmen(j) = hb hmep(j) = hb dhmen(j) = 0. dhmep(j) = 0. d2hmen (j) = 0. d2hmep (j) = 0. mesmn(j) = 1. mesmp(j) = 1. 94 continue kx = 1 e = 0. e22 = 1.43986* h f = 0. rc = r*1.09/r0 sym = (ma-2*jz)/ma z = 0. do 96 k = 1,nnn x = rm(k) dp(k) = 1./(1. + exp((x-rc)/.55)) 96 z = z + x*x*dp(k) y = float(jz)/(qp*z*h) do 97 k = 1,nnn 97 dp(k) = y*dp(k) do 98 k=1,nnn y = qp*rm(k)*rm(k) e = e + y*dp(k) f = f + y*dp(k)/rm(k) 98 vc(k) = e/rm(k) - f ff = f do 101 j=1,nnn vc(j) = e22*(vc(j) + ff) x = rm(j) pe = exp((x-r)/al) fpe = -1./(1.+pe) vn(j) = fpe vs(j) = - pe*fpe*fpe/al do 101 i=1,nos n = nn(i) l = ll(i) jm = lj(i) fl = (jm*(jm+2) - 4*l*(l+1) - 3)/8. isig = 2*lt(i) - 1 v = a1 + isig*a2*sym avp(j,i)=v*vn(j)+a3*cw2*vs(j)*2.*fl if (i.le.npo) avp(j,i) = avp(j,i) + vc(j) 101 continue do 102 i=1,nos ehf(i) = .3 * avp(1,i) 102 if(i.le.npo) ehf(i) = ehf(i) +.3*vc(1) 601 continue 107 continue itmax = ni ccc ITERATIONS xmu=0.3 104 iter = iter + 1 if(iter.eq.10) xmu=0.5 if(iter.eq.15) xmu=0.7 if(iter.eq.itmax) write(1,200),iter,xmu ccc SOLUTIONS OF SCHROEDINGER EQUATION epote=0. do 110 i=1,nos n = nn(i) l = ll(i) ll1 = l*(l+1) ei = ehf(i) k = 0 s2 = 0. do 300 j=1,nnn if( lt(i).eq.0) vqre(j) = mesmn(j) * (avp(j,i) - .25 * dhmen(j) * 1 dhmen(j) / hmen(j) + .5*d2hmen(j) ) + (1.-mesmn(j))*ei if( lt(i).eq.1) vqre(j) = mesmp(j) * (avp(j,i) - .25 * dhmep (j)* 1 dhmep(j) / hmep(j) + .5*d2hmep(j) ) + (1.-mesmp(j))*ei vqre(j) = vqre(j)/hb + ll1/(rm(j)*rm(j)) 300 continue go to 112 113 k=k+2 if ( k- 200) 114,114,115 115 write(1,204),st(i),lt(i) go to 111 114 ei = -.5*float(k) 112 e = ei * dmshb no = n - 1 call num1l(nnn,h,e,s2,vqre,u,no,.1e-03) if (no.lt.0) go to 113 111 continue do 12 j=1,nnn y = sqrt( mesmn(j)*(1-lt(i)) + mesmp(j)*lt(i)) u(j) = u(j) * y 12 fct(j) = u(j) call deriv s = 0. do 303 j= 1,nnn du(j) = df(j) 303 s = s + u(j)*u(j) s = sqrt(s*h) sig = u(1)/abs(u(1)) do 304 j=1,nnn du(j) = sig * du(j)/s 304 u(j) = sig*u(j)/s e = e/dmshb ehf(i) = e do 122 j=1,nnn dunl(j,i) = du(j) 122 unl(j,i) = u(j) if(iter.ne.itmax) go to 110 xkz=0. do 2222 j=1,nnn zy=(-unl(j,i)/rm(j)+dunl(j,i))/rm(j) zy1=unl(j,i)/rm2(j) 2222 xkz=xkz+(zy*zy+ll1*zy1*zy1)*rm2(j) xkz=sqrt(h*xkz) x = (e-ei)/e x = abs(x) write(1,2221),st(i),lt(i),e,deg(i),x,xkz 2221 format(2x,a6,i2,5x,'e=',e11.5,5x,'deg=',e11.4,5x,' d=', 1e8.2,5x,' k=',e12.5) epote=epote+e*deg(i) 110 continue ccc DENSITIES 352 do 143 j=1,nnn x = rm(j) r2 = rm2(j) r3 = rm3(j) e = 0. f = 0. e2= 0. f2 = 0. te = 0. tf= 0. do 142 i= 1,nos l = ll(i) n = nn(i) jm = lj(i) ll1 = l*(l+1) fl = (jm*(jm+2)-4 *l*(l+1) - 3.)/8. evso (i) = fl y = unl(j,i)*unl(j,i)*deg(i) e = e + y*(1-lt(i)) f = f + y*lt(i) ee = 2*fl*y e2 = e2 + ee*(1-lt(i)) f2 = f2 + ee*lt(i) y = (-unl(j,i)/x + dunl(j,i))/x y1 = unl(j,i)/r2 y2 = deg(i) *(y*y + ll1*y1*y1) te = te + y2*(1-lt(i)) tf = tf + y2 *lt(i) 142 continue dp(j) = f*uns4pi/r2 dn(j) = e* uns4pi/r2 dt(j) = dn(j) + dp(j) djn(j) = e2*uns4pi/r3 djp(j) = f2*uns4pi/r3 djt(j) = djn(j) + djp(j) taun(j) = te*uns4pi taup(j) = tf * uns4pi taut(j) = taun(j) + taup(j) 143 continue ccc DERIVATIVES OF DENSITES do 13 j=1,nnn 13 fct(j) = dn(j) call deriv do 14 j=1,nnn dn1(j) = df(j) 14 fct(j) = dp(j) call deriv do 15 j=1,nnn dp1(j) = df(j) dt1(j) = dn1(j) + dp1(j) 15 fct(j) = dn1(j) call deriv do 16 j=1,nnn dn2(j) = df(j) 16 fct(j) = dp1(j) call deriv do 17 j=1,nnn dp2(j) = df(j) dt2(j) = dp2(j) + dn2(j) 17 fct(j) = djn(j) call deriv do 18 j=1,nnn djn1(j) = df(j) 18 fct(j) = djp(j) call deriv do 19 j=1,nnn djp1(j) = df(j) 19 djt1(j) = djp1(j) + djn1(j) ccc CALCULATION OF FIELDS ccc COULOMB TERM e=0. f = 0. do 144 j=1,nnn x = rm(j) y = qp*x*x e = e + y*dp(j) f = f + y*dp(j)/x 144 vc(j) = e/x - f e222=((12./qp)**0.333333)*e22/h do 11 j=1,nnn 11 vc(j)=e22*(vc(j)+f)-e222*(dp(j)**0.333333) do 167 j = 1,nnn hn=hb+cro*dt(j)+croq*dn(j)+c13*(dt(j)**2-dn(j)**2) hp=hb+cro*dt(j)+croq*dp(j)+c13*(dt(j)**2-dp(j)**2) hmen(j)=hn*(1.-xmu)+hmen(j)*xmu hmep(j)=hp*(1.-xmu)+hmep(j)*xmu mesmn(j) = hb/hmen(j) mesmp(j) = hb/hmep(j) hn1=cro*dt1(j)+croq*dn1(j)+c13*2.*(dt(j)*dt1(j)-dn(j)*dn1(j)) hp1=cro*dt1(j)+croq*dp1(j)+c13*2.*(dt(j)*dt1(j)-dp(j)*dp1(j)) hn2=cro*dt2(j)+croq*dn2(j)+c13*2.*(dt1(j)**2+dt(j)*dt2(j) 1-dn1(j)**2-dn(j)*dn2(j)) hp2=cro*dt2(j)+croq*dp2(j)+c13*2.*(dt1(j)**2+dt(j)*dt2(j) 1-dp1(j)**2-dp(j)*dp2(j)) dhmen(j)=hn1*(1.-xmu)+dhmen(j)*xmu dhmep(j)=hp1*(1.-xmu)+dhmep(j)*xmu d2hmen(j)=hn2*(1.-xmu)+d2hmen(j)*xmu d2hmep(j)=hp2*(1.-xmu)+d2hmep(j)*xmu 167 continue do 1500 j=1,nnn vnt33=(1.-x3)*(alfe+2.)*dt(j)*dt(j) vnt3=(vnt33+(2.+4.*x3)*dp(j)*(dt(j)+alfe*dn(j)))/24. vpt3=(vnt33+(2.+4.*x3)*dn(j)*(dt(j)+alfe*dp(j)))/24. if(abs(alfe-1.).lt.1.e-06) go to 8888 vnt=dt(j)**(alfe-1.) vnt3=vnt3*vnt vpt3=vpt3*vnt 8888 continue vnt3=-vnt3 vpt3=-vpt3 vn(j) = txro*dt(j) + txdro *dt1(j)/rm(j) + txd2ro*dt2(j) + 1txtau*taut(j)-w3*(djt(j)/rm(j)+djt1(j)/2.) vn(j)=vn(j)-c13*(2.5*dt(j)*(dt2(j)+2.*dt1(j)/rm(j)) 1+1.25*(dt1(j)**2)-2.*dt(j)*taut(j)) vson(j) = w0* dt1(j)/(2.* rm(j)) - .125*( t1*x1 + t2*x2)*djt(j) 1 /rm(j) vp (j) = vn(j) + txroq*dp(j) + txdroq* dp1(j)/rm(j) +txd2rq * 1 dp2(j) +txtauq *taup(j) - w3*(djp(j)/rm(j) +.5*djp1(j)) - .25* 2t3*4.*vpt3 vp(j)=vp(j)-c13*(-2.5*dp(j)*(dp2(j)+2.*dp1(j)/rm(j))-1.25*(dp1(j) 1**2)+2.*dp(j)*taup(j)-0.25*(djp(j)**2)-0.50*(djn(j)**2)) vsop(j) = vson(j) + w0*dp1(j)/(2.*rm(j)) do 1502 i=1,npo 1502 avp(j,i) = (vp(j) +vsop(j)*evso(i) +vc(j))*(1. -xmu) + xmu*avp(j, 1 i) vn(j) = vn(j) + txroq*dn(j) + txdroq*dn1(j)/rm(j) + txd2rq* dn2(j 1) + txtauq * taun(j) - w3*( djn(j)/rm(j) + .5*djn1(j)) -.25*t3 2*4.*vnt3 vn(j)=vn(j)-c13*(-2.5*dn(j)*(dn2(j)+2.*dn1(j)/rm(j))-1.25*(dn1(j) 1**2)+2.*dn(j)*taun(j)-0.25*(djn(j)**2)-0.50*(djp(j)**2)) vson(j) = vson(j) + w0*dn1(j)/(2.*rm(j)) nneut = nos - npo do 1504 i2=1,nneut i = i2 + npo 1504 avp(j,i) = (vn(j) + vson(j)*evso(i))*(1.-xmu) + xmu*avp(j,i) 1500 continue if (iter. ge.itmax) go to 116 go to 104 ccc END OF ITERATION 116 continue ecbd=0. ecbe=0. do 5555 j=1,nnn ecbd=ecbd+vc(j)*qp*rm2(j)*dp(j)*h*0.5 5555 ecbe=ecbe-e22*rm2(j)*(dp(j)**1.333333) ecbe=ecbe*0.75*((12./qp)**0.333333)*qp write(1,5556),ecbd,ecbe 5556 format(//1x,'ecb(direct)=',e12.5,5x,'ecb(exchange)=',e12.5) if(ibox.ne.1) go to 461 450 format(i5,5e12.6) 451 format(10i3) 452 format(120('*')/2x,'l=',i3,3x,'j=',i3,'/2',3x,'nt=',i3,f10.5//) 453 format(2x,'n=',i4,5x,'not converged') 454 format(2x,'n=',i3,5x,10e10.4) 455 format(8f10.2) 456 format(50(1h*)) read 450,nmax,eps,em em1=em*dmshb 460 read 451,l,j,nt,ind,inf if(l) 461,462,462 462 inn=inf-ind+1 n=ind-1 jk=0 write(1,452),l,j,nt,em 463 n=n+1 jk=jk+1 if(n.gt.inf) go to 4460 ll1=l*(l+1) 474 k=0 s2=0. eso=(j*(j+2)-4*l*(l+1)-3.)/8. nmax1=nmax+1 do 464 kk=1,nmax1 if(kk-nnn) 465,465,466 465 if(nt.eq.0) vq(kk)=mesmn(kk)*(vn(kk)+eso*vson(kk)-0.25*dhmen(kk)* 1dhmen(kk)/hmen(kk)+0.5*d2hmen(kk)) if(nt.eq.1) vq(kk)=mesmp(kk)*(vp(kk)+vc(kk)+eso*vsop(kk)-0.25* 1dhmep(kk)*dhmep(kk)/hmep(kk)+0.5*d2hmep(kk)) vq(kk)=vq(kk)/hb+ll1/(rm(kk)*rm(kk)) rmsm(kk)=(1.-mesmn(kk))*(1.-nt)+(1.-mesmp(kk))*nt go to 464 466 vq(kk)=ll1/(h*h*float(kk*kk)) rmsm(kk)=0. 464 continue go to 480 470 write(1,453),n go to 463 480 no=n-1 call num2b(nmax,h,e,s2,vq,rmsm,uq,no,eps,em1,hb) if(no.lt.0) go to 470 e=e/dmshb nol=no+1 do 4461 kx=1,nmax vq(kx)=vq(kx)/dmshb 4461 wf(nol,kx)=uq(kx) write(1,454),nol,e write(1,211),(uq(ik),ik=1,nmax) 5006 format(2x,'l=',i2,2x,'j=',i2,'/2',2x,'lh=',i2,2x,'jh=',i2,'/2') 5007 format(10e12.5) 5003 continue write(1,456) write(1,1001) 1001 format(2x,'potential'/) write(1,211),(vq(ik),ik=1,nmax) write(1,456) write(1,456) 9998 format(6e12.6) go to 463 4460 write(1,456) do 4464 k1=ind,inf do 4462 k2=k1,inf ov(k1,k2)=0. do 4463 kx=1,nmax 4463 ov(k1,k2)=ov(k1,k2)+wf(k1,kx)*wf(k2,kx)*h ov(k2,k1)=ov(k1,k2) 4462 continue write(1,211),(ov(k1,k2),k2=ind,inf) 4464 continue go to 460 461 continue n64=64 1002 format(a6,4i5,e12.6) mash(1)=12 mash(2)=24 mash(3)=36 nx1=2 nx2=2 nx3=2 dal(1)=0.20 dal(2)=0.20 dal(3)=0.20 if(ipchpo.eq.0) go to 1789 open(unit=10,status='new',form='formatted') write(10,9997) t0,t1,t2,t3,t13,x0 write(10,9997) alfe,x3,x1,x2 fjz=float(jz) write(10,9997) ma,fjz,h,(dal(i),i=1,3) write(10,9996) (mash(i),i=1,3) nxxx=1 nd=1 write(10,9996) nxxx,nd,nx1,nx2,nx3 write(10,9998) (hmen(j),j=1,nnn),(hmep(j),j=1,nnn) write(10,9998) (vn(j),j=1,nnn),(vp(j),j=1,nnn) write(10,9998) (dhmen(j),j=1,nnn),(dhmep(j),j=1,nnn) write(10,9998) (vson(j),j=1,nnn),(vsop(j),j=1,nnn) write(10,9998) (vc(j),j=1,nnn) write(10,9998) (dt(j),j=1,nnn),(taut(j),j=1,nnn) write(10,9998) (dt1(j),j=1,nnn),(dt2(j),j=1,nnn) write(10,9998) (dp(j),j=1,nnn),(dn(j),j=1,nnn) 9997 format(6e12.6) 9996 format(15i3) 1789 if(ipchwf.eq.0) go to 1790 do 9999 i=1,nos if(deg(i).lt.0.1) go to 9999 write(10,1002) st(i),nn(i),ll(i),lj(i),lt(i),ehf(i) i1i=nd-1+nx1 i1f=mash(1)*nx1 i2i=i1f+nx2 i2f=i1f+(mash(2)-mash(1))*nx2 i3i=i2f+nx3 i3f=i2f+(mash(3)-mash(2))*nx3 9999 continue 1790 continue erea=0. do 1501 j=1,nnn ere=-0.5*c13*1.25*(dp(j)*(dp1(j)**2) 1+dn(j)*(dn1(j)**2)-dt(j)*(dt1(j)**2)) ere=ere-0.5*c13*(taun(j)*(dp(j)**2)+taup(j)*(dn(j)**2)+2.*dt(j) 1*(taup(j)*dp(j)+taun(j)*dn(j)-taut(j)*dt(j))) ere=ere-c13*0.125*((dp(j)-2.*dt(j))*(djp(j)**2) 1+(dn(j)-2.*dt(j))*(djn(j)**2)) erea=erea-rm2(j)*ere 1501 erea=erea-rm2(j)*(0.5*(1.-x3)*dt(j)*dt(j)+(2.*x3+1.)*dn(j)* 2dp(j))*(dt(j)**alfe)*t3*alfe/24. erea=erea*h*qp erea=erea/ma ectot=0. do 1510 j=1,nnn 1510 ectot=ectot+taut(j)*rm2(j) ectot=ectot*qp*h*hb epote=epote/ma ectot=ectot/ma h0=0.5*(epote+ectot) etot=h0+erea write(1,217),ectot,h0,erea,etot ccc CHARGE DENSITY do 1503 j=1,nnn x=rm(j) 1503 dc(j)=dsc(x) ccc RADII rn=0. rt = 0. rp=0. rc=0. do 105 j=1,nnn x4=rm2(j)*rm2(j) rn=rn+dn(j)*x4 rt = rt + dt(j)*x4 rp=rp+dp(j)*x4 rc=rc+dc(j)*x4 105 continue rn=sqrt(qp*(h*rn/(ma-jz))) rt = sqrt(qp*(h*rt/ma)) rp=sqrt(qp*(h*rp/jz)) rc=sqrt(qp*(h*rc/jz)) write(1,221),rn,rp,rc write(1,1201) write(1,213) write(1,212),(dn(j),j=1,nnn) write(1,214) write(1,212),(dp(j),j=1,nnn) write(1,216) write(1,212),(dc(j),j=1,nnn) write(1,1200) write(1,212),(dt(j),j=1,nnn) if(idcc.eq.1) write(10,1209)(dc(j),j=1,nnn) close(10) rewind(10) 1209 format(5e16.9) lj1=1 7001 rp=0. rpn=0. do 7000 j=1,nnn x4=rm2(j)**lj1 rp=rp+dt(j)*x4 rpn=rpn+dn(j)*dp(j)*x4 7000 continue rp=qp*h*rp/ma write(1,7010),lj1,rp rpn=qp*h*rpn write(1,7011),rpn 7011 format(2x,'integral of ro(n)*ro(p)=',e13.5) 7010 format(/2x,'r**(2*',i2,'-2)=',e13.5) lj1=lj1+1 if(lj1.le.12) go to 7001 sum = 0. do 150 j=1,nos 150 sum=sum+deg(j)*(2.*nn(j)+ll(j)-.5) sum = sqrt(sum/ma) b = rt/sum ecm = .75*41.47/(b*b*ma) write(1,231),ecm,b,rt 231 format(/,' ecm/a=', e12.5, ' b=',e12.5, ' rt=',e12.5/) do 153 j=1,nos l = ll(j) do 152 jj=1,8 152 rec(jj) = 0. e = 0. do 155 k = 1,8 n = k - 1 do 154 nrec = 1,nnn x = h * nrec unlosc(nrec,j) = hwf(n,l,b,x) 154 rec(k) = rec(k) + h*unl(nrec,j)*unlosc(nrec,j) 155 e = e + rec(k)*rec(k) write(1,229),st(j) , lt(j) write(1,211),(rec(m),m=1,8),e 153 continue 229 format(5x,a6,i2,' rec=',e12.5) do 156 i = 1,nos l = ll(i) n = nn(i) - 1 do 156 j=1,nnn x = h*j unlosc(j,i) = hwf(n,l,b,x) 156 continue do 4500 j=1,nnn xr=h*float(j) 4500 df(j)=(3.*dt(j)+xr*dt1(j))*xr*xr write(1,4501) 4501 format(//2x,' Tassie monopole transition density times r**2'/) write(1,212),(df(j),j=1,nnn) ccc CALCULATION OF (D**2) dsq=rt*rt*ma/12. do 1003 ia=1,nos lta=lt(ia) la=ll(ia) ja=lj(ia) do 1004 ib=1,nos if(lt(ib)-lta) 1004,1008,1004 1008 lab=ll(ib)+la+1 if(mod(lab,2)) 1004,1007,1004 1007 jb=lj(ib) rab=0. do 1005 ix=1,nnn 1005 rab=rab+unl(ix,ia)*unl(ix,ib)*ix*h*h ans=0. dsq=dsq-(ja+1.)*(jb+1.)*rab*rab*ans*ans/12. 1004 continue 1003 continue write(1,1006),dsq 1006 format(120('*')/2x,'dsquare=',e12.5/120('*')//) do 7700 l1=1,10 7701 format(2i10) 7702 format(//2x,'l=',i2/20(1h*)) 7706 format(/4x,'n',6x,'sum-rule'/) do 7703 no=1,nope 7704 format(5e15.8) sumr=0. call deriv do 7705 j=1,nnn x4=df(j)**2+jope*(jope+1.)*(fct(j)**2)/rm2(j) 7705 sumr=sumr+dt(j)*x4*rm2(j) sumr=sumr*qp*h/(qp*0.04823) 7703 continue 7707 format(3x,i2,5x,e12.5) 7700 continue go to 1 1000 continue return end SUBROUTINE NUM1L(N,H,E,S2,U,S,NO,EPS) ccc integration of the Schroedinger equation with Noumerov algorithm ccc for negative E ccc search of eigenenergy by Raphson-Newton method dimension u(1),s(1) data rap1,rap2/0.,0./ h12=h*h/12. ccc CHECK OF ASYMPTOTIC CONDITIONS if(e.gt..0) e=.0 dei=.0 epss=.1e-10 if(u(n-1).gt.epss) go to 10 dei=u(n-1)-epss do 8 k=1,n 8 u(k)=u(k)-dei 10 u(n)=u(n-1) ccc CALCULATION OF BOUND STATES NUMBER BY INTEGRATION AT ZERO ENERGY s(1)=1.e-10 b0=0. aa=h12*u(1) if (s2) 16,18,16 16 b0=-s(1)*aa 18 b1=s(1)*(1.-aa) do 38 k=2,n b2=12.*s(k-1)-10.*b1-b0 if (abs(b2).lt.1.e+10) go to 22 b2=b2*1.e-20 b1=b1*1.e-20 22 aa=h12*u(k) s(k)=b2/(1.-aa) b0=b1 38 b1=b2 do 42 k=5,n n0=k if(u(k).lt.0.) go to 44 42 continue 44 nel=0 do 52 k=n0,n if (s(k-1)*s(k)) 46,50,52 46 nel=nel+2 go to 52 50 nel=nel+1 52 continue nel=nel/2 if(nel.gt.no) go to 64 if(nel.eq.no) go to 60 62 no=-1 return 60 rap1=s(n-1)/s(n) rap2=exp(h*sqrt(u(n-1)-e)) if(rap1.lt.rap2) go to 62 ccc CALCULATION OF EMIN ET EMAX (EMIN < EIGENENERGY < EMAX) 64 umin=u(1) do 70 k=2,n if(u(k).lt.umin) umin=u(k) 70 continue emin=umin emax=0. ccc START OF EIGENENERGY SEARCH te=emax-emin ccc REJECTION OF THE TRIAL ENERGY if((e.lt.emin).or.(e.gt.emax)) e=emin+.5*te e1=emin e2=emax j=2 i=1 go to 102 ccc REDUCTION OF EMIN AND EMAX 90 emin=e1 emax=e2 te=emax-emin j=2 98 i=1 100 e=emin+te*float(i)/float(j) 102 de=0. 104 e=e+de if(e.gt.0.) go to 204 s(n)=1.e-10 n1=n-1 rap2=exp(h*sqrt((u(n-1)+u(n))/2.-e)) s(n1)=s(n)*rap2 aa=h12*(u(n1)-e) b0=s(n)*(1.-aa) b1=s(n1)*(1.-aa) n1=n-2 do 138 kaux=1,n1 k=n1-kaux+1 b2=12.*s(k+1)-10.*b1-b0 aa=h12*(u(k)-e) s(k)=b2/(1.-aa) b0=b1 b1=b2 if(u(k).lt.e) go to 140 138 continue 140 n1=k ccc WAVEFUNCTION NORMALIZATION AT S(N1) do 146 kaux=n1,n k=n-kaux+n1 146 s(k)=s(k)/s(n1) ccc START OF INTEGRATION UP TO N1 s(1)=1.e-10 b0=0. aa=h12*(u(1)-e) if(s2) 156,158,156 156 b0=-s(1)*aa 158 b1=s(1)*(1.-aa) do 170 k=2,n1 b2=12.*s(k-1)-10.*b1-b0 aa=h12*(u(k)-e) s(k)=b2/(1.-aa) b0=b1 170 b1=b2 ccc NORMALIZATION do 174 k=1,n1 174 s(k)=s(k)/s(n1) ccc CALCULATION OF ENERGY CORRECTION som=0. do 180 k=1,n 180 som=som+s(k)*s(k) de=((-s(n1-1)+2.-s(n1+1))/(h*h)+u(n1)-e)/som if(abs(de).gt.eps) go to 104 ccc CALCULATION OF NUMBER OF NODES do 182 k=5,n if(u(k).lt.e) go to 184 182 continue 184 n0=k nel=0 do 192 k=n0,n1 if(s(k-1)*s(k)) 186,190,192 186 nel=nel+2 go to 192 190 nel=nel+1 192 continue nel=nel/2 if(nel-no) 198,214,202 198 if(e.gt.e1) e1=e go to 204 202 if(e.lt.e2) e2=e 204 i=i+2 if (i.le.j) go to 100 j=2*j if(abs(e1-emin).gt.eps.or.abs(emax-e2).gt.eps) go to 90 go to 98 ccc NORMALIZATION 214 som=1./sqrt(som*h) do 218 k=1,n 218 s(k)=s(k)*som e=e+dei return 2000 format(/,35x,57hThe required state is not bound. End of num1l with 1 no=-1,/) end SUBROUTINE DERIV ccc five point differentation formula (Abramowitz, pag 914) parameter nnn=200,neoc=38 common/bpas/h,cof common/der/fct(nnn),df(nnn) dimension a(5,5) data a(1,1),a(1,2),a(1,3),a(1,4),a(1,5) + /-50.,96.,-72.,32.,-6./ data a(2,1),a(2,2),a(2,3),a(2,4),a(2,5) + /-6.,-20.,36.,-12.,2./ data a(3,1),a(3,2),a(3,3),a(3,4),a(3,5) + /2.,-16.,0.,16.,-2./ data a(4,1),a(4,2),a(4,3),a(4,4),a(4,5) + /-2.,12.,-36.,20.,6./ data a(5,1),a(5,2),a(5,3),a(5,4),a(5,5) + /6.,-32.,72.,-96.,50./ data emfact/24./ nmx=nnn-2 do 1 j=1,nnn k=3 if(j.lt.3)k=j if(j.gt.nmx)k=j-nnn+5 sum=0. do 2 i=1,5 jj=j+i-k 2 sum=sum+a(k,i)*fct(jj) 1 df(j)=sum/(h*emfact) return end SUBROUTINE NUM2B(N,H,E,S2,U,RMS,S,NO,EPS,EM,HB) dimension u(150),rms(150),s(150),p(150),t(150) n1=n+1 h12=h*h/12. go to 1 100 no=-1 return 1 e0=-2.8938 de=0.24115 2 e0=e0+de do3 k=1,n1 3 p(k)=u(k)+rms(k)*e0 umin=p(1) do 4 k=2,n1 if(p(k).lt.umin) umin=p(k) 4 continue do 63 ie=1,2 kk=0 go to(64,65),ie 64 e0=umin-de go to 5 65 de=0.24115 e0=e1-de 5 e0=e0+de kk=kk+1 if(kk.gt.50) go to 100 do 6 k=1,n1 6 p(k)=u(k)+rms(k)*e0 s(1)=1.e-10 b0=0. aa=h12*(p(1)-e0) if(s2) 7,8,7 7 b0=-s(1)*aa 8 b1=s(1)*(1.-aa) do 9 k=2,n1 b2=12.*s(k-1)-10.*b1-b0 aa=h12*(p(k)-e0) s(k)=b2/(1.-aa) b0=b1 9 b1=b2 nel=0 do 10 k=8,n if(s(k-1)*s(k)) 11,12,10 11 nel=nel+2 go to 10 12 nel=nel+1 10 continue nel=nel/2 go to(66,67),ie 66 if(nel-no) 5,13,14 13 e1=e0 se1=s(n) nel1=nel go to 63 14 e0=e0-de de=0.5*de go to 5 67 if(nel-no-1) 5,15,14 15 e2=e0 se2=s(n) nel2=nel 63 continue kk=0 do 25 k=1,n1 25 t(k)=s(k) 36 if(nel2-nel1-1) 18,17,18 18 et1=e1*hb et2=e2*hb write(1,19),et1,se1,nel1,et2,se2,nel2 19 format(2x,2e12.5,i5,2e12.5,i5) go to 100 17 if(se1*se2) 21,20,18 20 if(se1.eq.0.) go to 18 23 e=e2 nel=nel2 go to 101 21 de=0.5*(e2-e1) deb=de*hb if(deb.gt.1.e-04) go to 22 go to 23 22 e3=e1+de kk=kk+1 if(kk.gt.50) go to 18 do 24 k=1,n1 24 p(k)=u(k)+rms(k)*e3 s(1)=1.e-10 b0=0. aa=h12*(p(1)-e3) if(s2) 26,27,26 26 b0=-s(1)*aa 27 b1=s(1)*(1.-aa) do 28 k=2,n1 b2=12.*s(k-1)-10.*b1-b0 aa=h12*(p(k)-e3) s(k)=b2/(1.-aa) b0=b1 28 b1=b2 se3=s(n) nel3=0 do 29 k=8,n if(s(k-1)*s(k)) 31,30,29 31 nel3=nel3+2 go to 29 30 nel3=nel3+1 29 continue nel3=nel3/2 if(se1*se3) 32,33,34 32 nel2=nel3 se2=se3 e2=e3 do 35 k=1,n1 35 t(k)=s(k) go to 36 33 e=e3 nel=nel3 do 37 k=1,n1 37 t(k)=s(k) go to 101 34 nel1=nel3 se1=se3 e1=e3 go to 36 101 no=nel-1 do 43 k=n,4,-1 if(t(k-1)*t(k)) 44,45,43 45 kk=k-2 go to 46 44 kk=k-1 go to 46 43 continue 46 if(e) 47,48,48 47 do 53 k=kk,4,-1 if(abs(t(k-1))-abs(t(k))) 54,54,53 53 continue 54 kl=(k+kk)/2 k1=n-1 k2=n-2 s(n)=0. s(k1)=1.e-10 aa=h12*(p(k1)-e) b0=0. b1=s(k1)*(1.-aa) k4=k1-kl do 55 k3=1,k4 k=k2-k3+1 b2=12.*s(k+1)-10.*b1-b0 aa=h12*(p(k)-e) s(k)=b2/(1.-aa) b0=b1 b1=b2 55 continue fac=t(kl)/s(kl) do 56 k=kl,n 56 t(k)=s(k)*fac go to 50 48 do 51 k=kk,n 51 t(k)=t(k)-float(k-kk)*t(n)/float(n-kk) 50 write(1,52),kk 52 format(2x,'last node at kk=',i5/) som=0. do 38 k=1,n y=1.-rms(k) t(k)=t(k)*sqrt(y) 38 som=som+t(k)*t(k) som=sqrt(som*h) sig=t(10)/abs(t(10)) do 39 k=1,n 39 s(k)=sig*t(k)/som return end FUNCTION DSC(Z1) dimension yk(10),pi(10),sig(2) data (yk(i),i=1,10)/0.24534071,0.73747373,1.2340762,1.7385377,2.25 149740,2.7888061,3.3478546,3.9447640,4.6036824,5.3874809/ data(pi(i),i=1,10)/0.49092150,0.49384339,0.49992087,0.50967903, 10.52408035,0.54485174,0.57526244,0.62227870,0.70433296,0.89859196/ data sig/1.,-1./ x=z1 xmu=0.65 zw=1./(0.4431125*xmu**3) sg=0. xsm=x/xmu do 100 i=1,10 do 100 k=1,2 y=x+xmu*yk(i)*sig(k) if(y)100,100,101 101 f=zw*gugus(y)*y*y ysm=y/xmu sg=sg+xmu*pi(i)*f*vkg(0,xsm,ysm) 100 continue dsc=sg return end FUNCTION GUGUS(X) parameter nnn=200,neoc=38 common/der/fct(nnn),df(nnn) common/bpas/h,cof common rm(nnn),du(nnn),vc(nnn),u(nnn),vs(nnn),taun(nnn), +taup(nnn),taut(nnn),djn1(nnn),djp1(nnn),djt1(nnn),vn(nnn),vp(nnn), +vson(nnn),vsop(nnn),ecoulb(nnn),vqre(nnn),rm2(nnn),rm3(nnn), +hmen(nnn),hmep(nnn),dhmen(nnn),dhmep(nnn),d2hmen(nnn),d2hmep(nnn), +dn(nnn),dp(nnn),dt(nnn),dn1(nnn),dp1(nnn),dt1(nnn),dn2(nnn), +dp2(nnn),dt2(nnn),djn(nnn),djp(nnn),djt(nnn),dc(nnn),avp(nnn, +neoc),dunl(nnn,neoc),unl(nnn,neoc) x2=2.*h if(x.gt.x2) go to 2 gugus=((4.*dp(1)-dp(2))+(dp(2)-dp(1))*x*x)/3. return 2 i=x/h if(i.ge.nnn) go to 4 y=(x-i*h)/h s=dp(i)+.5*y*(dp(i+1)-dp(i-1)) gugus=s+.5*y*y*(dp(i+1)+dp(i-1)-2.*dp(i)) return 4 gugus=.0 return end FUNCTION VKG(I1,Z1,Z2) dimension fj(50) l=i1 x=z1 y=z2 delta=y-x a=2.*x*y if(a-15.)101,101,113 101 if(a-0.01)102,102,104 102 arg=exp(-(x*x+y*y)) fj0=arg*(1.+a*a/6.) if(l-1)103,103,105 103 fj(1)=fj0 fj(2)=arg*a*(10.+a*a)/30. go to 107 104 u=delta*delta v=(x+y)*(x+y) u=exp(-u) v=exp(-v) fj0=(u-v)/(2.*a) 105 l2=l+5 l2=l+10 fj(l2)=1.e-10*float(2*l2+1)/a fj(l2+1)=1.e-10 l3=l2-1 do 106 ll=1,l3 l1=l2-ll fj(l1)=float(2*l1+1)*fj(l1+1)/a+fj(l1+2) if(fj(l1)-1.e+30)106,106,111 111 do 112 l4=l1,l2 112 fj(l4)=1.e-10*fj(l4) 106 continue zz=fj0/fj(1) l2=l2-9 do 109 l1=1,l2 109 fj(l1)=zz*fj(l1) go to 107 113 u=delta*delta v=(x+y)*(x+y) u=exp(-u) v=exp(-v) fj(1)=(u-v)/(2.*a) fj(2)=((a-1.)*u+(a+1.)*v)/(2.*a*a) if(l-1)107,107,114 114 do 115 l1=2,l 115 fj(l1+1)=-float(2*l1-1)*fj(l1)/a+fj(l1-1) 107 vkg=float(l+l+1)*fj(l+1) return end FUNCTION HWF(I1,I2,Z1,Z2) n=i1 l=i2 b=z1 r=z2 x=r/b xsq=-2.*x*x y=-x*x/2. b3=b*b*b rho=1. vnl=1. if(n)101,101,102 102 xa=1. k=2*l+1 ns=n+1 ni=0 do 103 i=1,n k=k+2 ni=ni+1 ns=ns-1 xa=xa*xsq*ns/(k*ni) vnl=vnl+xa 103 rho=(k*rho)/(2.*ni) 101 dpl=1. xsl=1. if(l)104,104,105 105 k=1 do 106 i=1,l k=k+2 xk=k dpl=2.*dpl/xk 106 xsl=xsl*x 104 arg=rho*dpl/(b3*1.772454) xa=2.*xsl*r*sqrt(arg) hwf=xa*vnl*exp(y) return end SUBROUTINE LECPOT implicit real*8 (a-h,o-z) parameter nnn=200 common/bpot/xmn(nnn),xmp(nnn),vn(nnn),vp(nnn),yn(nnn),yp(nnn), 1vsn(nnn),vsp(nnn),vc(nnn),wnd(nnn),wpd(nnn) common/bpas/h,cof common/blecp1/t0,t1,t2,t3,t13,alfe,x0,x1,x2,x3,anucl,znucl common/blecp2/dt(nnn),taut(nnn),dt1(nnn),dt2(nnn),dp(nnn),dn(nnn) dimension dal(3),mash(3) 208 format(6e12.5) 209 format(15i3) open(unit=10,status='old',form='formatted') read(10,208) t0,t1,t2,t3,t13,x0 read(10,208) alfe,x3,x1,x2 read(10,208) anucl,znucl,h,(dal(i),i=1,3) read(10,209) (mash(i),i=1,3) read(10,209) i1,i2,i3,i4,i5 read(10,208) (xmn(j),j=1,nnn),(xmp(j),j=1,nnn) read(10,208) ( vn(j),j=1,nnn),( vp(j),j=1,nnn) do 400 j=1,nnn vn(j)=cof*vn(j) 400 vp(j)=cof*vp(j) read(10,208) (yn(j),j=1,nnn),(yp(j),j=1,nnn) read(10,208) (vsn(j),j=1,nnn),(vsp(j),j=1,nnn) read(10,208) (vc(j),j=1,nnn) read(10,208) (dt(j),j=1,nnn),(taut(j),j=1,nnn) read(10,208) (dt1(j),j=1,nnn),(dt2(j),j=1,nnn) read(10,208) (dp(j),j=1,nnn),(dn(j),j=1,nnn) close(10) n2=nnn-2 wnd(1)=(-11.D0*yn(1)+18.D0*yn(2)-9.D0*yn(3)+2.D0*yn(4))/ 1(6.D0*H) wpd(1)=(-11.D0*yp(1)+18.D0*yp(2)-9.D0*yp(3)+2.D0*yp(4))/ 1(6.D0*H) do1 i=2,n2 wnd(i)=(-2.D0*yn(i-1)-3.D0*yn(i)+6.D0*yn(i+1)-yn(i+2))/ 1(6.D0*H) wpd(i)=(-2.D0*yp(i-1)-3.D0*yp(i)+6.D0*yp(i+1)-yp(i+2))/ 1(6.D0*H) 1 continue wnd(nnn-1)=wnd(n2) wnd(nnn)=wnd(n2) wpd(nnn-1)=wpd(n2) wpd(nnn)=wpd(n2) return end SUBROUTINE QWG(B,NS,LS,JS,IQ) implicit real*8 (a-h,o-z) parameter nnn=200,nmx=16 dimension u(nnn,nmx),up(nnn,nmx),d(5,5),vu(nnn),wu(nnn) dimension a(nmx,nmx),v(nmx),t1(nmx),t2(nmx),t3(nmx) dimension e(nmx),tw(nmx),a1(nmx,nmx),vec(nmx),nord(nmx) common/bwg/wg(nnn,nmx),ewg(nmx),xor(nmx),dwg(nnn,nmx) common/bpas/h,cof common/bnri/nri common/hsqr1/t1/hsqr2/t2/hsqr3/t3 common/bpot/xmn(nnn),xmp(nnn),vn(nnn),vp(nnn),yn(nnn),yp(nnn), 1vsn(nnn),vsp(nnn),vc(nnn),wnd(nnn),wpd(nnn) data d(1,1),d(1,2),d(1,3),d(1,4),d(1,5)/-50.D0,96.D0,-72.D0, 132.D0,-6.D0/ data d(2,1),d(2,2),d(2,3),d(2,4),d(2,5)/-6.D0,-20.D0,36.D0, 1-12.D0,2.D0/ data d(3,1),d(3,2),d(3,3),d(3,4),d(3,5)/2.D0,-16.D0,0.D0, 116.D0,-2.D0/ data d(4,1),d(4,2),d(4,3),d(4,4),d(4,5)/-2.D0,12.D0,-36.D0, 120.D0,6.D0/ data d(5,1),d(5,2),d(5,3),d(5,4),d(5,5)/6.D0,-32.D0,72.D0, 1-96.D0,50.D0/ eps=(1.0d-015)*ns*ns nn2=nri-2 l=ls b3=b*b*b*1.772454D0 db3=b3 do 1 n1=1,ns n=n1-1 do 2 j=1,nri r=h*j x=r/b dx=x xsq=-2.D0*x*x y=-x*x/2.D0 dy=y dxsq=xsq dbr=r bh=24.d0*h dbh=bh rho=1.D0 vnl=1.D0 if(n) 101,101,102 102 xa=1.D0 k=2*l+1 nss=n1 ni=0 do 103 ii=1,n k=k+2 ni=ni+1 nss=nss-1 sni=float(nss)/(ni*float(k)) dni=sni xa=xa*dni*dxsq sk=float(k)/(2.d0*ni) dk=sk vnl=vnl+xa 103 rho=dk*rho 101 dpl=1.D0 xsl=1.D0 if(l) 104,104,105 105 k=1 do 106 ii=1,l k=k+2 xk=k dk=xk dpl=2.D0*dpl/dk 106 xsl=xsl*dx 104 arg=rho*dpl/db3 xa=2.d0*xsl*dbr*dsqrt(arg) u(j,n1)=xa*vnl*Dexp(Dy) 2 continue do 3 j=1,nri k=3 if(j.lt.3) k=j if(j.gt.nn2) k=j-nRI+5 sum=0.D0 do 4 ii=1,5 jj=j+ii-k dbd=d(k,ii) 4 sum=sum+dbd*u(jj,n1) 3 up(j,n1)=sum/(DBh) 1 continue xll=l*(l+1.D0) xjj=0.5D0*(0.25D0*js *(js +2.D0)-xll-0.75D0) go to (5,6),iq 5 do 7 j=1,nri x2=h*j*h*j wu(j)=xmn(j) 7 vu(j)=xmn(j)*xll/x2+vn(j)+xjj*vsn(j) go to 8 6 do 9 j=1,nri x2=h*j*h*j wu(j)=xmp(j) 9 vu(j)=xmp(j)*xll/x2+vp(j)+xjj*vsp(j)+vc(j) 8 do 10 i1=1,ns do 11 i2=i1,ns aa=0.D0 do 12 j=1,nri dwu=wu(j) dvu=vu(j) aa=aa+dwu*up(j,i1)*up(j,i2)+dvu*u(j,i1)*u(j,i2) 12 continue aa=aa*h a(i1,i2)=(aa) a(i2,i1)=a(i1,i2) 11 continue 10 continue call diasym(a,nmx,ns,v,e,vec,nord) do 20 n=1,ns xnorm=0.D0 sgn=0.D0 dy=0.d0 ewg(n)=v(n) do 50 ik=1,ns dy=dy+a(ik,n)**2 as=a(ik,n) 50 sgn=sgn+as*u(4,ik) xor(n)=dy sign=1.d0 if(sgn.lt.0.d0) sign=-1.D0 do 21 j=1,nri xx=0.d0 xyx=0.d0 do 22 ii=1,ns xyx=xyx+a(ii,n)*up(j,ii) 22 xx=xx+a(ii,n)*u(J,II) wg(j,n)=xx*sign dwg(j,n)=sign*xyx 21 xnorm=xnorm+(wg(j,n)**2)*h xnorm=dsqrt(xnorm) do 23 j=1,nri dwg(j,n)=dwg(j,n)/xnorm 23 wg(j,n)=wg(j,n)/xnorm 20 continue return end SUBROUTINE DIASYM(A,NP,N,D,E,VEC,NORD) ccc diagonalization of a real, symmetric matrix A IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(NP,NP),D(NP),E(NP),NORD(NP),VEC(NP) CALL TRED2(A,N,NP,D,E) CALL TQLI(D,E,N,NP,A) CALL REORDERING(D,A,N,NP,VEC,NORD) RETURN END SUBROUTINE TRED2(A,N,NP,D,E) ccc This performes a Householder reduction of a real, symmetric, ccc N*N matrix A, stored in an NP*NP physical array. IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(NP,NP),D(NP),E(NP) IF(N.GT.1)THEN DO 18 I=N,2,-1 L=I-1 H=0.D0 SCALE=0.D0 IF(L.GT.1)THEN DO 11 K=1,L SCALE=SCALE+DABS(A(I,K)) 11 CONTINUE IF(SCALE.EQ.0.D0)THEN E(I)=A(I,L) ELSE DO 12 K=1,L A(I,K)=A(I,K)/SCALE H=H+A(I,K)**2 12 CONTINUE F=A(I,L) G=-DSIGN(DSQRT(H),F) E(I)=SCALE*G H=H-F*G A(I,L)=F-G F=0.D0 DO 15 J=1,L A(J,I)=A(I,J)/H G=0.D0 DO 13 K=1,J G=G+A(J,K)*A(I,K) 13 CONTINUE IF(L.GT.J)THEN DO 14 K=J+1,L G=G+A(K,J)*A(I,K) 14 CONTINUE ENDIF E(J)=G/H F=F+E(J)*A(I,J) 15 CONTINUE HH=F/(H+H) DO 17 J=1,L F=A(I,J) G=E(J)-HH*F E(J)=G DO 16 K=1,J A(J,K)=A(J,K)-F*E(K)-G*A(I,K) 16 CONTINUE 17 CONTINUE ENDIF ELSE E(I)=A(I,L) ENDIF D(I)=H 18 CONTINUE ENDIF D(1)=0.D0 E(1)=0.D0 DO 23 I=1,N L=I-1 IF(D(I).NE.0.D0)THEN DO 21 J=1,L G=0.D0 DO 19 K=1,L G=G+A(I,K)*A(K,J) 19 CONTINUE DO 20 K=1,L A(K,J)=A(K,J)-G*A(K,I) 20 CONTINUE 21 CONTINUE ENDIF D(I)=A(I,I) A(I,I)=1.D0 IF(L.GE.1)THEN DO 22 J=1,L A(I,J)=0.D0 A(J,I)=0.D0 22 CONTINUE ENDIF 23 CONTINUE RETURN END SUBROUTINE TQLI(D,E,N,NP,Z) ccc Eigenvalues and eigenvectors of a real symmetric matrix reduced ccc by TRED2. Z is input as the matrix output of TRED2. Its K-th ccc column returns the normalized eigenvector corresponding to ccc D(K). IMPLICIT REAL*8(A-H,O-Z) DIMENSION D(NP),E(NP),Z(NP,NP) IF (N.GT.1) THEN DO 11 I=2,N E(I-1)=E(I) 11 CONTINUE E(N)=0.D0 DO 15 L=1,N ITER=0 1 DO 12 M=L,N-1 DD=DABS(D(M))+DABS(D(M+1)) IF (DABS(E(M))+DD.EQ.DD) GO TO 2 12 CONTINUE M=N 2 IF(M.NE.L)THEN IF(ITER.EQ.30)PAUSE 'too many iterations' ITER=ITER+1 G=(D(L+1)-D(L))/(2.*E(L)) R=DSQRT(G**2+1.D0) G=D(M)-D(L)+E(L)/(G+DSIGN(R,G)) S=1.D0 C=1.D0 P=0.D0 DO 14 I=M-1,L,-1 F=S*E(I) B=C*E(I) IF(DABS(F).GE.DABS(G))THEN C=G/F R=DSQRT(C**2+1.D0) E(I+1)=F*R S=1./R C=C*S ELSE S=F/G R=DSQRT(S**2+1.D0) E(I+1)=G*R C=1./R S=S*C ENDIF G=D(I+1)-P R=(D(I)-G)*S+2.D0*C*B P=S*R D(I+1)=G+P G=C*R-B DO 13 K=1,N F=Z(K,I+1) Z(K,I+1)=S*Z(K,I)+C*F Z(K,I)=C*Z(K,I)-S*F 13 CONTINUE 14 CONTINUE D(L)=D(L)-P E(L)=G E(M)=0.D0 GO TO 1 ENDIF 15 CONTINUE ENDIF RETURN END SUBROUTINE REORDERING(E,A,N,NP,VEC,NORD) IMPLICIT REAL*8(A-H,O-Z) DIMENSION NORD(NP),E(NP),VEC(NP),A(NP,NP) CCC ORDERING OF STATES DO 224 I=1,N 224 NORD(I)=I DO 225 I=1,N K=I X=E(I) DO 226 J=1,N 226 VEC(J)=A(J,I) DO 227 II=I,N IF(X-E(II)) 227,227,228 228 K=II X=E(II) DO 229 J=1,N 229 VEC(J)=A(J,II) 227 CONTINUE KK=NORD(K) NORD(K)=NORD(I) NORD(I)=KK E(K)=E(I) DO 230 J=1,N 230 A(J,K)=A(J,I) E(I)=X DO 231 J=1,N 231 A(J,I)=VEC(J) 225 CONTINUE RETURN END SUBROUTINE PHME(IA,IB,IC,ID,JT,NR,DR,RES) implicit real*8 (a-h,o-z) real*4 sla,slb,slc,sld,sja,sjb,sjc,sjd,sj,c1,c2,c3,c4, 1sl,slk1,slk2,sll,slka,slkb,slkc,slkd parameter nsp=250 common/betalf/calf(7),cbet(7) common/brint/ri0,ris,rx1,rx2(2,2) common/bqwf/nn(nsp),ll(nsp),lj(nsp),lt(nsp),ehf(nsp),xnor(nsp) dimension hh(9),vv(7,3) res=0.d0 do 1 i=1,9 1 hh(i)=0.d0 call rdint1(ia,ib,ic,id,nr,dr) djp1=2.d0*jt+1.d0 la=ll(ia) lb=ll(ib) lc=ll(ic) ld=ll(id) ja=lj(ia) jb=lj(ib) jc=lj(ic) jd=lj(id) sla=float(la) slb=float(lb) slc=float(lc) sld=float(ld) sja=float(ja)/2.d0 sjb=float(jb)/2.d0 sjc=float(jc)/2.d0 sjd=float(jd)/2.d0 sjt=float(jt) hh(8)=ri0*yl(la,ja,ld,jd,jt)*yl(lc,jc,lb,jb,jt)/djp1 ladm=iabs(la-ld) ladp=la+ld lbcm=iabs(lb-lc) lbcp=lb+lc j1m=iabs(jt-1) j1p=jt+1 l1=max0(ladm,lbcm,j1m) l2=min0(ladp,lbcp,j1p) x9=0.d0 if(l1.gt.l2) go to 2 do 3 l=l1,l2 3 x9=x9+tjl(la,ja,ld,jd,jt,l)*tjl(lc,jc,lb,jb,jt,l) 2 hh(9)=ris*x9/djp1 l1=max0(ladm,lbcm) l2=min0(ladp,lbcp) if(l1.gt.l2) go to 4 do 5 i1=1,7 do 6 i2=1,3 vv(i1,i2)=0.d0 6 continue 5 continue cofs=(ja+1.d0)*(jb+1.d0)*(jc+1.d0)*(jd+1.d0) cofs=dsqrt(cofs)*6.d0 is=la+lc+(jb+jd)/2+1 si=1.d0-2.d0*mod(is,2) call sixj(sja,sla,0.5,sld,sjd,sjt,c1) call sixj(sjc,slc,0.5,slb,sjb,sjt,c2) c1=c1*c2 cof0=si*cofs*dble(c1)/6.d0 ilpt=1 ilgd=3 do 7 il=ilpt,ilgd l=jt-2+il if(l.lt.0) go to 7 if((l-l1)*(l-l2).gt.0) go to 7 vv(1,il)=rx1*redy(la,l,ld)*redy(lc,l,lb)/(2.d0*l+1.d0) 7 continue do 8 iup=2,5 go to (9,9,10,11,12), iup 9 ka=ia kb=ib kc=ic kd=id go to 13 10 ka=id kb=ic kc=ib kd=ia go to 13 11 ka=ia kb=ic kc=ib kd=id go to 13 12 ka=ib kb=id kc=ia kd=ic 13 lka=ll(ka) lkb=ll(kb) lkc=ll(kc) lkd=ll(kd) slka=float(lka) slkb=float(lkb) slkc=float(lkc) slkd=float(lkd) call rdint2(ka,kb,kc,kd,nr,dr) xlkab=(2.d0*lka+1.d0)*(2.d0*lkb+1.d0) xlkab=dsqrt(xlkab) do 14 il=ilpt,ilgd l=jt-2+il sl=float(l) if(l.lt.0) go to 14 if((l-l1)*(l-l2).gt.0) go to 14 sgn1=1.d0-2.d0*mod(l,2) if(iup.eq.4.or.iup.eq.5) sgn1=1.d0 do 15 ik1=1,2 k1=lka-3+2*ik1 if(k1.lt.0) go to 15 slk1=float(k1) call troisj(slk1,slka,1.,0.,0.,0.,c1) cofk1=(1.d0-2.d0*mod(k1,2))*dsqrt(2.d0*k1+1.d0)*dble(c1) do 16 ik2=1,2 k2=lkb-3+2*ik2 if(k2.lt.0) go to 16 slk2=float(k2) call troisj(slk2,slkb,1.,0.,0.,0.,c2) cofk2=(1.d0-2.d0*mod(k2,2))*dsqrt(2.d0*k2+1.d0)*dble(c2) ll1m=iabs(l-1) ll2m=iabs(lkd-k1) ll3m=iabs(lkc-k2) ll1p=l+1 ll2p=lkd+k1 ll3p=lkc+k2 llm=max0(ll1m,ll2m,ll3m) llp=min0(ll1p,ll2p,ll3p) if(llm.gt.llp) go to 16 do 17 l1l=llm,llp sll=float(l1l) call sixj(sll,1.,sl,slka,slkd,slk1,c1) call sixj(sll,1.,sl,slkb,slkc,slk2,c2) c3=c1*c2 vv(iup,il)=vv(iup,il)+sgn1*xlkab*cofk1*cofk2*dble(c3)* 1redy(lkd,l1l,k1)*redy(lkc,l1l,k2)*rx2(ik1,ik2) 17 continue 16 continue 15 continue 14 continue 8 continue do 18 iup=6,7 ka=ia kb=ib kc=ic kd=id if(iup.eq.6) go to 19 ka=ib kb=ia kc=id kd=ic 19 lka=ll(ka) lkb=ll(kb) lkc=ll(kc) lkd=ll(kd) slka=float(lka) slkb=float(lkb) slkc=float(lkc) slkd=float(lkd) call rdint2(ka,kd,kb,kc,nr,dr) xlkad=(2.d0*lka+1.d0)*(2.d0*lkd+1.d0) xlkad=dsqrt(xlkad) do 20 il=ilpt,ilgd l=jt-2+il sl=float(l) if(l.lt.0) go to 20 if((l-l1)*(l-l2).gt.0) go to 20 cofl=redy(lkb,l,lkc)/(2.d0*l+1.d0) do 21 ik1=1,2 k1=lka-3+2*ik1 if(k1.lt.0) go to 21 slk1=float(k1) call troisj(slk1,slka,1.,0.,0.,0.,c1) cofk1=dsqrt(2.d0*k1+1.d0)*dble(c1) do 22 ik2=1,2 k2=lkd-3+2*ik2 if(k2.lt.0) go to 22 slk2=float(k2) call troisj(slk2,slkd,1.,0.,0.,0.,c2) cofk2=dsqrt(2.d0*k2+1.d0)*dble(c2) call sixj(slk2,slk1,sl,slka,slkd,1.,c3) cofk12=redy(k1,l,k2)*dble(c3) vv(iup,il)=vv(iup,il)+xlkad*cofl*cofk1*cofk2*cofk12*rx2 1(ik1,ik2) 22 continue 21 continue 20 continue 18 continue do 23 il=ilpt,ilgd l=jt-2+il if(l.lt.0) go to 23 sl=float(l) dlp1=2.d0*l+1.d0 call neufj(sja,sla,0.5,sjd,sld,0.5,sjt,sl,1.,c3) call neufj(sjc,slc,0.5,sjb,slb,0.5,sjt,sl,1.,c4) c4=c3*c4 do 24 iup=1,7 24 hh(iup)=hh(iup)+cofs*dlp1*dble(c4)*vv(iup,il) 23 continue do 25 iup=1,7 hh(iup)=cbet(iup)*hh(iup)+calf(iup)*cof0*vv(iup,2) 25 continue 4 continue do 26 iup=1,9 26 res=res+hh(iup) return end SUBROUTINE RDINT1(IA,IB,IC,ID,NR,DR) implicit real*8 (a-h,o-z) parameter nnn=200,nsp=250 common/brint/ri0,ris,rx1,rx2(2,2) common/bwf/wf(nnn,nsp),dwf(nnn,nsp) common/bqwf/nn(nsp),ll(nsp),lj(nsp),lt(nsp),ehf(nsp),xnor(nsp) common/bfg/f0(nnn),g0(nnn) ri0=0.D0 ris=0.D0 rx1=0.d0 xl=ll(ia)*(ll(ia)+1)+ll(ib)*(ll(ib)+1)+ll(ic)*(ll(ic)+1) xl=xl+ll(id)*(ll(id)+1) do 1 j=1,nr x1=j*dr x2=x1*x1 c=wf(j,ia)*wf(j,ib)*wf(j,ic)*wf(j,id)*dr/x2 ri0=ri0+f0(j)*c ris=ris+g0(j)*c c1=-xl*c/x2+2.d0*dr*((dwf(j,ia)*wf(j,ib)+wf(j,ia)*dwf(j,ib)) 1*wf(j,ic)*wf(j,id)+(dwf(j,ic)*wf(j,id)+wf(j,ic)*dwf(j,id))* 2wf(j,ia)*wf(j,ib))/(x1*x2) c1=c1-2.d0*dr*(dwf(j,ia)*(dwf(j,ib)*wf(j,ic)*wf(j,id)+wf(j, 1ib)*(dwf(j,ic)*wf(j,id)+wf(j,ic)*dwf(j,id)))+wf(j,ia)*(dwf( 2j,ib)*(dwf(j,ic)*wf(j,id)+wf(j,ic)*dwf(j,id))+wf(j,ib)*dwf( 3j,ic)*dwf(j,id)))/x2 rx1=rx1+c1 1 continue return end SUBROUTINE SIXJ(XJ1,XJ2,XJ3,XL1,XL2,XL3,C6J) DIMENSION FLOG(301) DATA(FLOG(I),I=2,31)/0.,.69314718,1.7917595,3.1780538,4.7874917,6. 15792511,8.5251613,10.604603,12.801827,15.104413,17.502307,19.98721 24,22.552163,25.191221,27.899271,30.671860,33.505072,36.395445,39.3 339884,42.335616,45.380139,48.471180,51.606674,54.784729,58.003604, 461.261702,64.557537,67.889743,71.257038,74.658235/ DATA(FLOG(I),I=32,61)/78.092223,81.557959,85.054466,88.580827,92.1 136175,95.719694,99.330612,102.96820,106.63176,110.32064,114.03421, 2117.77188,121.53308,125.31727,129.12393,132.95257,136.80272,140.67 3392,144.56574,148.47776,152.40959,156.36083,160.33112,164.32011,16 48.32744,172.35279,176.39584,180.45629,184.53383,188.62817/ DATA(FLOG(I),I=62,91)/192.73904,196.86618,201.00931,205.16820,209. 134258,213.53224,217.73693,221.95644,226.19054,230.43904,234.70172, 2238.97839,243.26885,247.57291,251.89040,256.22113,260.56494,264.92 3164,269.29110,273.67312,278.06757,282.47429,286.89313,291.32394,29 45.76659,300.22094,304.68685,309.16419,313.65283,318.15264/ DATA(FLOG(I),I=92,121)/322.66349,327.18529,331.71788,336.26118,340 1.81505,345.37940,349.95411,354.53908,359.13420,363.73937,368.35449 2,372.97946,377.61419,382.25859,386.91255,391.57598,396.24881,400.9 33094,405.62230,410.32277,415.03230,419.75080,424.47819,429.21439,4 433.95932,438.71291,443.47508,448.24576,453.02489,457.81238/ DATA(FLOG(I),I=122,151)/462.60817,467.41220,472.22438,477.04466,48 11.87298,486.70926,491.55345,496.40547,501.26529,506.13282,511.0080 22,515.89082,520.78117,525.67901,530.58428,535.49694,540.41692,545. 334417,550.27865,555.22029,560.16905,565.12488,570.08772,575.05753, 4580.03427,585.01787,590.00830,595.00552,600.00946,605.02010/ DATA(FLOG(I),I=152,181)/610.03738,615.06126,620.09170,625.12866,63 10.17208,635.22193,640.27818,645.34077,650.40968,655.48486,660.5662 26,665.65385,670.74760,675.84747,680.95341,686.06541,691.18340,696. 330735,701.43726,706.57306,711.71472,716.86221,722.01551,727.17456, 4732.33934,737.50983,742.68598,747.86776,753.05516,758.24811/ DATA(FLOG(I),I=182,211)/763.44661,768.65061,773.86010,779.07503,78 14.29539,789.52114,794.75224,799.98869,805.23044,810.47747,815.7297 23,820.98722,826.24991,831.51778,836.79078,842.06890,847.35209,852. 364036,857.93366,863.23199,868.53529,873.84356,879.15676,884.47488, 4889.79789,895.12577,900.45848,905.79603,911.13836,916.48547/ DATA(FLOG(I),I=212,241)/921.83732,927.19391,932.55521,937.92118,94 13.29181,948.66710,954.04699,959.43148,964.82056,970.21419,975.6123 25,981.01503,986.42220,991.83385,997.24995,1002.6705,1008.0954,1013 3.5248,1018.9585,1024.3966,1029.8389,1035.2857,1040.7367,1046.1920, 41051.6516,1057.1155,1062.5836,1068.0558,1073.5323,1079.0129/ DATA(FLOG(I),I=242,271)/1084.4977,1089.9866,1095.4797,1100.9768,11 106.4781,1111.9834,1117.4928,1123.0063,1128.5237,1134.0452,1139.570 26,1145.1001,1150.6335,1156.1708,1161.7120,1167.2573,1172.8063,1178 3.3593,1183.9161,1189.4768,1195.0413,1200.6097,1206.1818,1211.7577, 41217.3375,1222.9209,1228.5082,1234.0992,1239.6939,1245.2924/ DATA(FLOG(I),I=272,301)/1250.8944,1256.5003,1262.1097,1267.7228,12 173.3396,1278.9600,1284.5840,1290.2117,1295.8429,1301.4777,1307.116 20,1312.7580,1318.4034,1324.0524,1329.7048,1335.3609,1341.0203,1346 3.6833,1352.3497,1358.0196,1363.6929,1369.3697,1375.0499,1380.7334, 41386.4204,1392.1107,1397.8045,1403.5016,1409.2020,1414.9058/ DATA EPS1,EPS2/.1,-.2/ XN=-XJ1+XJ2+XJ3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N1=XN XN=-XL1+XL2+XJ3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N2=XN XN=-XL1+XJ2+XL3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N3=XN XN=-XJ1+XL2+XL3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N4=XN XN=XJ1-XJ2+XJ3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N5=XN XN=XL1-XL2+XJ3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N6=XN XN=XL1-XJ2+XL3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N7=XN XN=XJ1-XL2+XL3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N8=XN XN=XJ1+XJ2-XJ3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N9=XN XN=XL1+XL2-XJ3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N10=XN XN=XL1+XJ2-XL3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N11=XN XN=XJ1+XL2-XL3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N12=XN XN=-XJ1-XL1+XJ3+XL3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N13=XN XN=-XJ2-XL2+XJ3+XL3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N14=XN XN=XJ1+XL1+XJ2+XL2+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N15=XN N15=N15+1 XN=XJ1+XJ2+XJ3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N16=XN N16=N16+1 XN=XL1+XL2+XJ3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N17=XN N17=N17+1 XN=XL1+XJ2+XL3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N18=XN N18=N18+1 XN=XJ1+XL2+XL3+EPS1 IF (XN.LT.0.) XN=XN+EPS2 N19=XN N19=N19+1 IF(N9.LT.0) GO TO 50 IF(N5.LT.0) GO TO 50 IF(N1.LT.0) GO TO 50 IF(N10.LT.0) GO TO 50 IF(N6.LT.0) GO TO 50 IF(N2.LT.0) GO TO 50 IF(N11.LT.0) GO TO 50 IF(N7.LT.0) GO TO 50 IF(N3.LT.0) GO TO 50 IF(N12.LT.0) GO TO 50 IF(N8.LT.0) GO TO 50 IF(N4.LT.0) GO TO 50 K=N17+N18+N19-3 IF(K.GT.0) GO TO 54 C6J=1. RETURN 50 C6J=0. RETURN 54 K=0 L=-N13 IF(L.GT.K) K=L L=-N14 IF(L.GT.K) K=L L=N9 IF(N10.LT.L) L=N10 IF(N11.LT.L) L=N11 IF(N12.LT.L) L=N12 IF(N15.LT.L) L=N15 F=1. S=1. I=K+1 62 IF(I.GT.L) GO TO 80 IM1=I-1 NN=(N9-IM1)*(N10-IM1)*(N11-IM1)*(N12-IM1) ND=I*(N13+I)*(N14+I)*(N15-IM1) F=-F*FLOAT(NN)/FLOAT(ND) S=S+F I=I+1 GO TO 62 80 C2N=FLOG(N1+1)+FLOG(N2+1)+FLOG(N3+1)+FLOG(N4+1)+FLOG(N5+1)+FLOG(N6 1+1)+FLOG(N7+1)+FLOG(N8+1)+FLOG(N9+1)+FLOG(N10+1)+FLOG(N11+1)+FLOG( 2N12+1) C2N=.5*C2N C2D=FLOG(N16+1)+FLOG(N17+1)+FLOG(N18+1)+FLOG(N19+1) C2D=.5*C2D KM1=K-1 KP1=K+1 C2N=C2N+FLOG(N15-KM1) C2D=C2D+FLOG(KP1)+FLOG(N13+KP1)+FLOG(N14+KP1)+FLOG(N9-KM1)+FLOG(N1 10-KM1)+FLOG(N11-KM1)+FLOG(N12-KM1) F=C2D-C2N IF(F.GT.80.) GO TO 98 F=C2N/C2D IF((F.LT.1.01).AND.(F.GT.0.98)) GO TO 98 C6J=S*EXP(C2N-C2D) GO TO 106 98 IF(S) 100,50,102 100 S=ALOG(-S) C6J=-EXP(S+C2N-C2D) GO TO 106 102 S=ALOG(S) C6J=EXP(S+C2N-C2D) 106 L=N15+KM1 K=L/2 K=2*K IF(L.NE.K) C6J=-C6J RETURN END SUBROUTINE RDINT2(IA,IB,IC,ID,NR,DR) implicit real*8 (a-h,o-z) parameter nnn=200,nsp=250 dimension ua(2),ub(2) common/brint/ri0,ris,rx1,rx2(2,2) common/bwf/wf(nnn,nsp),dwf(nnn,nsp) common/bqwf/nn(nsp),ll(nsp),lj(nsp),lt(nsp),ehf(nsp),xnor(nsp) common/bfg/f0(nnn),g0(nnn) do 2 k1=1,2 do 6 k2=1,2 rx2(k1,k2)=0.d0 6 continue 2 continue do 1 j=1,nr x1=j*dr x2=x1*x1 do 3 k=1,2 sk=3.d0-2.d0*k ua(k)=dwf(j,ia)+sk*(ll(ia)+k-1)*wf(j,ia)/x1 ub(k)=dwf(j,ib)+sk*(ll(ib)+k-1)*wf(j,ib)/x1 3 continue do 4 k1=1,2 do 5 k2=1,2 c23=ua(k1)*ub(k2)*wf(j,ic)*wf(j,id) rx2(k1,k2)=rx2(k1,k2)+c23*dr/x2 5 continue 4 continue 1 continue return end SUBROUTINE TROISJ(XJ1,XJ2,XJ3,XM1,XM2,XM3,C3J) DIMENSION FLOG(301) DATA(FLOG(I),I=2,31)/0.,.69314718,1.7917595,3.1780538,4.7874917,6. 15792511,8.5251613,10.604603,12.801827,15.104413,17.502307,19.98721 24,22.552163,25.191221,27.899271,30.671860,33.505072,36.395445,39.3 339884,42.335616,45.380139,48.471180,51.606674,54.784729,58.003604, 461.261702,64.557537,67.889743,71.257038,74.658235/ DATA(FLOG(I),I=32,61)/78.092223,81.557959,85.054466,88.580827,92.1 136175,95.719694,99.330612,102.96820,106.63176,110.32064,114.03421, 2117.77188,121.53308,125.31727,129.12393,132.95257,136.80272,140.67 3392,144.56574,148.47776,152.40959,156.36083,160.33112,164.32011,16 48.32744,172.35279,176.39584,180.45629,184.53383,188.62817/ DATA(FLOG(I),I=62,91)/192.73904,196.86618,201.00931,205.16820,209. 134258,213.53224,217.73693,221.95644,226.19054,230.43904,234.70172, 2238.97839,243.26885,247.57291,251.89040,256.22113,260.56494,264.92 3164,269.29110,273.67312,278.06757,282.47429,286.89313,291.32394,29 45.76659,300.22094,304.68685,309.16419,313.65283,318.15264/ DATA(FLOG(I),I=92,121)/322.66349,327.18529,331.71788,336.26118,340 1.81505,345.37940,349.95411,354.53908,359.13420,363.73937,368.35449 2,372.97946,377.61419,382.25859,386.91255,391.57598,396.24881,400.9 33094,405.62230,410.32277,415.03230,419.75080,424.47819,429.21439,4 433.95932,438.71291,443.47508,448.24576,453.02489,457.81238/ DATA(FLOG(I),I=122,151)/462.60817,467.41220,472.22438,477.04466,48 11.87298,486.70926,491.55345,496.40547,501.26529,506.13282,511.0080 22,515.89082,520.78117,525.67901,530.58428,535.49694,540.41692,545. 334417,550.27865,555.22029,560.16905,565.12488,570.08772,575.05753, 4580.03427,585.01787,590.00830,595.00552,600.00946,605.02010/ DATA(FLOG(I),I=152,181)/610.03738,615.06126,620.09170,625.12866,63 10.17208,635.22193,640.27818,645.34077,650.40968,655.48486,660.5662 26,665.65385,670.74760,675.84747,680.95341,686.06541,691.18340,696. 330735,701.43726,706.57306,711.71472,716.86221,722.01551,727.17456, 4732.33934,737.50983,742.68598,747.86776,753.05516,758.24811/ DATA(FLOG(I),I=182,211)/763.44661,768.65061,773.86010,779.07503,78 14.29539,789.52114,794.75224,799.98869,805.23044,810.47747,815.7297 23,820.98722,826.24991,831.51778,836.79078,842.06890,847.35209,852. 364036,857.93366,863.23199,868.53529,873.84356,879.15676,884.47488, 4889.79789,895.12577,900.45848,905.79603,911.13836,916.48547/ DATA(FLOG(I),I=212,241)/921.83732,927.19391,932.55521,937.92118,94 13.29181,948.66710,954.04699,959.43148,964.82056,970.21419,975.6123 25,981.01503,986.42220,991.83385,997.24995,1002.6705,1008.0954,1013 3.5248,1018.9585,1024.3966,1029.8389,1035.2857,1040.7367,1046.1920, 41051.6516,1057.1155,1062.5836,1068.0558,1073.5323,1079.0129/ DATA(FLOG(I),I=242,271)/1084.4977,1089.9866,1095.4797,1100.9768,11 106.4781,1111.9834,1117.4928,1123.0063,1128.5237,1134.0452,1139.570 26,1145.1001,1150.6335,1156.1708,1161.7120,1167.2573,1172.8063,1178 3.3593,1183.9161,1189.4768,1195.0413,1200.6097,1206.1818,1211.7577, 41217.3375,1222.9209,1228.5082,1234.0992,1239.6939,1245.2924/ DATA(FLOG(I),I=272,301)/1250.8944,1256.5003,1262.1097,1267.7228,12 173.3396,1278.9600,1284.5840,1290.2117,1295.8429,1301.4777,1307.116 20,1312.7580,1318.4034,1324.0524,1329.7048,1335.3609,1341.0203,1346 3.6833,1352.3497,1358.0196,1363.6929,1369.3697,1375.0499,1380.7334, 41386.4204,1392.1107,1397.8045,1403.5016,1409.2020,1414.9058/ DATA EPS1,EPS2/.1,-.2/ XN=XJ2-XM2+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N1=XN XN=XJ3+XM3+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N2=XN IF(XN.LT.0.) XN=XN+EPS2 XN=XJ3-XM3+EPS1 N3=XN XN=XJ1+XM1+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N4=XN XN=XJ2+XJ3-XJ1+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N5=XN XN=XJ1+XJ3-XJ2+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N6=XN XN=XJ1+XJ2-XJ3+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N7=XN XN=XJ1-XM1+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N8=XN XN=XJ2+XM2+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N9=XN XN=XJ3-XJ1-XM2+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N10=XN XN=XJ3-XJ2+XM1+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N11=XN XN=XJ1+XJ2+XJ3+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N12=XN N12=N12+1 K=N4*N8 IF(K.LT.0) GO TO 59 K=N1*N9 IF(K.LT.0) GO TO 59 K=N2*N3 IF(K.LT.0) GO TO 59 IF(N5.LT.0) GO TO 59 IF(N6.LT.0) GO TO 59 IF(N7.LT.0) GO TO 59 L=N1-N2+N3-N4+N8-N9 IF(L.NE.0) GO TO 59 K=N12-1 IF(K.GT.0) GO TO 60 C3J=1. RETURN 59 C3J=0. RETURN 60 K=0 L=-N10 IF(L.GT.K) K=L L=-N11 IF(L.GT.K) K=L L=N7 IF(N8.GT.L) L=N8 IF(N9.GT.L) L=N9 F=1. S=1. I=K+1 62 IF(I.GT.L) GO TO 80 IM1=I-1 NN=(N7-IM1)*(N8-IM1)*(N9-IM1) ND=I*(N10+I)*(N11+I) F=-F*FLOAT(NN)/FLOAT(ND) S=S+F I=I+1 GO TO 62 80 C2N=FLOG(N1+1)+FLOG(N2+1)+FLOG(N3+1)+FLOG(N4+1)+FLOG(N5+1)+FLOG(N6 1+1)+FLOG(N7+1)+FLOG(N8+1)+FLOG(N9+1) C2N=.5*C2N KM1=K-1 KP1=K+1 C2D=FLOG(KP1)+FLOG(N7-KM1)+FLOG(N8-KM1)+FLOG(N9-KM1)+FLOG(N10+KP1) 1+FLOG(N11+KP1)+.5*FLOG(N12+1) F=C2D-C2N IF(F.GT.80.) GO TO 98 F=C2N/C2D IF((F.LT.1.01).AND.(F.GT.0.98)) GO TO 98 C3J=S*EXP(C2N-C2D) GO TO 106 98 IF(S) 100,59,102 100 S=ALOG(-S) C3J=-EXP(S+C2N-C2D) GO TO 106 102 S=ALOG(S) C3J=EXP(S+C2N-C2D) 106 XN=XJ1-XJ2-XM3+EPS1 IF(XN.LT.0.) XN=XN+EPS2 L=XN L=L+K K=L/2 K=2*K IF(L.NE.K) C3J=-C3J RETURN END SUBROUTINE NEUFJ(XJ11,XJ12,XJ13,XJ21,XJ22,XJ23,XJ31,XJ32,XJ33,C9J) DIMENSION FLOG(301) DATA(FLOG(I),I=2,31)/0.,.69314718,1.7917595,3.1780538,4.7874917,6. 15792511,8.5251613,10.604603,12.801827,15.104413,17.502307,19.98721 24,22.552163,25.191221,27.899271,30.671860,33.505072,36.395445,39.3 339884,42.335616,45.380139,48.471180,51.606674,54.784729,58.003604, 461.261702,64.557537,67.889743,71.257038,74.658235/ DATA(FLOG(I),I=32,61)/78.092223,81.557959,85.054466,88.580827,92.1 136175,95.719694,99.330612,102.96820,106.63176,110.32064,114.03421, 2117.77188,121.53308,125.31727,129.12393,132.95257,136.80272,140.67 3392,144.56574,148.47776,152.40959,156.36083,160.33112,164.32011,16 48.32744,172.35279,176.39584,180.45629,184.53383,188.62817/ DATA(FLOG(I),I=62,91)/192.73904,196.86618,201.00931,205.16820,209. 134258,213.53224,217.73693,221.95644,226.19054,230.43904,234.70172, 2238.97839,243.26885,247.57291,251.89040,256.22113,260.56494,264.92 3164,269.29110,273.67312,278.06757,282.47429,286.89313,291.32394,29 45.76659,300.22094,304.68685,309.16419,313.65283,318.15264/ DATA(FLOG(I),I=92,121)/322.66349,327.18529,331.71788,336.26118,340 1.81505,345.37940,349.95411,354.53908,359.13420,363.73937,368.35449 2,372.97946,377.61419,382.25859,386.91255,391.57598,396.24881,400.9 33094,405.62230,410.32277,415.03230,419.75080,424.47819,429.21439,4 433.95932,438.71291,443.47508,448.24576,453.02489,457.81238/ DATA(FLOG(I),I=122,151)/462.60817,467.41220,472.22438,477.04466,48 11.87298,486.70926,491.55345,496.40547,501.26529,506.13282,511.0080 22,515.89082,520.78117,525.67901,530.58428,535.49694,540.41692,545. 334417,550.27865,555.22029,560.16905,565.12488,570.08772,575.05753, 4580.03427,585.01787,590.00830,595.00552,600.00946,605.02010/ DATA(FLOG(I),I=152,181)/610.03738,615.06126,620.09170,625.12866,63 10.17208,635.22193,640.27818,645.34077,650.40968,655.48486,660.5662 26,665.65385,670.74760,675.84747,680.95341,686.06541,691.18340,696. 330735,701.43726,706.57306,711.71472,716.86221,722.01551,727.17456, 4732.33934,737.50983,742.68598,747.86776,753.05516,758.24811/ DATA(FLOG(I),I=182,211)/763.44661,768.65061,773.86010,779.07503,78 14.29539,789.52114,794.75224,799.98869,805.23044,810.47747,815.7297 23,820.98722,826.24991,831.51778,836.79078,842.06890,847.35209,852. 364036,857.93366,863.23199,868.53529,873.84356,879.15676,884.47488, 4889.79789,895.12577,900.45848,905.79603,911.13836,916.48547/ DATA(FLOG(I),I=212,241)/921.83732,927.19391,932.55521,937.92118,94 13.29181,948.66710,954.04699,959.43148,964.82056,970.21419,975.6123 25,981.01503,986.42220,991.83385,997.24995,1002.6705,1008.0954,1013 3.5248,1018.9585,1024.3966,1029.8389,1035.2857,1040.7367,1046.1920, 41051.6516,1057.1155,1062.5836,1068.0558,1073.5323,1079.0129/ DATA(FLOG(I),I=242,271)/1084.4977,1089.9866,1095.4797,1100.9768,11 106.4781,1111.9834,1117.4928,1123.0063,1128.5237,1134.0452,1139.570 26,1145.1001,1150.6335,1156.1708,1161.7120,1167.2573,1172.8063,1178 3.3593,1183.9161,1189.4768,1195.0413,1200.6097,1206.1818,1211.7577, 41217.3375,1222.9209,1228.5082,1234.0992,1239.6939,1245.2924/ DATA(FLOG(I),I=272,301)/1250.8944,1256.5003,1262.1097,1267.7228,12 173.3396,1278.9600,1284.5840,1290.2117,1295.8429,1301.4777,1307.116 20,1312.7580,1318.4034,1324.0524,1329.7048,1335.3609,1341.0203,1346 3.6833,1352.3497,1358.0196,1363.6929,1369.3697,1375.0499,1380.7334, 41386.4204,1392.1107,1397.8045,1403.5016,1409.2020,1414.9058/ DATA EPS1,EPS2/.1,-.2/ XN=-XJ11+XJ21+XJ31+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N1=XN XN=-XJ32+XJ33+XJ31+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N2=XN XN=XJ11-XJ21+XJ31+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N5=XN XN=XJ32-XJ33+XJ31+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N6=XN XN=XJ11+XJ21-XJ31+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N9=XN XN=XJ32+XJ33-XJ31+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N10=XN XN=XJ11+XJ32+XJ21+XJ33+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N15=XN N15=N15+1 XN=XJ11+XJ21+XJ31+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N16=XN N16=N16+1 XN=XJ32+XJ33+XJ31+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N17=XN N17=N17+1 XN=-XJ12+XJ22+XJ32+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N21=XN XN=-XJ21+XJ22+XJ23+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N23=XN XN=XJ12-XJ22+XJ32+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N25=XN XN=XJ21-XJ22+XJ23+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N27=XN XN=XJ12+XJ22-XJ32+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N29=XN XN=XJ21+XJ22-XJ23+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N31=XN XN=-XJ12-XJ21+XJ32+XJ23+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N33=XN XN=XJ12+XJ22+XJ32+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N36=XN N36=N36+1 XN=XJ21+XJ22+XJ23+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N38=XN N38=N38+1 XN=-XJ13+XJ23+XJ33+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N41=XN XN=-XJ13+XJ11+XJ12+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N44=XN XN=XJ13-XJ23+XJ33+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N45=XN XN=XJ13-XJ11+XJ12+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N48=XN XN=XJ13+XJ23-XJ33+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N49=XN XN=XJ13+XJ11-XJ12+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N52=XN XN=-XJ23-XJ11+XJ33+XJ12+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N54=XN XN=XJ13+XJ23+XJ33+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N56=XN N56=N56+1 XN=XJ13+XJ11+XJ12+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N59=XN N59=N59+1 IF(N9.LT.0) GO TO 50 IF(N5.LT.0) GO TO 50 IF(N1.LT.0) GO TO 50 IF(N10.LT.0) GO TO 50 IF(N6.LT.0) GO TO 50 IF(N2.LT.0) GO TO 50 IF(N29.LT.0) GO TO 50 IF(N25.LT.0) GO TO 50 IF(N21.LT.0) GO TO 50 IF(N31.LT.0) GO TO 50 IF(N27.LT.0) GO TO 50 IF(N23.LT.0) GO TO 50 IF(N49.LT.0) GO TO 50 IF(N45.LT.0) GO TO 50 IF(N41.LT.0) GO TO 50 IF(N52.LT.0) GO TO 50 IF(N48.LT.0) GO TO 50 IF(N44.LT.0) GO TO 50 K=N1+N2+N5+N6+N9+N10+N21+N23+N25+N27+N29+N31+N41+N44+N45+N48+N49+N 152 IF(K.GT.0) GO TO 54 C9J=1. RETURN 50 C9J=0. RETURN 54 XN=2.*(XJ21-XJ32)+EPS1 IF(XN.LT.0.) XN=-(XN+EPS2) JMIN=XN XN=2.*(XJ11-XJ33)+EPS1 IF(XN.LT.0.) XN=-(XN+EPS2) N=XN IF(N.GT.JMIN) JMIN=N XN=2.*(XJ12-XJ23)+EPS1 IF(XN.LT.0.) XN=-(XN+EPS2) N=XN IF(N.GT.JMIN) JMIN=N XN=2.*(XJ21+XJ32)+EPS1 JMAX=XN XN=2.*(XJ11+XJ33)+EPS1 N=XN IF(N.LT.JMAX) JMAX=N XN=2.*(XJ12+XJ23)+EPS1 N=XN IF(N.LT.JMAX) JMAX=N XJMIN=FLOAT(JMIN)/2. XJMAX=FLOAT(JMAX)/2. S=0. XJ=XJMIN XN=-XJ32+XJ21+XJ+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N3=XN XN=-XJ11+XJ33+XJ+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N4=XN XN=XJ32-XJ21+XJ+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N7=XN XN=XJ11-XJ33+XJ+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N8=XN XN=XJ32+XJ21-XJ+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N11=XN XN=XJ11+XJ33-XJ+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N12=XN XN=-XJ11-XJ32+XJ31+XJ+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N13=XN XN=-XJ21-XJ33+XJ31+XJ+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N14=XN XN=XJ32+XJ21+XJ+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N18=XN N18=N18+1 XN=XJ11+XJ33+XJ+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N19=XN N19=N19+1 XN=-XJ21+XJ+XJ32+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N22=XN XN=-XJ12+XJ23+XJ+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N24=XN XN=XJ21-XJ+XJ32+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N26=XN XN=XJ12-XJ+XJ23+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N28=XN XN=XJ21+XJ-XJ32+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N30=XN XN=XJ12+XJ-XJ23+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N32=XN XN=-XJ22-XJ+XJ32+XJ23+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N34=XN XN=XJ12+XJ21+XJ22+XJ+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N35=XN N35=N35+1 XN=XJ21+XJ+XJ32+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N37=XN N37=N37+1 XN=XJ12+XJ+XJ23+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N39=XN N39=N39+1 XN=-XJ+XJ11+XJ33+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N42=XN XN=-XJ+XJ23+XJ12+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N43=XN XN=XJ-XJ11+XJ33+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N46=XN XN=XJ-XJ23+XJ12+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N47=XN XN=XJ+XJ11-XJ33+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N50=XN XN=XJ+XJ23-XJ12+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N51=XN XN=-XJ13-XJ+XJ33+XJ12+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N53=XN XN=XJ13+XJ+XJ23+XJ11+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N55=XN N55=N55+1 XN=XJ+XJ11+XJ33+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N57=XN N57=N57+1 XN=XJ+XJ23+XJ12+EPS1 IF(XN.LT.0.) XN=XN+EPS2 N58=XN N58=N58+1 GO TO 10 52 IF(XJ.GT.XJMAX) GO TO 120 N3=N3+1 N4=N4+1 N7=N7+1 N8=N8+1 N11=N11-1 N12=N12-1 N13=N13+1 N14=N14+1 N18=N18+1 N19=N19+1 N22=N22+1 N24=N24+1 N26=N26-1 N28=N28-1 N30=N30+1 N32=N32+1 N34=N34-1 N35=N35+1 N37=N37+1 N39=N39+1 N42=N42-1 N43=N43-1 N46=N46+1 N47=N47+1 N50=N50+1 N51=N51+1 N53=N53-1 N55=N55+1 N57=N57+1 N58=N58+1 10 K1=0. L1=-N13 IF(L1.GT.K1) K1=L1 L1=-N14 IF(L1.GT.K1) K1=L1 L1=N9 IF(N10.LT.L1) L1=N10 IF(N11.LT.L1) L1=N11 IF(N12.LT.L1) L1=N12 IF(N15.LT.L1) L1=N15 F1=1. S1=1. I1=K1+1 62 IF(I1.GT.L1) GO TO 64 I1M1=I1-1 NN1=(N9-I1M1)*(N10-I1M1)*(N11-I1M1)*(N12-I1M1) ND1=I1*(N13+I1)*(N14+I1)*(N15-I1M1) F1=-F1*FLOAT(NN1)/FLOAT(ND1) S1=S1+F1 I1=I1+1 GO TO 62 64 K2=0 L2=-N33 IF(L2.GT.K2) K2=L2 L2=-N34 IF(L2.GT.K2) K2=L2 L2=N29 IF(N30.LT.L2) L2=N30 IF(N31.LT.L2) L2=N31 IF(N32.LT.L2) L2=N32 IF(N35.LT.L2) L2=N35 F2=1. S2=1. I2=K2+1 70 IF(I2.GT.L2) GO TO 80 I2M2=I2-1 NN2=(N29-I2M2)*(N30-I2M2)*(N31-I2M2)*(N32-I2M2) ND2=I2*(N33+I2)*(N34+I2)*(N35-I2M2) F2=-F2*FLOAT(NN2)/FLOAT(ND2) S2=S2+F2 I2=I2+1 GO TO 70 80 K3=0 L3=-N53 IF(L3.GT.K3) K3=L3 L3=-N54 IF(L3.GT.K3) K3=L3 L3=N49 IF(N50.LT.L3) L3=N50 IF(N51.LT.L3) L3=N51 IF(N52.LT.L3) L3=N52 IF(N55.LT.L3) L3=N55 F3=1. S3=1. I3=K3+1 84 IF(I3.GT.L3) GO TO 90 I3M3=I3-1 NN3=(N49-I3M3)*(N50-I3M3)*(N51-I3M3)*(N52-I3M3) ND3=I3*(N53+I3)*(N54+I3)*(N55-I3M3) F3=-F3*FLOAT(NN3)/FLOAT(ND3) S3=S3+F3 I3=I3+1 GO TO 84 90 S2N=FLOG(N3+1)+FLOG(N4+1)+FLOG(N7+1)+FLOG(N8+1)+FLOG(N11+1)+FLOG(N 112+1)+FLOG(N22+1)+FLOG(N24+1)+FLOG(N26+1)+FLOG(N28+1)+FLOG(N30+1)+ 2FLOG(N32+1)+FLOG(N42+1)+FLOG(N43+1)+FLOG(N46+1)+FLOG(N47+1)+FLOG(N 350+1)+FLOG(N51+1) S2N=.5*S2N S2D=FLOG(N18+1)+FLOG(N19+1)+FLOG(N37+1)+FLOG(N39+1)+FLOG(N57+1)+FL 1OG(N58+1) S2D=.5*S2D KM1=K1-1 KP1=K1+1 KM2=K2-1 KP2=K2+1 KM3=K3-1 KP3=K3+1 S2N=S2N+FLOG(N15-KM1)+FLOG(N35-KM2)+FLOG(N55-KM3) S2D=S2D+FLOG(KP1)+FLOG(KP2)+FLOG(KP3)+FLOG(N9-KM1)+FLOG(N10-KM1)+F 1LOG(N11-KM1)+FLOG(N12-KM1)+FLOG(N13+KP1)+FLOG(N14+KP1)+FLOG(N29-KM 22)+FLOG(N30-KM2)+FLOG(N31-KM2)+FLOG(N32-KM2)+FLOG(N33+KP2)+FLOG(N3 34+KP2)+FLOG(N49-KM3)+FLOG(N50-KM3)+FLOG(N51-KM3) 4+FLOG(N52-KM3)+FLOG(N53+KP3)+FLOG(N54+KP3) F=S2D-S2N IF (F.GT.80) GO TO 100 F=S2N/S2D IF((F.LT.1.01).AND.(F.GT.0.98)) GO TO 100 F=S1*S2*S3*EXP(S2N-S2D) GO TO 110 100 F=S1*S2*S3 IF(F) 102,112,104 102 F=ALOG(-F) F=-EXP(F+S2N-S2D) GO TO 110 104 F=ALOG(F) F=EXP(F+S2N-S2D) 110 F=F*(2.*XJ+1.) L1=K1+K2+K3 K1=L1/2 K1=2*K1 IF(L1.NE.K1) F=-F S=S+F 112 XJ=XJ+1. GO TO 52 120 C2N=FLOG(N1+1)+FLOG(N2+1)+FLOG(N5+1)+FLOG(N6+1)+FLOG(N9+1)+FLOG(N1 10+1)+FLOG(N21+1)+FLOG(N23+1)+FLOG(N25+1)+FLOG(N27+1)+FLOG(N29+1)+F 2LOG(N31+1)+FLOG(N41+1)+FLOG(N44+1)+FLOG(N45+1)+FLOG(N48+1)+FLOG(N4 39+1)+FLOG(N52+1) C2N=.5*C2N C2D=FLOG(N16+1)+FLOG(N17+1)+FLOG(N36+1)+FLOG(N38+1)+FLOG(N56+1)+FL 1OG(N59+1) C2D=.5*C2D F=C2D-C2N IF(F.GT.80.) GO TO 122 F=C2N/C2D IF((F.LT.1.01).AND.(F.GT.0.98)) GO TO 122 C9J=S*EXP(C2N-C2D) GO TO 130 122 IF(S) 124,50,126 124 S=ALOG(-S) C9J=-EXP(S+C2N-C2D) GO TO 130 126 S=ALOG(S) C9J=EXP(S+C2N-C2D) 130 K=N9+N16+N36+N56-1 L=K/2 L=2*L IF (L.NE.K) C9J=-C9J RETURN END SUBROUTINE DIAGMA(NCONF,IRPA) ccc method of Ullah-Rowe, Nucl.Phys. 163,p.257 implicit real*8 (a-h,o-z) parameter ncf=380,ncf2=2*ncf common/bdiam/aa(ncf,ncf),bb(ncf,ncf),evr(ncf),vecr(ncf2,ncf) common/hsqr1/t1/hsqr2/t2/hsqr3/t3 dimension x(ncf,ncf),v(ncf),t1(ncf),t2(ncf),t3(ncf),e(ncf) dimension cec(ncf),tep(ncf,ncf),alp(ncf),valp(ncf),vec(ncf) dimension nord(ncf) ndim=nconf n=ndim eps=(1.0d-015)*n*n go to (401,400),irpa 400 continue do 1 i=1,n do 1 j=1,i aa(i,j)=aa(i,j)+bb(i,j) bb(i,j)=aa(i,j)-2.d0*bb(i,j) aa(j,i)=aa(i,j) bb(j,i)=bb(i,j) x(i,j)=bb(i,j) x(j,i)=x(i,j) 1 continue call diasym(x,ncf,n,v,e,vec,nord) do 2 i=1,n cec(i)=v(i) do 2 j=1,n tep(i,j)=x(i,j) 2 continue do 1000 i=1,n do 1001 j=1,n xnew=0.d0 do 1002 k=1,n xnew=xnew+tep(k,i)*aa(k,j) 1002 continue x(i,j)=xnew 1001 continue 1000 continue do 1003 i=1,n do 1004 j=1,i bbij=0.d0 do 1005 k=1,n bbij=bbij+x(i,k)*tep(k,j) 1005 continue bb(i,j)=bbij bb(j,i)=bbij 1004 continue 1003 continue do 5 i=1,n do 5 j=1,i arpms=dsqrt(cec(i))*bb(i,j)*dsqrt(cec(j)) x(i,j)=arpms x(j,i)=x(i,j) 5 continue call diasym(x,ncf,n,v,e,vec,nord) do 6 i=1,n alp(i)=v(i) if (alp(i).gt.0.D0) go to 150 write(1,151),i 151 format(2x,'state number',i3,'is imaginary') go to 6 150 valp(i)=dsqrt(alp(i)) 6 continue do 60 i=1,n do 61 j=1,n x(i,j)=dsqrt(valp(j))*x(i,j) 61 continue 60 continue do 70 i=1,n do 7 j=1,n aa(i,j)=0.d0 bb(i,j)=0.d0 do 8 k=1,n yy1=0.5d0*tep(i,k)*dsqrt(cec(k))*x(k,j) yy2=(0.5d0*tep(i,k)/dsqrt(cec(k)))*x(k,j) aa(i,j)=aa(i,j)+yy1 bb(i,j)=bb(i,j)+yy2 8 continue aa(i,j)=aa(i,j)/valp(j) if(alp(i)) 25,25,24 24 aa(i,j)=aa(i,j)+bb(i,j) bb(i,j)=aa(i,j)-2.d0*bb(i,j) 25 continue 7 continue 70 continue nd2=2*n do 130 j=1,n evr(j)=valp(j) do 131 i=1,nd2 if(i.gt.n) go to 140 ev=aa(i,j) go to 141 140 ev=bb(i-n,j) 141 vecr(i,j)=ev 131 continue 130 continue return 401 do 402 i=1,n do 403 j=1,n x(i,j)=aa(i,j) 403 continue 402 continue call diasym(x,ncf,n,v,e,vec,nord) do 404 j=1,n evr(j)=v(j) do 405 i=1,n 405 vecr(i,j)=x(i,j) 404 continue return end FUNCTION REDY(L1,L2,L3) implicit real*8 (a-h,o-z) real*4 x1,x2,x3,xm,c3j redy=0.d0 l13m=iabs(l1-l3) l13p=l1+l3 if((l2-l13m)*(l2-l13p).gt.0) go to 1 l123=l13p+l2 l123=mod(l123,2) if(l123.ne.0) go to 1 qp=4.d0*3.14159265d0 x4=(2.d0*l1+1.d0)*(2.d0*l2+1.d0)*(2.d0*l3+1.d0) x4=x4/qp sg=1.d0-2.d0*mod(l1,2) x1=float(l1) x2=float(l2) x3=float(l3) xm=0. call troisj(x1,x2,x3,xm,xm,xm,c3j) redy=sg*dsqrt(x4)*dble(c3j) 1 return end FUNCTION YL(LA,JA,LB,JB,L) ccc these are the reduced matrix elements of spherical harmonics: ccc ja and jb are twice the true values; ccc other variables have their true values implicit real*8 (a-h,o-z) qp=4.D0*3.14159265D0 yl=0.D0 l1p=la+lb l1m=iabs(la-lb) l2p=(ja+jb)/2 l2m=iabs(ja-jb)/2 lp=min0(l1p,l2p) lm=max0(l1m,l2m) if((lp-lm).lt.0) go to 1 l3=(l-lp)*(l-lm) if(l3.gt.0) go to 1 l4=la+lb+l-2*((la+lb+l)/2) if(l4.ne.0) go to 1 ld=2*l ig=l+(jb-1)/2 sg=1.D0-2.D0*mod(ig,2) xja=ja/2.d0 xjb=jb/2.d0 xl=l jpha=(ja+jb)/2+1 jpha=1-2*mod(jpha,2) c=(ja+1.D0)*(jb+1.D0)/qp sg=sg*jpha yl=sg*dsqrt(c)*cofcg(xja,xjb,xl,0.5d0,-0.5d0,0.d0) 1 return end FUNCTION TJL(LA,JA,LB,JB,J,L) ccc ja and jb are twice the true values; ccc other variables have their true values implicit real*8 (a-h,o-z) tjl=0.D0 z=0.D0 qp=4.D0*3.14159265D0 fj=float(j) fl=float(l) jd=2*j lp=la+lb lm=iabs(la-lb) ll=(l-lp)*(l-lm) if(ll.gt.0) go to 1 ll=la+lb+l-2*((la+lb+l)/2) if(ll.ne.0) go to 1 jp=(ja+jb)/2 jm=iabs(ja-jb)/2 jj=(j-jp)*(j-jm) if(jj.gt.0) go to 1 ljp=l+j ljm=iabs(l-j) lj=(1-ljp)*(1-ljm) if(lj.gt.0) go to 1 if(l.ne.j) go to 2 ig=(ja+jb)/2+j cc=(jd+1.D0)/(fj*(fj+1.D0)) z=-0.5D0*dsqrt(cc)*(jb+1.D0+(1.D0-2.D0*mod(ig,2))*(ja 1+1.D0)) go to 3 2 ig=lb+(jb+1)/2 cc=(la-0.5D0*ja)*(ja+1.D0)+(lb-0.5D0*jb)*(jb+1.D0) sg=1.D0-2.D0*mod(ig,2) if(l.eq.(j+1)) z=sg*(cc+fl)/dsqrt(fl) if(l.eq.(j-1)) z=sg*(cc-fl-1.d0)/dsqrt(fl+1.d0) 3 c=(ja+1.D0)*(jb+1.d0)/((jd+1.d0)*qp) xja=ja/2.d0 xjb=jb/2.d0 xjd=j jpha=(ja+jb)/2+1 jpha=1-2*mod(jpha,2) tjl=(1.d0-2.d0*mod(la,2))*dsqrt(c)*z*jpha 1*cofcg(xja,xjb,xjd,0.5d0,-0.5d0,0.d0) 1 return end FUNCTION COFCG(A,B,C,X,Y,Z) implicit real*8 (a-h,o-z) real*4 xj1,xj2,xj3,xm1,xm2,xm3,c3j xj1=sngl(a) xj2=sngl(b) xj3=sngl(c) xm1=sngl(x) xm2=sngl(y) xm3=-sngl(z) call troisj(xj1,xj2,xj3,xm1,xm2,xm3,c3j) small=1.d-06 cof=a-b+z cof=dabs(cof)+small icof=int(cof) cof=1.d0-2.d0*mod(icof,2) cofcg=dble(c3j)*cof*dsqrt(2.d0*c+1.d0) return end