********************************************************************** PROGRAM EIGEN * finds the eigenvalue and eigenfunction of a particle with charge Z1, * in a WS + spin-orbit potential * * WS = V_0 [ 1 - F_so (l.s) (r_0/r) d/dr ] [ 1 + exp((r-R_0)/a) ]^(-1) * * with a nucleus with charge Z2. * * One chooses the nodes, angular momentum, etc., of the wanted state. *--------------- * DR = mesh size, NPNTS = no. of points in the radial mesh * RMAX = max. radius size, N_0 = nodes of the Wav.Func. (exclude origin) * J_0 = 2*total ang. mom., L_0 = orbital ang. mom. * V0 = depth of central pot., FSO = spin-orbit length, R0 = radius of * the pot., AA = diffuseness of the pot. Z1,Z2 = charges of the nuclei, * A1,A2 = masses of the nuclei (in nucleon mass units) ********************************************************************* * CALL SCHRO END * SUBROUTINE SCHRO IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(NLG=10000) COMMON /POTENTIAL/ U(NLG),WSO(NLG),UE(NLG) & ,VC(NLG) COMMON /WAVE_FUNCTIONS/ EIGWF(NLG),NRMAX COMMON /PAR/ DR,DR2,PI,FOUR_PI,NPNTS COMMON /POT/ V0,FSO,R0,AA,Z1Z2,ARED OPEN(12,FILE='EIGEN.IN', STATUS='UNKNOWN') OPEN(16,FILE='GSWF.IN', STATUS='UNKNOWN') OPEN(13,FILE='EIGEN.TXT', STATUS='UNKNOWN') WRITE(6,*) 'Bound-state in a WS+Coulomb+spin-orbit well' WRITE(13,*) 'Bound-state in a WS+Coulomb+spin-orbit well' * DR = mesh size, NPNTS = no. of points in the radial mesh * RMAX = max. radius size WRITE(13,*)' ' CALL RSHARPS(12) READ(12,*) NPNTS,RMAX NRMAX=INT(RMAX) DR=RMAX/NLG WRITE(13,*) 'NPNTS,RMAX' WRITE(13,*) NPNTS,RMAX * N_0 = nodes of the Wav.Func., J_0 = 2*total ang. mom. * L_0=orbital ang. mom. CALL RSHARPS(12) 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 * V0 = depth of central pot., FSO = spin-orbit length, * R0 = radius of the pot., AA = diffuseness of the pot. * Z1, Z2 = charges of the nuclei * A1, A2 = masses of the nuclei (in nucleon mass units) * WRITE(13,*)' ' CALL RSHARPS(12) READ(12,*)V0,FSO,R0,AA,Z1,Z2,A1,A2 WRITE(13,*) 'V0,FSO,R0,AA,Z1,Z2,A1,A2' WRITE(13,*) V0,FSO,R0,AA,Z1,Z2,A1,A2 Z1Z2=Z1*Z2 ARED=A1*A2/(A1+A2) * see if you are out of bounds IF(NPNTS.GT.NLG .OR. NRMAX.GT.250) THEN WRITE(6,*) ' THE INPUT VALUES FOR NPNTS OR NRMAX & WERE INADEQUATE AND HAVE BEEN CHANGED' IF(NPNTS.GT.NLG) NPNTS=NLG IF(NRMAX.GT.250) NRMAX=250 ENDIF * 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,' MeV') WRITE(13,*)' ' WRITE(13,*)' r (fm) x WAVEFUNCTION' WRITE(13,*)' ' DO IR=1,NPNTS XR=IR*DR WRITE(13,*)XR,EIGWF(IR) WRITE(16,*)XR,EIGWF(IR) ENDDO PRINT*,'- Output of r x wavefunction in EIGEN.TXT and GSWF.IN' PRINT*,'- GSWF.IN is prepared for use as input of wavefunction by' PRINT*,' another program' RETURN END * *===================================================================== SUBROUTINE POTENT PARAMETER(NLG=10000) * * Construct the potentials * IMPLICIT REAL*8 (A-H,O-Z) COMMON /POTENTIAL/ U(NLG),WSO(NLG),UE(NLG) & ,VC(NLG) COMMON /PAR/ DR,DR2,PI,FOUR_PI,NPNTS COMMON /POT/ V0,FSO,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 ) * Woods-Saxon U(IR) = V0/( 1.0D0+EXT ) & /B_MASS*DR2 * Spin-orbit R0R=1.25 WSO(IR) = V0*FSO*R0R*EXT/((1.0D0+EXT)**2*XR*AA) & /B_MASS*DR2/2. UE(IR) = DR2/B_MASS * 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) * add everything ENDDO RETURN END * *===================================================================== SUBROUTINE DIFF2(NODE_0,L_0,J_0,E) PARAMETER(NLG=10000) * * This subroutine solves the Schroedinger equation * IMPLICIT REAL*8 (A-H,O-Z) COMMON /WAVE_FUNCTIONS/ EIGWF(NLG),NRMAX COMMON /POTENTIAL/ U(NLG),WSO(NLG),UE(NLG) & ,VC(NLG) DIMENSION WF0(NLG), & VR(NLG),VE(NLG) 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 * WRITE(6,935) E,XR * 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 * *==================================================================== * subroutine rsharps(iunit) * * allows comments in input files * character*1 isharp 1 read(iunit,'(A1)')isharp if(isharp.eq."#")goto 1 backspace iunit return end *********************************************************************