c********************************************************************* PROGRAM EIGEN c finds the eigenvalue and eigenfunction of a particle with charge Z1, c in a WS + spin-orbit potential with a nucleus with charge Z2. You c the choose the nodes, ang. momentum, etc., of the state you want. C C.A. Bertulani - version 22/jun/96 C c********************************************************************* CALL SCHRO END c SUBROUTINE SCHRO IMPLICIT REAL*8 (A-H,O-Z) COMMON /POTENTIAL/ U(2000),WSO(2000),UE(2000) & ,VC(2000) COMMON /WAVE_FUNCTIONS/ EIGWF(2000),NRMAX COMMON /PAR/ DR,DR2,PI,FOUR_PI,NPNTS COMMON /POT/ V0,V_SO,R0,AA,Z1Z2,ARED OPEN(12,FILE='EIGEN.IN', STATUS='UNKNOWN') OPEN(13,FILE='EIGEN.OUT', STATUS='UNKNOWN') WRITE(6,*) 'Bound-state in a WS+Coulomb+spin-orbit well' WRITE(13,*) 'Bound-state in a WS+Coulomb+spin-orbit well' c DR=mesh size, NPNTS=no. of points in the radial mesh, NRMAX=max. radius size WRITE(13,*)' ' READ(12,*) DR,NPNTS,NRMAX WRITE(13,*) 'DR,NPNTS,NRMAX' WRITE(13,*) DR,NPNTS,NRMAX c N_0=nodes of the Wav.Func., J_0=2*total ang. mom., L_0=orbital ang. mom. READ(12,*) N_0,J_0,L_0 WRITE(13,*)' ' WRITE(13,*) 'Nodes, Total Ang. Mom. J, Orbital Ang. Mom. L' WRITE(13,*) N_0,J_0,'/2',L_0 c V0=depth of central pot., V_SO=depth of spin-orbit pot., R0=radius of the c pot., AA=diffuseness of the pot., Z1Z2=product of charges of the nuclei, c ARED=reduced mass of the nuclei WRITE(13,*)' ' READ(12,*)V0,V_SO,R0,AA,Z1Z2,ARED WRITE(13,*) 'V0,V_SO,R0,AA,Z1Z2,ARED' WRITE(13,*) V0,V_SO,R0,AA,Z1Z2,ARED c see if you are out of bounds IF(NPNTS.GT.2000 .OR. NRMAX.GT.250) THEN WRITE(6,*) ' THE INPUT VALUES FOR NPNTS OR NRMAX & WERE INADEQUATE AND HAVE BEEN CHANGED' IF(NPNTS.GT.2000) NPNTS=2000 IF(NRMAX.GT.250) NRMAX=250 ENDIF c famous constants PI = 3.14159265359 FOUR_PI = 4.0D0 *PI * * Setup nuclear + coulomb + spin-orbit potential CALL POTENT * * Integrate Schroedinger equation CALL DIFF2(N_0,L_0,J_0,ENER) * * Output WRITE(6,*)' ' WRITE(13,*)' ' WRITE(6,104) N_0,L_0,J_0,ENER WRITE(13,104) N_0,L_0,J_0,ENER 104 FORMAT(1X,' N = ',I2,' L = ',I2,' J =',I2,'/2', & ' EIGEN-ENERGY = ',F6.2) WRITE(13,*)' ' WRITE(13,*)' r (fm) x WAVEFUNCTION' WRITE(13,*)' ' DO IR=1,NPNTS XR=IR*DR WRITE(13,*)XR,EIGWF(IR) ENDDO PRINT*,'Output of r x wavefunction in EIGEN.OUT' RETURN END c c==================================================================== SUBROUTINE POTENT * * Construct the potentials * IMPLICIT REAL*8 (A-H,O-Z) COMMON /POTENTIAL/ U(2000),WSO(2000),UE(2000) & ,VC(2000) COMMON /PAR/ DR,DR2,PI,FOUR_PI,NPNTS COMMON /POT/ V0,V_SO,R0,AA,Z1Z2,ARED * * Construct the Coulomb potential * DR2 = DR**2 B_N = 938.9D0*ARED B_MASS= (197.329D0)**2/(2.0D0*B_N) CHARGE = 1.439976D0*Z1Z2 DO IR=1,NPNTS XR = IR*DR XRI=(XR-R0)/AA IF(XRI.GT.30.0D0)XRI=30.0D0 EXT=DEXP( XRI ) c Woods-Saxon U(IR) = V0/( 1.0D0+EXT ) & /B_MASS*DR2 c Spin-orbit WSO(IR) = V_SO*EXT/((1.0D0+EXT)**2*XR*AA) & /B_MASS*DR2/2. UE(IR) = DR2/B_MASS c Coulomb potential of a uniformly distributed charge VC(IR)=CHARGE/XR IF(XR.LT.R0) VC(IR)=CHARGE/R0* & (1.5D0-0.5D0*(XR/R0)**2) VC(IR)=VC(IR)*DR2/B_MASS U(IR)=U(IR)+VC(IR) c add everything ENDDO RETURN END c c==================================================================== SUBROUTINE DIFF2(NODE_0,L_0,J_0,E) * * This subroutine solves the Schroedinger equation * IMPLICIT REAL*8 (A-H,O-Z) COMMON /WAVE_FUNCTIONS/ EIGWF(2000),NRMAX COMMON /POTENTIAL/ U(2000),WSO(2000),UE(2000) & ,VC(2000) DIMENSION WF0(2000), & VR(2000),VE(2000) COMMON /PAR/ DR,DR2,PI,FOUR_PI,NPNTS NODE = NODE_0 L = L_0 J = J_0 * * Find a good initial estimate of the eigenvalue * IITER = 0 NITER = 0 INODE = NODE DEM = 0.5D0 COR0 = 0.5D0 IF(L.EQ.0) COR0 =0.75D0 COR = COR0+1.0D0/(DFLOAT(INODE)+1.0D0)/24.0D0 * -1.0D0/(DFLOAT(INODE)+1.0D0)**2/48.0D0 A12 = 1.0D0/12.0D0 ACCUR = 0.1D0 ACCUR0= 0.1D-8 * * find bottom of the well * AL = DFLOAT( L*(L+1) ) AJ = DFLOAT(J)/2.0D0 ALS = AJ*(AJ+1.0D0)-AL-0.75D0 EMIN = 0.0D0 DO IR=1,NPNTS XR = IR VE(IR) = UE(IR) VR(IR) = U(IR) + AL/IR**2 + ALS*WSO(IR) IF( EMIN .GT. ( VR(IR)/VE(IR) ) ) EMIN = VR(IR)/VE(IR) ENDDO IF( NITER. GT. 6 ) GO TO 100 IF( E .LE. EMIN .OR. E .GE. 0.0D0 ) E = EMIN*0.5D0 * * find turning point and WKB eigen-energy 9 NITER=NITER+1 S=0.0D0 S1=0.0D0 DO 1 IR=1,NPNTS P2 = E*VE(IR)-VR(IR) IF ( P2.LE. 0.0D0) GO TO 1 P = DSQRT( P2 ) S = S+P S1 = S1+VE(IR)/P ITURN = IR 1 CONTINUE S = S-PI*( DFLOAT(INODE)+COR ) DE=-S/S1*2.0D0 108 IF( (E+DE) .GT. EMIN .AND. (E+DE) .LT. 0.0D0 ) GO TO 109 DE=DE*0.5D0 GO TO 108 109 E=E+DE IF( NITER .GT. 50 ) GO TO 1000 IF( DABS(DE) .GT. ACCUR ) GO TO 9 XR=ITURN*DR c WRITE(6,935) E,XR c 935 FORMAT(1X,' E(WKB) = ', F8.2,' Turning Point = ',F5.2) 100 IITER=IITER+1 110 CONTINUE * * Integrate outward from origin to the turning point * A = 0.0D0 B = 1.0D-10 WF0(1) = B D = VR(1)-VE(1)*E B = B*(1.0D0-A12*D) SEI = WF0(1)**2*VE(1) ICHECK=0 DO 101 IR=2,ITURN C = 2.0D0*B-A+WF0(IR-1)*D A = B B = C D = VR(IR)-VE(IR)*E C = C/(1.0D0-A12*D) SEI = SEI+C**2*VE(IR) WF0(IR) = C IF(IR.LT.5) GO TO 101 IF( WF0(IR)*WF0(IR-1) .LT. 0.0D0 .OR. WF0(IR) .EQ. 0.0D0) & ICHECK=ICHECK+1 101 CONTINUE IF(ICHECK.GT.INODE) THEN E=E*1.2D0 DO 1200 IR=1,NPNTS IF( VE(IR)*E .GT. VR(IR) ) ITURN=IR 1200 CONTINUE GO TO 110 ENDIF IF(ICHECK.LT.INODE) THEN E=E*0.9D0 DO 1201 IR=1,NPNTS IF( VE(IR)*E .GT. VR(IR)) ITURN=IR 1201 CONTINUE GO TO 110 ENDIF SEI = SEI/C**2-0.5D0*VE(ITURN) WFI = C WF1I=(274.0D0*WF0(ITURN)-600.0D0*WF0(ITURN-1) * +600.0D0*WF0(ITURN-2) * -400.0D0*WF0(ITURN-3)+150.0D0*WF0(ITURN-4) * -24.0D0*WF0(ITURN-5))/(120.0D0*DR) DWFI = WF1I/WFI SEI = SEI/DR+A12*DWFI*DR**2 * * Integrate inward from infinity to the turning point * D = VR(NPNTS)-VE(NPNTS)*E A = DEXP(-DFLOAT(NPNTS)*DSQRT(DABS(D))) IF ( A .LT. 1.0D-12) A =0.0D0 WF0(NPNTS) = A A = A*(1.0D0-A12*D) D = VR(NPNTS-1)-VE(NPNTS-1)*E B = DEXP(-DFLOAT(NPNTS-1)*DSQRT(DABS(D))) IF ( B .LE. 1.0D-12 ) B = 1.0D-12 WF0(NPNTS-1) = B SEO = B**2*VE(NPNTS-1) B = B*(1.0D0-A12*D) DO 201 IR=NPNTS-2,ITURN,-1 C = 2.0D0*B-A+WF0(IR+1)*D A = B B = C D = VR(IR)-VE(IR)*E C = C/(1.0D0-A12*D) SEO = SEO+C**2*VE(IR) WF0(IR)=C 201 CONTINUE WFO = C WF1O=(-274.0D0*WF0(ITURN)+600.0D0*WF0(ITURN+1) * -600.0D0*WF0(ITURN+2) * +400.0D0*WF0(ITURN+3)-150.0D0*WF0(ITURN+4) * +24.0D0*WF0(ITURN+5))/(120.0D0*DR) DWFO = WF1O/WFO SEO = SEO/WFO**2-0.50D0*VE(ITURN) SEO = SEO/DR-A12*DWFO*DR**2 SE = DWFI-DWFO * * Compute the correction to the eigenvlue * DE = SE/(SEI+SEO) 221 IF( DABS(DE) .LT. DEM ) GO TO 222 DE = DE*0.8D0 GO TO 221 222 E1 = E+DE IF( E1 .GE. 0.0D0 .OR. E1 .LE. EMIN ) THEN DE = DE*0.9D0 GO TO 222 ENDIF E = E1 IF( IITER .GT. 50 ) GO TO 1000 XTEST = DABS(DE)/(DABS(E)+1.D-25) IF( XTEST .GT. ACCUR0 ) GO TO 100 * * Match the wavefunction * SE = WFI/WFO DO 303 IR=ITURN,NPNTS 303 WF0(IR)=WF0(IR)*SE * * Normalize th waveunction * SN = 0.0D0 DX = DR/15.0D0 AX = DX*7.0D0 NUP= NPNTS IF( MOD(NPNTS,2) .NE.0 ) NUP=NPNTS-1 DO 301 IR=NUP,1,-1 SN = SN+WF0(IR)**2*VE(IR)*AX AX = DR+DX DX =-DX 301 CONTINUE SN=1.0D0/DSQRT(SN) DO 302 IR=1,NPNTS AX = WF0(IR)*SN*DSQRT(VE(IR)) XR = IR*DR WF0(IR) = AX/XR 302 EIGWF(IR) = AX RETURN 1000 WRITE(6,*) 'THE ITERATION PROCCES FAILED NITER=',NITER,' & IITER=' * ,IITER,' NODE=',NODE,' E=',E STOP END