c VERSION with 1P1H +pairing cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccc c VERSION TESTED WITH HOUDA FOR 1P1H+2P2H c c'est impeccable! c sauf le calcul des elements la difference est 1p et 1n c(mais ca donne presque rien) c c +Lanczos from Etienne Caurier c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c the time-reversal symmetry is not conserved for some 4p4h terms c there are minor changes in the densities and connected quantities ccc c==================================================================== c==================================================================== c==================================================================== c==================================================================== c PN BLOCK PN BLOCK PN BLOCK c==================================================================== c==================================================================== c==================================================================== c==================================================================== C======================================================================= c Calculate matrix element between the single states C for pn part; c N1,N3 MUST be chosen neutron and N2,N4 proton states C C======================================================================= DOUBLE PRECISION FUNCTION EMADPN(N1,N2,N3,N4) INCLUDE 'All.f' c IF (N1.LT.0) THEN JZA = SS2K(-N1,1,1) IPA = SSIP(-N1,1,1) ELSE JZA = SS2K(N1,2,1) IPA = SSIP(N1,2,1) ENDIF IF (N2.LT.0) THEN JZB = SS2K(-N2,1,2) IPBB = SSIP(-N2,1,2) ELSE JZB = SS2K(N2,2,2) IPBB = SSIP(N2,2,2) ENDIF IF (N3.LT.0) THEN JZC = SS2K(-N3,1,1) IPCC = SSIP(-N3,1,1) ELSE JZC = SS2K(N3,2,1) IPCC = SSIP(N3,2,1) ENDIF IF (N4.LT.0) THEN JZD = SS2K(-N4,1,2) IPD = SSIP(-N4,1,2) ELSE JZD = SS2K(N4,2,2) IPD = SSIP(N4,2,2) ENDIf c ETMP1 = 0D0 ETMP0 = 0D0 c IF (JZA + JZB - JZC - JZD.EQ.0) THEN C In fact, we don't need this check because the many-body basis C is the constructed in such a way that this condition is C verified automatically and the same for the CONDITION OF PARITY c c we calculate overlap of any two states. It should be positive--assumption c we multiply proton functions by sign. c important for off diagonal elements only c OVERAB=0.D0 OVERCD=0.D0 DO IH=1, NGHP DO IL=1,NGL OVERAB=OVERAB+(SPF(IH,IL,N1,1)*SPF(IH,IL,N2,2)+ $ SMF(IH,IL,N1,1)*SMF(IH,IL,N2,2))* $ GH(IH)*GL(IL) OVERCD=OVERCD+(SPF(IH,IL,N3,1)*SPF(IH,IL,N4,2)+ $ SMF(IH,IL,N3,1)*SMF(IH,IL,N4,2))* $ GH(IH)*GL(IL) END DO END DO ZNAKB=DSIGN(1.d0,OVERAB) ZNAKD=DSIGN(1.D0,OVERCD) C DO IH = 1, NGHP DO IL = 1, NGL PHIPA = SPF(IH,IL,N1,1) PHIPB = SPF(IH,IL,N2,2)*ZNAKB PHIPC = SPF(IH,IL,N3,1) PHIPD = SPF(IH,IL,N4,2)*ZNAKD PHIMA = SMF(IH,IL,N1,1) PHIMB = SMF(IH,IL,N2,2)*ZNAKB PHIMC = SMF(IH,IL,N3,1) PHIMD = SMF(IH,IL,N4,2)*ZNAKD C C ETMP1-T=1 ISOSPIN; ETMP0-T=0 ISOSPIN; C ETMP1 = ETMP1+(PHIPA*PHIMB-PHIMA*PHIPB)* & (PHIPC*PHIMD-PHIMC*PHIPD) & *GH(IH)*GL(IL) C ETMP0 = ETMP0+((PHIPA*PHIMB+PHIMA*PHIPB)* & (PHIPC*PHIMD+PHIMC*PHIPD) + & 2*PHIPA*PHIPB*PHIPC*PHIPD+ & 2*PHIMA*PHIMB*PHIMC*PHIMD) & *GH(IH)*GL(IL) C ENDDO ENDDO c ETMP1 = ETMP1*FHTDA1/BZP2/(2D0*PI)/2D0 !additional 1/2 as compared to pp and nn ETMP0 = ETMP0*FHTDA0/BZP2/(2D0*PI)/2D0 C c C 1/(2*Pi) = 2*Pi/(SQRT(2Pi)**4) C Pi comme from the right formule C SQRT(2Pi) comme from the angle part that we "forget" to include in C the function (theta or phi in polar coordinate, axial symmetry) C 4 function A B C D (N1 N2 N3 N4) give 4Pi**2 C 2 comme from the symmetry on axe z: we have taken only z>0 C C Smooth on two-body state ETMP1 = ETMP1*WFU(N1,1)*WFU(N2,2)*WFU(N3,1)*WFU(N4,2) ETMP0 = ETMP0*WFU(N1,1)*WFU(N2,2)*WFU(N3,1)*WFU(N4,2) END IF EMADPN = ETMP0+ETMP1 RETURN END C======================================================================= c Determining the pn pairs by the calcualtion of overlaps C Done only for K>0 c======================================================================= SUBROUTINE PARTNER () INCLUDE 'All.f' C C first hole levels C DO IHN=1,2*NHLD(1)-1,2 SOVER=0.D0 DO IHP=1,2*NHLD(2)-1,2 AOVER=0.D0 OVER=0.D0 DO IH=1,NGHP DO IL=1,NGL OVER=OVER+(SPF(IH,IL,-IHN,1)*SPF(IH,IL,-IHP,2)+ $ SMF(IH,IL,-IHN,1)*SMF(IH,IL,-IHP,2))* $ GH(IH)*GL(IL) END DO END DO AOVER=DABS(OVER) IF (AOVER.GT.SOVER) THEN SOVER=AOVER IMPN=IHP END IF END DO IHPN(IHN)=IMPN IHPN(IHN+1)=IMPN+1 IHPNTR(IHN)=IMPN+1 IHPNTR(IHN+1)=IMPN c END DO C C particle levels C DO IPN=1,2*NPLD(1)-1,2 SOVER=0.D0 DO IPP=1,2*NPLD(2)-1,2 AOVER=0.D0 OVER=0.D0 DO IH=1,NGHP DO IL=1,NGL OVER=OVER+(SPF(IH,IL,IPN,1)*SPF(IH,IL,IPP,2)+ $ SMF(IH,IL,IPN,1)*SMF(IH,IL,IPP,2))* $ GH(IH)*GL(IL) END DO END DO AOVER=DABS(OVER) IF (AOVER.GT.SOVER) THEN SOVER=AOVER IMPN=IPP END IF C END DO IPPN(IPN)=IMPN IPPN(IPN+1)=IMPN+1 IPPNTR(IPN)=IMPN+1 IPPNTR(IPN+1)=IMPN END DO C RETURN END C======================================================================= C Determine the excitation window for the treatment by HTDA method c dans le cas delta C======================================================================= SUBROUTINE NEWEXCWIN1(EDEL1, EDEL2, EDELD1, EDELD2) INCLUDE 'All.f' C LOGICAL NEWBL INTEGER TIME(3) CHARACTER * 80 FMT C CALL ITIME(TIME) WRITE(20,99) TIME PRINT 99,TIME C IF (INP.EQ.1) THEN OPEN(UNIT=30,FILE='HTDAnwindow') ELSEIF (INP.EQ.2) THEN OPEN(UNIT=30,FILE='HTDApwindow') ENDIF NBW = 0 DO I = 1, NOB NSIB(I) = 0 ENDDO NM = NUN(1) + NUN(2) WRITE(30,100) CSYM,NM,NAME,NUN,NM C======================================================================= C Determine the excitation windows C======================================================================= ITMP = NUN(INP)/2 C Determine Energie of the Fermi level FERMIE = (EHF(ITMP+1) + EHF(ITMP)) / 2 ELO = FERMIE - EDEL1 EHI = FERMIE + EDEL2 ELOD = FERMIE - EDELD1 EHID = FERMIE + EDELD2 C The first level in the windows (PETITE) ILOD = ITMP do while (ehf(ilod-1).ge.elod.and.ilod.gt.1) ilod=ilod-1 enddo C The last level in the windows (PETITE) IHID = ITMP do while (ehf(ihid+1).le.ehid.and.ihid.lt.nov) ihid=ihid+1 enddo c NHLD(inp) = ITMP - ILOD + 1 NPLD(inp) = IHID - ITMP c C The first level in the window (GRANDE) ILO = ILOD do while (ehf(ilo-1).ge.elo.and.ilo.gt.1) ilo=ilo-1 enddo C The last level in the window (GRANDE) IHI = IHID do while (ehf(ihi+1).le.ehi.and.ihi.lt.nov) ihi=ihi+1 enddo c NHL(inp) = ITMP - ILO + 1 NPL(inp) = IHI - ITMP c cbk c c IF (INP.EQ.2) THEN c c NPL(inp) = NPL(inp) -1 c NPLD(inp) = NPLD(inp)-1 c c NFERMI = NUN(INP)/2 c DO I=NFERMI+1,NOV-1 c EHF(I)=EHF(I+1) c NBRE(I)=NBRE(I+1) c END DO c END IF cbk C Integral Constant for (r,theta,phi)=>(eta,xi) C X = Pi/BZP2 * 2 /(2Pi) : 2 from symmetry on axe z C 2Pi= (sqrt(2pi)**2) : sqrt(2pi) comme from angle part of function C that we have "forgoten" X = 1D0/BZP2 C WRITE(20,106) NHL(inp),NPL(inp),NHLD(inp),NPLD(inp) c ICT1 = 0 WFU(0,INP) = 1D0 C======================================================================= C Determination of states in the window and their wave function C======================================================================= DO ICT = ILO, ITMP ICT1 = ICT1 + 1 ISW(-ICT1,INP) = ICT C IBLOCK : Block that contain this single state IBLOCK = MB(NBRE(ICT)) SSE(ICT1,1,INP) = EHF(ICT) SS2K(ICT1,1,INP) = K2B(IBLOCK) SSIP(ICT1,1,INP) = IPB(IBLOCK) SSSZ(ICT1,1,INP) = 0D0 SSR2(ICT1,1,INP) = 0D0 TWFU = (1 + DEXP(-WFX/WFMU))/ & (1 + DEXP((DABS(SSE(ICT1,1,INP)-FERMIE)-WFX)/WFMU)) TWFU = DSQRT(TWFU) WFU(-ICT1,INP) = TWFU C NVTB: Number of Vectors in This Block C INB : Position of this vector in this block NHBM1 = NHB(IBLOCK) - 1 NVTB = NVB(IBLOCK) INB = NBRE(ICT) - NHBM1 NEWBL = .TRUE. DO I = 1, NBW IF (IOB(I).EQ.IBLOCK) THEN NEWBL = .FALSE. NSIB(I) = NSIB(I) + 1 ISB(I,NSIB(I)) =-ICT1 IRBB(I,NSIB(I)) = INB IS2W(ICT) = NSIB(I) ENDIF ENDDO IF (NEWBL) THEN NBW = NBW + 1 NSIB(NBW) = 1 ISB(NBW,1) =-ICT1 IRBB(NBW,1) = INB IOB(NBW) = IBLOCK IS2W(ICT) = 1 ENDIF DO IH = 1, NGHP XI = XH(IH) COEF1 = DEXP(-XI*XI/2D0) DO IL = 1, NGL ETA = XL(IL) COEF2 = DEXP(-ETA/2D0) SPF(IH,IL,-ICT1,INP) = 0D0 SMF(IH,IL,-ICT1,INP) = 0D0 DO I = 1, NVTB NI = NHBM1 + I NSI = MS(NI) NLI = ML(NI) NRI = MR(NI) NZI = MZ(NI) COEF3 = 1D0 IF (NLI.NE.0) COEF3 = DSQRT(ETA**NLI) QHI = QH(IH,NZI) QLI = QL(IL,NRI,NLI) C DECF: Decomposition facteur DECF= COEF(I,INB,IBLOCK,INP) C PFUN: Partial Function PFUN= QHI*QLI*DECF*COEF3*COEF2*COEF1*BP*DSQRT(BZ)*DSQRT(2D0) IF (NSI.EQ.1) THEN IF (IH.EQ.1.AND.IL.EQ.1) & SSSZ(ICT1,1,INP) = SSSZ(ICT1,1,INP) + DECF**2 SPF(IH,IL,-ICT1,INP) = SPF(IH,IL,-ICT1,INP) + PFUN ELSE IF (IH.EQ.1.AND.IL.EQ.1) & SSSZ(ICT1,1,INP) = SSSZ(ICT1,1,INP) - DECF**2 SMF(IH,IL,-ICT1,INP) = SMF(IH,IL,-ICT1,INP) + PFUN ENDIF ENDDO SSR2(ICT1,1,INP) = SSR2(ICT1,1,INP) + ((XI/BZ)**2+ETA/(BP**2))* &(SPF(IH,IL,-ICT1,INP)**2+SMF(IH,IL,-ICT1,INP)**2) &*GH(IH)*GL(IL) ENDDO ENDDO SSSZ(ICT1,1,INP) = SSSZ(ICT1,1,INP)/2D0 SSR2(ICT1,1,INP) = SSR2(ICT1,1,INP) * X C ICT1N: ICT1 New ICT1N = ICT1 + 1 ISW(-ICT1N,INP) =-ICT SSE(ICT1N,1,INP) = EHF(ICT) SS2K(ICT1N,1,INP)=-K2B(IBLOCK) SSIP(ICT1N,1,INP)= IPB(IBLOCK) SSSZ(ICT1N,1,INP)=-SSSZ(ICT1,1,INP) SSR2(ICT1N,1,INP)= SSR2(ICT1,1,INP) WFU(-ICT1N,INP) = TWFU cc-->> xphase=1.d0 c if (k2b(iblock).eq.1) xphase=-1.d0 cc<<-- DO IH = 1, NGHP DO IL = 1, NGL SPF(IH,IL,-ICT1N,INP) = SMF(IH,IL,-ICT1,INP) SMF(IH,IL,-ICT1N,INP) =-SPF(IH,IL,-ICT1,INP) cc-->> c spf(ih,il,ict1n) = xphase*smf(ih,il,ict1) c smf(ih,il,ict1n) =-xphase*spf(ih,il,ict1) cc<<-- ENDDO ENDDO ICT1 = ICT1N ENDDO C ICT1 = 0 ITMP = ITMP + 1 DO ICT = ITMP, IHI ICT1 = ICT1 + 1 ISW(ICT1,INP) = ICT IBLOCK = MB(NBRE(ICT)) SSE(ICT1,2,INP) = EHF(ICT) SS2K(ICT1,2,INP) = K2B(IBLOCK) SSIP(ICT1,2,INP) = IPB(IBLOCK) SSSZ(ICT1,2,INP) = 0D0 SSR2(ICT1,2,INP) = 0D0 TWFU = (1 + DEXP(-WFX/WFMU))/ & (1 + DEXP((DABS(SSE(ICT1,2,INP)-FERMIE)-WFX)/WFMU)) TWFU = DSQRT(TWFU) WFU(ICT1,INP) = TWFU NHBM1 = NHB(IBLOCK) - 1 NVTB = NVB(IBLOCK) INB = NBRE(ICT) - NHBM1 NEWBL = .TRUE. DO I = 1, NBW IF (IOB(I).EQ.IBLOCK) THEN NEWBL = .FALSE. NSIB(I) = NSIB(I) + 1 ISB(I,NSIB(I)) = ICT1 IRBB(I,NSIB(I)) = INB IS2W(ICT) = NSIB(I) ENDIF ENDDO IF (NEWBL) THEN NBW = NBW + 1 NSIB(NBW) = 1 ISB(NBW,1) = ICT1 IRBB(NBW,1) = INB IOB(NBW) = IBLOCK IS2W(ICT) = 1 ENDIF DO IH = 1, NGHP XI = XH(IH) COEF1 = DEXP(-XI*XI/2D0) DO IL = 1, NGL ETA = XL(IL) COEF2 = DEXP(-ETA/2D0) SPF(IH,IL,ICT1,INP) = 0D0 SMF(IH,IL,ICT1,INP) = 0D0 DO I = 1, NVTB NI = NHBM1 + I NSI = MS(NI) NLI = ML(NI) NRI = MR(NI) NZI = MZ(NI) COEF3 = 1D0 IF (NLI.NE.0) COEF3= DSQRT(ETA**NLI) QHI = QH(IH,NZI) QLI = QL(IL,NRI,NLI) DECF= COEF(I,INB,IBLOCK,INP) PFUN= QHI*QLI*DECF*COEF3*COEF2*COEF1*BP*DSQRT(BZ)*DSQRT(2D0) IF (NSI.EQ.1) THEN IF (IH.EQ.1.AND.IL.EQ.1) & SSSZ(ICT1,2,INP) = SSSZ(ICT1,2,INP) + DECF**2 SPF(IH,IL,ICT1,INP) = SPF(IH,IL,ICT1,INP) + PFUN ELSE IF (IH.EQ.1.AND.IL.EQ.1) & SSSZ(ICT1,2,INP) = SSSZ(ICT1,2,INP) - DECF**2 SMF(IH,IL,ICT1,INP) = SMF(IH,IL,ICT1,INP) + PFUN ENDIF ENDDO SSR2(ICT1,2,INP) = SSR2(ICT1,2,INP)+((XI/BZ)**2+ETA/(BP**2))* &(SPF(IH,IL,ICT1,INP)**2+SMF(IH,IL,ICT1,INP)**2)*GH(IH)*GL(IL) ENDDO ENDDO SSSZ(ICT1,2,INP) = SSSZ(ICT1,2,INP)/2D0 SSR2(ICT1,2,INP) = SSR2(ICT1,2,INP) * X C ICT1N = ICT1 + 1 ISW(ICT1N,INP) =-ICT SSE(ICT1N,2,INP) = EHF(ICT) SS2K(ICT1N,2,INP)=-K2B(IBLOCK) SSIP(ICT1N,2,INP)= IPB(IBLOCK) SSSZ(ICT1N,2,INP)=-SSSZ(ICT1,2,INP) SSR2(ICT1N,2,INP)= SSR2(ICT1,2,INP) WFU(ICT1N,INP) = TWFU cc-->> xphase=1.d0 c if (k2b(iblock).eq.1) xphase=-1.d0 cc<<-- DO IH = 1, NGHP DO IL = 1, NGL SPF(IH,IL,ICT1N,INP) = SMF(IH,IL,ICT1,INP) SMF(IH,IL,ICT1N,INP) =-SPF(IH,IL,ICT1,INP) cc-->> c spf(ih,il,ict1n) = xphase*smf(ih,il,ict1) c smf(ih,il,ict1n) =-xphase*spf(ih,il,ict1) cc<<-- ENDDO ENDDO ICT1 = ICT1N ENDDO C WRITE(30,101) NCL(INP), NHL(inp), NPL(inp) C WRITE(30,102) EHI if (NPL(inp).EQ.NPLD(inp)) WRITE(30,202) EHID C Particles states DO ICT = 1, NPL(inp) ICT1 = 2*(NPL(inp)-ICT+1) - 1 WRITE(30,109) ICT1,ICT1+2*NHL(inp), & SS2K(ICT1,2,INP),PTY(SSIP(ICT1,2,INP)), & SSE(ICT1,2,INP), & SSE(ICT1+1,2,INP),SS2K(ICT1+1,2,INP),PTY(SSIP(ICT1+1,2,INP)), & ICT1+1, & ICT1+1+2*NHL(inp),WFU(ICT1,INP),SSSZ(ICT1,2,INP), & SSR2(ICT1,2,INP),IHI-ICT+1 IF (ICT.EQ.NPL(inp)-NPLD(inp)) WRITE(30,202) EHID ENDDO C WRITE(30,104) FERMIE C Hole states DO ICT = 1, NHL(inp) ICT1 = 2*(NHL(inp)-ICT+1) - 1 WRITE(30,103) ICT1,SS2K(ICT1,1,INP), & PTY(SSIP(ICT1,1,INP)),SSE(ICT1,1,INP), & SSE(ICT1+1,1,INP),SS2K(ICT1+1,1,INP), &PTY(SSIP(ICT1+1,1,INP)),ICT1+1, & WFU(-ICT1,INP),SSSZ(ICT1,1,INP),SSR2(ICT1,1,INP), &ILO+NHL(inp)-ICT IF (ICT.EQ.NHLD(inp)) WRITE(30,202) ELOD ENDDO C WRITE(30,105) ELO C WRITE(30,107) NBW ICH0 = ICHAR('0') DO I = 1,NBW FMT = '(I5,I11,8X,1H|,I3)' IF (NSIB(I).GT.1) THEN FMT = '(I5,I11,8X,1H|,I3, (1H,,I3))' FMT(19:19) = CHAR(ICH0+NSIB(I)/10) FMT(20:20) = CHAR(ICH0+NSIB(I)-NSIB(I)/10*10-1) ENDIF C WRITE(30,FMT) IOB(I),NSIB(I),(ISB(I,J),J=1,NSIB(I)) WRITE(30,*) IOB(I),NSIB(I),(ISB(I,J),J=1,NSIB(I)) ENDDO CLOSE(30) C Inforamation used to verify if there's any level oscillation around C Fermi level ICT1 = 1 WRITE(80+INP,*) 'Iteration : ',ITHF+ITHTDA WRITE(80+INP,103) ICT1,SS2K(ICT1,2,INP),PTY(SSIP(ICT1,2,INP)), & SSE(ICT1,2,INP), & SSE(ICT1+1,2,INP),SS2K(ICT1+1,2,INP), &PTY(SSIP(ICT1+1,2,INP)),ICT1+1, & WFU(ICT1,INP),SSSZ(ICT1,2,INP),SSR2(ICT1,1,INP),IHI-ICT+1 WRITE(80+INP,104) FERMIE ICT1 = 2*NHL(inp) - 1 WRITE(80+INP,103) ICT1,SS2K(ICT1,1,INP),PTY(SSIP(ICT1,1,INP)), & SSE(ICT1,1,INP), & SSE(ICT1+1,1,INP),SS2K(ICT1+1,1,INP), &PTY(SSIP(ICT1+1,1,INP)),ICT1+1, & WFU(-ICT1,INP),SSSZ(ICT1,1,INP), &SSR2(ICT1,1,INP),ILO+NHL(inp)-ICT WRITE(80+INP,*) C "HTDAnpwindow" is used by another program to find the C K/pi many body configuration. This file is created only C when we do HTDA method with 0 HF iteration, 1 HTDA iteration C Write out (in binary) the information of excitation window C in the file HTDAnpwindow which will be used by another code C (findNPconf) to find autmatically configurations of low energy p-h. IF (MNHF.EQ.0.AND.MNHTDA.EQ.1) THEN IF (INP.EQ.1) THEN OPEN(UNIT=23,FORM='UNFORMATTED',FILE='HTDAnpwindow') REWIND 23 ENDIF WRITE(23) NHL(inp),NPL(inp) WRITE(23) &(SS2K(I,1,INP),SSE(I,1,INP),SSIP(I,1,INP),I = 1, 2*NHL(inp)-1, 2) WRITE(23) &(SS2K(I,2,INP),SSE(I,2,INP),SSIP(I,2,INP),I = 1, 2*NPL(inp)-1, 2) IF (INP.EQ.2) CLOSE(23) C Find blocks from the deepest sp state to the excitation window C NSIB --> NSIBP = NSIB' C NBW --> NBWP = NBW' NBWP = NBW DO I = 1, NBW NSIBP(I) = NSIB(I) ENDDO DO I = NBW + 1, NOB NSIBP(I) = 0 ENDDO C DO I = 1, ILO-1 NEWBL = .TRUE. IBLOCK = MB(NBRE(I)) NHBM1 = NHB(IBLOCK) - 1 INB = NBRE(I) - NHBM1 DO I1 = 1, NBWP IF (IBLOCK.EQ.IOB(I1)) THEN NEWBL = .FALSE. NSIBP(I1) = NSIBP(I1) + 1 !!! Sign attention ISB(I1,NSIBP(I1)) =-I IRBB(I1,NSIBP(I1)) = INB IS2W(I) = NSIBP(I1) ENDIF ENDDO IF (NEWBL) THEN NBWP = NBWP + 1 NSIBP(NBWP) = 1 !!! Sign attention ISB(NBWP,1) =-I IRBB(NBWP,1)= INB IOB(NBWP) = IBLOCK IS2W(I) = 1 ENDIF ENDDO ENDIF C 99 FORMAT(4X,36H*Determine the excitation window at ,25X, &I2.2,1H:,I2.2,1X,I2.2,1Hs) 100 FORMAT(A,1H-,I3,2H :,5X,A,8X,4HN = ,I3,5X,4HZ = ,I3,10X, & 4HA = ,I3) 101 FORMAT(/,22HExcitation window for ,A,//, &I3,12H hole levels,16X,I3,16H particle levels,/) 102 FORMAT(79(1H=),/, &54H n K,pi Energy Energy K,pi n Weight, &8H ,16H Level,/,79(1H=),/, &15(1H-),F8.3,3X,6H(MeV) ,28(1H-),19H Hight limit Energy) 202 FORMAT(15(1H-),F8.3,3X,6H(MeV) ,21(1H-), $ 26H Small window limit Energy) 103 FORMAT(4x,I3,1X,I3,2H/2,A1,2X,F7.3,1X, & F7.3,1X,I3,2H/2,A1,5X,I3,1X,F8.4,F8.3,F10.3,I6) 104 FORMAT(15(1H-),F8.3,3X,6H(MeV) ,35(1H-),12H Fermi Level) 105 FORMAT(14(1H-),F8.3,4X,6H(MeV) ,29(1H-),17H Low limit Energy, &/,79(1H=)) 106 FORMAT(6X,I2,12H hole levels,16X,I2,16H particle levels, $ ' in large window',/,6X,I2,12H hole levels,16X,I2, & 16H particle levels,' in small window') 107 FORMAT(/,43HNumber of Blocks in the excitation window :,I3,/, &79(1H=),/,32HBlock Number of States | States,/,79(1H=)) 108 FORMAT(I5,I11,8X,1H|,12(I3,1H,)) 109 FORMAT(I3,1H(,I3,1H),I3,2H/2,A1,2X,F7.3,1X, & F7.3,1X,I3,2H/2,A1,1X,I3,1H(,I3,1H),F8.4,F8.3,F10.3,I6) C RETURN END CC======================================================================= C Determine the matrix for HTDA calculation+tr correction c in the presence of neiutron-proton coupling C======================================================================= SUBROUTINE HTDAMAPN INCLUDE 'All.f' C c replace |Psi> with (|Psi>+\bar|Psi>)/sqrt(2) c c C VABL : in the line C ICOL : Column of the matrix element DIMENSION VABL(NE),ICOL(NE),DIAGEI(NE) INTEGER TIME(3) C c first we construct the table of SD pair due to tr symmetry c c NT1=1 NPLUS=0 DO I=1,NTOTPN IBARPN=MBREVERSEDPN(I) IF (I.LE.IABS(IBARPN)) THEN PHIPLUS(NT1,1)=I IF(IBARPN.EQ.I) THEN PHIPLUS(NT1,2)=0 !STATE + NPLUS=NPLUS+1 ELSE PHIPLUS(NT1,2)=IBARPN !STATE - END IF NT1=NT1+1 END IF END DO NTOT1PN=NT1-1 write(*,*) 'Number of tr even mb states',NT1-1 C CALL ITIME(TIME) WRITE(20,100) TIME PRINT 100,TIME C C NDIFF1N = 0 NDIF0HN = 0 NDIF0PN = 0 NDIFF1P = 0 NDIF0HP = 0 NDIF0PP = 0 cc<<-- C C EIDEN: Identity Matrice * EIDEN c DO I = 1, NTOT1PN C C |i+>=|k>+\bar|k> C c K =PHIPLUS(I,1) c KBAR=PHIPLUS(I,2) c K1 =IABS(KBAR) C c KMAX = 0 c DO J = I, NTOT1PN C c L =PHIPLUS(J,1) c LBAR=PHIPLUS(J,2) c L1 =IABS(LBAR) C C matrix elements 1/2 C c E1=0D0 c E2=0D0 c E3=0D0 c E4=0D0 c EMTMP=0D0 C c E1=RINEWPN(K,L) ! c IF(L1.NE.0) E2=RINEWPN(K,L1) ! c IF(K1.NE.0) E3=RINEWPN(K1,L) !<\bar psi|V|psi> c IF(K1.NE.0.AND.L1.NE.0) E4=RINEWPN(K1,L1) !<\bar psi|V|\bar psi> c c IF (L1.NE.0.AND.K1.EQ.0) THEN c EMTMP=(E1+(LBAR/L1)*E2)/DSQRT(2.D0) c ELSEIF (L1.EQ.0.AND.K1.NE.0) THEN c EMTMP=(E1+(KBAR/K1)*E3)/DSQRT(2.D0) c ELSEIF (L1.NE.0.AND.K1.NE.0) THEN c EMTMP=(E1+E4+(LBAR/L1)*E2+(KBAR/K1)*E3)/2d0 c ELSEIF (L1.EQ.0.AND.K1.EQ.0) THEN c EMTMP=E1 c END IF C c c if (i.eq.1.or.j.eq.1) emtmp=0d0 c IF (I.EQ.J) THEN c DIAGE = EMTMP +EPHPN(K)+ EIDEN c EMTMP = DIAGE c ENDIF c IF (DABS(DIAGE).GT.1D-30) THEN c IF (DABS(EMTMP/DIAGE).GT.EPHTDA) THEN c KMAX = KMAX + 1 c VABL(KMAX) = EMTMP c ICOL(KMAX) = J c ENDIF c ELSE c IF (DABS(EMTMP).GT.EPHTDA) THEN c KMAX = KMAX + 1 c VABL(KMAX) = EMTMP c ICOL(KMAX) = J c ENDIF c ENDIF c ENDDO c NF=29 c WRITE(NF) KMAX, (VABL(J),ICOL(J), J = 1, KMAX) c ENDDO cc<<-- c*********************************************************** open (unit=33, & file='diage',form='unformatted',status='unknown') open (unit=35, & file='nondiag',form='unformatted',status='unknown') c this two files contain what newlan will read. It has to be c in this order. c************************************************************ C EIDEN: Identity Matrice * EIDEN DO I = 1, NTOT1PN C C |i+>=|k>+\bar|k> EE1=0d0 EE2=0d0 EE3=0d0 EE4=0d0 EMTMP=0d0 c K =PHIPLUS(I,1) KBAR=PHIPLUS(I,2) K1 =IABS(KBAR) C EE1=RINEWPN(K,K) EE2=RINEWPN(K,K1) EE3=RINEWPN(K1,K) EE4=RINEWPN(K1,K1) c IF (K1.NE.0) THEN EMTMP=(EE1+EE4+(KBAR/K1)*EE2+(KBAR/K1)*EE3)/2d0 ELSEIF (K1.EQ.0) THEN EMTMP=EE1 END IF DIAGEI(I)=EPHPN(K)+EMTMP write(65,*) i, diagei(i) c KMAX = 0 c DO J = I+1, NTOT1PN C L =PHIPLUS(J,1) LBAR=PHIPLUS(J,2) L1 =IABS(LBAR) C C matrix elements 1/2 C E1=0D0 E2=0D0 E3=0D0 E4=0D0 EMTMP=0D0 C E1=RINEWPN(K,L) ! IF(L1.NE.0) E2=RINEWPN(K,L1) ! IF(K1.NE.0) E3=RINEWPN(K1,L) !<\bar psi|V|psi> IF(K1.NE.0.AND.L1.NE.0) E4=RINEWPN(K1,L1) !<\bar psi|V|\bar psi> c IF (L1.NE.0.AND.K1.EQ.0) THEN EMTMP=(E1+(LBAR/L1)*E2)/DSQRT(2.D0) ELSEIF (L1.EQ.0.AND.K1.NE.0) THEN EMTMP=(E1+(KBAR/K1)*E3)/DSQRT(2.D0) ELSEIF (L1.NE.0.AND.K1.NE.0) THEN EMTMP=(E1+E4+(LBAR/L1)*E2+(KBAR/K1)*E3)/2d0 ELSEIF (L1.EQ.0.AND.K1.EQ.0) THEN EMTMP=E1 END IF c IF (DABS(EMTMP).GT.EPHTDA) THEN KMAX = KMAX + 1 VABL(KMAX) = EMTMP ICOL(KMAX) = J ENDIF C ENDDO c WRITE(35) I,I,KMAX WRITE(35) (VABL(J),J=1,KMAX) WRITE(35) (ICOL(J),J=1,KMAX) WRITE(35) KMAX,I,I ENDDO WRITE(33) (DIAGEI(I),I=1,NTOT1PN) CLOSE(33) CLOSE(35) c PRINT 101, EPHTDA WRITE(20,101) EPHTDA PRINT 102, NDIFF1N,NDIF0HN,NDIF0PN WRITE(20,102) NDIFF1N,NDIF0HN,NDIF0PN PRINT 102, NDIFF1P,NDIF0HP,NDIF0PP WRITE(20,102) NDIFF1P,NDIF0HP,NDIF0PP c REWIND(25) REWIND(26) REWIND(27) WRITE(24) NDIFF1N,NDIF0HN,NDIF0PN DO I=1,NDIFF1N READ(25) N1,N2,N3,N4,N5 WRITE(24) N1,N2,N3,N4,N5 ENDDO DO I=1,NDIF0HN READ(26) N1,N2 WRITE(24) N1,N2 ENDDO DO I=1,NDIF0PN READ(27) N1,N2 WRITE(24) N1,N2 ENDDO c REWIND(45) REWIND(46) REWIND(47) WRITE(44) NDIFF1P,NDIF0HP,NDIF0PP DO I=1,NDIFF1P READ(45) N1,N2,N3,N4,N5 WRITE(44) N1,N2,N3,N4,N5 ENDDO DO I=1,NDIF0HP READ(46) N1,N2 WRITE(44) N1,N2 ENDDO DO I=1,NDIF0PP READ(47) N1,N2 WRITE(44) N1,N2 ENDDO 100 FORMAT(4X,39H*Calculate many-body matrix elements at,22X, &I2.2,1H:,I2.2,1X,I2.2,1Hs) 101 FORMAT(6X,39HCut-off for value of matrix elements : ,D7.1) 102 FORMAT(6X,19HDensity supplement :,/,10X,7HDIFF1 =,I10,10X, &10HDIFF0 H =,I10,5X,3HP =,I10) c RETURN END C C======================================================================= C Calculate many-body matrix elements of residual interaction C======================================================================= FUNCTION RINEWPN(IA,IB) INCLUDE 'All.f' DIMENSION IEXOPN(0:4),IEXOPP(0:4) DOUBLE PRECISION EMADPN,EMAD 11 format(13x,'+',/,'|',i3,'> = ',i2,' * a a |',i3,'>',/, # 13x,i3,2x,i3) 12 format(13x,'+ +',/,'|',i3,'> = ',i2,' * a a a a', # 5x,'|',i3,'>',/,13x,i3,2x,i3,2x,i3,2x,i3) c ETMP=0.D0 ETMP1=0.D0 ETMP2=0.D0 ETMP3=0.D0 c CALL MBCOMPPN(IA,IB,ITYPE,NPH,IEXOPN,IEXOPP) IF (NPH.LT.0.OR.NPH.GT.2) THEN RINEWPN=ETMP RETURN ENDIF c IF(NPH.EQ.0) THEN C C |A>=|B> DIAGONAL MATRIX ELEMENTS C C***************************************************************** C !SUMS FOR NEUTRONS c***************************************************************** IPH=1 !SUM OVER PARTICLE STATES DO WHILE ( (MBSTN(IPH,2,IA).NE.0).AND.(IPH.LE.NHPMAX) ) !SUM OVER K JPH=IPH+1 NP1N=MBSTN(IPH,2,IA) WRITE(27) ISW(NP1N,1),IA NDIF0PN=NDIF0PN+1 DO WHILE ( (MBSTN(JPH,2,IA).NE.0).AND.(JPH.LE.NHPMAX) )! SUM OVER L NP2N=MBSTN(JPH,2,IA) ETMP1=ETMP1 + EMAD(NP1N,NP2N,NP1N,NP2N,1) JPH=JPH+1 END DO IPH=IPH+1 END DO C IPH=1 !SUM OVER HOLE STATES DO WHILE ( (MBSTN(IPH,1,IA).NE.0).AND.(IPH.LE.NHPMAX) ) !SUM OVER K JPH=IPH+1 NT1N=-MBSTN(IPH,1,IA) WRITE(26) ISW(NT1N,1),IA NDIF0HN=NDIF0HN+1 DO WHILE ( (MBSTN(JPH,1,IA).NE.0).AND.(JPH.LE.NHPMAX) )! SUM OVER L NT2N=-MBSTN(JPH,1,IA) ETMP1=ETMP1 + EMAD(NT1N,NT2N,NT1N,NT2N,1) JPH=JPH+1 END DO IPH=IPH+1 END DO C IPH=1 !SUM OVER HOLE AND PARTICLE STATES DO WHILE ( (MBSTN(IPH,2,IA).NE.0).AND.(IPH.LE.nhpmax) ) !SUM OVER K JPH=1 NPN=MBSTN(IPH,2,IA) DO WHILE ( (MBSTN(JPH,1,IA).NE.0).AND.(JPH.LE.nhpmax) )! SUM OVER L NTN=-MBSTN(JPH,1,IA) ETMP1=ETMP1 - EMAD(NTN,NPN,NTN,NPN,1) JPH=JPH+1 END DO IPH=IPH+1 END DO c*********************************************************************** C !SUMS FOR PROTONS c*********************************************************************** IPH=1 !SUM OVER PARTICLE STATES DO WHILE ( (MBSTP(IPH,2,IA).NE.0).AND.(IPH.LE.nhpmax) ) !SUM OVER K JPH=IPH+1 NP1P=MBSTP(IPH,2,IA) WRITE(47) ISW(NP1P,2),IA NDIF0PP=NDIF0PP+1 DO WHILE ( (MBSTP(JPH,2,IA).NE.0).AND.(JPH.LE.nhpmax) )! SUM OVER L NP2P=MBSTP(JPH,2,IA) ETMP2=ETMP2 + EMAD(NP1P,NP2P,NP1P,NP2P,2) JPH=JPH+1 END DO IPH=IPH+1 END DO C IPH=1 !SUM OVER HOLE STATES DO WHILE ( (MBSTP(IPH,1,IA).NE.0).AND.(IPH.LE.nhpmax) ) !SUM OVER K JPH=IPH+1 NT1P=-MBSTP(IPH,1,IA) WRITE(46) ISW(NT1P,2),IA NDIF0HP=NDIF0HP+1 DO WHILE ( (MBSTP(JPH,1,IA).NE.0).AND.(JPH.LE.nhpmax) )! SUM OVER L NT2P=-MBSTP(JPH,1,IA) ETMP2=ETMP2 + EMAD(NT1P,NT2P,NT1P,NT2P,2) JPH=JPH+1 END DO IPH=IPH+1 END DO C IPH=1 !SUM OVER HOLE AND PARTICLE STATES DO WHILE ( (MBSTP(IPH,2,IA).NE.0).AND.(IPH.LE.nhpmax) ) !SUM OVER K JPH=1 NPP=MBSTP(IPH,2,IA) DO WHILE ( (MBSTP(JPH,1,IA).NE.0).AND.(JPH.LE.nhpmax) )! SUM OVER L NTP=-MBSTP(JPH,1,IA) ETMP2=ETMP2 - EMAD(NTP,NPP,NTP,NPP,2) JPH=JPH+1 END DO IPH=IPH+1 END DO C******************************************************************** C SUMS FOR PROTON-NEUTRON PART (SUMS WITH RESPECT TO QVACUUM) C******************************************************************** C IPH=1 !SUM OVER NEUTRON PARTICLE STATES DO WHILE ( (MBSTN(IPH,2,IA).NE.0).AND.(IPH.LE.nhpmax) ) !SUM OVER K JPH=1 NPN=MBSTN(IPH,2,IA) DO WHILE ( (MBSTP(JPH,2,IA).NE.0).AND.(JPH.LE.nhpmax) )! SUM OVER L PARTICLES NPP=MBSTP(JPH,2,IA) ETMP3=ETMP3 + EMADPN(NPN,NPP,NPN,NPP) ! JPH=JPH+1 END DO IPH=IPH+1 END DO c IPH=1 DO WHILE ( (MBSTN(IPH,2,IA).NE.0).AND.(IPH.LE.nhpmax) ) !SUM OVER K JPH1=1 NPN=MBSTN(IPH,2,IA) DO WHILE ((MBSTP(JPH1,1,IA).NE.0).AND.(JPH1.LE.nhpmax))!SUM OVER L HOLES NTP=-MBSTP(JPH1,1,IA) ETMP3=ETMP3 - EMADPN(NPN,NTP,NPN,NTP) ! JPH1=JPH1+1 END DO IPH=IPH+1 END DO C c c IPH=1 !SUM OVER NEUTRON HOLE STATES DO WHILE ( (MBSTN(IPH,1,IA).NE.0).AND.(IPH.LE.nhpmax) ) !SUM OVER K JPH=1 NTN=-MBSTN(IPH,1,IA) DO WHILE ( (MBSTP(JPH,2,IA).NE.0).AND.(JPH.LE.nhpmax) )! SUM OVER L PARTICLES NPP=MBSTP(JPH,2,IA) ETMP3=ETMP3 - EMADPN(NTN,NPP,NTN,NPP) ! JPH=JPH+1 END DO IPH=IPH+1 END DO c IPH=1 DO WHILE ( (MBSTN(IPH,1,IA).NE.0).AND.(IPH.LE.nhpmax) ) !SUM OVER K JPH1=1 NTN=-MBSTN(IPH,1,IA) DO WHILE ((MBSTP(JPH1,1,IA).NE.0).AND.(JPH1.LE.nhpmax))!SUM OVER L HOLES NTP=-MBSTP(JPH1,1,IA) ETMP3=ETMP3 + EMADPN(NTN,NTP,NTN,NTP) JPH1=JPH1+1 END DO IPH=IPH+1 END DO c RINEWPN=ETMP1+ETMP2+ETMP3 c RETURN END IF c*************************************************************************************** c non-diagonal matrix elements c*************************************************************************************** c |A> differ from |B> by one nucleon c*************************************************************************************** IF (NPH.EQ.1) THEN IF (ITYPE.EQ.1) THEN !!!one neutron difference IPN=IEXOPN(1) ITN=IEXOPN(3) WRITE(25) IEXOPN(0),ISW(IPN,1),ISW(ITN,1),IA,IB NDIFF1N=NDIFF1N+1 c c DIAGONAL IN PROTONS c IPH=1 !SUM OVER PARTICLE STATES DO WHILE ( (MBSTP(IPH,2,IA).NE.0).AND.(IPH.LE.nhpmax) ) !SUM OVER K JPH=IPH+1 NP1P=MBSTP(IPH,2,IA) DO WHILE ( (MBSTP(JPH,2,IA).NE.0).AND.(JPH.LE.nhpmax) )! SUM OVER L NP2P=MBSTP(JPH,2,IA) ETMP2=ETMP2 + EMAD(NP1P,NP2P,NP1P,NP2P,2) JPH=JPH+1 END DO IPH=IPH+1 END DO C IPH=1 !SUM OVER HOLE STATES DO WHILE ( (MBSTP(IPH,1,IA).NE.0).AND.(IPH.LE.nhpmax) ) !SUM OVER K JPH=IPH+1 NT1P=-MBSTP(IPH,1,IA) DO WHILE ( (MBSTP(JPH,1,IA).NE.0).AND.(JPH.LE.nhpmax) )! SUM OVER L NT2P=-MBSTP(JPH,1,IA) ETMP2=ETMP2 + EMAD(NT1P,NT2P,NT1P,NT2P,2) JPH=JPH+1 END DO IPH=IPH+1 END DO C IPH=1 !SUM OVER HOLE AND PARTICLE STATES DO WHILE ( (MBSTP(IPH,2,IA).NE.0).AND.(IPH.LE.nhpmax) ) !SUM OVER K JPH=1 NPP=MBSTP(IPH,2,IA) DO WHILE ( (MBSTP(JPH,1,IA).NE.0).AND.(JPH.LE.nhpmax) )! SUM OVER L NTP=-MBSTP(JPH,1,IA) ETMP2=ETMP2 - EMAD(NTP,NPP,NTP,NPP,2) JPH=JPH+1 END DO IPH=IPH+1 END DO c c 1P1H FOR NEUTRONS c IPH=1 DO WHILE ( (MBSTN(IPH,2,IA).NE.0).AND.(IPH.LE.nhpmax)) IPN1=MBSTN(IPH,2,IA) ETMP1=ETMP1+EMAD(IPN1,IPN,IPN1,ITN,1) !!! IPH=IPH+1 END DO IPH=1 DO WHILE ( (MBSTN(IPH,1,IA).NE.0).AND.(IPH.LE.nhpmax)) ITN1=-MBSTN(IPH,1,IA) ETMP1=ETMP1-EMAD(ITN1,IPN,ITN1,ITN,1) !! IPH=IPH+1 END DO c c AND PROTON-NEUTRON INETRACTION c IPH=1 DO WHILE ( (MBSTP(IPH,2,IA).NE.0).AND.(IPH.LE.nhpmax))!SUM OVER PARTICLE STATES NPP=MBSTP(IPH,2,IA) ETMP3=ETMP3+EMADPN(ITN,NPP,IPN,NPP) IPH=IPH+1 END DO IPH=1 DO WHILE ( (MBSTP(IPH,1,IA).NE.0).AND.(IPH.LE.NHPMAX))!SUM OVER HOLE STATES NTP=-MBSTP(IPH,1,IA) ETMP3=ETMP3-EMADPN(ITN,NTP,IPN,NTP) IPH=IPH+1 END DO c RINEWPN=ETMP1*IEXOPN(0)+ETMP2*IEXOPN(0) & +ETMP3*IEXOPN(0)*IEXOPP(0) RETURN ELSE !!!ONE PROTON DIFFERENCE IPP=IEXOPP(1) ITP=IEXOPP(3) c WRITE(45) IEXOPP(0),ISW(IPP,2),ISW(ITP,2),IA,IB NDIFF1P=NDIFF1P+1 c c DIAGONAL IN NEUTRONS c IPH=1 !SUM OVER PARTICLE STATES DO WHILE ( (MBSTN(IPH,2,IA).NE.0).AND.(IPH.LE.NHPMAX) ) !SUM OVER K JPH=IPH+1 NP1N=MBSTN(IPH,2,IA) DO WHILE ( (MBSTN(JPH,2,IA).NE.0).AND.(JPH.LE.NHPMAX) )! SUM OVER L NP2N=MBSTN(JPH,2,IA) ETMP1=ETMP1 + EMAD(NP1N,NP2N,NP1N,NP2N,1) JPH=JPH+1 END DO IPH=IPH+1 END DO C IPH=1 !SUM OVER HOLE STATES DO WHILE ( (MBSTN(IPH,1,IA).NE.0).AND.(IPH.LE.NHPMAX) ) !SUM OVER K JPH=IPH+1 NT1N=-MBSTN(IPH,1,IA) DO WHILE ( (MBSTN(JPH,1,IA).NE.0).AND.(JPH.LE.NHPMAX) )! SUM OVER L NT2N=-MBSTN(JPH,1,IA) ETMP1=ETMP1 + EMAD(NT1N,NT2N,NT1N,NT2N,1) JPH=JPH+1 END DO IPH=IPH+1 END DO C IPH=1 !SUM OVER HOLE AND PARTICLE STATES DO WHILE ( (MBSTN(IPH,2,IA).NE.0).AND.(IPH.LE.nhpmax) ) !SUM OVER K JPH=1 NPN=MBSTN(IPH,2,IA) DO WHILE ( (MBSTN(JPH,1,IA).NE.0).AND.(JPH.LE.nhpmax) )! SUM OVER L NTN=-MBSTN(JPH,1,IA) ETMP1=ETMP1 - EMAD(NTN,NPN,NTN,NPN,1) JPH=JPH+1 END DO IPH=IPH+1 END DO C C 1p1H IN PROTONS C IPH=1 !SUM OVER PARTICLE STATES DO WHILE ( (MBSTP(IPH,2,IA).NE.0).AND.(IPH.LE.nhpmax)) NPP=MBSTP(IPH,2,IA) ETMP2=ETMP2+EMAD(NPP,IPP,NPP,ITP,2) IPH=IPH+1 END DO IPH=1 !SUM OVER HOLE STATES DO WHILE ( (MBSTP(IPH,1,IA).NE.0).AND.(IPH.LE.nhpmax)) NTP=-MBSTP(IPH,1,IA) ETMP2=ETMP2-EMAD(NTP,IPP,NTP,ITP,2) IPH=IPH+1 END DO C c SUMS FOR PN PART IPH=1 DO WHILE ( (MBSTN(IPH,2,IA).NE.0).AND.(IPH.LE.nhpmax))!SUM OVER PARTICLE STATES NPN=MBSTN(IPH,2,IA) ETMP3=ETMP3+EMADPN(NPN,ITP,NPN,IPP) IPH=IPH+1 END DO IPH=1 DO WHILE ( (MBSTN(IPH,1,IA).NE.0).AND.(IPH.LE.NHPMAX))!SUM OVER HOLE STATES NTN=-MBSTN(IPH,1,IA) ETMP3=ETMP3-EMADPN(NTN,ITP,NTN,IPP) IPH=IPH+1 END DO c RINEWPN=ETMP2*IEXOPP(0)+ETMP1*IEXOPN(0) &+ETMP3*IEXOPN(0)*IEXOPP(0) c RETURN END IF END IF c****************************************************************************************** c c 2p2h excitations c c****************************************************************************************** IF (NPH.EQ.2) THEN rinewpn=0 etmp1=0 etmp2=0 etmp3=0 IF (ITYPE.EQ.2) THEN !!! |A> DIFFER FROM |B> BY ONE PROTON AND ONE NEUTRON IPN=IEXOPN(1) ITN=IEXOPN(3) IPP=IEXOPP(1) ITP=IEXOPP(3) C c c SUMS FOR NEUTRONS IPH=1 DO WHILE ( (MBSTN(IPH,2,IA).NE.0).AND.(IPH.LE.NHPMAX)) IPN1=MBSTN(IPH,2,IA) ETMP1=ETMP1+EMAD(IPN1,ITN,IPN1,IPN,1) IPH=IPH+1 END DO IPH=1 DO WHILE ( (MBSTN(IPH,1,IA).NE.0).AND.(IPH.LE.NHPMAX)) IPH1=-MBSTN(IPH,1,IA) ETMP1=ETMP1-EMAD(IPH1,ITN,IPH1,IPN,1) !! IPH=IPH+1 END DO c C SUMS FOR PROTONS IPH=1 !SUM OVER PARTICLE STATES DO WHILE ( (MBSTP(IPH,2,IA).NE.0).AND.(IPH.LE.NHPMAX)) NPP=MBSTP(IPH,2,IA) ETMP2=ETMP2+EMAD(NPP,IPP,NPP,ITP,2) IPH=IPH+1 END DO IPH=1 !SUM OVER HOLE STATES DO WHILE ( (MBSTP(IPH,1,IA).NE.0).AND.(IPH.LE.NHPMAX)) NTP=-MBSTP(IPH,1,IA) ETMP2=ETMP2-EMAD(NTP,IPP,NTP,ITP,2) IPH=IPH+1 END DO C PROTON-NEUTRON PART ETMP3=-EMADPN(IPN,IPP,ITN,ITP) !!! c RINEWPN=ETMP1*IEXOPN(0)+ETMP2*IEXOPP(0)+ & ETMP3*IEXOPP(0)*IEXOPN(0) c RETURN END IF IF (ITYPE.EQ.3) THEN !!!|A> DIFFER FROM |B> BY TWO NEUTRONS IPN= IEXOPN(1) IPN1=IEXOPN(2) ITN= IEXOPN(3) ITN1=IEXOPN(4) ETMP1=-EMAD(IPN,IPN1,ITN,ITN1,1) RINEWPN=ETMP1*IEXOPN(0) RETURN END IF IF (ITYPE.EQ.4) THEN !!!|A> DIFFER FROM |B> BY TWO PROTONS IPP= IEXOPP(1) IPP1=IEXOPP(2) ITP= IEXOPP(3) ITP1=IEXOPP(4) ETMP2=-EMAD(IPP,IPP1,ITP,ITP1,2) RINEWPN=ETMP2*IEXOPP(0) RETURN END IF END IF c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc RETURN c END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE MBCOMPPN(IA,IB,ITYPE,NPH,IEXOPN,IEXOPP) INCLUDE 'All.f' ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c |A> = a^+_pa a_ha |0> et |B> = a^+_pb a_hb |0> c ==> |B> = a^+_pb a_hb a^+_ha a a_pa |A> c = (-)^m0 a_hb a^+_ha a a^+_pb a_pa |A> c avec m0=npib(nhib+nhia) c nh et np contractions -> phase m1 c ==> |B> = (-)^m0+m1+m2 a^+_ha a^+_pb a_pa a_hb |A> c avec m2=(nhib-nh)(nhia-nh+npia+npib-2np) c meme nbe de part => m0+m2=nhib+(npia-1)npib+nh c nhib=npib => m0+m2=npib*npia+nh c DIMENSION IEXOPN(0:4),IEXOPP(0:4) dimension nab(2),nhpab(2,2,2) dimension nwick(-nhpmax:nhpmax,2,2),ncontr(2,2) cc-->> dimension mwick(-nhpmax:nhpmax,2,2) cc<<-- c IEXOPN(0)=0 IEXOPN(1)=0 IEXOPN(2)=0 IEXOPN(3)=0 IEXOPN(4)=0 IEXOPP(0)=0 IEXOPP(1)=0 IEXOPP(2)=0 IEXOPP(3)=0 IEXOPP(4)=0 c NAB(1)=IA NAB(2)=IB NPH=0 ITYPE=0 c do iab=1,2 j=3-2*iab do ihp=1,2 in=NHPMAX do while (MBSTN(in,ihp,nab(iab)).eq.0.and.in.ge.1) nwick(in*j,ihp,1)=0 in=in-1 enddo ip=nhpmax do while (MBSTP(ip,ihp,nab(iab)).eq.0.and.ip.ge.1) nwick(ip*j,ihp,2)=0 ip=ip-1 enddo nhpab(ihp,iab,1)=in nhpab(ihp,iab,2)=ip do in=1,nhpab(ihp,iab,1) nwick(in*j,ihp,1)=MBSTN(nhpab(ihp,iab,1)-in+1,ihp,nab(iab)) enddo c do ip=1,nhpab(ihp,iab,2) nwick(ip*j,ihp,2)=MBSTP(nhpab(ihp,iab,2)-ip+1,ihp,nab(iab)) enddo c enddo ! ihp c if (NAB(1).eq.NAB(2)) RETURN c ! états many-body identiques enddo ! iab c c c NPH=-1 if (nhpab(1,1,1)-nhpab(2,1,1).ne.nhpab(1,2,1)-nhpab(2,2,1)) RETURN ! differents nombres de particles and holes in neutron SD if (nhpab(1,1,2)-nhpab(2,1,2).ne.nhpab(1,2,2)-nhpab(2,2,2)) RETURN ! differents nombres de particles and holes in proton SD if (iabs(nhpab(1,1,1)-nhpab(2,1,1)).gt.2) RETURN if (iabs(nhpab(1,1,2)-nhpab(2,1,2)).gt.2) RETURN ! avec |n-m| > 2 c IEXOP(0)=nhpab(2,2)*(nhpab(1,1)+nhpab(1,2)) c MPHN=0 do ihp=1,2 ncontr(ihp,1)=0 i=-nhpab(ihp,2,1) do while (i.le.-1) j=nhpab(ihp,1,1)-ncontr(ihp,1) do while (j.ge.1) if (nwick(j,ihp,1).eq.nwick(i,ihp,1)) then ncontr(ihp,1)=ncontr(ihp,1)+1 IEXOPN(0)=IEXOPN(0)-i+j do k=i,-nhpab(ihp,2,1)+ncontr(ihp,1),-1 nwick(k,ihp,1)=nwick(k-1,ihp,1) enddo nwick(-nhpab(ihp,2,1)+ncontr(ihp,1)-1,ihp,1)=0 do k=j,nhpab(ihp,1,1)-ncontr(ihp,1) nwick(k,ihp,1)=nwick(k+1,ihp,1) enddo nwick(nhpab(ihp,1,1)-ncontr(ihp,1)+1,ihp,1)=0 j=1 elseif (nwick(j,ihp,1).lt.nwick(i,ihp,1)) then j=1 endif j=j-1 enddo i=i+1 enddo MPHN=MPHN+nhpab(ihp,ihp,1)-ncontr(ihp,1) if (MPHN.gt.2) RETURN ! difference greater than 2p-2h !!!neutrons enddo MPHP=0 do ihp=1,2 ncontr(ihp,2)=0 i=-nhpab(ihp,2,2) do while (i.le.-1) j=nhpab(ihp,1,2)-ncontr(ihp,2) do while (j.ge.1) if (nwick(j,ihp,2).eq.nwick(i,ihp,2)) then ncontr(ihp,2)=ncontr(ihp,2)+1 IEXOPP(0)=IEXOPP(0)-i+j do k=i,-nhpab(ihp,2,2)+ncontr(ihp,2),-1 nwick(k,ihp,2)=nwick(k-1,ihp,2) enddo nwick(-nhpab(ihp,2,2)+ncontr(ihp,2)-1,ihp,2)=0 do k=j,nhpab(ihp,1,2)-ncontr(ihp,2) nwick(k,ihp,2)=nwick(k+1,ihp,2) enddo nwick(nhpab(ihp,1,2)-ncontr(ihp,2)+1,ihp,2)=0 j=1 elseif (nwick(j,ihp,2).lt.nwick(i,ihp,2)) then j=1 endif j=j-1 enddo i=i+1 enddo MPHP=MPHP+nhpab(ihp,ihp,2)-ncontr(ihp,2) if (MPHP.gt.2) RETURN ! difference greater than 2p-2h !!! protons enddo c IEXOPN(0)=IEXOPN(0)+ncontr(1,1)+nhpab(2,2,1)*nhpab(2,1,1) IEXOPN(0)=1-2*mod(IEXOPN(0),2) do ihp=1,2 j=3-2*ihp m=2*(ihp-1) do i=1,nhpab(ihp,1,1)-ncontr(ihp,1) m=m+1 IEXOPN(m)=-j*nwick(i,ihp,1) enddo do i=-nhpab(3-ihp,2,1)+ncontr(3-ihp,1),-1 m=m+1 IEXOPN(m)=j*nwick(i,3-ihp,1) enddo enddo C IEXOPP(0)=IEXOPP(0)+ncontr(1,2)+nhpab(2,2,2)*nhpab(2,1,2) IEXOPP(0)=1-2*mod(IEXOPP(0),2) do ihp=1,2 j=3-2*ihp m=2*(ihp-1) do i=1,nhpab(ihp,1,2)-ncontr(ihp,2) m=m+1 IEXOPP(m)=-j*nwick(i,ihp,2) enddo do i=-nhpab(3-ihp,2,2)+ncontr(3-ihp,2),-1 m=m+1 IEXOPP(m)=j*nwick(i,3-ihp,2) enddo enddo NPH=MPHN+MPHP IF (NPH.GT.2) RETURN !Difference more than 2 nucleons IF (MPHN.EQ.1.AND.MPHP.EQ.0) ITYPE=1 IF (MPHN.EQ.1.AND.MPHP.EQ.1) ITYPE=2 IF (MPHN.EQ.2.AND.MPHP.EQ.0) ITYPE=3 IF (MPHN.EQ.0.AND.MPHP.EQ.2) ITYPE=4 c write(53,*) ia,ib c write(53,12) (IEXOPN(I)) c write(54,*) ia,ib c write(54,12) (iexopp(I)) 5 format('hole states : (',3i3,'|',i1,'|',3i3,')',/, # 'part states : (',3i3,'|',i1,'|',3i3,')') 6 format(2('(',3i3,'=>',3i3,')',5x)) 11 format(13x,'+',/,8x,i2,' * a a',/, # 13x,i3,2x,i3) 12 format(13x,'+ +',/,8x,i2,' * a a a a', #/,13x,i3,2x,i3,2x,i3,2x,i3) c RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE PNBASE INCLUDE 'All.f' c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc IF (INP.EQ.1) CALL CONFIG2(NUK2(INP),NUIP(INP),0,0,ECHTDA) c CALL CONFIG1(NUK2(INP),NUIP(INP),1,1,ECHTDA) CALL CONFIG1(NUK2(INP),NUIP(INP),2,2,ECHTDA) cccccccccccccccccccc test protons x neutrons ccccccccccccccccccccc IF (inp.eq.2) CALL CONFIG1(NUK2(INP),NUIP(INP),8,8,ECHTDA) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc IF (INP.EQ.2) CALL PARTNER() c IF (INP.EQ.2) CALL CONFIG2(NUK2(INP),NUIP(INP),2,2,ECHTDA) c CALL CONFIG1(NUK2(INP),NUIP(INP),4,4,ECHTDA) c c IF (INP.EQ.2) CALL CONFIG2(NUK2(INP),NUIP(INP),4,4,ECHTDA) c RETURN END C======================================================================= C To check C Calculate the correlated ground state C======================================================================= SUBROUTINE PNSTATES INCLUDE 'All.f' INTEGER TIME(3), MAXPIS(200) DOUBLE PRECISION CONTRI(0:10),CONTRP(0:10) INTEGER NCONTR(0:10),NCONTP(0:10) DIMENSION ZV(NVALMAX) C C************ C Determine how to define the order: by the number of holes C or by the number of particles: NORDER = 1 !Hole C NDPH:The difference between the number of particles and the number C of holes NDPH = 0 C IORMA: Maximal order IORDMA = NHPMAX - NDPH C DO I=1,200 MAXPIS(I)=0 END DO C Diagonalise the matrix C 10 CALL ITIME(TIME) PRINT 100, TIME WRITE(20,100) TIME c CALL newDIAG c &(NVALMAX,NE,NTOT1PN,-neigen,INP,nf,EVEC1,EVALUE) c c********************************************************************* c call Lanczos: c ntot:our matrix dimension, c neigen:number of eigenvalues, c nloop:number of Lanczos c iterations (at least as large as number of eigenvalues required) c zprec:precision of the Lanczos calculation. c filv: file in which vectors are stored c zv-eigenvalues c I put zprec and nloop and ifil here for a moment, but it would be nice c to have them in the data file c REWIND(IFIL) CALL NEWLAN(ZV,ZPREC,NTOT1PN,NEIGEN,NLOOP,IFIL) write(*,*) (zv(i),i=1,neigen),ntot1PN,'finished' c************ reading eigenvectors from ifil ********************* REWIND(IFIL) DO I=1,NEIGEN READ(IFIL)(EVEC1(K,I),K=1,NTOT1PN) END DO c***************************************************************** c c rearragement due to time-reversal symmetry conservation c and probabilities of various components c DO L=1,NEIGEN c DO K=1,NTOTPN EVEC(K) =0D0 VAB(K,L)=0D0 END DO c DO I=1,NTOT1PN K =PHIPLUS(I,1) KBAR=PHIPLUS(I,2) K1 =IABS(KBAR) C EVEC(K) =0D0 VAB(K1,L)=0D0 EVEC(K1) =0D0 VAB(K1,L)=0D0 C IF(K1.EQ.0) THEN EVEC(K) =EVEC1(I,1) VAB(K,L)=EVEC1(I,L) ELSE EVEC(K) =EVEC1(I,1)/DSQRT(2D0) EVEC(K1) =DFLOAT(KBAR/K1)*EVEC1(I,1)/DSQRT(2D0) VAB(K,L) =EVEC1(I,L)/DSQRT(2D0) VAB(K1,L)=DFLOAT(KBAR/K1)*EVEC1(I,L)/DSQRT(2D0) END IF END DO C END DO sumeph=0d0 sumchi2=0d0 do i=1,ntotpn sumeph=sumeph+ephpn(i)*evec(i)**2d0 sumchi2=sumchi2+evec(I)**2d0 end do write(20,*)'=', SUMEPH ,' sxi2=',sumchi2 C C MAXPI : Index of the configuration of Maximal probability C PROBM : Maximal Probability C PROBL = 1D-3 CALL VMAXC(NTOTPN,EVEC,MAXPI) CALL VMAXCS(NTOTPN,EVEC,20,MAXPIS) c c EPHMPN = EPHPN(MAXPI) DO K=1,NEIGEN EVALUEPN(K)=zv(k)!evalue(k,1)+evalue(k,2)-eiden!ZV(K) ENDDO PRINT 101, EVALUEPN(1) WRITE(20,101) EVALUEPN(1) c xi00=0D0 xi02=0d0 xi20=0d0 xi11=0d0 xi22=0d0 xi31=0d0 xi13=0d0 xi40=0d0 xi04=0d0 c n2040=nss(1)+nfour(1)+1 n0204=n2040+nss(2)+nfour(2) c do i=1,ntotpn if (i.eq.1) xi00=xi00+evec(i)**2d0 c if (I.gt.1.and.i.le.(nss(1)+1)) & xi20=xi20+evec(i)**2D0 c if (i.gt.(nss(1)+1).and.i.le.n2040) & xi40=xi40+evec(i)**2d0 c if (I.gt.n2040.and.i.le.(n2040+nss(2))) & xi02=xi02+evec(i)**2d0 c if (i.gt.(n2040+nss(2)).and.i.le.(n0204)) & xi04=xi04+evec(i)**2d0 c if (i.gt.(n0204).and.i.le.(n0204+nsspn(2))) & xi11=xi11+evec(i)**2d0 c if (i.gt.(n0204+nsspn(2)).and.i.le. & (n0204+nsspn(1)+nsspn(2))) & xi31=xi31+evec(i)**2D0 c if (i.gt.(n0204+nsspn(1)+nsspn(2)).and.i.le. & (n0204+nsspn(1)+nsspn(2)+nsspn(3)) ) & xi13=xi13+evec(i)**2d0 c if (i.gt.(n0204+nsspn(2)+nsspn(1)+nsspn(3))) & xi22=xi22+evec(i)**2d0 end do c x2=xi11+xi20+xi02 x4=xi40+xi04+xi22+xi13+xi31 c write(*,*)'nall',nallpn(1),nss(1),nfour(1),nss(2),nfour(2), &nsspn(2),nsspn(1),nsspn(3),nsspn(4) c write(20,*) write(20,*) 'WAVE FUNCTION DECOMPOSITION' write(20,*) '0p0h=',xi00, ' 2p2h=',x2, ' 4p4h=',x4 write(20,*) write(20,*) '(20=',xi20,' 11=',xi11,' 02=',xi02,')' write(20,*) &'(40=',xi40,' 31=',xi31,' 22=',xi22,' 13=',xi13,' 04=',xi04,')' C* PROBM = EVEC(MAXPI)**2 PRINT 102, MAXPI,((MBSTN(I,J,MAXPI),I = 1, NHPMAX), J = 1, 2), &((MBSTP(I,J,MAXPI),I = 1, NHPMAX), J = 1, 2), &EPHPN(MAXPI),PROBM WRITE(20,102) MAXPI,((MBSTN(I,J,MAXPI),I = 1, NHPMAX), J = 1, 2), &((MBSTP(i,j,MAXPI),I = 1, NHPMAX), J = 1, 2), &EPHPN(MAXPI),PROBM WRITE(20,103) 50,PROBL, &(K,MAXPIS(K),((MBSTN(I,J,MAXPIS(K)),I = 1, NHPMAX), J = 1, 2), &((MBSTP(I,J,MAXPIS(K)),I = 1, NHPMAX), J = 1, 2), & EPHPN(MAXPIS(K)),EVEC(MAXPIS(K))**2,K=1,20) c C 100 FORMAT(4X,43H*Diagonalise the matrix for Ground State at,18X, &I2.2,1H:,I2.2,1X,I2.2,1Hs) 101 FORMAT(6X,18HCorrelation energy,F10.4,5H(MeV)) 102 FORMAT(6X,31HConfiguration the most probable,I9,2H :,/,15X, &1H(,4(I3),3H =>,4(I3),1H),1H(,4(I3),3H =>,4(I3),1H),/, &10X,11HEnergy p-h:,F7.3,4H MeV,8X,13HProbability :,2P F6.2,2H %) 103 FORMAT(6X,5HNext ,I3,21H greatest components:, &15H (probability >,2P F7.3,2H%),/, &200(I3,I9,4H (,4(I3),3H =>,4(I3),1H),1H(,4(I3),3H =>,4(I3),1H), &0 P F12.3,2P F12.6,2H %,/)) 106 FORMAT(6X,41HStates of probability greater than 0,1% :,/,6X, &6HNumber,3X,13HConfiguration,22X,12H Energy(MeV),12H Probability) 107 FORMAT(21HContribution of order,I2,2H :,2P F7.3,1H%,I10,7H states, & 2X,10H>= 1 pair:,2P F7.3,1H%,I9) 108 FORMAT(79(1H-),/,25HNon fondamental contri. :,2P F7.3,1H%, &I10,7H states,2X,10H>= 1 pair:,2P F7.3,1H%,I9) C END C======================================================================== SUBROUTINE HTDAMEPN C======================================================================== INCLUDE 'All.f' INTEGER TIME(3),MAXPIS(200) DIMENSION IEXOP(0:4) C** C WRITE(20,100) PRINT 100 C EXCITATION WINDOW CALL NEWEXCWIN1(DELE,DELE,DELED,DELED) c C Determine THE MANY-BODY BASIS {|npmh>} do i=1,nhpmax nall(i,inp)=0 nss(inp)=0 nfour(inp)=0 end do if (inp.eq.1) then do i=1,nhpmax nallpn(i)=0d0 nsspn(i)=0d0 end do end if nallpn(1)=0 CALL PNBASE c c test with neutron spectra only c c do i=1,nw c sse(i,1,2)=sse(i,1,1) c sse(i,2,2)=sse(i,2,1) c ss2k(i,1,2)=ss2k(i,1,1) c ss2k(i,2,2)=ss2k(i,2,1) c ssip(i,1,2)=ssip(i,1,1) c ssip(i,2,2)=ssip(i,2,1) c end do c do i=-nw,nw c do j=1,nh/2 c do k=1,nl c spf(j,k,i,2)=spf(j,k,i,1) c smf(j,k,i,2)=smf(j,k,i,1) c enddo c end do c end do c IF (INP.EQ.2) THEN CALL HTDAMAPN CALL PNSTATES END IF C 100 FORMAT(2X,34H+Pairing correction by HTDA method) 101 FORMAT(4X, &55H*Calculated the density matrix in single state basis at, &6X,I2.2,1H:,I2.2,1X,I2.2,1Hs) 102 FORMAT(10X,A,26H radius calculated by HTDA,F11.5) END C C======================================================================= C Find all possible configurations (NHOLE holes - NPAR particles) C of good K,pi and of excited energy less than CUTOFF c for pn configurations C======================================================================= SUBROUTINE CONFIG2(KCONF2,IPCONF,NHOLE,NPAR,CUTOFF) c2*Nuclid's quantum number K,nuclid parity,number of holes, number of particles, echtda c cutoff energy INCLUDE 'All.f' C INTEGER NCONFS LOGICAL NEWCON C C EREF: Energy_Reference C NPS: Number of Particle States C NHS: Number of Hole States cc-->> EREF = 0D0 IF (NHOLE-NPAR.EQ.1) THEN EREF=-SSE(NHS,1,INP) ELSEIF (NHOLE-NPAR.EQ.-1) THEN !!!PO CO TO EREF= SSE(1,2,inp) ENDIF C I1 = IHC(1) I2 = IHC(2) I3 = IHC(3) J1 = IPC(1) J2 = IPC(2) J3 = IPC(3) cc-->> c T |Phi0> = + |Phi0> c car T a^+i a^+ibar = a^+ibar (- a^+i) = + a^+i a^+ibar cc<<-- C======================================================================= C |0p-0h> C======================================================================= IF (NHOLE.EQ.0) THEN IF (NPAR.EQ.0) THEN ntotpn=0 c nallpn(1)=0 IF (KCONF2.EQ.0) THEN IF (IPCONF.EQ.1) THEN NALLPN(1) = NALLPN(1) + 1 NSSPN(0) = NSSPN(0) + 1 NTOTPN = NTOTPN + 1 IF (NTOTPN.GT.NE) GOTO 999 DO I = 1, NHPMAX DO J = 1, 2 MBSTP(I,J,NTOTPN)=0 MBSTN(I,J,NTOTPN)=0 ENDDO ENDDO EPHPN(NTOTPN) = 0D0 cc-->> MBREVERSEDPN(NTOTPN)=NTOTPN cc<<-- ENDIF ENDIF GOTO 5 ENDIF ENDIF C======================================================================= C |1p-1h> |1P-1H> ONE PROTON AND ONE NEUTRON C C======================================================================= c c mbmax=0 IF (NHOLE.EQ.2) THEN IF (NPAR.EQ.2) THEN NPSDN = 2*NPLD(1) NHSDN = 2*NHLD(1) NPSDP = 2*NPLD(2) NHSDP = 2*NHLD(2) c DO ICTN = 1,NHSDN !!!neutron hole state if( (mod(ictn,2).eq.1)) mbmin=ntotpn+1 if( (mod(ictn,2).eq.0)) mbmax=ntotpn c DO ICTP = 1, NHSDP !!!proton hole state c KTPN1 = SS2K(ICTN,1,1) + SS2K(ICTP,1,2) !!!hole pn omega DO ICTN1 = 1, NPSDN !!!neutron particle state DO ICTP1 = 1, NPSDP !!!proton particle state C KTPN2 = SS2K(ICTN1,2,1) + SS2K(ICTP1,2,2) !!!particle pn omega IF (KTPN1-KTPN2.EQ.KCONF2) THEN !!!omega condition IPPN1 = SSIP(ICTN,1,1) * SSIP(ICTP,1,2) IPPN2 = SSIP(ICTN1,2,1)* SSIP(ICTP1,2,2) !!!parity of pn configs c IF (IPPN1*IPPN2.EQ.IPCONF) THEN c NALLPN(2) = NALLPN(2) + 1 !!!number of pn configs of type (1p1h)_n (1p1h)_p EEXCI = SSE(ICTN1,2,1) + SSE(ICTP1,2,2) & - SSE(ICTN,1,1) - SSE(ICTP,1,2) !!!excitation energy of this configuration c IF (EEXCI.LE.CUTOFF+EREF) THEN C Check if limit at one pair transfers IF (KOPAIR.EQ.1) THEN IF (.NOT. & (ICTP.EQ.IHPN(ICTN).OR.ICTP.EQ.IHPNTR(ICTN)) ) GO TO 2 IF (.NOT. & (ICTP1.EQ.IPPN(ICTN1).OR.(ICTP1.EQ.IPPNTR(ICTN1)))) GO TO 2 END IF ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c this is the condition to avoid transfering a particle from normal to c time-reversal state if( (mod(ictp+ictp1,2).eq.0).and.(mod(ictn+ictn1,2).eq.0)) then ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc NSSPN(2) = NSSPN(2) + 1 NTOTPN = NTOTPN + 1 IF (NTOTPN.GT.NE) GOTO 999 DO I = 1, NHPMAX do j=1,2 MBSTP(I,J,NTOTPN)=0 MBSTN(I,J,NTOTPN)=0 end do ENDDO C MBSTN(1,1,NTOTPN) = ICTN !!!many body state(particle or hole number, particle(2) or hole(1), config number) MBSTP(1,1,NTOTPN) = ICTP MBSTN(1,2,NTOTPN) = ICTN1 MBSTP(1,2,NTOTPN) = ICTP1 C EPHPN(NTOTPN) = EEXCI c if (mod(ictn,2).eq.0) then jctn1=ictn1+2*mod(ictn1,2)-1 jctp1=ictp1+2*mod(ictp1,2)-1 jctn=ictn+2*mod(ictn,2)-1 jctp=ictp+2*mod(ictp,2)-1 isign=1-2*mod(ictn+ictn1,2) isigp=1-2*mod(ictp+ictp1,2) isig=isigp*isign mbrev=mbmin do while(mbrev.le.mbmax) if((mbstn(1,2,mbrev).eq.jctn1).and. & (mbstp(1,2,mbrev).eq.jctp1).and. & (mbstp(1,1,mbrev).eq.jctp).and. & (mbstn(1,1,mbrev).eq.jctn)) then mbreversedpn(ntotpn)=isig*mbrev mbreversedpn(mbrev)=isig*ntotpn mbrev=mbmax end if mbrev=mbrev+1 end do end if c c c ENDIF ! extra condition for transfer normal to tr states ENDIF !eexci.le.cutoff 2 CONTINUE ENDIF !parity ENDIF !omega ENDDO !particle state p ENDDO !particle state n cc<<-- ENDDO !hole state p ENDDO !hole state n GOTO 5 ENDIF !1p1h ENDIF !1p1h C======================================================================= C |1p-1h> |1P-1H> ONE PROTON AND ONE NEUTRON + 1 particle like pair C +1 neutron pair C======================================================================= c c mbmax=0 IF (NHOLE.EQ.4) THEN IF (NPAR.EQ.4) THEN NPSDN = 2*NPLD(1) NHSDN = 2*NHLD(1) NPSDP = 2*NPLD(2) NHSDP = 2*NHLD(2) c DO ICTN1 = 1,NHSDN !!!neutron hole state DO ICTP1 = 1, NHSDP !!!proton hole state c KTPN1 = SS2K(ICTN1,1,1) + SS2K(ICTP1,1,2) !!!hole pn omega DO ICTN2 = 1, NPSDN !!!neutron particle state DO ICTP2 = 1, NPSDP !!!proton particle state C KTPN2 = SS2K(ICTN2,2,1) + SS2K(ICTP2,2,2) !!!particle pn omega c IPPN1 = SSIP(ICTN1,1,1)* SSIP(ICTP1,1,2) IPPN2 = SSIP(ICTN2,2,1)* SSIP(ICTP2,2,2) !!!parity of pn configs c EEXCI1 = SSE(ICTN2,2,1) + SSE(ICTP2,2,2) & - SSE(ICTN1,1,1) - SSE(ICTP1,1,2) !!!excitation energy of this configuration DO ICTN3=ICTN1+1,NHSDN-1 DO ICTN4=ICTN3+1,NHSDN C KTPN3=SS2K(ICTN3,1,1)+SS2K(ICTN4,1,1) IPPN3=SSIP(ICTN3,1,1)*SSIP(ICTN4,1,1) C DO ICTN5=ICTN2+1,NPSDN-1 DO ICTN6=ICTN5+1,NPSDN C KTPN4=SS2K(ICTN5,2,1)+SS2K(ICTN6,2,1) IPPN4=SSIP(ICTN5,2,1)*SSIP(ICTN6,2,1) C EEXCI2 = SSE(ICTN6,2,1) + SSE(ICTN5,2,2) & - SSE(ICTN4,1,1) - SSE(ICTN3,1,2) !!!excitation energy C EEXCI=EEXCI1+EEXCI2 IF ((KTPN1-KTPN2).EQ.KCONF2.AND.(KTPN3-KTPN4).EQ.KCONF2) THEN !!!omega condition c IF ((IPPN1*IPPN2).EQ.IPCONF.AND.(IPPN3*IPPN4).EQ.IPCONF) THEN !!!parity condition NALLPN(3) = NALLPN(3) + 1 !!!number of pn configs of type (1p1h)_n (1p1h)_p c IF (EEXCI.LE.CUTOFF+EREF) THEN C Check if limit at one pair transfers IF (KOPAIR.EQ.1) THEN IF (.NOT. &(ICTP1.EQ.IHPN(ICTN1).OR.ICTP1.EQ.IHPNTR(ICTN1)) ) GOTO 23 IF (.NOT. &(ICTP2.EQ.IPPN(ICTN2).OR.(ICTP2.EQ.IPPNTR(ICTN2))))GOTO 23 IF (.NOT.(ICTN3+1.EQ.ICTN4.AND.MOD(ICTN4,2).EQ.0)) GOTO 23 IF (.NOT.(ICTN5+1.EQ.ICTN6.AND.MOD(ICTN6,2).EQ.0)) GOTO 23 END IF ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c this is the condition to avoid transfering a particle from normal to c time-reversal state if( (mod(ictp1+ictp2,2).eq.0).and.(mod(ictn1+ictn2,2).eq.0)) then ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc NSSPN(1) = NSSPN(1) + 1 NTOTPN = NTOTPN + 1 IF (NTOTPN.GT.NE) GOTO 999 DO I = 1, NHPMAX do j=1,2 MBSTP(I,J,NTOTPN)=0 MBSTN(I,J,NTOTPN)=0 end do ENDDO C MBSTN(1,1,NTOTPN) = ICTN1 !!!many body state(particle or hole number, particle(2) or hole(1), config number) MBSTP(1,1,NTOTPN) = ICTP1 MBSTN(1,2,NTOTPN) = ICTN2 MBSTP(1,2,NTOTPN) = ICTP2 MBSTN(2,1,NTOTPN) = ICTN3 MBSTN(3,1,NTOTPN) = ICTN4 MBSTN(2,2,NTOTPN) = ICTN5 MBSTN(3,2,NTOTPN) = ICTN6 c EPHPN(NTOTPN) = EEXCI C MBREVERSEDPN(NTOTPN)=NTOTPN c ENDIF ! extra condition for transfer normal to tr states ENDIF !eexci.le.cutoff 23 CONTINUE ENDIF !parity ENDIF !omega ENDDO !particle state p ENDDO !particle state n cc<<-- ENDDO !hole state p ENDDO !hole state n c END DO END DO END DO END DO GOTO 50 ENDIF !1p1h ENDIF !1p1h 50 CONTINUE c WRITE(20,100) & NHOLE,NPAR,NSSPN(NHOLE-3),NALLPN(NHOLE-1) PRINT 100, NHOLE,NPAR,NSSPN(NHOLE-3),NALLPN(NHOLE-1) C======================================================================= C |1p-1h> |1P-1H> ONE PROTON AND ONE NEUTRON + 1 particle like pair C +1 proton pair C======================================================================= c c IF (NHOLE.EQ.4) THEN IF (NPAR.EQ.4) THEN NALLPN(3)=0 NPSDN = 2*NPLD(1) NHSDN = 2*NHLD(1) NPSDP = 2*NPLD(2) NHSDP = 2*NHLD(2) c DO ICTN1 = 1,NHSDN !!!neutron hole state DO ICTP1 = 1, NHSDP !!!proton hole state c KTPN1 = SS2K(ICTN1,1,1) + SS2K(ICTP1,1,2) !!!hole pn omega DO ICTN2 = 1, NPSDN !!!neutron particle state DO ICTP2 = 1, NPSDP !!!proton particle state C KTPN2 = SS2K(ICTN2,2,1) + SS2K(ICTP2,2,2) !!!particle pn omega c IPPN1 = SSIP(ICTN1,1,1)* SSIP(ICTP1,1,2) IPPN2 = SSIP(ICTN2,2,1)* SSIP(ICTP2,2,2) !!!parity of pn configs c EEXCI1 = SSE(ICTN2,2,1) + SSE(ICTP2,2,2) & - SSE(ICTN1,1,1) - SSE(ICTP1,1,2) !!!excitation energy of this configuration DO ICTP3=ICTP1+1,NHSDP-1 DO ICTP4=ICTP3+1,NHSDP C KTPN3=SS2K(ICTP3,1,2)+SS2K(ICTP4,1,2) IPPN3=SSIP(ICTP3,1,2)*SSIP(ICTP4,1,2) C DO ICTP5=ICTP2+1,NPSDP-1 DO ICTP6=ICTP5+1,NPSDP C KTPN4=SS2K(ICTP5,2,2)+SS2K(ICTP6,2,2) IPPN4=SSIP(ICTP5,2,2)*SSIP(ICTP6,2,2) C EEXCI2 = SSE(ICTP6,2,2) + SSE(ICTP5,2,2) & - SSE(ICTP4,1,2) - SSE(ICTP3,1,2) !!!excitation energY C EEXCI=EEXCI1+EEXCI2 IF ((KTPN1-KTPN2).EQ.KCONF2.AND.(KTPN3-KTPN4).EQ.KCONF2) THEN !!!omega condition c IF ((IPPN1*IPPN2).EQ.IPCONF.AND.(IPPN3*IPPN4).EQ.IPCONF) THEN !!!parity condition NALLPN(3) = NALLPN(3) + 1 !!!number of pn configs of type (1p1h)_n (1p1h)_p c IF (EEXCI.LE.CUTOFF+EREF) THEN C Check if limit at one pair transfers IF (KOPAIR.EQ.1) THEN IF (.NOT. &(ICTP1.EQ.IHPN(ICTN1).OR.ICTP1.EQ.IHPNTR(ICTN1)) ) GO TO 25 IF (.NOT. &(ICTP2.EQ.IPPN(ICTN2).OR.(ICTP2.EQ.IPPNTR(ICTN2)))) GO TO 25 IF (.NOT.(ICTP3+1.EQ.ICTP4.AND.MOD(ICTP4,2).EQ.0)) GOTO 25 IF (.NOT.(ICTP5+1.EQ.ICTP6.AND.MOD(ICTP6,2).EQ.0)) GOTO 25 END IF ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c this is the condition to avoid transfering a particle from normal to c time-reversal state if( (mod(ictp1+ictp2,2).eq.0).and.(mod(ictn1+ictn2,2).eq.0)) then ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc NSSPN(3) = NSSPN(3) + 1 NTOTPN = NTOTPN + 1 IF (NTOTPN.GT.NE) GOTO 999 DO I = 1, NHPMAX do j=1,2 MBSTP(I,J,NTOTPN)=0 MBSTN(I,J,NTOTPN)=0 end do ENDDO C MBSTN(1,1,NTOTPN) = ICTN1 !!!many body state(particle or hole number, particle(2) or hole(1), config number) MBSTP(1,1,NTOTPN) = ICTP1 MBSTN(1,2,NTOTPN) = ICTN2 MBSTP(1,2,NTOTPN) = ICTP2 MBSTP(2,1,NTOTPN) = ICTP3 MBSTP(3,1,NTOTPN) = ICTP4 MBSTP(2,2,NTOTPN) = ICTP5 MBSTP(3,2,NTOTPN) = ICTP6 C C EPHPN(NTOTPN) = EEXCI c MBREVERSEDPN(NTOTPN)=NTOTPN c ENDIF ! extra condition for transfer normal to tr states ENDIF !eexci.le.cutoff 25 CONTINUE ENDIF !parity ENDIF !omega ENDDO !particle state p ENDDO !particle state n cc<<-- ENDDO !hole state p ENDDO !hole state n c END DO END DO END DO END DO GOTO 55 ENDIF !1p1h ENDIF !1p1h 55 CONTINUE c WRITE(20,100) & NHOLE,NPAR,NSSPN(NHOLE-1),NALLPN(NHOLE-1) PRINT 100, NHOLE,NPAR,NSSPN(NHOLE-1),NALLPN(NHOLE-1) C======================================================================= C |1p-1h> |1P-1H> ONE PROTON AND ONE NEUTRON !!!4p4h C X|1p-1h> |1P-1H> ONE PROTON AND ONE NEUTRON C======================================================================= c c IF (NHOLE.EQ.4) THEN IF (NPAR.EQ.4) THEN NPSDN = 2*NPLD(1) NHSDN = 2*NHLD(1) NPSDP = 2*NPLD(2) NHSDP = 2*NHLD(2) c DO ICTN1 = 1,NHSDN !!!neutron hole state DO ICTP1 = 1, NHSDP !!!proton hole state KTPN1 = SS2K(ICTN1,1,1)+SS2K(ICTP1,1,2) !!!hole 1 pn omega IPPN1 = SSIP(ICTN1,1,1)*SSIP(ICTP1,1,2) c DO ICTN3 = ICTN1+1,NHSDN !!!neutron 2 hole state DO ICTP3 = ICTP1+1, NHSDP !!!proton 2 hole state KTPN3 = SS2K(ICTN3,1,1)+SS2K(ICTP3,1,2) !!!hole 2 pn omega IPPN3 = SSIP(ICTN3,1,1)*SSIP(ICTP3,1,2) c DO ICTN2 = 1, NPSDN !!!neutron particle state DO ICTP2 = 1, NPSDP !!!proton particle state KTPN2 = SS2K(ICTN2,2,1)+SS2K(ICTP2,2,2) !!!particle 1 pn omega IPPN2 = SSIP(ICTN2,2,1)*SSIP(ICTP2,2,2) DO ICTN4 = ICTN2+1, NPSDN !!!neutron 2 particle state DO ICTP4 = ICTP2+1, NPSDP !!!proton 2 particle state KTPN4 = SS2K(ICTN4,2,1)+ SS2K(ICTP4,2,2) !!!particle 2 pn omega IPPN4 = SSIP(ICTN4,2,1)* SSIP(ICTP4,2,2) c IF ((KTPN1-KTPN2).EQ.KCONF2.AND.(KTPN3-KTPN4).EQ.KCONF2) THEN !!!omega condition c c IF ((IPPN1*IPPN2).EQ.IPCONF.AND.(IPPN3*IPPN4).EQ.IPCONF) THEN EEXCI1 = SSE(ICTN2,2,1) + SSE(ICTP2,2,2) & - SSE(ICTN1,1,1) - SSE(ICTP1,1,2) !!!excitation energy of this configuration c EEXCI2 = SSE(ICTN4,2,1) + SSE(ICTP4,2,2) & - SSE(ICTN3,1,1) - SSE(ICTP3,1,2) !!!excitation energy of this configuration c EEXCI=EEXCI1+EEXCI2 NALLPN(4)=NALLPN(4)+1 c IF (EEXCI.LE.CUTOFF+EREF) THEN C Check if limit at one pair transfers IF (KOPAIR.EQ.1) THEN IF (.NOT. &(ICTP1.EQ.IHPN(ICTN1).OR.ICTP1.EQ.IHPNTR(ICTN1)) ) GO TO 21 IF (.NOT. &(ICTP2.EQ.IPPN(ICTN2).OR.(ICTP2.EQ.IPPNTR(ICTN2)))) GO TO 21 IF (.NOT. &(ICTP3.EQ.IHPN(ICTN3).OR.ICTP3.EQ.IHPNTR(ICTN3)) ) GO TO 21 IF (.NOT. &(ICTP4.EQ.IPPN(ICTN4).OR.(ICTP4.EQ.IPPNTR(ICTN4)))) GO TO 21 END IF ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c this is the condition to avoid transfering a particle from normal to c time-reversal state if( (mod(ictp1+ictp2,2).eq.0).and.(mod(ictn1+ictn2,2).eq.0) &.and.(mod(ictp3+ictp4,2).eq.0) &.and.(mod(ictn3+ictn4,2).eq.0) )then ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc NSSPN(4)=NSSPN(4)+1 NTOTPN = NTOTPN + 1 IF (NTOTPN.GT.NE) GOTO 999 DO I = 1, NHPMAX do j=1,2 MBSTP(I,J,NTOTPN)=0 MBSTN(I,J,NTOTPN)=0 end do ENDDO C MBSTN(1,1,NTOTPN) = ICTN1 !!!many body state MBSTP(1,1,NTOTPN) = ICTP1 MBSTN(1,2,NTOTPN) = ICTN2 MBSTP(1,2,NTOTPN) = ICTP2 C MBSTN(2,1,NTOTPN) = ICTN3 !!!many body state MBSTP(2,1,NTOTPN) = ICTP3 MBSTN(2,2,NTOTPN) = ICTN4 MBSTP(2,2,NTOTPN) = ICTP4 C EPHPN(NTOTPN) = EEXCI c MBREVERSEDPN(NTOTPN)=NTOTPN c ENDIF ! extra condition for transfer normal to tr states ENDIF !eexci.le.cutoff 21 CONTINUE ENDIF !parity ENDIF !omega ENDDO !particle state p ENDDO !particle state n cc<<-- ENDDO !hole state p ENDDO !hole state n CCC ENDDO !particle state p ENDDO !particle state n cc<<-- ENDDO !hole state p ENDDO !hole state n GOTO 5 ENDIF !1p1h ENDIF !1p1h 5 CONTINUE c WRITE(20,100) & NHOLE,NPAR,NSSPN(NHOLE),NALLPN(NHOLE) PRINT 100, NHOLE,NPAR,NSSPN(NHOLE),NALLPN(NHOLE) RETURN 999 WRITE(20,*) 'Number of configuration is too big',NTOTPN,NE CLOSE(20) STOP 'Number of configuration is too big' c 100 FORMAT $(6X,28HNumber of pn configuration |,I1,1Hh,I1,3Hp>:,I7,1H/,I8) END C======================================================================= C Find all possible configurations (NHOLE holes - NPAR particles) C of good K,pi and of excited energy less than CUTOFF C======================================================================= SUBROUTINE CONFIG1(KCONF2,IPCONF,NHOLE,NPAR,CUTOFF) c2*Nuclid's quantum number K,nuclid parity,number of holes, number of particles, echtda c cutoff energy in HTDA must be quite big for precision of Laso method c for Ge nuclei, cutoff=100 gives the same result as cutoff=1000, sufficient 100 INCLUDE 'All.f' C INTEGER NCONFS LOGICAL NEWCON C C EREF: Energy_Reference C NPS: Number of Particle State C NHS: Number of Hole State c NSS(inp)=0 c NCONFS = 0 c NTOTPN=0 c NFOUR(inp)=0 NPS = 2*NPL(inp) NHS = 2*NHL(inp) NPSD = 2*NPLD(inp) NHSD = 2*NHLD(inp) c EREF = 0D0 IF (NHOLE-NPAR.EQ.1) THEN EREF=-SSE(NHS,1,INP) ELSEIF (NHOLE-NPAR.EQ.-1) THEN EREF= SSE(1,2,inp) ENDIF I1 = IHC(1) I2 = IHC(2) I3 = IHC(3) J1 = IPC(1) J2 = IPC(2) J3 = IPC(3) cc-->> c T |Phi0> = + |Phi0> c car T a^+i a^+ibar = a^+ibar (- a^+i) = + a^+i a^+ibar cc<<-- C======================================================================= C |1p-1h> C======================================================================= IF (NHOLE.EQ.1) THEN IF (NPAR.EQ.1) THEN DO ICT1 = 1, NHS cc-->> if (mod(ict1,2).eq.1) mbmin=ntotpn+1 if (mod(ict1,2).eq.0) mbmax=ntotpn cc<<-- DO ICT2 = 1, NPS IF (SS2K(ICT2,2,inp)-SS2K(ICT1,1,inp).EQ.KCONF2) THEN IF (SSIP(ICT2,2,inp)*SSIP(ICT1,1,inp).EQ.IPCONF) THEN NALL(1,inp) = NALL(1,inp) + 1 EEXCI = SSE(ICT2,2,inp) - SSE(ICT1,1,inp) IF (EEXCI.LE.CUTOFF+EREF) THEN NSS(inp) = NSS(inp) + 1 NTOTPN = NTOTPN+ 1 IF (NTOTPN.GT.NE) GOTO 333 DO I = 1, NHPMAX DO J = 1, 2 if (inp.eq.1) then MBSTN(I,J,NTOTPN) = 0 else MBSTP(i,j,ntotpn) = 0 end if ENDDO ENDDO if (inp.eq.1) then MBSTn(1,1,NTOTPN) = ICT1 MBSTn(1,2,NTOTPN) = ICT2 else mbstp(1,1,ntotpn) = ict1 mbstp(1,2,ntotpn) = ict2 end if EPHPN(ntotpn) = EEXCI cc-->> if (mod(ict1,2).eq.0) then jct2=ict2 + 2*mod(ict2,2)-1 isig=1-2*mod(ict1+ict2,2) mbrev=mbmin do while (mbrev.le.mbmax) if (inp.eq.1) then if (mbstn(1,2,mbrev).eq.jct2) then mbreversedpn(ntotpn)=isig*mbrev mbreversedpn(mbrev)=isig*ntotpn mbrev=mbmax endif else if (mbstp(1,2,mbrev).eq.jct2) then mbreversedpn(ntotpn)=isig*mbrev mbreversedpn(mbrev)=isig*ntotpn mbrev=mbmax endif end if mbrev=mbrev+1 enddo endif cc<<-- ENDIF !eexci.le.cutoff ENDIF !parity ENDIF !omega ENDDO !particle state ENDDO !hole state GOTO 35 ENDIF !1p ENDIF !1h 35 continue if (nhole.eq.1) then if (inp.eq.1) then WRITE(20,101) & NHOLE,NPAR,NSS(inp),NALL(1,inp) PRINT 101, NHOLE,NPAR,NSS(inp),NALL(1,inp) else WRITE(20,100) & NHOLE,NPAR,NSS(inp),NALL(1,inp) PRINT 100, NHOLE,NPAR,NSS(inp),NALL(1,inp) end if end if C======================================================================= C |2p-2h> C !!!! uniquement dans la petite fenetre !!!! C======================================================================= IF (NHOLE.EQ.2) THEN IF (NPAR.EQ.2) THEN DO ICT1 = NHS-NHSD+1, NHS- 1 cc-->> if (mod(ict1,2).eq.1) mbmin=NTOTPN+1 cc<<-- DO ICT2 = ICT1+1, NHS KT12 = SS2K(ICT1,1,inp) + SS2K(ICT2,1,inp) DO ICT3 = 1, NPSD - 1 DO ICT4 = ICT3 + 1, NPSD KT34 = SS2K(ICT3,2,inp) + SS2K(ICT4,2,inp) IF (KT34-KT12.EQ.KCONF2) THEN IP12 = SSIP(ICT1,1,inp) * SSIP(ICT2,1,inp) IP34 = SSIP(ICT3,2,inp) * SSIP(ICT4,2,inp) IF (IP12*IP34.EQ.IPCONF) THEN NALL(2,Inp) = NALL(2,inp) + 1 EEXCI = SSE(ICT4,2,inp) + SSE(ICT3,2,inp) & - SSE(ICT2,1,inp) - SSE(ICT1,1,inp) IF (EEXCI.LE.CUTOFF+EREF) THEN C Check if limit at one pair transfers IF (KOPAIR.EQ.1) THEN IF (.NOT.(ICT1+1.EQ.ICT2.AND.MOD(ICT2,2).EQ.0)) GOTO 2 IF (.NOT.(ICT3+1.EQ.ICT4.AND.MOD(ICT4,2).EQ.0)) GOTO 2 ENDIF NSS(inp) = NSS(inp) + 1 NTOTPN = NTOTPN + 1 IF (NTOTPN.GT.NE) GOTO 333 DO I = 1, NHPMAX DO J = 1, 2 MBSTN(I,J,NTOTPN) = 0 MBSTP(I,J,NTOTPN) = 0 ENDDO ENDDO IF (INP.EQ.1) THEN MBSTn(2,1,NTOTPN ) = ICT1 MBSTN(1,1,NTOTPN ) = ICT2 MBSTN(2,2,NTOTPN ) = ICT3 MBSTN(1,2,NTOTPN ) = ICT4 ELSE MBSTP(2,1,NTOTPN ) = ICT1 MBSTP(1,1,NTOTPN ) = ICT2 MBSTP(2,2,NTOTPN ) = ICT3 MBSTP(1,2,NTOTPN ) = ICT4 END IF C EPHPN(NTOTPN) = EEXCI cc-->> jct2=ict2 + 2*mod(ict2,2)-1 jct4=ict4 + 2*mod(ict4,2)-1 if (mod(ict1,2).eq.1) then if (ict1.eq.jct2.and.ict3.eq.jct4) then mbreversedpn(ntotPN)=ntotPN ! 1paire_h 1paire_p else mbreversedpn(ntotPN)=0 ! unknown yet endif else jct3=ict3 + 2*mod(ict3,2)-1 if (jct3.gt.jct4) then jct3=jct4 jct4=ict3 + 2*mod(ict3,2)-1 end if mbrev=mbmin if (inp.eq.1) then do while (mbrev.le.mbmax) if ((mbstn(1,1,mbrev).eq.jct2).and. #(mbstn(2,2,mbrev).eq.jct3.and.mbstn(1,2,mbrev).eq.jct4)) then isig=1-2*mod(ict1+ict2,2) ! 0paire_h 1paire_p if ((mod(jct3,2).ne.1).or.(jct4.ne.jct3+1)) # isig=isig*(1-2*mod(ict3+ict4,2))! 0paire_h 0paire_p mbreversedpn(ntotpn)=isig*mbrev mbreversedpn(mbrev)=isig*ntotPN mbrev=mbmax endif mbrev=mbrev+1 enddo end if c if (inp.eq.2) then do while (mbrev.le.mbmax) if ((mbstp(1,1,mbrev).eq.jct2).and. #(mbstp(2,2,mbrev).eq.jct3.and.mbstp(1,2,mbrev).eq.jct4)) then isig=1-2*mod(ict1+ict2,2) ! 0paire_h 1paire_p if ((mod(jct3,2).ne.1).or.(jct4.ne.jct3+1)) # isig=isig*(1-2*mod(ict3+ict4,2))! 0paire_h 0paire_p mbreversedpn(ntotpn)=isig*mbrev mbreversedpn(mbrev)=isig*ntotPN mbrev=mbmax endif mbrev=mbrev+1 enddo end if endif cc<<-- ENDIF 2 CONTINUE ENDIF ENDIF ENDDO ENDDO cc-->> if (inp.eq.1) then if ((mod(ict1,2).eq.1).and.(ict2.eq.ict1+1)) then do mb1=mbmin,ntotPN if (mbreversedpn(mb1).eq.0) then jct3=mbstn(1,2,mb1)+2*mod(mbstn(1,2,mb1),2)-1 jct4=mbstn(2,2,mb1)+2*mod(mbstn(2,2,mb1),2)-1 isig=(1-2*mod(ict3+ict4,2)) ! 1paire_h 0paire_p mbrev=mbmin do while (mbrev.le.ntotPN) if ((mbstn(1,2,mbrev).eq.jct3).and. # (mbstn(2,2,mbrev).eq.jct4)) then mbreversedpn(mb1)=isig*mbrev mbreversedpn(mbrev)=isig*mb1 mbrev=ntotPN endif mbrev=mbrev+1 enddo endif enddo endif end if if (inp.eq.2) then if ((mod(ict1,2).eq.1).and.(ict2.eq.ict1+1)) then do mb1=mbmin,ntotPN if (mbreversedpn(mb1).eq.0) then jct3=mbstp(1,2,mb1)+2*mod(mbstp(1,2,mb1),2)-1 jct4=mbstp(2,2,mb1)+2*mod(mbstp(2,2,mb1),2)-1 isig=(1-2*mod(ict3+ict4,2)) ! 1paire_h 0paire_p mbrev=mbmin do while (mbrev.le.ntotPN) if ((mbstp(1,2,mbrev).eq.jct3).and. # (mbstp(2,2,mbrev).eq.jct4)) then mbreversedpn(mb1)=isig*mbrev mbreversedpn(mbrev)=isig*mb1 mbrev=ntotPN endif mbrev=mbrev+1 enddo endif enddo endif end if cc<<-- ENDDO cc-->> if (mod(ict1,2).eq.1) then mbmax=ntotPN endif cc<<-- ENDDO GOTO 5 ENDIF ENDIF 5 CONTINUE if (nhole.eq.2) then if (inp.eq.1) then WRITE(20,101) & NHOLE,NPAR,NSS(inp),NALL(2,inp) PRINT 101, NHOLE,NPAR,NSS(inp),NALL(2,inp) else WRITE(20,100) & NHOLE,NPAR,NSS(inp),NALL(2,inp) PRINT 100, NHOLE,NPAR,NSS(inp),NALL(2,inp) end if end if 100 FORMAT $(6X,28HNumber of pp configuration |,I1,1Hh,I1,3Hp>:,I7,1H/,I8) 101 FORMAT $(6X,28HNumber of nn configuration |,I1,1Hh,I1,3Hp>:,I7,1H/,I8) C======================================================================= C |2p-2h> \times |2p2h> two pairs C !!!! uniquement dans la petite fenetre !!!! C======================================================================= IF (NHOLE.EQ.4) THEN IF (NPAR.EQ.4) THEN c DO ICT1 = NHS-NHSD+1, NHS- 1 DO ICT2 = ICT1+1, NHS KT12 = SS2K(ICT1,1,inp) + SS2K(ICT2,1,inp) DO ICT3 = 1, NPSD - 1 DO ICT4 = ICT3+1, NPSD KT34 = SS2K(ICT3,2,inp) + SS2K(ICT4,2,inp) IF (KT34-KT12.EQ.KCONF2) THEN IP12 = SSIP(ICT1,1,inp) * SSIP(ICT2,1,inp) IP34 = SSIP(ICT3,2,inp) * SSIP(ICT4,2,inp) IF (IP12*IP34.EQ.IPCONF) THEN EEXCI1 = SSE(ICT4,2,inp) + SSE(ICT3,2,inp) & - SSE(ICT2,1,inp) - SSE(ICT1,1,inp) c DO ICT5 = ICT1+2, NHS-1 DO ICT6 = ICT5+1, NHS KT56 = SS2K(ICT5,1,inp) + SS2K(ICT6,1,inp) DO ICT7 = ICT3+2, NPSD - 1 DO ICT8 = ICT7+1, NPSD KT78 = SS2K(ICT7,2,inp) + SS2K(ICT8,2,inp) IF (KT78-KT56.EQ.KCONF2) THEN IP56 = SSIP(ICT5,1,inp) * SSIP(ICT6,1,inp) IP78 = SSIP(ICT7,2,inp) * SSIP(ICT8,2,inp) IF (IP56*IP78.EQ.IPCONF) THEN EEXCI2 = SSE(ICT8,2,inp) + SSE(ICT7,2,inp) & - SSE(ICT6,1,inp) - SSE(ICT5,1,inp) C if (eexci1.eq.eexci2) go to 20 !not the same pair twice NALL(4,inp) = NALL(4,inp) + 1 EEXCI=EEXCI1+EEXCI2 C IF (EEXCI.LE.CUTOFF+EREF) THEN C Check if limit at one pair transfers IF (KOPAIR.EQ.1) THEN IF (.NOT.(ICT1+1.EQ.ICT2.AND.MOD(ICT2,2).EQ.0)) GOTO 20 IF (.NOT.(ICT3+1.EQ.ICT4.AND.MOD(ICT4,2).EQ.0)) GOTO 20 IF (.NOT.(ICT5+1.EQ.ICT6.AND.MOD(ICT6,2).EQ.0)) GOTO 20 IF (.NOT.(ICT7+1.EQ.ICT8.AND.MOD(ICT8,2).EQ.0)) GOTO 20 ENDIF Nfour(inp) = Nfour(inp) + 1 NTOTPN = NTOTPN + 1 IF (NTOTPN.GT.NE) GOTO 333 DO I = 1, NHPMAX DO J = 1, 2 if (inp.eq.1) then MBSTn(I,J,NTOTPN) = 0 else mbstp(i,j,NTOTPN) = 0 end if ENDDO ENDDO if (inp.eq.1) then MBSTN(2,1,NTOTPN) = ICT1 MBSTN(1,1,NTOTPN) = ICT2 MBSTN(2,2,NTOTPN) = ICT3 MBSTN(1,2,NTOTPN) = ICT4 mbstn(4,1,ntotpn) = ict5 mbstn(3,1,ntotpn) = ict6 mbstn(4,2,ntotpn) = ict7 mbstn(3,2,ntotpn) = ict8 ELSE MBSTP(2,1,NTOTPN) = ICT1 MBSTP(1,1,NTOTPN) = ICT2 MBSTP(2,2,NTOTPN) = ICT3 MBSTP(1,2,NTOTPN) = ICT4 mbstp(4,1,ntotpn) = ict5 mbstp(3,1,ntotpn) = ict6 mbstp(4,2,ntotpn) = ict7 mbstp(3,2,ntotpn) = ict8 END IF C EPHPN(NTOTPN) = EEXCI c mbreversedpn(ntotpn)=ntotpn cc<<-- ENDIF !eexci 20 CONTINUE END IF !SPIN END IF !PARITY END DO !ICT8 END DO !ICT7 END DO !ICT6 END DO !ICT5 C ENDIF !spin ENDIF !parity ENDDO !ict4 ENDDO !ict3 ENDDO !ict2 ENDDO ! ict1 GOTO 50 ENDIF !4h ENDIF !4p c 50 CONTINUE c c RETURN c 333 WRITE(*,*) 'Number of Configuration is too big',NTOTPN,NE c STOP 'Number of configuration is too big' c END C======================================================================= C |2p-2h> X |2p-2h> only two pair transfers case C C======================================================================= IF (NHOLE.EQ.8) THEN IF (NPAR.EQ.8) THEN NPSN = 2*NPL(1) NHSN = 2*NHL(1) NPSDN = 2*NPLD(1) NHSDN = 2*NHLD(1) c NPSP = 2*npl(2) nhsp = 2*nhl(2) npsdp=2*npld(2) nhsdp=2*nhld(2) c 2p2h-neutrons DO ICT1N = NHSN-NHSDN+1, NHSN- 1 DO ICT2N = ICT1N+1, NHSN KT12 = SS2K(ICT1N,1,1) + SS2K(ICT2N,1,1) DO ICT3N = 1, NPSDN - 1 DO ICT4N = ICT3N + 1, NPSDN KT34 = SS2K(ICT3N,2,1) + SS2K(ICT4N,2,1) IF (KT34-KT12.EQ.KCONF2) THEN IP12 = SSIP(ICT1N,1,1) * SSIP(ICT2N,1,1) IP34 = SSIP(ICT3N,2,1) * SSIP(ICT4N,2,1) IF (IP12*IP34.EQ.IPCONF) THEN EEXCIN= SSE(ICT4N,2,1) + SSE(ICT3N,2,1) & - SSE(ICT2N,1,1) - SSE(ICT1N,1,1) c 2p2h-protons DO ICT1P = NHSP-NHSDP+1, NHSP- 1 DO ICT2P = ICT1P+1, NHSP KT12 = SS2K(ICT1P,1,2) + SS2K(ICT2P,1,2) DO ICT3P = 1, NPSDP - 1 DO ICT4P = ICT3P + 1, NPSDP KT34 = SS2K(ICT3P,2,2) + SS2K(ICT4P,2,2) IF (KT34-KT12.EQ.KCONF2) THEN IP12 = SSIP(ICT1P,1,2) * SSIP(ICT2P,1,2) IP34 = SSIP(ICT3P,2,2) * SSIP(ICT4P,2,2) IF (IP12*IP34.EQ.IPCONF) THEN EEXCIP = SSE(ICT4P,2,2) + SSE(ICT3P,2,2) & - SSE(ICT2P,1,2) - SSE(ICT1P,1,2) EEXCI=EEXCIN+EEXCIP IF (EEXCI.LE.CUTOFF+EREF) THEN c Check if limit at one pair transfers IF (KOPAIR.EQ.1) THEN IF (.NOT.(ICT1N+1.EQ.ICT2N.AND.MOD(ICT2N,2).EQ.0)) GOTO 21 IF (.NOT.(ICT3N+1.EQ.ICT4N.AND.MOD(ICT4N,2).EQ.0)) GOTO 21 IF (.NOT.(ICT1P+1.EQ.ICT2P.AND.MOD(ICT2P,2).EQ.0)) GOTO 21 IF (.NOT.(ICT3P+1.EQ.ICT4P.AND.MOD(ICT4P,2).EQ.0)) GOTO 21 ENDIF NTOTPN = NTOTPN + 1 NFOUR(2)=NFOUR(2)+1 IF (NTOTPN.GT.NE) GOTO 333 DO I = 1, NHPMAX DO J = 1, 2 MBSTN(I,J,NTOTPN) = 0 MBSTP(I,J,NTOTPN) = 0 ENDDO ENDDO MBSTn(2,1,NTOTPN ) = ICT1N MBSTN(1,1,NTOTPN ) = ICT2N MBSTN(2,2,NTOTPN ) = ICT3N MBSTN(1,2,NTOTPN ) = ICT4N MBSTP(2,1,NTOTPN ) = ICT1P MBSTP(1,1,NTOTPN ) = ICT2P MBSTP(2,2,NTOTPN ) = ICT3P MBSTP(1,2,NTOTPN ) = ICT4P c EPHPN(NTOTPN) = EEXCI c mbreversedpn(ntotPN)=ntotPN ! 1paire_h 1paire_p cc<<-- ENDIF !cutoff for EPH 21 CONTINUE endif !parity N endif !parity P end do !ictp4 end do !ictp3 end do !ictp2 end do !ictp1 ENDIF !omega N ENDIF !parity N ENDDO !ictn4 ENDDO !ictn3 ENDDO !ictn2 ENDDO !ictn1 GOTO 51 ENDIF !4 holes ENDIF !4part ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 51 continue if (nhole.eq.8) print*, 'Number of nn x pp configs ',nfour(2) if (nhole.eq.4) then if (inp.eq.1) then PRINT 105, NHOLE,NPAR,NFOUR(INP),NALL(NHOLE,INP) write(20,105) nhole,npar,NFOUR(INP),NALL(NHOLE,INP) else PRINT 107, NHOLE,NPAR,NFOUR(INP),NALL(NHOLE,INP) write(20,107) nhole,npar,NFOUR(INP),NALL(NHOLE,INP) end if end if c 105 FORMAT $(6X,28HNumber of nn configuration |,I1,1Hh,I1,3Hp>:,I7,1H/,I8) 107 FORMAT $(6X,28HNumber of nn configuration |,I1,1Hh,I1,3Hp>:,I7,1H/,I8) RETURN 333 WRITE(*,*) 'Number of Configuration is too big',NTOTPN,NE STOP 'Number of configuration is too big' END