This file contains the fortran program "Jellyrpa" followed by a data file to calculate the response of the cluster Na_20^+. PROGRAM JELLYRPA C Written by G. Bertsch, based on a nuclear physics program C described in Nuclear Physics A243 (1975), p. 507 C Corrected Sept. 13, 1990; thanks to Sasumo Saito and Stuart Pollack C July 5, 1995; thanks to CC. Chen C Parameters: C NRMAX - maximum number of points in spatial grid C NHMAX - maximum number of hole wave functions C MATMAX - maximum dimensionality of matrix C PARAMETER (NRMAX=100,NHMAX=50,MATMAX=100) C units of length are a.u. (1 a.u.= 0.53 A) C units of energy are eV externally and Ry internally. C Principal Variables: C DEL - step size in radial mesh C NGRID - number of points in wave function C N=NGRID/NX - number of spatial points in response matrix C UH - radial wave function of orbital C LH,NODE,NOCC- angular momentum, radial nodes, and C occupation number of orbital C L,EX - multipolarity and energy (eV) of response C G - independent particle response matrix C B - rpa polarizability matrix (1+VG)**(-1) C COMPLEX B(MATMAX,MATMAX),PIFREE,PIRPA DIMENSION VEXT(MATMAX) COMMON/DENSITY/RHO(0:NRMAX) COMMON/PARAM/DEL,DEL2,N,NX,NGRID,DD COMMON/WFN/LH(NHMAX),EH(NHMAX),NODE(NHMAX),OCC(NHMAX), & UH(NHMAX,NRMAX),NH COMMON/POTEN/V(NRMAX),EL COMMON/POTN2/VRES(MATMAX,MATMAX) C Number of mesh points and approximate spacing for response c Initialize clebsch routine CALL FACTOR C Set up wave functions for ground state READ(5,*) NGRID,N,NX WRITE(6,81) NGRID,N,NX if(nx.ne.1.and.ngrid.ne.n*nx) then print*,' NGRID does not equal NX*N' stop endif 81 FORMAT(1X,I4,' grid points in wfn;',I5,' in matrix;',I3, & ' relative step size') CALL STATIC(Z) DEL2=NX*DEL c Multipolarity,starting energy, ending energy, and energy step c for response function READ(5,*) L,EX,EXM,DEX,GAM CALL VRESID(L) PRINT*,' L= ',L,' GAM= ', GAM C Make multipole external field CALL EXTF(VEXT,L) C Loop on energy WRITE(*,'('' Energy Free response'', & '' RPA Response'')') WRITE(*,'('' (eV) Real Imaginary'', & '' Real Imaginary'')') 11 IF(EX.GT.EXM) GO TO 40 CALL RESPONSE(L,EX/13.6,GAM/13.6,B,VEXT,PIFREE,PIRPA) IF(EX.EQ.0.0.AND.L.EQ.1) THEN ALPHA = REAL(PIRPA)*8*3.1416/3.0 WRITE(*,87) ALPHA 87 FORMAT(' Polarizability is ',e12.4,' a.u.**3') ENDIF WRITE(*,'(F11.5,4E12.4)') EX,PIFREE,PIRPA SUMF=SUMF+AIMAG(PIFREE)*EX*DEX/3.1416 SUMRPA=SUMRPA+AIMAG(PIRPA)*EX*DEX/3.1416 EX=EX+DEX GO TO 11 40 CALL SUMRULE(VEXT,VALUE,L) PRINT*,' TOTAL STRENGTH IN FREE RESPONSE, RPA RESPONSE, & AND SUM RULE' WRITE(*,'(21X,E12.4,E14.4,E15.4)') SUMF,SUMRPA,VALUE*(13.6)**2 IF(L.EQ.1) THEN SUMRPA=SUMRPA*32*3.1416**3*0.529**2/3.0/137.0/13.6 WRITE(*,85) SUMRPA 85 FORMAT(' energy-integrated photon cross section, (ev-A^2) ', & E12.4) ENDIF END SUBROUTINE STATIC (Z) C spherical jellium model PARAMETER (NRMAX=100,NHMAX=50) REAL*8 U,SX,ETRIAL,EMAX,EMIN DIMENSION U(NRMAX),VIN(0:NRMAX),VOUT(0:NRMAX) DIMENSION V1(NRMAX),V2(NRMAX) COMMON/DENSITY/RHO(0:NRMAX) COMMON/PARAM/DEL,DEL2,N,NX,NGRID,DD COMMON/WFN/LH(NHMAX),EH(NHMAX),NODE(NHMAX),OCC(NHMAX), & UH(NHMAX,NRMAX),NH COMMON/POTEN/V(NRMAX),EL DATA V0,A0,ALPH/-0.40,2.0,0.2/ C RS =(3*density/ 4 * pi )^(1/3), in A.U.=0.529 A. C Z is charge of jellium READ(5,*) RS,Z,ITER WRITE(*,81)RS,Z,ITER 81 FORMAT(' RS,Z,ITER= ',F5.2,F5.1,I3) RR=RS*Z**(0.33333) PRINT*,' R= ',RR DEL=RR*2.0/NGRID DD=DEL**2 C initial trial potential is Woods-Saxon DO 6 I = 1, NGRID R=DEL*I TEMP=(R-RR)/A0 IF(TEMP.GT.30) TEMP=30 EX=EXP(TEMP) V(I)=V0*(1+ALPH*(R/RR)**2)/(1+EX) 6 CONTINUE NH=0 ZE=0.0 C Read in quantum numbers of occupied states 1 NH=NH+1 READ(5,*) LH(NH),NODE(NH),NOCC OCC(NH)=NOCC IF (NOCC.LT.0) OCC(NH)=2.0*(2.0*LH(NH)+1) ZE=ZE+OCC(NH) WRITE(6,*)LH(NH),NODE(NH),OCC(NH) IF (LH(NH).GE.0) GO TO 1 NH=NH-1 PRINT*,'NET CHARGE = ' ,Z-ZE C iterate Hartree loop DO 33 KK=1,ITER PRINT*,'iteration number ',KK DO 3 I=1,NGRID 3 RHO(I)=0.0 DO 32 IH=1,NH ALL = LH(IH)*(LH(IH)+1) C Energy of state is found by bisection EMIN=-1.2 EMAX=0.0 DO 7 K = 1,40 C Integrate schroedinger equation ETRIAL=(EMIN+EMAX)/2.0 U(1)=(0.1)**(LH(IH)+1) U(2)=(0.2)**(LH(IH)+1) ND=0 DO 8 I =2,NGRID-1 R=DEL*I VV=V(I)+ALL/R**2 SX=DD*(VV-ETRIAL) U(I+1)=(2+SX)*U(I)-U(I-1) IF(U(I+1)*U(I).LT.0)ND=ND+1 8 CONTINUE IF(ND.GT.NODE(IH)) THEN EMAX=ETRIAL ELSE EMIN=ETRIAL END IF 7 CONTINUE SUM=0.0 DO 9 I=1,NGRID SUM=SUM+U(I)**2*DEL 9 CONTINUE SUM=SQRT(SUM) DO 10 I=1,NGRID R=DEL*I UH(IH,I)=U(I)/SUM RHO(I)=RHO(I)+OCC(IH)*(UH(IH,I)/R)**2/12.56 10 CONTINUE EH(IH)=ETRIAL IF (KK.LE.ITER) THEN WRITE(*,'('' L,NODE,E='',2I3,F10.5,''eV'')') & LH(IH),NODE(IH),ETRIAL*13.6 END IF 32 CONTINUE IF (KK.EQ.ITER) GO TO 33 C construct electron-electron potential C by the Green's function 1/r_> VIN(0)=0.0 VOUT(NGRID)=0.0 T1=0.0 T3=0.0 DO 27 I=1,NGRID T2=RHO(I)*12.56*(DEL*I)**2 J=NGRID-I VIN(I)=VIN(I-1)+(T1+T2)/2.0 T1=T2 T4=RHO(J)*12.56*(DEL*J) VOUT(J)=VOUT(J+1)+(T3+T4)/2.0 T3=T4 27 CONTINUE DO 28 I=1,NGRID RSS=0.620/RHO(I)**0.33333 R=DEL*I C jellium potential is vcore VCORE=-2*Z/R IF(I.LT.NGRID/2) VCORE=-2*Z/RR*(1.5-0.5*(R/RR)**2) VEE=2.0*DEL*(VIN(I)/R+VOUT(I)) C vexc is from Gunnarsson and Lundqvist, P.R. B13 4274 (1976). VEC=-1.222/RSS-0.0666*ALOG(1.0+11.4/RSS) V1(I)=VCORE+VEE+VEC-V(I) c convergence parameter 0.28 must be decreased for larger clusters V(I)=V(I)+0.28*V1(I) 28 CONTINUE 33 CONTINUE END C SUBROUTINE RESPONSE(L,EX,GAM,B,VEXT,PIFREE,PIRPA) C Construct independent particle response PARAMETER (NRMAX=100,NHMAX=50,MATMAX=100) COMPLEX A(MATMAX,MATMAX),G(MATMAX,MATMAX), & B(MATMAX,MATMAX),UP(NRMAX),WP(NRMAX),UP1(NRMAX), & WP1(NRMAX),DRHOF(MATMAX),DRHORPA(MATMAX),PIFREE,PIRPA & ,VNEW(NRMAX) DIMENSION VNN(MATMAX),VNP(MATMAX),VEXT(MATMAX) COMMON/DENSITY/RHO(0:NRMAX) COMMON/PARAM/DEL,DEL2,N,NX,NGRID,DD COMMON/WFN/LH(NHMAX),EH(NHMAX),NODE(NHMAX),OCC(NHMAX), & UH(NHMAX,NRMAX),NH COMMON/POTEN/V(NRMAX),EL COMMON/POTN2/VRES(MATMAX,MATMAX) DO 13 M1=1,N DO 12 M2=1,N 12 G(M1,M2)=0.0 13 CONTINUE C LOOP OVER HOLES DO 15 I2=1,NH C Loop over particle angular momentum (twice j) LMIN=IABS(LH(I2)-L) LMAX=LH(I2)+L C Angular momentum triangle condition DO 14 LP=LMIN,LMAX CALL CLEBSCH(2*LH(I2),2*LP,L*2,0,0,0,T1) C factor of 2 below is for two spin states TEMP=OCC(I2)*(2*LP+1)*T1**2/(12.566*(2*L+1)) IF (TEMP.EQ.0.0) go to 14 C Particle Green's function CALL GREEN(UP,WP,EX+EH(I2),GAM,LP,i2) CALL GREEN(UP1,WP1,-EX+EH(I2),-GAM,LP,i2) C Independent particle response DO 17 M1=1,N MA=M1*NX DO 16 M2 = M1,N MB=M2*NX 16 G(M1,M2)=G(M1,M2)+TEMP*UH(I2,MA) & *UH(I2,MB)*(UP(MA)*WP(MB)+UP1(MA)*WP1(MB)) 17 CONTINUE 14 CONTINUE 15 CONTINUE C Multiply G by the residual interaction V DO 22 M1=1,N DO 27 M2=1,N G(M2,M1)=G(M1,M2) A(M1,M2)=0.0 DO 28 K=1,N 28 A(M1,M2)=A(M1,M2)+G(M1,K)*VRES(K,M2)*DEL2 IF(M1.EQ.M2) A(M1,M2)=1.0+A(M1,M2) 27 CONTINUE 22 CONTINUE C Invert 1+VG CALL MATR(N,A,B) C Evaluate response to field VEXT PIFREE=0.0 DO 35 M1=1,N DRHOF(M1)=0.0 DO 34 M2=1,N 34 DRHOF(M1)=DRHOF(M1)+G(M1,M2)*VEXT(M2)*DEL2 C Free response to field VEXT 35 PIFREE=PIFREE+VEXT(M1)*DRHOF(M1)*DEL2 PIRPA=0.0 DO 36 M1=1,N DRHORPA(M1)=0.0 DO 37 M2 = 1,N 37 DRHORPA(M1)=DRHORPA(M1)+B(M1,M2)*DRHOF(M2) C RPA response to field VEXT 36 PIRPA=PIRPA+VEXT(M1)*DRHORPA(M1)*DEL2 81 FORMAT(1X,10F8.4) END c SUBROUTINE VRESID(L) PARAMETER (NRMAX=100,MATMAX=100) COMMON/DENSITY/RHO(0:NRMAX) COMMON/PARAM/DEL,DEL2,N,NX,NGRID,DD COMMON/POTN2/VRES(MATMAX,MATMAX) DO 1 I=1,N RI=I*DEL2 DO 2 J=1,I RJ=J*DEL2 C Note: 4*pi*e*e=25.13 VRES(I,J)=DEL2*25.13*RJ**L/(RI**(L+1)*(2.0*L+1.0)) VRES(J,I)=VRES(I,J) 2 CONTINUE RH=AMAX1(RHO(I*nx),0.000001) RSS=0.6204/RH**0.3333 VRES(I,I)=VRES(I,I)-(0.407/RSS+0.253/(RSS+11.4))/RH/RI**2 1 CONTINUE RETURN END SUBROUTINE EXTF(F,L) C pure multipole field DIMENSION F(1) COMMON/PARAM/DEL,DEL2,N,NX,NGRID,DD IF(L.EQ.0) THEN K=2 ELSE K=L END IF DO 1 M1=1,N R=DEL2*M1 F(M1)=R**K 1 CONTINUE END SUBROUTINE FACTOR COMMON/BL1/FACLOG(130) FACLOG(1)=0 FACLOG(2)=0 FN=1 DO 95 M=3,130 FN=FN+1 95 FACLOG(M)=FACLOG(M-1)+ALOG(FN) RETURN END SUBROUTINE CLEBSCH(J1,J2,J3,M1,M2,M3,ANS) C calculates Clebsch-Gordon coefficients C Definition found in: "Elementary Theory of Angular Momentum", M.E. Ros C Equation 3.19 (John Wiley) COMMON/BL1/FACLOG(130) IF(M1+M2.NE.M3)GO TO 130 IF(J3.LE.J1+J2.AND.J3.GE.ABS (J1-J2))GO TO 110 GO TO 130 110 IA=(J1+J2-J3+2)/2.0 IB=(J3+J1-J2+2)/2.0 IC=(J3+J2-J1+2)/2.0 ID=(J1+J2+J3+4)/2.0 IE=(J1+M1+2)/2.0 IF=(J1-M1+2)/2.0 IG=(J2+M2+2)/2.0 IH=(J2-M2+2)/2.0 II=(J3+M3+2)/2.0 IJ=(J3-M3+2)/2.0 FIRST=0.5*(ALOG(J3+1.0)+FACLOG(IA)+FACLOG(IB)+FACLOG(IC) & +FACLOG(IE)+FACLOG(IF)+FACLOG(IG)+FACLOG(IH)+FACLOG(II) & +FACLOG(IJ)-FACLOG(ID)) PART1=EXP (FIRST) NUMIN1= ABS(AMIN1(((J3-J2+M1)/2.0),((J3-J1-M2)/2.0),0.0))+1.0 NUMAX1=MIN1(((J1+J2-J3)/2.0),((J1-M1)/2.0),((J2+M2)/2.0))+1 KB=(J3-J2+M1)/2.0+1.0 KC=(J3-J1-M2)/2.0+1.0 FF=0 DO 120 NUM1=NUMIN1,NUMAX1 NU=NUM1-1 SECOND=-(FACLOG(NUM1)+FACLOG(IA-NU)+FACLOG(IF-NU) & +FACLOG(IG-NU)+FACLOG(KB+NU)+FACLOG(KC+NU)) CON=((-1)**NU)*EXP (SECOND) FF=FF+CON 120 CONTINUE ANS=PART1*FF GO TO 140 130 ANS=0 140 CONTINUE RETURN END SUBROUTINE GREEN(UP,WP,E,GAM,L,i2) PARAMETER (NRMAX=100,NHMAX=50) REAL ALL,R,Z,PK INTEGER I,K,NMAX1,NP COMPLEX UP(1),WP(1),ONE,EI,WW,ECOMP,S(NRMAX) COMMON/DENSITY/RHO(0:NRMAX) COMMON/WFN/LH(NHMAX),EH(NHMAX),NODE(NHMAX), & UH(NHMAX,NRMAX),NH COMMON/POTEN/V(NRMAX),EL COMMON/PARAM/DEL,DEL2,N,NX,NMAX,DD ONE=(1.,0.) EI=(0.,1.) ECOMP=E*ONE+GAM*EI/2. ALL=L*(L+1) NMAX1=NMAX+1 C Initialize regular solution UP(1)=(0.1)**(L+1) UP(2)=(0.2)**(L+1) DO 10 I=2,NMAX1 K=I-1 R=DEL*(I-1) S(I-1)=DD*(V(K)+ALL/R**2-ECOMP) IF(I.LT.3.OR.I.GT.NMAX) GO TO 10 UP(I)=(2.+S(I-1))*UP(I-1)-UP(I-2) 10 CONTINUE Z=-S(NMAX-1) IF(Z.GT.0.0) THEN C Initialize irregular solution to outgoing wave PK=SQRT(Z) WP(NMAX-1)=(1.0,0.0) WP(NMAX)=ONE*(1.-PK**2./2.)+EI*(PK) ELSE IF (Z.GT.EH(I2)) THEN PK=-SQRT(-Z) WP(NMAX-1)=1.0 WP(NMAX)=EXP(PK) ELSE WP(NMAX)=(0.0,0.0) WP(NMAX-1)=0.001 ENDIF DO 30 K=3,NMAX I=NMAX1-K WP(I)=(2.+S(I+1) )*WP(I+1)-WP(I+2) 30 CONTINUE C Wronskian NP=NMAX/2 WW=-(UP(NP)*WP(NP+1)-UP(NP+1)*WP(NP))/DEL DO 50 I=1,NMAX WP(I)=WP(I)/WW 50 CONTINUE RETURN END SUBROUTINE SUMRULE(F,VALUE,J) C calculates energy-weighted sum rule PARAMETER (NRMAX=100,MATMAX=100) DIMENSION F(1),FP(MATMAX) COMMON/DENSITY/RHO(0:NRMAX) COMMON/PARAM/DEL,DEL2,N,NX,NGRID,DD DO 3 I=1,N IF(I.EQ.1) THEN FL=0.0 ELSE FL=F(I-1) ENDIF IF(I.EQ.N) THEN FG=(2*F(I)-F(I-1)) ELSE FG=F(I+1) ENDIF FP(I)=(FG-FL)/2.0/DEL2 3 CONTINUE DO 5 I=1,N R=I*DEL2 SS=SS+((FP(I)**2+(F(I)/R)**2*J*(J+1))*RHO(I*NX))* & R**2*12.56*DEL2 5 CONTINUE VALUE=SS/4./3.1416 RETURN END SUBROUTINE MATR(NMAX,C,D) PARAMETER (MATMAX=100) IMPLICIT COMPLEX (A-H,O-Z) DIMENSION U(MATMAX),V(MATMAX),C(MATMAX,MATMAX),D(MATMAX,MATMAX) DET=1.0 DO 1 N=1,NMAX DO 1 M=1,NMAX D(N,M)=0.0 IF(N.EQ.M) D(N,M)=1.0 1 CONTINUE DO 10 N=1,NMAX T=C(N,N) IF(CABS(T).LE.1.E-10) GO TO 3 GO TO 6 3 J=N 4 IF(J.GT.NMAX) GO TO 11 J=J+1 T=C(N,J) IF(CABS(T).LE.1.E-10) GO TO 4 DO 5 K=1,NMAX U(K)=C(N,K) V(K)=D(N,K) C(N,K)=C(J,K) D(N,K)=D(J,K) C(J,K)=U(K) 5 D(J,K)=V(K) 6 DO 8 K=1,NMAX IF(K.EQ.N) GO TO 8 A=C(K,N)/C(N,N) DO 7 L=1,NMAX C(K,L)=C(K,L)-A*C(N,L) D(K,L)=D(K,L)-A*D(N,L) 7 CONTINUE 8 CONTINUE B=C(N,N) DET=DET*B DO 10 M=1,NMAX C(N,M)=C(N,M)/B D(N,M)=D(N,M)/B 10 CONTINUE RETURN 11 WRITE (6,'('' Matrix not invertable '')') RETURN END 50 45 1 3.93 21 15 0 0 -1 1 0 -1 2 0 -1 0 1 -1 -1 0 0 1 0.00 15.0 0.2 0.4