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) TRO01770 DIMENSION FLOG(301) TRO01780 DATA(FLOG(I),I=2,31)/0.,.69314718,1.7917595,3.1780538,4.7874917,6.TRO01820 15792511,8.5251613,10.604603,12.801827,15.104413,17.502307,19.98721TRO01830 24,22.552163,25.191221,27.899271,30.671860,33.505072,36.395445,39.3TRO01840 339884,42.335616,45.380139,48.471180,51.606674,54.784729,58.003604,TRO01850 461.261702,64.557537,67.889743,71.257038,74.658235/ TRO01860 DATA(FLOG(I),I=32,61)/78.092223,81.557959,85.054466,88.580827,92.1TRO01870 136175,95.719694,99.330612,102.96820,106.63176,110.32064,114.03421,TRO01880 2117.77188,121.53308,125.31727,129.12393,132.95257,136.80272,140.67TRO01890 3392,144.56574,148.47776,152.40959,156.36083,160.33112,164.32011,16TRO01900 48.32744,172.35279,176.39584,180.45629,184.53383,188.62817/ TRO01910 DATA(FLOG(I),I=62,91)/192.73904,196.86618,201.00931,205.16820,209.TRO01920 134258,213.53224,217.73693,221.95644,226.19054,230.43904,234.70172,TRO01930 2238.97839,243.26885,247.57291,251.89040,256.22113,260.56494,264.92TRO01940 3164,269.29110,273.67312,278.06757,282.47429,286.89313,291.32394,29TRO01950 45.76659,300.22094,304.68685,309.16419,313.65283,318.15264/ TRO01960 DATA(FLOG(I),I=92,121)/322.66349,327.18529,331.71788,336.26118,340TRO01970 1.81505,345.37940,349.95411,354.53908,359.13420,363.73937,368.35449TRO01980 2,372.97946,377.61419,382.25859,386.91255,391.57598,396.24881,400.9TRO01990 33094,405.62230,410.32277,415.03230,419.75080,424.47819,429.21439,4TRO02000 433.95932,438.71291,443.47508,448.24576,453.02489,457.81238/ TRO02010 DATA(FLOG(I),I=122,151)/462.60817,467.41220,472.22438,477.04466,48TRO02020 11.87298,486.70926,491.55345,496.40547,501.26529,506.13282,511.0080TRO02030 22,515.89082,520.78117,525.67901,530.58428,535.49694,540.41692,545.TRO02040 334417,550.27865,555.22029,560.16905,565.12488,570.08772,575.05753,TRO02050 4580.03427,585.01787,590.00830,595.00552,600.00946,605.02010/ TRO02060 DATA(FLOG(I),I=152,181)/610.03738,615.06126,620.09170,625.12866,63TRO02070 10.17208,635.22193,640.27818,645.34077,650.40968,655.48486,660.5662TRO02080 26,665.65385,670.74760,675.84747,680.95341,686.06541,691.18340,696.TRO02090 330735,701.43726,706.57306,711.71472,716.86221,722.01551,727.17456,TRO02100 4732.33934,737.50983,742.68598,747.86776,753.05516,758.24811/ TRO02110 DATA(FLOG(I),I=182,211)/763.44661,768.65061,773.86010,779.07503,78TRO02120 14.29539,789.52114,794.75224,799.98869,805.23044,810.47747,815.7297TRO02130 23,820.98722,826.24991,831.51778,836.79078,842.06890,847.35209,852.TRO02140 364036,857.93366,863.23199,868.53529,873.84356,879.15676,884.47488,TRO02150 4889.79789,895.12577,900.45848,905.79603,911.13836,916.48547/ TRO02160 DATA(FLOG(I),I=212,241)/921.83732,927.19391,932.55521,937.92118,94TRO02170 13.29181,948.66710,954.04699,959.43148,964.82056,970.21419,975.6123TRO02180 25,981.01503,986.42220,991.83385,997.24995,1002.6705,1008.0954,1013TRO02190 3.5248,1018.9585,1024.3966,1029.8389,1035.2857,1040.7367,1046.1920,TRO02200 41051.6516,1057.1155,1062.5836,1068.0558,1073.5323,1079.0129/ TRO02210 DATA(FLOG(I),I=242,271)/1084.4977,1089.9866,1095.4797,1100.9768,11TRO02220 106.4781,1111.9834,1117.4928,1123.0063,1128.5237,1134.0452,1139.570TRO02230 26,1145.1001,1150.6335,1156.1708,1161.7120,1167.2573,1172.8063,1178TRO02240 3.3593,1183.9161,1189.4768,1195.0413,1200.6097,1206.1818,1211.7577,TRO02250 41217.3375,1222.9209,1228.5082,1234.0992,1239.6939,1245.2924/ TRO02260 DATA(FLOG(I),I=272,301)/1250.8944,1256.5003,1262.1097,1267.7228,12TRO02270 173.3396,1278.9600,1284.5840,1290.2117,1295.8429,1301.4777,1307.116TRO02280 20,1312.7580,1318.4034,1324.0524,1329.7048,1335.3609,1341.0203,1346TRO02290 3.6833,1352.3497,1358.0196,1363.6929,1369.3697,1375.0499,1380.7334,TRO02300 41386.4204,1392.1107,1397.8045,1403.5016,1409.2020,1414.9058/ TRO02310 DATA EPS1,EPS2/.1,-.2/ TRO02320 XN=-XJ1+XJ2+XJ3+EPS1 TRO02360 IF (XN.LT.0.) XN=XN+EPS2 TRO02370 N1=XN TRO02380 XN=-XL1+XL2+XJ3+EPS1 TRO02390 IF (XN.LT.0.) XN=XN+EPS2 TRO02400 N2=XN TRO02410 XN=-XL1+XJ2+XL3+EPS1 TRO02420 IF (XN.LT.0.) XN=XN+EPS2 TRO02430 N3=XN TRO02440 XN=-XJ1+XL2+XL3+EPS1 TRO02450 IF (XN.LT.0.) XN=XN+EPS2 TRO02460 N4=XN TRO02470 XN=XJ1-XJ2+XJ3+EPS1 TRO02480 IF (XN.LT.0.) XN=XN+EPS2 TRO02490 N5=XN TRO02500 XN=XL1-XL2+XJ3+EPS1 TRO02510 IF (XN.LT.0.) XN=XN+EPS2 TRO02520 N6=XN TRO02530 XN=XL1-XJ2+XL3+EPS1 TRO02540 IF (XN.LT.0.) XN=XN+EPS2 TRO02550 N7=XN TRO02560 XN=XJ1-XL2+XL3+EPS1 TRO02570 IF (XN.LT.0.) XN=XN+EPS2 TRO02580 N8=XN TRO02590 XN=XJ1+XJ2-XJ3+EPS1 TRO02600 IF (XN.LT.0.) XN=XN+EPS2 TRO02610 N9=XN TRO02620 XN=XL1+XL2-XJ3+EPS1 TRO02630 IF (XN.LT.0.) XN=XN+EPS2 TRO02640 N10=XN TRO02650 XN=XL1+XJ2-XL3+EPS1 TRO02660 IF (XN.LT.0.) XN=XN+EPS2 TRO02670 N11=XN TRO02680 XN=XJ1+XL2-XL3+EPS1 TRO02690 IF (XN.LT.0.) XN=XN+EPS2 TRO02700 N12=XN TRO02710 XN=-XJ1-XL1+XJ3+XL3+EPS1 TRO02720 IF (XN.LT.0.) XN=XN+EPS2 TRO02730 N13=XN TRO02740 XN=-XJ2-XL2+XJ3+XL3+EPS1 TRO02750 IF (XN.LT.0.) XN=XN+EPS2 TRO02760 N14=XN TRO02770 XN=XJ1+XL1+XJ2+XL2+EPS1 TRO02780 IF (XN.LT.0.) XN=XN+EPS2 TRO02790 N15=XN TRO02800 N15=N15+1 TRO02810 XN=XJ1+XJ2+XJ3+EPS1 TRO02820 IF (XN.LT.0.) XN=XN+EPS2 TRO02830 N16=XN TRO02840 N16=N16+1 TRO02850 XN=XL1+XL2+XJ3+EPS1 TRO02860 IF (XN.LT.0.) XN=XN+EPS2 TRO02870 N17=XN TRO02880 N17=N17+1 TRO02890 XN=XL1+XJ2+XL3+EPS1 TRO02900 IF (XN.LT.0.) XN=XN+EPS2 TRO02910 N18=XN TRO02920 N18=N18+1 TRO02930 XN=XJ1+XL2+XL3+EPS1 TRO02940 IF (XN.LT.0.) XN=XN+EPS2 TRO02950 N19=XN TRO02960 N19=N19+1 TRO02970 IF(N9.LT.0) GO TO 50 TRO03010 IF(N5.LT.0) GO TO 50 TRO03020 IF(N1.LT.0) GO TO 50 TRO03030 IF(N10.LT.0) GO TO 50 TRO03040 IF(N6.LT.0) GO TO 50 TRO03050 IF(N2.LT.0) GO TO 50 TRO03060 IF(N11.LT.0) GO TO 50 TRO03070 IF(N7.LT.0) GO TO 50 TRO03080 IF(N3.LT.0) GO TO 50 TRO03090 IF(N12.LT.0) GO TO 50 TRO03100 IF(N8.LT.0) GO TO 50 TRO03110 IF(N4.LT.0) GO TO 50 TRO03120 K=N17+N18+N19-3 TRO03130 IF(K.GT.0) GO TO 54 TRO03140 C6J=1. TRO03150 RETURN TRO03160 50 C6J=0. TRO03170 RETURN TRO03180 54 K=0 TRO03220 L=-N13 TRO03230 IF(L.GT.K) K=L TRO03240 L=-N14 TRO03250 IF(L.GT.K) K=L TRO03260 L=N9 TRO03270 IF(N10.LT.L) L=N10 TRO03280 IF(N11.LT.L) L=N11 TRO03290 IF(N12.LT.L) L=N12 TRO03300 IF(N15.LT.L) L=N15 TRO03310 F=1. TRO03320 S=1. TRO03330 I=K+1 TRO03340 62 IF(I.GT.L) GO TO 80 TRO03350 IM1=I-1 TRO03360 NN=(N9-IM1)*(N10-IM1)*(N11-IM1)*(N12-IM1) TRO03370 ND=I*(N13+I)*(N14+I)*(N15-IM1) TRO03380 F=-F*FLOAT(NN)/FLOAT(ND) TRO03390 S=S+F TRO03400 I=I+1 TRO03410 GO TO 62 TRO03420 80 C2N=FLOG(N1+1)+FLOG(N2+1)+FLOG(N3+1)+FLOG(N4+1)+FLOG(N5+1)+FLOG(N6TRO03460 1+1)+FLOG(N7+1)+FLOG(N8+1)+FLOG(N9+1)+FLOG(N10+1)+FLOG(N11+1)+FLOG(TRO03470 2N12+1) TRO03480 C2N=.5*C2N TRO03490 C2D=FLOG(N16+1)+FLOG(N17+1)+FLOG(N18+1)+FLOG(N19+1) TRO03500 C2D=.5*C2D TRO03510 KM1=K-1 TRO03520 KP1=K+1 TRO03530 C2N=C2N+FLOG(N15-KM1) TRO03540 C2D=C2D+FLOG(KP1)+FLOG(N13+KP1)+FLOG(N14+KP1)+FLOG(N9-KM1)+FLOG(N1TRO03550 10-KM1)+FLOG(N11-KM1)+FLOG(N12-KM1) TRO03560 F=C2D-C2N TRO03600 IF(F.GT.80.) GO TO 98 TRO03610 F=C2N/C2D TRO03620 IF((F.LT.1.01).AND.(F.GT.0.98)) GO TO 98 TRO03630 C6J=S*EXP(C2N-C2D) TRO03640 GO TO 106 TRO03650 98 IF(S) 100,50,102 TRO03660 100 S=ALOG(-S) TRO03670 C6J=-EXP(S+C2N-C2D) TRO03680 GO TO 106 TRO03690 102 S=ALOG(S) TRO03700 C6J=EXP(S+C2N-C2D) TRO03710 106 L=N15+KM1 TRO03750 K=L/2 TRO03760 K=2*K TRO03770 IF(L.NE.K) C6J=-C6J TRO03780 RETURN TRO03790 END TRO03800 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) TRO00010 DIMENSION FLOG(301) TRO00020 DATA(FLOG(I),I=2,31)/0.,.69314718,1.7917595,3.1780538,4.7874917,6.TRO00060 15792511,8.5251613,10.604603,12.801827,15.104413,17.502307,19.98721TRO00070 24,22.552163,25.191221,27.899271,30.671860,33.505072,36.395445,39.3TRO00080 339884,42.335616,45.380139,48.471180,51.606674,54.784729,58.003604,TRO00090 461.261702,64.557537,67.889743,71.257038,74.658235/ TRO00100 DATA(FLOG(I),I=32,61)/78.092223,81.557959,85.054466,88.580827,92.1TRO00110 136175,95.719694,99.330612,102.96820,106.63176,110.32064,114.03421,TRO00120 2117.77188,121.53308,125.31727,129.12393,132.95257,136.80272,140.67TRO00130 3392,144.56574,148.47776,152.40959,156.36083,160.33112,164.32011,16TRO00140 48.32744,172.35279,176.39584,180.45629,184.53383,188.62817/ TRO00150 DATA(FLOG(I),I=62,91)/192.73904,196.86618,201.00931,205.16820,209.TRO00160 134258,213.53224,217.73693,221.95644,226.19054,230.43904,234.70172,TRO00170 2238.97839,243.26885,247.57291,251.89040,256.22113,260.56494,264.92TRO00180 3164,269.29110,273.67312,278.06757,282.47429,286.89313,291.32394,29TRO00190 45.76659,300.22094,304.68685,309.16419,313.65283,318.15264/ TRO00200 DATA(FLOG(I),I=92,121)/322.66349,327.18529,331.71788,336.26118,340TRO00210 1.81505,345.37940,349.95411,354.53908,359.13420,363.73937,368.35449TRO00220 2,372.97946,377.61419,382.25859,386.91255,391.57598,396.24881,400.9TRO00230 33094,405.62230,410.32277,415.03230,419.75080,424.47819,429.21439,4TRO00240 433.95932,438.71291,443.47508,448.24576,453.02489,457.81238/ TRO00250 DATA(FLOG(I),I=122,151)/462.60817,467.41220,472.22438,477.04466,48TRO00260 11.87298,486.70926,491.55345,496.40547,501.26529,506.13282,511.0080TRO00270 22,515.89082,520.78117,525.67901,530.58428,535.49694,540.41692,545.TRO00280 334417,550.27865,555.22029,560.16905,565.12488,570.08772,575.05753,TRO00290 4580.03427,585.01787,590.00830,595.00552,600.00946,605.02010/ TRO00300 DATA(FLOG(I),I=152,181)/610.03738,615.06126,620.09170,625.12866,63TRO00310 10.17208,635.22193,640.27818,645.34077,650.40968,655.48486,660.5662TRO00320 26,665.65385,670.74760,675.84747,680.95341,686.06541,691.18340,696.TRO00330 330735,701.43726,706.57306,711.71472,716.86221,722.01551,727.17456,TRO00340 4732.33934,737.50983,742.68598,747.86776,753.05516,758.24811/ TRO00350 DATA(FLOG(I),I=182,211)/763.44661,768.65061,773.86010,779.07503,78TRO00360 14.29539,789.52114,794.75224,799.98869,805.23044,810.47747,815.7297TRO00370 23,820.98722,826.24991,831.51778,836.79078,842.06890,847.35209,852.TRO00380 364036,857.93366,863.23199,868.53529,873.84356,879.15676,884.47488,TRO00390 4889.79789,895.12577,900.45848,905.79603,911.13836,916.48547/ TRO00400 DATA(FLOG(I),I=212,241)/921.83732,927.19391,932.55521,937.92118,94TRO00410 13.29181,948.66710,954.04699,959.43148,964.82056,970.21419,975.6123TRO00420 25,981.01503,986.42220,991.83385,997.24995,1002.6705,1008.0954,1013TRO00430 3.5248,1018.9585,1024.3966,1029.8389,1035.2857,1040.7367,1046.1920,TRO00440 41051.6516,1057.1155,1062.5836,1068.0558,1073.5323,1079.0129/ TRO00450 DATA(FLOG(I),I=242,271)/1084.4977,1089.9866,1095.4797,1100.9768,11TRO00460 106.4781,1111.9834,1117.4928,1123.0063,1128.5237,1134.0452,1139.570TRO00470 26,1145.1001,1150.6335,1156.1708,1161.7120,1167.2573,1172.8063,1178TRO00480 3.3593,1183.9161,1189.4768,1195.0413,1200.6097,1206.1818,1211.7577,TRO00490 41217.3375,1222.9209,1228.5082,1234.0992,1239.6939,1245.2924/ TRO00500 DATA(FLOG(I),I=272,301)/1250.8944,1256.5003,1262.1097,1267.7228,12TRO00510 173.3396,1278.9600,1284.5840,1290.2117,1295.8429,1301.4777,1307.116TRO00520 20,1312.7580,1318.4034,1324.0524,1329.7048,1335.3609,1341.0203,1346TRO00530 3.6833,1352.3497,1358.0196,1363.6929,1369.3697,1375.0499,1380.7334,TRO00540 41386.4204,1392.1107,1397.8045,1403.5016,1409.2020,1414.9058/ TRO00550 DATA EPS1,EPS2/.1,-.2/ TRO00560 XN=XJ2-XM2+EPS1 TRO00600 IF(XN.LT.0.) XN=XN+EPS2 TRO00610 N1=XN TRO00620 XN=XJ3+XM3+EPS1 TRO00630 IF(XN.LT.0.) XN=XN+EPS2 TRO00640 N2=XN TRO00650 IF(XN.LT.0.) XN=XN+EPS2 TRO00660 XN=XJ3-XM3+EPS1 TRO00670 N3=XN TRO00680 XN=XJ1+XM1+EPS1 TRO00690 IF(XN.LT.0.) XN=XN+EPS2 TRO00700 N4=XN TRO00710 XN=XJ2+XJ3-XJ1+EPS1 TRO00720 IF(XN.LT.0.) XN=XN+EPS2 TRO00730 N5=XN TRO00740 XN=XJ1+XJ3-XJ2+EPS1 TRO00750 IF(XN.LT.0.) XN=XN+EPS2 TRO00760 N6=XN TRO00770 XN=XJ1+XJ2-XJ3+EPS1 TRO00780 IF(XN.LT.0.) XN=XN+EPS2 TRO00790 N7=XN TRO00800 XN=XJ1-XM1+EPS1 TRO00810 IF(XN.LT.0.) XN=XN+EPS2 TRO00820 N8=XN TRO00830 XN=XJ2+XM2+EPS1 TRO00840 IF(XN.LT.0.) XN=XN+EPS2 TRO00850 N9=XN TRO00860 XN=XJ3-XJ1-XM2+EPS1 TRO00870 IF(XN.LT.0.) XN=XN+EPS2 TRO00880 N10=XN TRO00890 XN=XJ3-XJ2+XM1+EPS1 TRO00900 IF(XN.LT.0.) XN=XN+EPS2 TRO00910 N11=XN TRO00920 XN=XJ1+XJ2+XJ3+EPS1 TRO00930 IF(XN.LT.0.) XN=XN+EPS2 TRO00940 N12=XN TRO00950 N12=N12+1 TRO00960 K=N4*N8 TRO01000 IF(K.LT.0) GO TO 59 TRO01010 K=N1*N9 TRO01020 IF(K.LT.0) GO TO 59 TRO01030 K=N2*N3 TRO01040 IF(K.LT.0) GO TO 59 TRO01050 IF(N5.LT.0) GO TO 59 TRO01060 IF(N6.LT.0) GO TO 59 TRO01070 IF(N7.LT.0) GO TO 59 TRO01080 L=N1-N2+N3-N4+N8-N9 TRO01090 IF(L.NE.0) GO TO 59 TRO01100 K=N12-1 TRO01110 IF(K.GT.0) GO TO 60 TRO01120 C3J=1. TRO01130 RETURN TRO01140 59 C3J=0. TRO01150 RETURN TRO01160 60 K=0 TRO01200 L=-N10 TRO01210 IF(L.GT.K) K=L TRO01220 L=-N11 TRO01230 IF(L.GT.K) K=L TRO01240 L=N7 TRO01250 IF(N8.GT.L) L=N8 TRO01260 IF(N9.GT.L) L=N9 TRO01270 F=1. TRO01280 S=1. TRO01290 I=K+1 TRO01300 62 IF(I.GT.L) GO TO 80 TRO01310 IM1=I-1 TRO01320 NN=(N7-IM1)*(N8-IM1)*(N9-IM1) TRO01330 ND=I*(N10+I)*(N11+I) TRO01340 F=-F*FLOAT(NN)/FLOAT(ND) TRO01350 S=S+F TRO01360 I=I+1 TRO01370 GO TO 62 TRO01380 80 C2N=FLOG(N1+1)+FLOG(N2+1)+FLOG(N3+1)+FLOG(N4+1)+FLOG(N5+1)+FLOG(N6TRO01420 1+1)+FLOG(N7+1)+FLOG(N8+1)+FLOG(N9+1) TRO01430 C2N=.5*C2N TRO01440 KM1=K-1 TRO01450 KP1=K+1 TRO01460 C2D=FLOG(KP1)+FLOG(N7-KM1)+FLOG(N8-KM1)+FLOG(N9-KM1)+FLOG(N10+KP1)TRO01470 1+FLOG(N11+KP1)+.5*FLOG(N12+1) TRO01480 F=C2D-C2N TRO01520 IF(F.GT.80.) GO TO 98 TRO01530 F=C2N/C2D TRO01540 IF((F.LT.1.01).AND.(F.GT.0.98)) GO TO 98 TRO01550 C3J=S*EXP(C2N-C2D) TRO01560 GO TO 106 TRO01570 98 IF(S) 100,59,102 TRO01580 100 S=ALOG(-S) TRO01590 C3J=-EXP(S+C2N-C2D) TRO01600 GO TO 106 TRO01610 102 S=ALOG(S) TRO01620 C3J=EXP(S+C2N-C2D) TRO01630 106 XN=XJ1-XJ2-XM3+EPS1 TRO01670 IF(XN.LT.0.) XN=XN+EPS2 TRO01680 L=XN TRO01690 L=L+K TRO01700 K=L/2 TRO01710 K=2*K TRO01720 IF(L.NE.K) C3J=-C3J TRO01730 RETURN TRO01740 END TRO01750 SUBROUTINE NEUFJ(XJ11,XJ12,XJ13,XJ21,XJ22,XJ23,XJ31,XJ32,XJ33,C9J)TRO03820 DIMENSION FLOG(301) TRO03830 DATA(FLOG(I),I=2,31)/0.,.69314718,1.7917595,3.1780538,4.7874917,6.TRO03870 15792511,8.5251613,10.604603,12.801827,15.104413,17.502307,19.98721TRO03880 24,22.552163,25.191221,27.899271,30.671860,33.505072,36.395445,39.3TRO03890 339884,42.335616,45.380139,48.471180,51.606674,54.784729,58.003604,TRO03900 461.261702,64.557537,67.889743,71.257038,74.658235/ TRO03910 DATA(FLOG(I),I=32,61)/78.092223,81.557959,85.054466,88.580827,92.1TRO03920 136175,95.719694,99.330612,102.96820,106.63176,110.32064,114.03421,TRO03930 2117.77188,121.53308,125.31727,129.12393,132.95257,136.80272,140.67TRO03940 3392,144.56574,148.47776,152.40959,156.36083,160.33112,164.32011,16TRO03950 48.32744,172.35279,176.39584,180.45629,184.53383,188.62817/ TRO03960 DATA(FLOG(I),I=62,91)/192.73904,196.86618,201.00931,205.16820,209.TRO03970 134258,213.53224,217.73693,221.95644,226.19054,230.43904,234.70172,TRO03980 2238.97839,243.26885,247.57291,251.89040,256.22113,260.56494,264.92TRO03990 3164,269.29110,273.67312,278.06757,282.47429,286.89313,291.32394,29TRO04000 45.76659,300.22094,304.68685,309.16419,313.65283,318.15264/ TRO04010 DATA(FLOG(I),I=92,121)/322.66349,327.18529,331.71788,336.26118,340TRO04020 1.81505,345.37940,349.95411,354.53908,359.13420,363.73937,368.35449TRO04030 2,372.97946,377.61419,382.25859,386.91255,391.57598,396.24881,400.9TRO04040 33094,405.62230,410.32277,415.03230,419.75080,424.47819,429.21439,4TRO04050 433.95932,438.71291,443.47508,448.24576,453.02489,457.81238/ TRO04060 DATA(FLOG(I),I=122,151)/462.60817,467.41220,472.22438,477.04466,48TRO04070 11.87298,486.70926,491.55345,496.40547,501.26529,506.13282,511.0080TRO04080 22,515.89082,520.78117,525.67901,530.58428,535.49694,540.41692,545.TRO04090 334417,550.27865,555.22029,560.16905,565.12488,570.08772,575.05753,TRO04100 4580.03427,585.01787,590.00830,595.00552,600.00946,605.02010/ TRO04110 DATA(FLOG(I),I=152,181)/610.03738,615.06126,620.09170,625.12866,63TRO04120 10.17208,635.22193,640.27818,645.34077,650.40968,655.48486,660.5662TRO04130 26,665.65385,670.74760,675.84747,680.95341,686.06541,691.18340,696.TRO04140 330735,701.43726,706.57306,711.71472,716.86221,722.01551,727.17456,TRO04150 4732.33934,737.50983,742.68598,747.86776,753.05516,758.24811/ TRO04160 DATA(FLOG(I),I=182,211)/763.44661,768.65061,773.86010,779.07503,78TRO04170 14.29539,789.52114,794.75224,799.98869,805.23044,810.47747,815.7297TRO04180 23,820.98722,826.24991,831.51778,836.79078,842.06890,847.35209,852.TRO04190 364036,857.93366,863.23199,868.53529,873.84356,879.15676,884.47488,TRO04200 4889.79789,895.12577,900.45848,905.79603,911.13836,916.48547/ TRO04210 DATA(FLOG(I),I=212,241)/921.83732,927.19391,932.55521,937.92118,94TRO04220 13.29181,948.66710,954.04699,959.43148,964.82056,970.21419,975.6123TRO04230 25,981.01503,986.42220,991.83385,997.24995,1002.6705,1008.0954,1013TRO04240 3.5248,1018.9585,1024.3966,1029.8389,1035.2857,1040.7367,1046.1920,TRO04250 41051.6516,1057.1155,1062.5836,1068.0558,1073.5323,1079.0129/ TRO04260 DATA(FLOG(I),I=242,271)/1084.4977,1089.9866,1095.4797,1100.9768,11TRO04270 106.4781,1111.9834,1117.4928,1123.0063,1128.5237,1134.0452,1139.570TRO04280 26,1145.1001,1150.6335,1156.1708,1161.7120,1167.2573,1172.8063,1178TRO04290 3.3593,1183.9161,1189.4768,1195.0413,1200.6097,1206.1818,1211.7577,TRO04300 41217.3375,1222.9209,1228.5082,1234.0992,1239.6939,1245.2924/ TRO04310 DATA(FLOG(I),I=272,301)/1250.8944,1256.5003,1262.1097,1267.7228,12TRO04320 173.3396,1278.9600,1284.5840,1290.2117,1295.8429,1301.4777,1307.116TRO04330 20,1312.7580,1318.4034,1324.0524,1329.7048,1335.3609,1341.0203,1346TRO04340 3.6833,1352.3497,1358.0196,1363.6929,1369.3697,1375.0499,1380.7334,TRO04350 41386.4204,1392.1107,1397.8045,1403.5016,1409.2020,1414.9058/ TRO04360 DATA EPS1,EPS2/.1,-.2/ TRO04370 XN=-XJ11+XJ21+XJ31+EPS1 TRO04410 IF(XN.LT.0.) XN=XN+EPS2 TRO04420 N1=XN TRO04430 XN=-XJ32+XJ33+XJ31+EPS1 TRO04440 IF(XN.LT.0.) XN=XN+EPS2 TRO04450 N2=XN TRO04460 XN=XJ11-XJ21+XJ31+EPS1 TRO04470 IF(XN.LT.0.) XN=XN+EPS2 TRO04480 N5=XN TRO04490 XN=XJ32-XJ33+XJ31+EPS1 TRO04500 IF(XN.LT.0.) XN=XN+EPS2 TRO04510 N6=XN TRO04520 XN=XJ11+XJ21-XJ31+EPS1 TRO04530 IF(XN.LT.0.) XN=XN+EPS2 TRO04540 N9=XN TRO04550 XN=XJ32+XJ33-XJ31+EPS1 TRO04560 IF(XN.LT.0.) XN=XN+EPS2 TRO04570 N10=XN TRO04580 XN=XJ11+XJ32+XJ21+XJ33+EPS1 TRO04590 IF(XN.LT.0.) XN=XN+EPS2 TRO04600 N15=XN TRO04610 N15=N15+1 TRO04620 XN=XJ11+XJ21+XJ31+EPS1 TRO04630 IF(XN.LT.0.) XN=XN+EPS2 TRO04640 N16=XN TRO04650 N16=N16+1 TRO04660 XN=XJ32+XJ33+XJ31+EPS1 TRO04670 IF(XN.LT.0.) XN=XN+EPS2 TRO04680 N17=XN TRO04690 N17=N17+1 TRO04700 XN=-XJ12+XJ22+XJ32+EPS1 TRO04710 IF(XN.LT.0.) XN=XN+EPS2 TRO04720 N21=XN TRO04730 XN=-XJ21+XJ22+XJ23+EPS1 TRO04740 IF(XN.LT.0.) XN=XN+EPS2 TRO04750 N23=XN TRO04760 XN=XJ12-XJ22+XJ32+EPS1 TRO04770 IF(XN.LT.0.) XN=XN+EPS2 TRO04780 N25=XN TRO04790 XN=XJ21-XJ22+XJ23+EPS1 TRO04800 IF(XN.LT.0.) XN=XN+EPS2 TRO04810 N27=XN TRO04820 XN=XJ12+XJ22-XJ32+EPS1 TRO04830 IF(XN.LT.0.) XN=XN+EPS2 TRO04840 N29=XN TRO04850 XN=XJ21+XJ22-XJ23+EPS1 TRO04860 IF(XN.LT.0.) XN=XN+EPS2 TRO04870 N31=XN TRO04880 XN=-XJ12-XJ21+XJ32+XJ23+EPS1 TRO04890 IF(XN.LT.0.) XN=XN+EPS2 TRO04900 N33=XN TRO04910 XN=XJ12+XJ22+XJ32+EPS1 TRO04920 IF(XN.LT.0.) XN=XN+EPS2 TRO04930 N36=XN TRO04940 N36=N36+1 TRO04950 XN=XJ21+XJ22+XJ23+EPS1 TRO04960 IF(XN.LT.0.) XN=XN+EPS2 TRO04970 N38=XN TRO04980 N38=N38+1 TRO04990 XN=-XJ13+XJ23+XJ33+EPS1 TRO05000 IF(XN.LT.0.) XN=XN+EPS2 TRO05010 N41=XN TRO05020 XN=-XJ13+XJ11+XJ12+EPS1 TRO05030 IF(XN.LT.0.) XN=XN+EPS2 TRO05040 N44=XN TRO05050 XN=XJ13-XJ23+XJ33+EPS1 TRO05060 IF(XN.LT.0.) XN=XN+EPS2 TRO05070 N45=XN TRO05080 XN=XJ13-XJ11+XJ12+EPS1 TRO05090 IF(XN.LT.0.) XN=XN+EPS2 TRO05100 N48=XN TRO05110 XN=XJ13+XJ23-XJ33+EPS1 TRO05120 IF(XN.LT.0.) XN=XN+EPS2 TRO05130 N49=XN TRO05140 XN=XJ13+XJ11-XJ12+EPS1 TRO05150 IF(XN.LT.0.) XN=XN+EPS2 TRO05160 N52=XN TRO05170 XN=-XJ23-XJ11+XJ33+XJ12+EPS1 TRO05180 IF(XN.LT.0.) XN=XN+EPS2 TRO05190 N54=XN TRO05200 XN=XJ13+XJ23+XJ33+EPS1 TRO05210 IF(XN.LT.0.) XN=XN+EPS2 TRO05220 N56=XN TRO05230 N56=N56+1 TRO05240 XN=XJ13+XJ11+XJ12+EPS1 TRO05250 IF(XN.LT.0.) XN=XN+EPS2 TRO05260 N59=XN TRO05270 N59=N59+1 TRO05280 IF(N9.LT.0) GO TO 50 TRO05320 IF(N5.LT.0) GO TO 50 TRO05330 IF(N1.LT.0) GO TO 50 TRO05340 IF(N10.LT.0) GO TO 50 TRO05350 IF(N6.LT.0) GO TO 50 TRO05360 IF(N2.LT.0) GO TO 50 TRO05370 IF(N29.LT.0) GO TO 50 TRO05380 IF(N25.LT.0) GO TO 50 TRO05390 IF(N21.LT.0) GO TO 50 TRO05400 IF(N31.LT.0) GO TO 50 TRO05410 IF(N27.LT.0) GO TO 50 TRO05420 IF(N23.LT.0) GO TO 50 TRO05430 IF(N49.LT.0) GO TO 50 TRO05440 IF(N45.LT.0) GO TO 50 TRO05450 IF(N41.LT.0) GO TO 50 TRO05460 IF(N52.LT.0) GO TO 50 TRO05470 IF(N48.LT.0) GO TO 50 TRO05480 IF(N44.LT.0) GO TO 50 TRO05490 K=N1+N2+N5+N6+N9+N10+N21+N23+N25+N27+N29+N31+N41+N44+N45+N48+N49+NTRO05500 152 TRO05510 IF(K.GT.0) GO TO 54 TRO05520 C9J=1. TRO05530 RETURN TRO05540 50 C9J=0. TRO05550 RETURN TRO05560 54 XN=2.*(XJ21-XJ32)+EPS1 TRO05600 IF(XN.LT.0.) XN=-(XN+EPS2) TRO05610 JMIN=XN TRO05620 XN=2.*(XJ11-XJ33)+EPS1 TRO05630 IF(XN.LT.0.) XN=-(XN+EPS2) TRO05640 N=XN TRO05650 IF(N.GT.JMIN) JMIN=N TRO05660 XN=2.*(XJ12-XJ23)+EPS1 TRO05670 IF(XN.LT.0.) XN=-(XN+EPS2) TRO05680 N=XN TRO05690 IF(N.GT.JMIN) JMIN=N TRO05700 XN=2.*(XJ21+XJ32)+EPS1 TRO05710 JMAX=XN TRO05720 XN=2.*(XJ11+XJ33)+EPS1 TRO05730 N=XN TRO05740 IF(N.LT.JMAX) JMAX=N TRO05750 XN=2.*(XJ12+XJ23)+EPS1 TRO05760 N=XN TRO05770 IF(N.LT.JMAX) JMAX=N TRO05780 XJMIN=FLOAT(JMIN)/2. TRO05790 XJMAX=FLOAT(JMAX)/2. TRO05800 S=0. TRO05810 XJ=XJMIN TRO05820 XN=-XJ32+XJ21+XJ+EPS1 TRO05830 IF(XN.LT.0.) XN=XN+EPS2 TRO05840 N3=XN TRO05850 XN=-XJ11+XJ33+XJ+EPS1 TRO05860 IF(XN.LT.0.) XN=XN+EPS2 TRO05870 N4=XN TRO05880 XN=XJ32-XJ21+XJ+EPS1 TRO05890 IF(XN.LT.0.) XN=XN+EPS2 TRO05900 N7=XN TRO05910 XN=XJ11-XJ33+XJ+EPS1 TRO05920 IF(XN.LT.0.) XN=XN+EPS2 TRO05930 N8=XN TRO05940 XN=XJ32+XJ21-XJ+EPS1 TRO05950 IF(XN.LT.0.) XN=XN+EPS2 TRO05960 N11=XN TRO05970 XN=XJ11+XJ33-XJ+EPS1 TRO05980 IF(XN.LT.0.) XN=XN+EPS2 TRO05990 N12=XN TRO06000 XN=-XJ11-XJ32+XJ31+XJ+EPS1 TRO06010 IF(XN.LT.0.) XN=XN+EPS2 TRO06020 N13=XN TRO06030 XN=-XJ21-XJ33+XJ31+XJ+EPS1 TRO06040 IF(XN.LT.0.) XN=XN+EPS2 TRO06050 N14=XN TRO06060 XN=XJ32+XJ21+XJ+EPS1 TRO06070 IF(XN.LT.0.) XN=XN+EPS2 TRO06080 N18=XN TRO06090 N18=N18+1 TRO06100 XN=XJ11+XJ33+XJ+EPS1 TRO06110 IF(XN.LT.0.) XN=XN+EPS2 TRO06120 N19=XN TRO06130 N19=N19+1 TRO06140 XN=-XJ21+XJ+XJ32+EPS1 TRO06150 IF(XN.LT.0.) XN=XN+EPS2 TRO06160 N22=XN TRO06170 XN=-XJ12+XJ23+XJ+EPS1 TRO06180 IF(XN.LT.0.) XN=XN+EPS2 TRO06190 N24=XN TRO06200 XN=XJ21-XJ+XJ32+EPS1 TRO06210 IF(XN.LT.0.) XN=XN+EPS2 TRO06220 N26=XN TRO06230 XN=XJ12-XJ+XJ23+EPS1 TRO06240 IF(XN.LT.0.) XN=XN+EPS2 TRO06250 N28=XN TRO06260 XN=XJ21+XJ-XJ32+EPS1 TRO06270 IF(XN.LT.0.) XN=XN+EPS2 TRO06280 N30=XN TRO06290 XN=XJ12+XJ-XJ23+EPS1 TRO06300 IF(XN.LT.0.) XN=XN+EPS2 TRO06310 N32=XN TRO06320 XN=-XJ22-XJ+XJ32+XJ23+EPS1 TRO06330 IF(XN.LT.0.) XN=XN+EPS2 TRO06340 N34=XN TRO06350 XN=XJ12+XJ21+XJ22+XJ+EPS1 TRO06360 IF(XN.LT.0.) XN=XN+EPS2 TRO06370 N35=XN TRO06380 N35=N35+1 TRO06390 XN=XJ21+XJ+XJ32+EPS1 TRO06400 IF(XN.LT.0.) XN=XN+EPS2 TRO06410 N37=XN TRO06420 N37=N37+1 TRO06430 XN=XJ12+XJ+XJ23+EPS1 TRO06440 IF(XN.LT.0.) XN=XN+EPS2 TRO06450 N39=XN TRO06460 N39=N39+1 TRO06470 XN=-XJ+XJ11+XJ33+EPS1 TRO06480 IF(XN.LT.0.) XN=XN+EPS2 TRO06490 N42=XN TRO06500 XN=-XJ+XJ23+XJ12+EPS1 TRO06510 IF(XN.LT.0.) XN=XN+EPS2 TRO06520 N43=XN TRO06530 XN=XJ-XJ11+XJ33+EPS1 TRO06540 IF(XN.LT.0.) XN=XN+EPS2 TRO06550 N46=XN TRO06560 XN=XJ-XJ23+XJ12+EPS1 TRO06570 IF(XN.LT.0.) XN=XN+EPS2 TRO06580 N47=XN TRO06590 XN=XJ+XJ11-XJ33+EPS1 TRO06600 IF(XN.LT.0.) XN=XN+EPS2 TRO06610 N50=XN TRO06620 XN=XJ+XJ23-XJ12+EPS1 TRO06630 IF(XN.LT.0.) XN=XN+EPS2 TRO06640 N51=XN TRO06650 XN=-XJ13-XJ+XJ33+XJ12+EPS1 TRO06660 IF(XN.LT.0.) XN=XN+EPS2 TRO06670 N53=XN TRO06680 XN=XJ13+XJ+XJ23+XJ11+EPS1 TRO06690 IF(XN.LT.0.) XN=XN+EPS2 TRO06700 N55=XN TRO06710 N55=N55+1 TRO06720 XN=XJ+XJ11+XJ33+EPS1 TRO06730 IF(XN.LT.0.) XN=XN+EPS2 TRO06740 N57=XN TRO06750 N57=N57+1 TRO06760 XN=XJ+XJ23+XJ12+EPS1 TRO06770 IF(XN.LT.0.) XN=XN+EPS2 TRO06780 N58=XN TRO06790 N58=N58+1 TRO06800 GO TO 10 TRO06810 52 IF(XJ.GT.XJMAX) GO TO 120 TRO06820 N3=N3+1 TRO06830 N4=N4+1 TRO06840 N7=N7+1 TRO06850 N8=N8+1 TRO06860 N11=N11-1 TRO06870 N12=N12-1 TRO06880 N13=N13+1 TRO06890 N14=N14+1 TRO06900 N18=N18+1 TRO06910 N19=N19+1 TRO06920 N22=N22+1 TRO06930 N24=N24+1 TRO06940 N26=N26-1 TRO06950 N28=N28-1 TRO06960 N30=N30+1 TRO06970 N32=N32+1 TRO06980 N34=N34-1 TRO06990 N35=N35+1 TRO07000 N37=N37+1 TRO07010 N39=N39+1 TRO07020 N42=N42-1 TRO07030 N43=N43-1 TRO07040 N46=N46+1 TRO07050 N47=N47+1 TRO07060 N50=N50+1 TRO07070 N51=N51+1 TRO07080 N53=N53-1 TRO07090 N55=N55+1 TRO07100 N57=N57+1 TRO07110 N58=N58+1 TRO07120 10 K1=0. TRO07160 L1=-N13 TRO07170 IF(L1.GT.K1) K1=L1 TRO07180 L1=-N14 TRO07190 IF(L1.GT.K1) K1=L1 TRO07200 L1=N9 TRO07210 IF(N10.LT.L1) L1=N10 TRO07220 IF(N11.LT.L1) L1=N11 TRO07230 IF(N12.LT.L1) L1=N12 TRO07240 IF(N15.LT.L1) L1=N15 TRO07250 F1=1. TRO07260 S1=1. TRO07270 I1=K1+1 TRO07280 62 IF(I1.GT.L1) GO TO 64 TRO07290 I1M1=I1-1 TRO07300 NN1=(N9-I1M1)*(N10-I1M1)*(N11-I1M1)*(N12-I1M1) TRO07310 ND1=I1*(N13+I1)*(N14+I1)*(N15-I1M1) TRO07320 F1=-F1*FLOAT(NN1)/FLOAT(ND1) TRO07330 S1=S1+F1 TRO07340 I1=I1+1 TRO07350 GO TO 62 TRO07360 64 K2=0 TRO07370 L2=-N33 TRO07380 IF(L2.GT.K2) K2=L2 TRO07390 L2=-N34 TRO07400 IF(L2.GT.K2) K2=L2 TRO07410 L2=N29 TRO07420 IF(N30.LT.L2) L2=N30 TRO07430 IF(N31.LT.L2) L2=N31 TRO07440 IF(N32.LT.L2) L2=N32 TRO07450 IF(N35.LT.L2) L2=N35 TRO07460 F2=1. TRO07470 S2=1. TRO07480 I2=K2+1 TRO07490 70 IF(I2.GT.L2) GO TO 80 TRO07500 I2M2=I2-1 TRO07510 NN2=(N29-I2M2)*(N30-I2M2)*(N31-I2M2)*(N32-I2M2) TRO07520 ND2=I2*(N33+I2)*(N34+I2)*(N35-I2M2) TRO07530 F2=-F2*FLOAT(NN2)/FLOAT(ND2) TRO07540 S2=S2+F2 TRO07550 I2=I2+1 TRO07560 GO TO 70 TRO07570 80 K3=0 TRO07580 L3=-N53 TRO07590 IF(L3.GT.K3) K3=L3 TRO07600 L3=-N54 TRO07610 IF(L3.GT.K3) K3=L3 TRO07620 L3=N49 TRO07630 IF(N50.LT.L3) L3=N50 TRO07640 IF(N51.LT.L3) L3=N51 TRO07650 IF(N52.LT.L3) L3=N52 TRO07660 IF(N55.LT.L3) L3=N55 TRO07670 F3=1. TRO07680 S3=1. TRO07690 I3=K3+1 TRO07700 84 IF(I3.GT.L3) GO TO 90 TRO07710 I3M3=I3-1 TRO07720 NN3=(N49-I3M3)*(N50-I3M3)*(N51-I3M3)*(N52-I3M3) TRO07730 ND3=I3*(N53+I3)*(N54+I3)*(N55-I3M3) TRO07740 F3=-F3*FLOAT(NN3)/FLOAT(ND3) TRO07750 S3=S3+F3 TRO07760 I3=I3+1 TRO07770 GO TO 84 TRO07780 90 S2N=FLOG(N3+1)+FLOG(N4+1)+FLOG(N7+1)+FLOG(N8+1)+FLOG(N11+1)+FLOG(NTRO07820 112+1)+FLOG(N22+1)+FLOG(N24+1)+FLOG(N26+1)+FLOG(N28+1)+FLOG(N30+1)+TRO07830 2FLOG(N32+1)+FLOG(N42+1)+FLOG(N43+1)+FLOG(N46+1)+FLOG(N47+1)+FLOG(NTRO07840 350+1)+FLOG(N51+1) TRO07850 S2N=.5*S2N TRO07860 S2D=FLOG(N18+1)+FLOG(N19+1)+FLOG(N37+1)+FLOG(N39+1)+FLOG(N57+1)+FLTRO07870 1OG(N58+1) TRO07880 S2D=.5*S2D TRO07890 KM1=K1-1 TRO07900 KP1=K1+1 TRO07910 KM2=K2-1 TRO07920 KP2=K2+1 TRO07930 KM3=K3-1 TRO07940 KP3=K3+1 TRO07950 S2N=S2N+FLOG(N15-KM1)+FLOG(N35-KM2)+FLOG(N55-KM3) TRO07960 S2D=S2D+FLOG(KP1)+FLOG(KP2)+FLOG(KP3)+FLOG(N9-KM1)+FLOG(N10-KM1)+FTRO07970 1LOG(N11-KM1)+FLOG(N12-KM1)+FLOG(N13+KP1)+FLOG(N14+KP1)+FLOG(N29-KMTRO07980 22)+FLOG(N30-KM2)+FLOG(N31-KM2)+FLOG(N32-KM2)+FLOG(N33+KP2)+FLOG(N3TRO07990 34+KP2)+FLOG(N49-KM3)+FLOG(N50-KM3)+FLOG(N51-KM3) TRO08000 4+FLOG(N52-KM3)+FLOG(N53+KP3)+FLOG(N54+KP3) TRO08010 F=S2D-S2N TRO08020 IF (F.GT.80) GO TO 100 TRO08030 F=S2N/S2D TRO08040 IF((F.LT.1.01).AND.(F.GT.0.98)) GO TO 100 TRO08050 F=S1*S2*S3*EXP(S2N-S2D) TRO08060 GO TO 110 TRO08070 100 F=S1*S2*S3 TRO08080 IF(F) 102,112,104 TRO08090 102 F=ALOG(-F) TRO08100 F=-EXP(F+S2N-S2D) TRO08110 GO TO 110 TRO08120 104 F=ALOG(F) TRO08130 F=EXP(F+S2N-S2D) TRO08140 110 F=F*(2.*XJ+1.) TRO08150 L1=K1+K2+K3 TRO08190 K1=L1/2 TRO08200 K1=2*K1 TRO08210 IF(L1.NE.K1) F=-F TRO08220 S=S+F TRO08230 112 XJ=XJ+1. TRO08240 GO TO 52 TRO08250 120 C2N=FLOG(N1+1)+FLOG(N2+1)+FLOG(N5+1)+FLOG(N6+1)+FLOG(N9+1)+FLOG(N1TRO08290 10+1)+FLOG(N21+1)+FLOG(N23+1)+FLOG(N25+1)+FLOG(N27+1)+FLOG(N29+1)+FTRO08300 2LOG(N31+1)+FLOG(N41+1)+FLOG(N44+1)+FLOG(N45+1)+FLOG(N48+1)+FLOG(N4TRO08310 39+1)+FLOG(N52+1) TRO08320 C2N=.5*C2N TRO08330 C2D=FLOG(N16+1)+FLOG(N17+1)+FLOG(N36+1)+FLOG(N38+1)+FLOG(N56+1)+FLTRO08340 1OG(N59+1) TRO08350 C2D=.5*C2D TRO08360 F=C2D-C2N TRO08370 IF(F.GT.80.) GO TO 122 TRO08380 F=C2N/C2D TRO08390 IF((F.LT.1.01).AND.(F.GT.0.98)) GO TO 122 TRO08400 C9J=S*EXP(C2N-C2D) TRO08410 GO TO 130 TRO08420 122 IF(S) 124,50,126 TRO08430 124 S=ALOG(-S) TRO08440 C9J=-EXP(S+C2N-C2D) TRO08450 GO TO 130 TRO08460 126 S=ALOG(S) TRO08470 C9J=EXP(S+C2N-C2D) TRO08480 130 K=N9+N16+N36+N56-1 TRO08520 L=K/2 TRO08530 L=2*L TRO08540 IF (L.NE.K) C9J=-C9J TRO08550 RETURN TRO08560 END TRO08570 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