C============================================================================== C This is the main program for HTDA calculations. C code1 C - 25/9/2003 : Can stop program using STOP file 3.1.1 C - 25/9/2003 : partial energy Correction 4.0 C - 29/9/2003 : 5.0 like 3.x but the subroutine WRFIEL is add C - 8/10/2003 : 6.0 Calculate the spatial part of sp. wave function C in different direction C - 12/11/2003: Detect one error in the determination of C 3h-2p configuration (IP1 instead of IP12) C - 19/11/2003: 7.0 role of 2 paires C - 15/04/2004: 8.0 archive version before transfering the code C============================================================================== c -hfaxial.f( fenetre d'excitation+les etats excites avec laso +le temps de calcul c des elements de matrice) pour comparer avec hfaxial3.f C======================================================================= C Version à deux fenetres avec interactions delta, Q.Q et O.O C======================================================================= C Initialition of constants : Pi,4Pi,eV,2**x,n!,options messages C======================================================================= cbk c canonical basis-corrected; epair replaced; c LAnczos from Etienne CAurrier cbk c======================================================================= c version with pn pairing ; some options and routins for RPA calculations c removed or do not work; c pn pairing up to 2p2h excitations; c to check: rinewpn, emadpn, densities calculations c c======================================================================= SUBROUTINE INICTS INCLUDE 'All.f' DATA NCL,PTY /8Hneutrons,8Hprotons ,1H-,1H0,1H+/ DATA SKM /4HSIII,4HSKM*,4HSLY4/ C PI = 4D0*DATAN(1D0) QP = 4D0*PI EV = 1439965D-6 C FACT: Factorial C PO2: Power Of 2 FACT(0) = 1D0 PO2(0) = 1D0 DO I = 1, NZ FACT(I) = FACT(I-1)*DBLE(I) PO2(I) = PO2(I-1)*2D0 ENDDO C OPTS(0,1)='calculated' OPTS(1,1)='read from file' OPTS(0,2)='||Take all 2p2h configurations' OPTS(1,2)='||Limit at one pair transfers' OPTS(0,9)='delta' OPTS(1,9)='delta+QQ' OPTS(2,9)='delta+OO' OPTS(3,9)='delta+QQ+OO' opts(4,9)='none' opts(5,9)='QQ' opts(6,9)='OO' opts(7,9)='QQ+OO' OPTS(0,8)=' None' ! 0000 OPTS(1,8)=' R2' ! 0001 OPTS(2,8)=' Q20' ! 0010 OPTS(3,8)=' R2 and Q20' ! 0011 OPTS(4,8)=' Q40' ! 0100 OPTS(5,8)=' Q40 and R2' ! 0101 C END C======================================================================= C= Default values for variables C======================================================================= SUBROUTINE INITDA INCLUDE 'All.f' C C--------------------------------------------- C Nuclid being calculated by default is Hf-178 C--------------------------------------------- NUN(1) = 106 NUN(2) = 72 C----- C SIII C----- NSKY = 1 C------------------------ C Parameters of the basis C------------------------ N0 = 12 B0 = 4332D-4 Q = 11204D-4 C------------------------------------ C Degree of Gauss intergration method C------------------------------------ NGH = 20 NGL = 10 NGHP = NGH/2 C---------------------------------- C Do a HF calculation (without BCS) C---------------------------------- MNHF = 1000 MNHTDA = 1 MNBCS = 1000 C---------------------------------------------------- C No contraints on quadruple and hexadecapole moments C---------------------------------------------------- C2 = 0D0 Q2 = 0D0 C4 = 0D0 Q4 = 0D0 C------------------------------ C Parameters of BCS calculation C------------------------------ G(1) = 161D-1 G(2) = 133D-1 ECBCS = 65D-1 C------------------------------- C Parameters of HTDA calculation C------------------------------- C Intensities of Delta force FHTDA(1) = -535D0 FHTDA(2) = -690D0 FHTDA1=0D0 FHTDA0=0D0 C Excitation Window pour RPA (grande fenetre) DELE = 6.5D0 C fenetre réduite pour PAIRING (petite fenetre) DELED=6.5D0 ECHTDA = 10D0 WFX = 6D0 WFMU = 2D-1 C EIDEN, EPHTDA EIDEN = 1D0 EPHTDA = 1D-15 C------------------ C Mixing parameters C------------------ HFM = 5D-1 HTDAM = 5D-1 C------------------ C Execution options C------------------ C Initial Potential: Wood-Saxon KOIP = 0 C******** c interaction résiduelle: DELTA KORINT = 0 C------------------------------ C Take all 2p2h configuration KOPAIR = 0 c c Take pn coupling into account KOPN=1 c C-------------------- C Converged condition C-------------------- EPSE = 5D-6 EPQ20 = 1D-1 EPQ40 = 1D-1 C--------------------------------------------- C K/pi configuration condition of each isospin C--------------------------------------------- NUK2(1) = 0 NUIP(1) = 1 NUK2(2) = 0 NUIP(2) = 1 C----------------- C Don't do simplex C----------------- ICTRL = 0 C----------------------------------------------------------------- C Number of eigenvector to be calculated: for neutrons and protons C----------------------------------------------------------------- NEIGEN=1 C c---------------------------------------------------------------- c Calcul des matrices HTDA c---------------------------------------------------------------- KOMAT=0 C END C======================================================================= C= Input from the file "data" this will overwrite default values C======================================================================= SUBROUTINE DATAIN INCLUDE 'All.f' CHARACTER * 80 GLINE C OPEN(UNIT=30,FILE='data',STATUS='OLD',ERR=999) PRINT *, 'File "data" exist' WRITE(20,*) 'File "data" exist' C GLINE : Guide Line C------------------------------- C Determine the nuclid to calcul C------------------------------- GLINE = '' REWIND 30 1 READ(30,100,END=203) GLINE IF (INDEX(GLINE,'Nuclid').NE.0) GOTO 2 GOTO 1 2 READ(30,*,END=203,ERR=203) NUN C----------------------------------- C Determine the type of Skyrme force C----------------------------------- 203 REWIND 30 3 READ(30,100,END=205) GLINE IF (INDEX(GLINE,'Skyrme force').NE.0) GOTO 4 GOTO 3 4 READ(30,*,END=205,ERR=205) NSKY C------------------------------- C Determine the basis parameters C------------------------------- 205 REWIND 30 5 READ(30,100,END=207) GLINE IF (INDEX(GLINE,'Truncated').NE.0) GOTO 6 GOTO 5 6 READ(30,*,END=207,ERR=207) N0,B0,Q C-------------------------------------------------- C Determine the degree of Gauss intergration method C-------------------------------------------------- 207 REWIND 30 7 READ(30,100,END=209) GLINE IF (INDEX(GLINE,'Gauss').NE.0) GOTO 8 GOTO 7 8 READ(30,*,END=209,ERR=209) NGH,NGL C-------------------------------------------------------- C Determine the maximum number of each kind of iteration: C HF, BCS, HTDA C-------------------------------------------------------- 209 REWIND 30 9 READ(30,100,END=211) GLINE IF (INDEX(GLINE,'iterations').NE.0) GOTO 10 GOTO 9 10 READ(30,*,END=211,ERR=211) MNHF,MNHTDA,MNBCS C------------------------------------- C Determine the deformation contraints C quadruple and hexadecapole moments C------------------------------------- 211 REWIND 30 11 READ(30,100,END=213) GLINE IF (INDEX(GLINE,'contraints').NE.0) GOTO 12 GOTO 11 12 READ(30,*,END=213,ERR=213) C2,Q2,C4,Q4 C-------------------------------------------------------- C Determine the parameters of BCS calculation G0, CUT-OFF C We used BCS method inside a HF calculation C-------------------------------------------------------- 213 IF (MNHF.GT.0) THEN IF (MNBCS.GT.0) THEN REWIND 30 13 READ(30,100,END=215) GLINE IF (INDEX(GLINE,'G0').NE.0) GOTO 14 GOTO 13 14 READ(30,*,END=215,ERR=215) G,ECBCS ENDIF ENDIF C--------------------------------------------- C Determine the parameters of HTDA calculation C--------------------------------------------- 215 IF (MNHTDA.GT.0) THEN C Intensity of the Delta force REWIND 30 15 READ(30,100,END=217) GLINE IF (INDEX(GLINE,'Delta force').NE.0) GOTO 16 GOTO 15 16 READ(30,*,END=217,ERR=217) FHTDA,FHTDA1,FHTDA0 C Excitation Window dans le cas de delta 217 REWIND 30 17 READ(30,100,END=219) GLINE IF (INDEX(GLINE,'Window').NE.0) GOTO 18 GOTO 17 18 READ(30,*,END=219,ERR=219) DELE C Parameters for Many-body matrix element diagonalisation C (EIDEN, EPHTDA) 219 REWIND 30 19 READ(30,100,END=221) GLINE IF (INDEX(GLINE,'EPHTDA').NE.0) GOTO 20 GOTO 19 20 READ(30,*,END=221,ERR=221) EIDEN,EPHTDA ENDIF C--------------------------- C Determine mixing parameter C--------------------------- 221 REWIND 30 21 READ(30,100,END=223) GLINE IF (INDEX(GLINE,'Mixing').NE.0) GOTO 22 GOTO 21 22 READ(30,*,END=223,ERR=223) HFM,HTDAM C---------------------------- C Determine execution options C---------------------------- C Initial Potential 223 REWIND 30 23 READ(30,100,END=231) GLINE IF (INDEX(GLINE,'Woods-Saxon').NE.0) GOTO 24 GOTO 23 24 READ(30,*,END=231,ERR=231) KOIP C------------------------------------------ C Which kind of residual interaction to use C------------------------------------------ C KORINT = 0-3 we use delta C KORINT = 1,3,5,7 we use Q.Q C KORINT = 2,3,6,7 we use O.O 229 IF (MNHTDA.GT.0) THEN REWIND 30 29 READ(30,100,END=231) GLINE IF (INDEX(GLINE,'Residual').NE.0) GOTO 30 GOTO 29 30 READ(30,*,END=231,ERR=231) KORINT ENDIF C------------------------------ C Determine converged condition C------------------------------ 231 REWIND 30 31 READ(30,100,END=233) GLINE IF (INDEX(GLINE,'Converged').NE.0) GOTO 32 GOTO 31 32 READ(30,*,END=233,ERR=233) EPSE,EPQ20,EPQ40 C----------------------------------------------------------- C Determine the K/pi configuration condition of each isospin C For HTDA calculation only C----------------------------------------------------------- 233 REWIND 30 IF (MNHTDA.GT.0) THEN 33 READ(30,100,END=235) GLINE IF (INDEX(GLINE,'K/pi').NE.0) GOTO 34 GOTO 33 34 READ(30,*,END=235,ERR=235) NUK2(1), NUIP(1), NUK2(2), NUIP(2) ENDIF C--------------------------------------- C Check if one want to do simplex or not C Not for HTDA calculation C--------------------------------------- 235 REWIND 30 C IF (MNHTDA.EQ.0) THEN 35 READ(30,100,END=239) GLINE IF (INDEX(GLINE,'simplex').NE.0) GOTO 36 GOTO 35 36 READ(30,*,END=239,ERR=239) ICTRL C ENDIF C------------------------------------------------ C Determine number of eigenvalue to be calculated C------------------------------------------------ 239 IF (MNHTDA.GT.0) THEN REWIND 30 39 READ(30,100,END=247) GLINE IF (INDEX(GLINE,'eigen').NE.0) GOTO 40 GOTO 39 40 READ(30,*,END=247,ERR=247) NEIGEN,ZPREC,NLOOP,IFIL ENDIF C------------------------------------------------------------- C Limit at one pair transfers or not? C------------------------------------------------------------- 247 REWIND 30 47 READ(30,100,END=249) GLINE IF (INDEX(GLINE,'one pair').NE.0) GOTO 48 GOTO 47 48 READ(30,*,END=249,ERR=249) KOPAIR write(*,*) 'kopair=',kopair C------------------------------------------------------------- C------------------------------------------------------------------- 249 REWIND 30 49 READ(30,100,END=251) GLINE IF (INDEX(GLINE,'REDUITE').NE.0) GOTO 50 GOTO 49 50 READ(30,*,END=251,ERR=251) DELED,ECHTDA,WFX,WFMU c******************************************************************** c Calcul des matrices HTDA c******************************************************************* 251 REWIND 30 51 READ(30,100,END=253) GLINE IF (INDEX(GLINE,'Matrix').NE.0) GOTO 52 GOTO 51 52 READ(30,*,END=253,ERR=253) KOMAT c cbk Take pn coupling into account or not c 253 REWIND 30 53 READ(30,100,END=990) GLINE IF(INDEX(GLINE, 'PN').NE.0) GOTO 54 GO TO 53 54 READ(30,*,END=990,ERR=990) KOPN cbk cbk 990 CLOSE(30) RETURN 999 PRINT *,'File "data" does not exist' WRITE(20,*) 'File "data" does not exist' C 100 FORMAT(A) C if (DELE.LE.DELED) DELE=DELED END C======================================================================= C= Get informations of the nuclid whose proton number is NUN(2): C= Name, Chemical Symbol C= from the data file "nuclids" in "fixedata" directory C= C= Output: NAME, CSYM C= In the case the information is not found, return: C= NAME = 'Unknown', CSYM = ' X' C======================================================================= SUBROUTINE INFON INCLUDE 'All.f' C IF (NUN(1).LT.0) THEN PRINT *, 'Number of neutrons is less than 0' WRITE(20,*) 'Number of neutrons is less than 0 - Stop' CLOSE (20) STOP 'Fatal error' ELSEIF (NUN(2).LT.0) THEN PRINT *, 'Number of protons is less than 0' WRITE(20,*) 'Number of protons is less than 0 - Stop' CLOSE (20) STOP 'Fatal error' ELSEIF (NUN(1)+NUN(2).EQ.0) THEN PRINT *, 'Mass number is equal to 0' WRITE(20,*) 'Mass number is equal to 0 - Stop' CLOSE (20) STOP 'Fatal error' ENDIF C OPEN(UNIT=30,FILE = 'fixedata/nuclids',STATUS='OLD',ERR=999) C 1 READ(30, 100, END = 2,ERR = 2) NPTMP, CSYM, NAME IF (NPTMP.EQ.NUN(2)) THEN CLOSE(30) RETURN ENDIF GOTO 1 2 CSYM = ' X' NAME = 'Unknown' CLOSE(30) RETURN 999 PRINT *,'Error opening file "nuclids" in directory "fixedata"' PRINT *,'File is not existed' PRINT *,'' WRITE(20,*) 'Error opening file "nuclids" in directory "fixedata"' WRITE(20,*) 'File is not existed' WRITE(20,*) CSYM = ' X' NAME = 'Unknown' RETURN C 100 FORMAT(I3,2X,A,2X,A) C END C======================================================================= C= Initialize the parameters coresponding to the choosen Skyrme force C======================================================================= SUBROUTINE SKYRME(N) INCLUDE 'All.f' C IF (N.EQ.1) THEN CCC SIII (Nucl. Phys. A264, 197) T0 =-112875D-2 T1 = 395D0 T2 =-95D0 T3 = 14000D0 X0 = 45D-2 X1 = 0D0 X2 = 0D0 X3 = 1D0 W = 120D0 XBET3 = 1D0 ELSEIF (N.EQ.2) THEN CCC SKM* (Nucl. Phys. A635, 231) T0 =-2645D0 T1 = 410D0 T2 =-135D0 T3 = 15595D0 X0 = 9D-2 X1 = 0D0 X2 = 0D0 X3 = 0D0 W = 130D0 XBET3 = 1D0/6D0 ELSEIF (N.EQ.3) THEN CCC SLY4 (Nucl. Phys. A635, 231) T0 =-248891D-2 T1 = 48682D-2 T2 =-54639D-2 T3 = 13777D0 X0 = 834D-3 X1 =-344D-3 X2 =-1D0 X3 = 1354D-3 W = 123D0 XBET3 = 1D0/6D0 ELSEIF (N.EQ.4) THEN C Interaction delta for neutrons T0 =-193D0 T1 = 0D0 T2 = 0D0 T3 = 0D0 X0 =-1D0 X1 = 0D0 X2 = 0D0 X3 = 0D0 W = 0D0 XBET3 = 1D0 ELSEIF (N.EQ.5) THEN C Interaction delta for protons T0 =-249D0 T1 = 0D0 T2 = 0D0 T3 = 0D0 X0 =-1D0 X1 = 0D0 X2 = 0D0 X3 = 0D0 W = 0D0 XBET3 = 1D0 ELSE WRITE (20,*) 'Fatal error: Bad type of Skyrme force' CLOSE (20) STOP 'Fatal error: Bad type of Skyrme force' ENDIF C B1 = T0*(1D0+X0/2D0)/2D0 B2 =-T0*(X0+5D-1)/2D0 B3 = (T1*(1D0+X1/2D0)+T2*(1D0+X2/2D0))/4D0 B4 =-(T1*(5D-1+X1)-T2*(5D-1+X2))/4D0 B5 =-(3D0*T1*(1D0+X1/2D0)-T2*(1D0+X2/2D0))/16D0 B6 = (3D0*T1*(5D-1+X1)+T2*(5D-1+X2))/16D0 B7 = T3*(1D0+X3/2D0)/12D0 B8 =-T3*(5D-1+X3)/12D0 B9 =-W/2D0 C END C======================================================================= C= Construction of the truncated basis C======================================================================= SUBROUTINE BASE INCLUDE 'All.f' C IF (N0.GT.20) THEN WRITE(20,*)'N0 = ',N0,' is too grand. Limit at 20' CLOSE (20) STOP 'Fatal error: N0 too grand ' ENDIF C OPEN(UNIT=30,FILE='base') WRITE(30,*) N0,Q DSOQ = DSQRT(Q) OSDQ = 1D0/DSOQ OSQ = 1D0/Q QF(1) = DSOQ-OSDQ QF(2) = DSOQ+OSDQ QF(3) = Q + OSQ - 2D0 QF(4) = Q + OSQ + 2D0 QF(5) = Q - OSQ C NLMAX = IDINT(DBLE(N0+2)/Q**(1D0/3D0)-1D0-5D-1/Q) C NRMAX = NLMAX/2 NZMAX = IDINT(DBLE(N0+2)*Q**(2D0/3D0)-5D-1-Q) C NZMT: Nz Max True, NRMT: Nr, NMMT: Lambda NZMT = 0 NRMT = 0 NLMT = 0 C EC : ECUT EC = DBLE(N0+2) Q1S3= Q**(1D0/3D0) Q2S3= Q**(2D0/3D0) NOV = 0 NOB = 0 MBD = 0 DO IO = 1, 2*NLMAX + 1, 2 C LIL: Low IL LIL = (IO-1)/2 DO IP = 1, -1, -2 C NVBT : Number of Vectors in the Block temp. NVBT = 0 DO IR = 0, NLMAX DO IZ = 0, NZMAX DO IL = LIL, LIL+1 C Condition on K IS = IO - 2*IL C Condition on Parity IF (1-2*MOD(IZ+IL,2).EQ.IP) THEN C Condition on Energy E = Q1S3*DBLE(IR+IR+IL+1) + (DBLE(IZ)+5D-1)/Q2S3 IF (E.LE.EC) THEN NVBT = NVBT + 1 NOV = NOV + 1 MZ(NOV) = IZ MR(NOV) = IR ML(NOV) = IL MS(NOV) = IS MB(NOV) = NOB + 1 IF (NZMT.LT.IZ) NZMT = IZ IF (NRMT.LT.IR) NRMT = IR IF (NLMT.LT.IL) NLMT = IL ENDIF ENDIF ENDDO ENDDO ENDDO IF (NVBT.GT.0) THEN NOB = NOB + 1 NHB(NOB) = NOV - NVBT + 1 NVB(NOB) = NVBT IF (MBD.LT.NVBT) MBD = NVBT K2B(NOB) = IO IPB(NOB) = IP WRITE(30,*) NOB,NVBT,IO,PTY(IP) WRITE(30,*) WRITE(30,*) & (I-NHB(NOB)+1,I,MZ(I),MR(I),ML(I),MS(I), & 2*MR(I)+ML(I)+MZ(I),MZ(I),ML(I),MS(I),I=NHB(NOB),NOV) ENDIF ENDDO ENDDO NZMAX = NZMT NRMAX = NRMT NLMAX = NLMT WRITE(30,*) NOV,NOB,MBD,NZMAX,NRMAX,NLMAX CLOSE (30) C 100 FORMAT(29HStates of the truncated base , & 38H(axially deformed harmonic-oscillator),//, &I2,13H major shells,10X,25HDeformation parameter q =,F8.4) 101 FORMAT(/,79(1H=),/,8H|| Block,I4,2H :,I4,8H vectors, & I26,2H/2,A1,22X,2H||,/,79(1H=)) 102 FORMAT(11X,5HIndex,11X,20H|Nz Nr Lambda Sigma>, &6X,19H|N Nz Lambda Sigma>,/,79(1H-)) 103 FORMAT(I12,2H :,I3,10X,1H|,I3,I3,I3,I3,5H/2 >,8X, & 1H|,I3,I3,I3,I3,5H/2 >) 104 FORMAT(//,79(1H*),/,1H*,I4,11H vectors in,I3,8H blocks.,51X, &1H*,/,1H*,27H Maximum block dimension : ,I3,47X,1H*,/,1H*, &10H Max(Nz) =,I3,5X,9HMax(Nr) =,I3,5X,13HMax(Lambda) =,I3, &26X,1H*,/,79(1H*)) C END C============================================================================== C Find all the roots of Hermite polynomial of degree N C============================================================================== SUBROUTINE ROOTHE(X,N) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(*) PARAMETER (EPS=3D-14,PIM4=7511255444649425D-16,MAXIT=100) C IF (MOD(N,2).EQ.0) THEN M = N/2 ELSE M = (N-1)/2 X(M+1) = 0D0 ENDIF DO I = 1, M IF (I.EQ.1) THEN Z = DSQRT(DBLE(2*N+1)) - 185575D-5*(2*N+1)**(-1D0/6D0) ELSEIF (I.EQ.2) THEN Z = Z-114D-2*N**426-3/Z ELSE IF (I.EQ.3) THEN Z = 186D-2*Z - 86D-2*X(1) ELSE IF (I.EQ.4) THEN Z = 191D-2*Z - 91D-2*X(2) ELSE Z = 2D0*Z - X(I-2) ENDIF DO ITS = 1, MAXIT P1 = PIM4 P2 = 0D0 DO J = 1, N P3 = P2 P2 = P1 P1 = Z*DSQRT(2D0/J)*P2 - DSQRT(DBLE(J-1)/DBLE(J))*P3 ENDDO PP = DSQRT(2D0*N)*P2 Z1 = Z Z = Z1-P1/PP IF (DABS(Z-Z1).LE.EPS) GOTO 1 ENDDO PAUSE 'Too many iterations in gauher' 1 X(I) = Z ENDDO C END C============================================================================== C Find all the roots of Laguerre generalize polynomial of degree N (ALF) C============================================================================== SUBROUTINE ROOTLA(X,N,ALF) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(*) PARAMETER (EPS=3D-14, MAXIT=100) C DO I = 1, N IF (I.EQ.1) THEN Z = (1D0+ALF)*(3D0+92D-2*ALF)/(1D0+24D-1*N+18D-1*ALF) ELSEIF (I.EQ.2) THEN Z = Z + (15D0+625D-2*ALF)/(1D0+9D-1*ALF+25D-1*N) ELSE AI = I - 2 Z = Z + ((1D0+255D-2*AI)/(19D-1*AI) + & 126D-2*AI*ALF/(1D0+35D-1*AI))*(Z-X(I-2))/(1D0+3D-1*ALF) ENDIF DO ITS = 1, MAXIT P1 = 1D0 P2 = 0D0 DO J = 1, N P3 = P2 P2 = P1 P1 = ((DBLE(2*J-1)+ALF-Z)*P2 - (DBLE(J-1)+ALF)*P3)/J ENDDO PP = (N*P1-(DBLE(N)+ALF)*P2)/Z Z1 = Z Z = Z1 - P1/PP IF (DABS(Z-Z1).LE.EPS) GOTO 1 ENDDO PAUSE 'Too many iterations in GauLag' 1 X(I) = Z ENDDO C END C======================================================================= C= Enter the values of Mesh Points C======================================================================= SUBROUTINE FILLMP INCLUDE 'All.f' DIMENSION X(100) C XH(0) = 0D0 XL(0) = 0D0 CALL ROOTHE(X,NGH) DO I = 1, NGHP XH(I) = X(NGHP-I+1) ENDDO CALL ROOTLA(X,NGL,0D0) DO I = 1, NGL XL(I) = X(I) ENDDO C END C======================================================================= C Calculate Hermite and Laguerre polynomials using C the recurrence relation C======================================================================= DOUBLE PRECISION FUNCTION RECUR(X,FN,FNM1) INCLUDE 'All.f' C RECUR = ((A2N+A3N*X)*FN-A4N*FNM1)/A1N C END C======================================================================= C Calculate polynomials Hermite or Laguerre at the mesh point C XH(K) or XL(K) C IHL = 1 : Hermite C IHL = 2 : Laguerre C======================================================================= SUBROUTINE CALPLS(K,IHL) INCLUDE 'All.f' C C Hermite polynomials C IF (IHL.EQ.1) THEN C X = XH(K) QH(K,0) = 1D0 QH(K,1) = 2D0*X A1N = 1D0 A2N = 0D0 A3N = 2D0 QH1(K,0) =-X QH1(K,1) = 2D0 - 2D0*X*X DO I = 1, NZMAX-1 A4N = DBLE(2*I) QH(K,I+1) = RECUR(X,QH(K,I),QH(K,I-1)) QH1(K,I+1) = DBLE(I+I+2)*QH(K,I)-X*QH(K,I+1) ENDDO IF (K.NE.0) THEN IF (NZMAX.LT.NGH-1) THEN DO I = NZMAX, NGH-2 A4N = DBLE(2*I) QH(K,I+1) = RECUR(X,QH(K,I),QH(K,I-1)) ENDDO ENDIF ENDIF C ELSEIF (IHL.EQ.2) THEN C C Laguerre polynomials C X = XL(K) A3N =-1 DO L = 0, NLMAX QL(K,0,L) = 1D0 QL(K,1,L) = DBLE(1+L) - X QL1(K,0,L) = DBLE(L) - X QL1(K,1,L) = X*X - X*DBLE(L+L+3) + DBLE(L*L+L) DO I = 1, NRMAX - 1 A1N = DBLE(I + 1) A2N = DBLE(I + I + L + 1) A4N = DBLE(I + L) QL(K,I+1,L) = RECUR(X,QL(K,I,L),QL(K,I-1,L)) QL1(K,I+1,L) = (DBLE(I+I+2+L)-X)*QL(K,I+1,L) & - DBLE(I+I+L+L+2)*QL(K,I,L) ENDDO ENDDO IF (K.NE.0) THEN IF (NRMAX.LT.NGL) THEN DO I = NRMAX, NGL - 1 L = 0 A1N = DBLE(I + 1) A2N = DBLE(I + I + L + 1) A4N = DBLE(I + L) QL(K,I+1,L) = RECUR(X,QL(K,I,L),QL(K,I-1,L)) ENDDO ENDIF ENDIF C ENDIF C END C======================================================================= C= Calculate values at the mesh points of the polynoms used in C= integrating. C======================================================================= SUBROUTINE POLYHL INCLUDE 'All.f' C IF (NGH.GT.50.OR.NGL.GT.16) THEN WRITE(20,*) 'Fatal error: Too many points for Gauss methode' CLOSE(20) STOP 'Fatal error: Too many points for Gauss methode' ENDIF C---------------- CALL FILLMP C---------------- DO I = 1, NGHP CALL CALPLS(I,1) ENDDO DO I = 1, NGL CALL CALPLS(I,2) ENDDO C-------------------------- C Calculate weight function C-------------------------- DO I = 1, NGHP WH(I) = DSQRT(PI)*PO2(NGH-1)*FACT(NGH-1)/ & (DBLE(NGH)*QH(I,NGH-1)**2) GH(I) = DEXP(XH(I)*XH(I))*WH(I) ENDDO DO I = 1, NGL WL(I) =-XL(I)/(DBLE(NGL*NGL)*QL(I,NGL-1,0)* & (QL(I,NGL,0)-QL(I,NGL-1,0))) GL(I) = DEXP(XL(I))*WL(I) ENDDO C-------------- C Normalisation C-------------- DO N = 0, NZMAX X = DSQRT(DSQRT(PI)*PO2(N)*FACT(N)) DO I = 1, NGHP QH(I,N) = QH(I,N)/X QH1(I,N) = QH1(I,N)/X ENDDO ENDDO C DO N = 0, NRMAX X = 1D0 DO L = 1, NLMAX X = X * DSQRT(DBLE(N+L)) DO I = 1, NGL QL(I,N,L) = QL(I,N,L)/X QL1(I,N,L) = QL1(I,N,L)/X ENDDO ENDDO ENDDO C END C======================================================================= C= Show out the data being used and do some basic tasks C======================================================================= SUBROUTINE SHOWDA INCLUDE 'All.f' C C---------------------------------------- C Write out the constrains of the program C---------------------------------------- WRITE(20,100) NV,NB,NH,NL,NZ,NR NM = NUN(1) + NUN(2) C--------------- CALL INFON C--------------- PRINT 101,CSYM,NM,NAME,NUN,NM WRITE(20,101) CSYM,NM,NAME,NUN,NM C---------------------- CALL SKYRME(NSKY) C---------------------- PRINT 102,SKM(NSKY) WRITE(20,102) SKM(NSKY) C-------------- CALL BASE C-------------- PRINT 103,N0,B0,Q WRITE(20,103) N0,B0,Q NGHP = NGH/2 C---------------- CALL POLYHL C---------------- PRINT 104,NGH,NGL WRITE(20,104) NGH,NGL C IF (MNHF.GT.1000.OR.MNHTDA.GT.1000) THEN WRITE(20,*) &'Fatal error: Maximum number of iterations: HF <=1000,HTDA<=1000' CLOSE(20) STOP &'Fatal error: Maximum number of iterations: HF <=1000,HTDA<=1000' ENDIF IF (MNHF.LE.0.AND.MNHTDA.LE.0) THEN WRITE(20,*) 'Fatal error: No iteration. Program stopped' CLOSE(20) STOP 'Fatal error: No iteration. Program stopped' ENDIF C PRINT 105,MNHF,MNHTDA,MNBCS WRITE(20,105) MNHF,MNHTDA,MNBCS C PRINT 106,C2,Q2,C4,Q4 WRITE(20,106) C2,Q2,C4,Q4 C Change Q2 in fm^2 Q4 in fm^4 Q2 = Q2*1D2 Q4 = Q4*1D4 C IF (MNHF.GT.0.AND.MNBCS.GT.0) THEN PRINT 107, ECBCS,G WRITE(20,107) ECBCS,G ENDIF C IF (MNHTDA.GT.0) THEN PRINT *, &FHTDA,FHTDA1,FHTDA0,DELE,ECHTDA,WFX,WFMU, &EIDEN,EPHTDA,NEIGEN,DELED WRITE(20,*) &FHTDA,FHTDA1,FHTDA0,DELE,ECHTDA,WFX,WFMU, &EIDEN,EPHTDA,NEIGEN,DELED ENDIF C PRINT 109,HFM,HTDAM WRITE(20,109) HFM,HTDAM C IF (KOIP.LT.2.AND.KOIP.GE.0) THEN PRINT 110,OPTS(KOIP,1) WRITE(20,110) OPTS(KOIP,1) ELSE WRITE (20,*) 'Fatal error: Bad option for initial potential' CLOSE(20) STOP 'Fatal error: Bad option for initial potential' ENDIF C Possibly use O1HTDA only when MNHTDA > 0 IF (MNHTDA.GT.0) THEN IF (KOPAIR.EQ.1.OR.KOPAIR.EQ.0) THEN PRINT 118,OPTS(KOPAIR,2) WRITE(20,118) OPTS(KOPAIR,2) ELSE STOP 'Bad option for one pair transfers' ENDIF ENDIF C************** IF (MNHTDA.GT.0) THEN PRINT 114,OPTS(KORINT,9) WRITE(20,114) OPTS(KORINT,9) ENDIF C PRINT 116,EPSE,EPQ20,EPQ40 WRITE(20,116) ,EPSE,EPQ20,EPQ40 C IF (MNHTDA.GT.0) THEN PRINT 117, NUK2(1), PTY(NUIP(1)), NUK2(2), PTY(NUIP(2)) WRITE(20,117) NUK2(1), PTY(NUIP(1)), NUK2(2), PTY(NUIP(2)) ENDIF C IF (ICTRL.EQ.0) THEN PRINT 118,'Don''t do simplex' ELSE IF (MNHTDA.EQ.0) THEN PRINT 118,'Do simplex' ELSE WRITE(20,*) & 'Please check, this program does not do simplex for HTDA' CLOSE(20) STOP 'Please check, this program does not do simplex for HTDA' ENDIF ENDIF C 100 FORMAT(52HThis program is configured to calculated in a basis , &2Hof,I5,19H vectors at maximum,/,27Hand they are distributed in, &8H at most,I4,34H blocks. It can do integrals up to,I3,/, &24Hpoints by G-H method and,I3,22H points by G-L method., &29H Limits for degree of Hermite,/, &28Hand Laguerre polynomials are, &I3,4H and,I3,14H respectively.,/,79(1H=),/) 101 FORMAT(A,1H-,I3,2H :,5X,A,8X,4HN = ,I3,5X,4HZ = ,I3,10X,4HA = ,I3) 102 FORMAT(20HSkyrme force type : ,20X,A4) 103 FORMAT(24HNumber of major shells :,I18,/, &26HOscillator parameter : b =,F8.4,3H/fm,3X, &27HDeformation parameter : q =,F8.4) 104 FORMAT(22HQuadrature formulas : , &I2,21H points Gauss-Hermite,4X,I2,22H points Gauss-Laguerre) 105 FORMAT(30HMaximum number of iterations :,/, &5X,3HHF:,I4,13X,5HHTDA:,I4,11X,4HBCS:,I5) 106 FORMAT(24HDeformation contraints :,/, &18H,Quadruple moment:,12X,4HC2 =,F12.7,5X,4HQ2 =,F10.5,6H(barn),/, &21H,Hexadecapole moment:,9X,4HC4 =,F12.9,5X,4HQ4 =,F10.5, & 6H(barn)) 107 FORMAT(16HBCS parameters :,/, &16H*Cut-off Energy:,F19.2,5H(MeV),/, &17H*Seniority force:,F18.2,10X,F10.2) c 108 FORMAT(17H,HTDA parameters :,/,13H,Delta force:, c &F22.2,10X,F10.2,F22.2,10X,F10.2/, c &19H,Excitation window:,F16.2,5H(MeV),5X,/, c &30H,Cut-off in excitation energy:,F25.2,5H(MeV),/, c &17H,Weight function:,3X,3HX =,F10.2,10X,4HMu =,F9.2,/, c &31H,Diagonalisation: EIDEN = ,F6.2,12H EPHTDA =,D10.1,/, c &46H,Number of Eigenvalue : neutrons and protons = ,I2,/, c &25H,Fenêtre réduite dans RPA:,F8.2,3HMeV) 109 FORMAT(43HMixing parameter (HF,HTDA) X*NEW+(1-X)*OLD:,F8.2,F18.2) 110 FORMAT(24H||Initial potential is: ,A) 116 FORMAT(21HConverged condition :,/,8H*Energy:,F13.6,3HMeV,2X, &4H*Q20:,F10.2,3Hfm2,2X,4H*Q40:,F10.2,3Hfm4) 117 FORMAT(32HK/pi configuration of neutrons :,I8,2H/2,A1,/, &32HK/pi configuration of protons :,I8,2H/2,A1) 118 FORMAT(A) 114 FORMAT(24H||Residual interaction: ,A9) C END C======================================================================= C= Read input data C======================================================================= SUBROUTINE REDATA INCLUDE 'All.f' C C Default values, ones that will be used when there's no other source C of data CALL INITDA C Commons values for some calculations:read from "data" file C These ones mask initialized values in the program CALL DATAIN C Show value of all variables CAll SHOWDA C END C======================================================================= C Initialial potential C======================================================================= SUBROUTINE INIPOT INCLUDE 'All.f' C X = Q**(1D0/6D0) BP = X*B0 BP2 = BP*BP BZ = BP/DSQRT(Q) BZ2 = BZ*BZ BPZ = BZ*BP BZP2= BZ*BP2 BZP2PI = BZP2/PI AA = DBLE(NUN(1)+NUN(2)) HB0 = 20.735530D0*(AA-1D0)/AA IF (KOIP.EQ.0) THEN C C Woods-Saxon potential construction C GR: Grand R C GR = 135D-2*AA**(1D0/3D0) !!!! AL = 6D-1 !!!! V0 = 70D0 !!!! DO I = 1, NGHP DO J = 1, NGL R = DSQRT(Q*XL(J)+XH(I)*XH(I))/(X*BP) V =-V0/(1D0+DEXP((R-GR)/AL)) DO K = 1, 2 EM(I,J,K) = HB0 CE(I,J,K) = V SO(I,J,K) =-35D-2*V ENDDO ENDDO ENDDO ELSEIF (KOIP.EQ.1) THEN C Initial HF potential is read from disk C M : NGHP Read, N : NGL Read OPEN(UNIT=21,FORM='UNFORMATTED',FILE='HFfields.in') READ(21) M, N IF (M.NE.NGHP.OR.N.NE.NGL) THEN PRINT *,'-------------------------' PRINT *,'Read data NGH = ',2*M,' NGL = ',N PRINT *,'Desired data NGH = ',NGH,' NGL = ',NGL WRITE(20,*) 'Bad data file for HF fields' CLOSE(20) STOP 'Bad data file for HF fields' ENDIF READ(21) (((EM(I,J,K),I = 1, M),J = 1, N),K = 1, 2) READ(21) (((CE(I,J,K),I = 1, M),J = 1, N),K = 1, 2) READ(21) (((SO(I,J,K),I = 1, M),J = 1, N),K = 1, 2) CLOSE(21) ELSE WRITE(20,*) 'Bad choice of initial potential option' CLOSE(20) STOP 'Bad choice of initial potential option' ENDIF C END C======================================================================= C Calculation of HF Matrix Elements C======================================================================= SUBROUTINE CAHFME(IBLOCK) INCLUDE 'All.f' C C NHBM1 N Head Block Minus 1 N = NVB(IBLOCK) NHBM1 = NHB(IBLOCK)-1 DO I = 1, N IA = NHBM1 + I NZA = MZ(IA) NRA = MR(IA) LA = ML(IA) NSA = MS(IA) KA2 = LA + LA + NSA !2*Omega_z XLA = DBLE(LA) XX = BP2*DBLE(LA*NSA) DO J = 1, I IB = NHBM1 + J LB = ML(IB) NSB = MS(IB) KB2 = LB + LB + NSB X = 0D0 NZB = MZ(IB) NRB = MR(IB) XLB = DBLE(LB) XLAB= DBLE(LA*LB) IF (NSA.NE.NSB) THEN LC = (LA + LB - 1)/2 DO IL = 1, NGL QL00 = QL(IL,NRA,LA)*QL(IL,NRB,LB) QL01 = QL(IL,NRA,LA)*QL1(IL,NRB,LB) QL10 = QL1(IL,NRA,LA)*QL(IL,NRB,LB) X1 = XLB*(QL00+QL01) - XLA*QL01 X2 = XLA*(QL00+QL10) - XLB*QL10 ETA = XL(IL) ETALC= 1D0 IF (LC.NE.0) ETALC = ETA**LC ETAWL = ETALC*WL(IL) DO IH = 1, NGHP QH10 = QH1(IH,NZA)*QH(IH,NZB) QH01 = QH(IH,NZA)*QH1(IH,NZB) X3 = QH10*X1 + QH01*X2 X = X + ETAWL*X3*SO(IH,IL,INP)*WH(IH) ENDDO ENDDO X = (X + X)*BPZ ELSE C XSA2 = 2*Spin_A XSA2 = DBLE(NSA) DO IL = 1, NGL ETA = XL(IL) ETALA= 1D0 IF (LA.NE.0) ETALA = ETA**LA ETAWL=ETALA/ETA*WL(IL) QL00 = QL(IL,NRA,LA)*QL(IL,NRB,LB) QL01 = QL(IL,NRA,LA)*QL1(IL,NRB,LB) QL10 = QL1(IL,NRA,LA)*QL(IL,NRB,LB) QL11 = QL1(IL,NRA,LA)*QL1(IL,NRB,LB) QL00E= QL00*ETA X12 = BP2*(QL11 + XLAB*QL00) X3SUP= XX*(QL10+QL01) DO IH = 1, NGHP QH00 = QH(IH,NZA)*QH(IH,NZB) QH11 = QH1(IH,NZA)*QH1(IH,NZB) X11 = BZ2*QL00E X1 = (QH11*X11 + QH00*X12)*EM(IH,IL,INP) X2 = QH00*QL00E*CE(IH,IL,INP) X3 =-X3SUP*QH00*SO(IH,IL,INP) X = X + (X1+X2+X3)*ETAWL*WH(IH) ENDDO ENDDO X = X + X ENDIF A(I,J) = X A(J,I) = X ENDDO ENDDO C END C======================================================================= C Calculation of HF Matrix Elements C = C Sum_alpha,beta C======================================================================= SUBROUTINE CALEHF(IBLOCK) INCLUDE 'All.f' C N = NVB(IBLOCK) NHBM1 = NHB(IBLOCK) - 1 DO I = 1, N J = NHBM1 + I TMP = 0D0 DO IA = 1, N DO IB = 1, N TMP = TMP+COEF(IA,I,IBLOCK,INP)*A(IA,IB)*COEF(IB,I,IBLOCK,INP) ENDDO ENDDO EHF(J) = TMP NBRE(J) = J ENDDO C END C======================================================================= C Diagonalization les blocks C======================================================================= SUBROUTINE DIAGON(N,IBLOCK) INCLUDE 'All.f' DIMENSION D(ND),E(ND),Z(ND,ND) DATA EPS/9D-20/ DATA TOL/1D-36/ C Z(1,1) = 1D0 IF (N.LE.1) GOTO 10 DO I = 1, N DO J = 1, N Z(I,J) = A(I,J) ENDDO ENDDO DO II = 2, N L = N - II I = L + 2 F = Z(I,I-1) GG = 0D0 IF (L.NE.0) THEN DO K = 1, L GG = GG + Z(I,K)**2 ENDDO ENDIF H = GG + F*F IF (GG.LE.TOL) THEN E(I) = F H = 0D0 ELSE L = L + 1 GG =-DSIGN(DSQRT(H),F) E(I ) = GG H = H - F*GG Z(I,I-1) = F - GG F = 0D0 DO J = 1, L Z(J,I) = Z(I,J)/H GG = 0D0 DO K = 1, J GG = GG + Z(J,K)*Z(I,K) ENDDO J1 = J + 1 IF (J1.LE.L) THEN DO K = J1, L GG = GG + Z(K,J)*Z(I,K) ENDDO ENDIF E(J) = GG/H F = F + GG*Z(J,I) ENDDO HH = F/(H+H) DO J = 1, L F = Z(I,J) GG = E(J) - HH*F E(J) = GG DO K = 1, J Z(J,K) = Z(J,K) - F*E(K) - GG*Z(I,K) ENDDO ENDDO ENDIF D(I) = H ENDDO D(1) = 0D0 E(1) = 0D0 DO I = 1, N L = I - 1 IF (DABS(D(I)).GE.1D-36) THEN IF (L.NE.0) THEN DO J = 1, L GG = 0D0 DO K = 1, L GG = GG + Z(I,K)*Z(K,J) ENDDO DO K = 1, L Z(K,J) = Z(K,J) - GG*Z(K,I) ENDDO ENDDO ENDIF ENDIF D(I) = Z(I,I) Z(I,I) = 1D0 IF (L.NE.0) THEN DO J = 1, L Z(I,J)=0D0 Z(J,I)=0D0 ENDDO ENDIF ENDDO DO I = 2, N E(I-1) = E(I) ENDDO E(N) = 0D0 B = 0D0 F = 0D0 DO L = 1, N J = 0 H = EPS * (DABS(D(L))+DABS(E(L))) IF (B.LT.H) B = H DO M = L, N IF (DABS(E(M)).LE.B) GOTO 1 ENDDO 1 IF (M.NE.L) THEN 2 IF (J.EQ.30) THEN WRITE(20,*) 'Problem in Diagonalization a block number ',IBLOCK CLOSE(20) STOP 'Problem in Diagonalization a bloc ' ENDIF J = J + 1 P = (D(L+1)-D(L))/(2D0*E(L)) R = DSQRT(P*P + 1D0) H = D(L) - E(L)/(P+DSIGN(R,P)) DO I = L, N D(I) = D(I) - H ENDDO F = F + H P = D(M) C = 1D0 S = 0D0 M1 = M - 1 MLL = M1 + L DO II = L, M1 I = MLL - II GG = C*E(I) H = C*P IF (DABS(P).GE.DABS(E(I))) THEN C = E(I)/P R = DSQRT(C*C+1D0) E(I+1) = S*P*R S = C/R C = 1D0/R ELSE C = P/E(I) R = DSQRT(C*C+1D0) E(I+1) = S*E(I)*R S = 1D0/R C = C/R ENDIF P = C*D(I) - S*GG D(I+1) = H + S*(C*GG+S*D(I)) DO K = 1, N H = Z(K,I+1) Z(K,I+1) = S*Z(K,I) + C*H Z(K,I) = C*Z(K,I) - S*H ENDDO ENDDO E(L) = S*P D(L) = C*P IF (DABS(E(L)).GT.B) GOTO 2 ENDIF D(L) = D(L) + F ENDDO N1 = N - 1 DO I = 1, N1 K = I P = D(I) II = I + 1 DO J = II, N IF (D(J).LT.P) THEN K = J P = D(J) ENDIF ENDDO IF (K.NE.I) THEN D(K) = D(I) D(I) = P DO J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P ENDDO ENDIF ENDDO DO I = 1, N A(I,I) = D(I) ENDDO 10 DO I = 1, N DO J = 1, N COEF(I,J,IBLOCK,INP) = Z(I,J) ENDDO ENDDO C END C======================================================================= C BCS to deal the pairing C======================================================================= SUBROUTINE BCS INCLUDE 'All.f' DIMENSION DELTA(2),AMBDA(2),G1(2) C IF (INP.EQ.1) THEN OPEN(UNIT=30,FILE='BCSndis') OPEN(UNIT=31,FILE='BCSnKPi') ELSEIF (INP.EQ.2) THEN OPEN(UNIT=30,FILE='BCSpdis') OPEN(UNIT=31,FILE='BCSpKPi') ENDIF NM = NUN(1) + NUN(2) WRITE(20,99) C PRINT 99 XYX = 0D0 XX = 0D0 NF = NUN(INP)/2 AMBDA(INP) = (EHF(NF)+EHF(NF+1))/2D0 DELTA(INP) = 1D0 CUT = AMBDA(INP) + ECBCS c IF (CUT.LT.0D0) CUT = 0D0 IF (G(INP).GT.1D-3) THEN AMBDAI = AMBDA(INP) DELTAI = DELTA(INP) AA = DBLE(2*NF) G1(INP)= G(INP)/(11D0 + AA) PP1 = 1D0/DBLE(ITHF+1) C IF (ITHF.GT.9) PP1 = 0D0 PP = 1D0 - PP1 PRI = 1D-8 NIT = 0 C Calcul de Lambda et Delta 1 SOER = 0D0 SOR = 0D0 SOMM = 0D0 DO I = 1, NOV IF (EHF(I).LE.CUT) THEN RA = DSQRT((EHF(I)-AMBDAI)**2 + DELTAI**2) SOMM = SOMM + G1(INP)*DELTAI/RA SOER = SOER + 1D0-EHF(I)/RA SOR = SOR + 1D0/RA ENDIF ENDDO DELTA(INP) = 5D-1*SOMM AMBDA(INP) = (AA-SOER)/SOR C C Test de convergence C IF (PRI.LT.DABS(AMBDAI-AMBDA(INP))) THEN IF (MNBCS.GT.NIT) THEN AMBDAI = AMBDA(INP) DELTAI = DELTA(INP)*PP + DELTAI*PP1 NIT = NIT + 1 GOTO 1 ENDIF ENDIF IF (PRI.LT.DABS(DELTAI-DELTA(INP))) THEN IF (MNBCS.GT.NIT) THEN DELTAI = DELTA(INP)*PP + DELTAI*PP1 NIT = NIT + 1 GOTO 1 ENDIF ENDIF C PRINT 100, NIT,MNBCS,AMBDAI,AMBDA(INP),DELTAI,DELTA(INP) WRITE(20,100) NIT,MNBCS,AMBDAI,AMBDA(INP),DELTAI,DELTA(INP) C Calcul niveau occupation probabilities WRITE(30,*) CSYM,NM,NAME,NUN,NM,AMBDA(INP) WRITE(30,*) WRITE(31,*) CSYM,NM,NAME,NUN,NM,AMBDA(INP) WRITE(31,*) SUIVI = 0D0 DO I = 1, NOV K = NBRE(I) B = EHF(I) - AMBDA(INP) C = DSQRT(B*B + DELTA(INP)**2) IF (EHF(I).GT.CUT) C = B DENS(K,K) = 5D-1*(1D0-B/C) RHOII = DENS(K,K) VI = DSQRT(RHOII) UI = DSQRT(1-RHOII) SUIVI = SUIVI + UI*VI IF (EHF(I).LE.CUT) THEN WRITE(30,103) I,EHF(I),RHOII,VI,UI,UI*VI C Find the K and Parity of the level IB = MB(K) WRITE(31,106) I,EHF(I),RHOII,K2B(IB),PTY(IPB(IB)) ENDIF ENDDO XX =-DELTA(INP)**2*SOR ENDIF C DO I = 1, NOV K = NBRE(I) XYX = XYX + 2D0*DENS(K,K)*EHF(I) ENDDO EPAIR = EPAIR + 0.5d0*XX WRITE(20,101) XX,SUIVI C PRINT 101, XX,SUIVI ESOMIN = ESOMIN + XYX CC PRINT *,'ESOMIN = ',ESOMIN,'XYX = ',XYX CLOSE(30) CLOSE(31) C 99 FORMAT(2X,33H+Pairing correction by BCS method) 100 FORMAT(5X,6HAfter ,I5,1H/,I5,12H iterations:,/,9X, &7HLambda:,F30.10,3H =>,F30.10,/,9X,7HDelta :,F30.10,3H =>,F30.10) 101 FORMAT(7X,9HPairing :,F12.4,6H (MeV),13X,12HDiffusivity:,F12.4) 102 FORMAT(A,1H-,I3,2H :,5X,A,8X,4HN = ,I3,5X,4HZ = ,I3,10X, &4HA = ,I3,/,8HLambda =,F12.4) 103 FORMAT(I3,5(F15.4)) 104 FORMAT(3HInx,9X,6HEnergy,9X,6Hrho_ij,13X,2HVi,13X,2HUi,9X,5HUi*Vi) 105 FORMAT(5HIndex,9X,6HEnergy,9X,6Hrho_ij,9X,4HK/Pi) 106 FORMAT(I5,2(F15.4),I10,2H/2,A) C END c*************************************************************************** ****** ****** Calcul des elements de matrice de Q_2^mu dans la base OH ****** (methode analytique d'apres la these de P. Quentin) ****** DOUBLE PRECISION FUNCTION QOH(I1,I2,IRV) INCLUDE 'All.f' dsqflt(i)=dsqrt(dfloat(i)) c nsa=ms(i1) nsb=irv*ms(i2) QOH=0d0 if (nsa.ne.nsb) return q2mu=0d0 c nza=mz(i1) nra=mr(i1) la=ml(i1) lla=la c nzb=mz(i2) nrb=mr(i2) lb=ml(i2) llb=irv*lb c nzpz=nza-nzb mu=lla-llb naa=(2*nra+la+lla)/2 nba=naa-lla napa=naa-(2*nrb+lb+llb)/2 nbpb=napa-mu ndiff=iabs(napa)+iabs(nbpb)+iabs(nzpz) if (ndiff.gt.2) return c iz=nza+(iabs(nzpz)-nzpz)/2 ia=naa+(iabs(napa)-napa)/2 ib=nba+(iabs(nbpb)-nbpb)/2 if (ndiff.eq.0) then ! mu = 0 q2mu=dfloat(1+2*iz)-dfloat(1+ia+ib)/q elseif (iabs(nzpz).eq.2) then q2mu=dsqflt(iz*(iz-1)) elseif (napa*nbpb.eq.1) then q2mu=dsqflt(ia*ib)/q elseif (iabs(nzpz*napa).eq.1) then ! mu = +/- 1 q2mu=-dfloat(mu)*dsqflt(iz*ia) elseif (iabs(nzpz*nbpb).eq.1) then q2mu=dfloat(mu)*dsqflt(iz*ib) elseif (iabs(napa).eq.2) then ! mu = +/- 2 q2mu=dsqflt(ia*(ia-1)) elseif (iabs(nbpb).eq.2) then q2mu=dsqflt(ib*(ib-1)) elseif (napa*nbpb.eq.-1) then q2mu=-dsqflt(4*ia*ib) endif if (mu.eq.0) QOH = q2mu/(2.*bz*bz) if (iabs(mu).eq.1) QOH = q2mu*dsqrt(3d0)/(2.*bz*bp) if (iabs(mu).eq.2) QOH = q2mu*dsqrt(6d0)/(4.*bp*bp) jzap1=lla+(nsa+1)/2 jzbp1=llb+(nsb+1)/2 if (mod(jzap1,2).eq.-1) QOH = -QOH if (mod(jzbp1,2).eq.-1) QOH = -QOH c RETURN END C============================================================================== C To check (houda) C Calculate the Quadruple matrix element if necessary C============================================================================== SUBROUTINE CALQHF INCLUDE 'All.f' C CQ1 = DSQRT(3D0)/(2D0*BPZ) CQ2 = DSQRT(6D0)/(4D0*BP2) CQ3 = BZ2 + BZ2 CQ4 = BP2 + BP2 AA = DBLE(NUN(1)+NUN(2)) X = 110D0/AA*4D0*PI/5D0 X = X/(2073D-2*AA**(1D0/3D0)/41D0)**2 Q00(1) = X/4225D-2 Q00(2) = X/3025D-2 c***************************** c Q00(1) = 20D-3 c Q00(2) = 30D-3 C****************************************************** DO I = -2*NHL(inp),2*NPL(inp) DO J = -2*NHL(inp),2*NPL(inp) QHF(I,J) = 0D0 ENDDO ENDDO DO I = 1, NBW DO J = 1, NBW IBL1 = IOB(I) IBL2 = IOB(J) IF (IPB(IBL1).EQ.IPB(IBL2)) THEN K21 = K2B(IBL1) K22 = K2B(IBL2) IF (IABS(K21-K22).LE.4) THEN DO I1 = 1, NVB(IBL1) DO I2 = 1, NVB(IBL2) C IOH*: Index of Oscillator harmonic state IOH1 = NHB(IBL1) + I1 - 1 IOH2 = NHB(IBL2) + I2 - 1 TMP=QOH(IOH1,IOH2,1) DO II = 1, NSIB(I) DO IJ = 1, NSIB(J) IHF1 = ISB(I,II) IHF2 = ISB(J,IJ) COFP = COEF(I1,IRBB(I,II),IBL1,INP) * & COEF(I2,IRBB(J,IJ),IBL2,INP) QHF(IHF1,IHF2) = QHF(IHF1,IHF2) + COFP * TMP ENDDO ENDDO ENDDO ENDDO ICOEF=1 IF (IABS(K21-K22).EQ.2) ICOEF=-1 DO II = 1, NSIB(I) DO IJ = 1, NSIB(J) IHF1 = ISB(I,II) IHF2 = ISB(J,IJ) IHF1R = ISIGN(IABS(IHF1)+1,IHF1) IHF2R = ISIGN(IABS(IHF2)+1,IHF2) QHF(IHF1R,IHF2R) = ICOEF*QHF(IHF1,IHF2) ENDDO ENDDO ENDIF IF (IABS(K21+K22).LE.4) THEN DO I1 = 1, NVB(IBL1) DO I2 = 1, NVB(IBL2) C IOH*: Index of Oscillator harmonic state IOH1 = NHB(IBL1) + I1 - 1 IOH2 = NHB(IBL2) + I2 - 1 TMP=QOH(IOH1,IOH2,-1) DO II = 1, NSIB(I) DO IJ = 1, NSIB(J) IHF1 = ISB(I,II) IHF2 = ISB(J,IJ) IHF2R= ISIGN(IABS(IHF2)+1,IHF2) COFP = COEF(I1,IRBB(I,II),IBL1,INP) * & COEF(I2,IRBB(J,IJ),IBL2,INP) QHF(IHF1,IHF2R) = QHF(IHF1,IHF2R) + COFP * TMP ENDDO ENDDO ENDDO ENDDO ICOEF=-1 IF (IABS(K21+K22).EQ.2) ICOEF=1 DO II = 1, NSIB(I) DO IJ = 1, NSIB(J) IHF1 = ISB(I,II) IHF2 = ISB(J,IJ) IHF1R = ISIGN(IABS(IHF1)+1,IHF1) IHF2R = ISIGN(IABS(IHF2)+1,IHF2) QHF(IHF1R,IHF2) = ICOEF*QHF(IHF1,IHF2R) ENDDO ENDDO ENDIF ENDIF ENDDO ENDDO C END C============================================================================== C LE CALCUL DES ÉLÉMENTS DE LA MATRICE DE L'OPÉRATEUR OCTIPOLAIRE c :(i,j)les indices de l'oscillateur harmonique et mu=-3..3 C============================================================================== DOUBLE PRECISION FUNCTION OCOH(I,J,MU,IRV1,IRV2) INCLUDE 'All.f' C IF (IABS(IRV1).NE.1.OR.IABS(IRV2).NE.1) PRINT *,'Error' TMPC = 0D0 LS1 = MS(I)*IRV1 LS2 = MS(J)*IRV2 IF (LS1.EQ.LS2) THEN !Spin LZ1 = MZ(I) LZ2 = MZ(J) IF (IRV1.EQ.1) THEN LA1 = MR(I) + ML(I) LB1 = MR(I) ELSE LB1 = MR(I) + ML(I) LA1 = MR(I) ENDIF IF (IRV2.EQ.1) THEN LA2 = MR(J) + ML(J) LB2 = MR(J) ELSE LB2 = MR(J) + ML(J) LA2 = MR(J) ENDIF C MU=3******************************** IF (MU.EQ.3)THEN IF (LZ1.EQ.LZ2)THEN IF (LB1.EQ.LB2)THEN TMPC =-DSQRT(DBLE(LA1*(LA2+2)*(LA2+1)))*CC3 ENDIF IF(LA1.EQ.LA2)THEN TMPC = DSQRT(DBLE(LB2*(LB2-1)*(LB1-1)))*CC3 ENDIF IF (LA1.EQ.LA2+1)THEN TMPC =-3D0*DSQRT(DBLE(LB2*(LB2-1)*LA1))*CC3 ENDIF IF (LA1.EQ.LA2+2)THEN TMPC = 3D0*DSQRT(DBLE(LB2*LA1*(LA2+1)))*CC3 ENDIF ENDIF ENDIF C****MU=-3****************** IF (MU.EQ.-3)THEN IF(LZ1.EQ.LZ2)THEN IF (LA1.EQ.LA2)THEN TMPC =-DSQRT(DBLE(LB1*(LB2+1)*(LB2+2)))*CC3 ENDIF IF (LB1.EQ.LB2)THEN TMPC = DSQRT(DBLE(LA2*(LA2-1)*LA1))*CC3 ENDIF IF (LA1.EQ.LA2-1)THEN TMPC = 3D0*DSQRT(DBLE(LA2*LB1*(LB2+1)))*CC3 ENDIF ENDIF ENDIF C***MU=2***** IF (MU.EQ.2)THEN IF (LZ1.EQ.LZ2-1)THEN IF (LB1.EQ.LB2)THEN TMPC = DSQRT(DBLE(LZ2*LA1*(LA2+1)))*CC2 ENDIF IF (LA1.EQ.LA2)THEN TMPC = DSQRT(DBLE(LZ2*LB2*(LB2-1)))*CC2 ENDIF IF (LA1.EQ.LA2+1)THEN TMPC =-2D0*DSQRT(DBLE(LZ2*LB2*LA1))*CC2 ENDIF ENDIF IF (LZ1.EQ.LZ2+1)THEN IF (LB1.EQ.LB2)THEN TMPC = DSQRT(DBLE(LZ1*LA1*(LA2+1)))*CC2 ENDIF IF (LA1.EQ.LA2)THEN TMPC = DSQRT(DBLE(LZ1*LB2*(LB2-1)))*CC2 ENDIF IF (LA1.EQ.LA2+1)THEN TMPC =-2D0*DSQRT(DBLE(LZ1*LB2*LA1))*CC2 ENDIF ENDIF ENDIF C***MU=-2 IF (MU.EQ.-2)THEN IF (LZ1.EQ.LZ2-1)THEN IF (LA1.EQ.LA2)THEN TMPC = DSQRT(DBLE(LZ2*LB1*(LB2+1)))*CC2 ENDIF IF (LB1.EQ.LB2)THEN TMPC = DSQRT(DBLE(LZ2*LA2*(LA1+1)))*CC2 ENDIF IF (LB1.EQ.LB2+1)THEN TMPC =-2D0*DSQRT(DBLE(LZ2*LA2*LB1))*CC2 ENDIF ENDIF IF (LZ1.EQ.LZ2+1)THEN IF (LA1.EQ.LA2)THEN TMPC = DSQRT(DBLE(LZ1*LB1*(LB2+1)))*CC2 ENDIF IF (LB1.EQ.LB2)THEN TMPC = DSQRT(DBLE(LZ1*LA2*(LA1+1)))*CC2 ENDIF IF (LB1.EQ.LB2+1)THEN TMPC =-2D0*DSQRT(DBLE(LZ1*LA2*LB1))*CC2 ENDIF ENDIF ENDIF C****MU=1 IF (MU.EQ.1)THEN IF (LZ1.EQ.LZ2)THEN IF (LB1.EQ.LB2)THEN TMP11 = DBLE((4D0*LZ2+2)*CC11) TMP12 = DBLE((2D0*LB2+LA2+2)*CC12) TMP13 = DBLE((1+2*LA2+LB2)*CC12) TMP =-DSQRT(DBLE(LA1))*(TMP11-TMP12) ENDIF IF (LA1.EQ.LA2)THEN TMPC = DSQRT(DBLE(LB2))*(TMP11-TMP13) ENDIF IF (LA1.EQ.LA2+2)THEN TMPC =-DSQRT(DBLE((LA2+1)*LA1*LB1))*CC12 ENDIF IF (LA1.EQ.LA2-1)THEN TMPC = DSQRT(DBLE(LA2*LB2*(LB2-1)))*CC12 ENDIF ENDIF IF (LZ1.EQ.LZ2+2.AND.LB1.EQ.LB2)THEN TMPC =-2D0*SQRT(DBLE((LZ2+1)*LZ1*LA1))*CC11 ENDIF IF (LZ1.EQ.LZ2+2.AND.LA1.EQ.LA2)THEN TMPC = 2D0*SQRT(DBLE((LZ2+1)*LZ1*LB2))*CC11 ENDIF IF (LZ1.EQ.LZ2-2.AND.LB1.EQ.LB2)THEN TMPC =-2D0*DSQRT(DBLE((LZ1+1)*LZ2*LA1))*CC11 ENDIF IF (LZ1.EQ.LZ2-2.AND.LA1.EQ.LA2)THEN TMPC = 2D0*DSQRT(DBLE((LZ1+1)*LZ2*LA2))*CC11 ENDIF ENDIF C****MU=-1 IF (MU.EQ.-1)THEN IF (LZ1.EQ.LZ2)THEN IF (LA1.EQ.LA2)THEN TMP14 = DBLE((4D0*LZ2+2)*CC11) TMP15 = DBLE((2D0*LA2+LB2+2)*CC12) TMP16 = DBLE((1+LA2+2*LB2)*CC12) TMPC = DSQRT(DBLE(LB1))*(TMP15-TMP14) ENDIF IF (LB1.EQ.LB2)THEN TMPC = DSQRT(DBLE(LA2))*(TMP14-TMP16) ENDIF IF (LA1.EQ.LA2+1)THEN TMPC =-DSQRT(DBLE(LA1*LB1*(LB2+1)))*CC12 ENDIF IF (LB1.EQ.LB2-1)THEN TMPC = DSQRT(DBLE(LA2*LB2*(LA2-1)))*CC12 ENDIF ENDIF IF (LZ1.EQ.LZ2+2.AND.LA1.EQ.LA2)THEN TMPC =-2D0*SQRT(DBLE((LZ2+1)*LZ1*LB1))*CC11 ENDIF IF (LZ1.EQ.LZ2+2.AND.LB1.EQ.LB2)THEN TMPC = 2D0*SQRT(DBLE((LZ2+1)*LZ1*LA2))*CC11 ENDIF IF (LZ1.EQ.LZ2-2.AND.LA1.EQ.LA2)THEN TMPC =-2D0*DSQRT(DBLE((LZ1+1)*LZ2*LB1))*CC11 ENDIF IF (LZ1.EQ.LZ2-2.AND.LA1.EQ.LA2)THEN TMPC = 2D0*DSQRT(DBLE((LZ1+1)*LZ2*LA2))*CC11 ENDIF ENDIF C***MU=0 IF (MU.EQ.0)THEN IF (LA1.EQ.LA2.AND.LB1.EQ.LB2)THEN IF (LZ1.EQ.LZ2+1)THEN TMP17 = DBLE(LZ1*CQ01-(1+LA2+LB2)*CC02) TMPC = 3D0*SQRT(DBLE(LZ1))*TMP17 ENDIF IF (LZ1.EQ.LZ2-1)THEN TMP18 = DBLE(LZ2*CQ01-(1+LA2+LB2)*CC02) TMPC = 3D0*DSQRT(DBLE(LZ2))*TMP18 ENDIF IF (LZ1.EQ.LZ2-3)THEN TMPC= DSQRT(DBLE(LZ2*(LZ1+1)*(LZ1+2)))*CC01 ENDIF IF (LZ1.EQ.LZ2+3)THEN TMPC = DSQRT(DBLE(LZ1*(LZ2+1)*(LZ2+2)))*CC01 ENDIF ENDIF IF (LA1.EQ.LA2+1.AND.LZ1.EQ.LZ2+1)THEN TMPC = 3D0*DSQRT(DBLE(LA1*LB1*LZ1))*CC02 ENDIF IF (LA1.EQ.LA2+1.AND.LZ1.EQ.LZ2-1)THEN TMPC = 3D0*DSQRT(DBLE(LA1*LB1*LZ2))*CC02 ENDIF IF (LA1.EQ.LA2-1.AND.LZ1.EQ.LZ2+1)THEN TMPC = 3D0*DSQRT(DBLE(LA2*LB2*LZ1))*CC02 ENDIF IF (LA1.EQ.LA2-1.AND.LZ1.EQ.LZ2-1)THEN TMPC = 3D0*DSQRT(DBLE(LA2*LB2*LZ2))*CC02 ENDIF ENDIF ENDIF OCOH=TMPC C END C============================================================================== C To check C Calculate the octupole matrix element if necessary C============================================================================== SUBROUTINE CALOCHF INCLUDE 'All.f' C CC3 = DSQRT(5D0)/(4D0*BP3) CC2 = DSQRT(15D0)/(4D0*BZ*BP2) CC11 = DSQRT(3D0)/(4D0*BP*BZ2) CC12 = DSQRT(3D0)/(4D0*BP3) CC01 = DSQRT(2D0)/(4D0*BZ3) CC02 = DSQRT(2D0)/(4D0*BZ*BP2) AA = DBLE(NUN(1)+NUN(2)) X = 110D0/AA*4D0*PI/5D0 X = X/(2073D-2*AA**(1D0/3D0)/41D0)**2 OC00(1) = X/4225D-2 OC00(2) = X/3025D-2 DO I = -3*NHL(inp),3*NPL(inp) DO J = -3*NHL(inp),3*NPL(inp) OCHF(I,J) = 0D0 ENDDO ENDDO DO I = 1, NBW DO J = 1, NBW IBL1 = IOB(I) IBL2 = IOB(J) IF (IPB(IBL1).EQ.-IPB(IBL2)) THEN K21 = K2B(IBL1) K22 = K2B(IBL2) IF (IABS(K21-K22).LE.6) THEN MU = (K21-K22)/2 DO I1 = 1, NVB(IBL1) DO I2 = 1, NVB(IBL2) C IOH*: Index of Oscillator harmonic state IOH1 = NHB(IBL1) + I1 - 1 IOH2 = NHB(IBL2) + I2 - 1 TMP = OCOH(IOH1,IOH2,MU,1,1) DO II = 1, NSIB(I) DO IJ = 1, NSIB(J) IHF1 = ISB(I,II) IHF2 = ISB(J,IJ) COFP = COEF(I1,IRBB(I,II),IBL1,INP) * & COEF(I2,IRBB(J,IJ),IBL2,INP) OCHF(IHF1,IHF2) = OCHF(IHF1,IHF2) + COFP * TMPC ENDDO ENDDO ENDDO ENDDO ICOEF = 1 IF (MOD((K21+K22)/2,3).EQ.0) ICOEF = -1 DO II = 1, NSIB(I) DO IJ = 1, NSIB(J) IHF1 = ISB(I,II) IHF2 = ISB(J,IJ) IHF1R = ISIGN(IABS(IHF1)+1,IHF1) IHF2R = ISIGN(IABS(IHF2)+1,IHF2) OCHF(IHF1R,IHF2R) = OCHF(IHF1,IHF2)*ICOEF ENDDO ENDDO ENDIF IF (IABS(K21+K22).LE.6) THEN MU = (K21+K22)/2 DO I1 = 1, NVB(IBL1) DO I2 = 1, NVB(IBL2) C IOH*: Index of Oscillator harmonic state IOH1 = NHB(IBL1) + I1 - 1 IOH2 = NHB(IBL2) + I2 - 1 TMPC = OCOH(IOH1,IOH2,MU,1,-1) DO II = 1, NSIB(I) DO IJ = 1, NSIB(J) IHF1 = ISB(I,II) IHF2 = ISB(J,IJ) IHF2R= ISIGN(IABS(IHF2)+1,IHF2) COFP = COEF(I1,IRBB(I,II),IBL1,INP) * & COEF(I2,IRBB(J,IJ),IBL2,INP) OCHF(IHF1,IHF2R) = OCHF(IHF1,IHF2R) + COFP * TMPC ENDDO ENDDO ENDDO ENDDO IF (MOD((K22+1)/2,3).NE.0) THEN DO II = 1, NSIB(I) DO IJ = 1, NSIB(J) IHF1 = ISB(I,II) IHF2 = ISB(J,IJ) IHF2R= ISIGN(IABS(IHF2)+1,IHF2) OCHF(IHF1,IHF2R) =-OCHF(IHF1,IHF2R) ENDDO ENDDO ENDIF ICOEF = 1 IF (MOD((K21-K22)/2,3).NE.0) ICOEF = -1 DO II = 1, NSIB(I) DO IJ = 1, NSIB(J) IHF1 = ISB(I,II) IHF2 = ISB(J,IJ) IHF1R = ISIGN(IABS(IHF1)+1,IHF1) IHF2R = ISIGN(IABS(IHF2)+1,IHF2) OCHF(IHF1R,IHF2) = OCHF(IHF1,IHF2R)*ICOEF ENDDO ENDDO ENDIF ENDIF ENDDO ENDDO C END C======================================================================= C This subroutine compare the configuration number NTOT with ones C from number NCOM+1 to NTOT-1. If it differ from all the confs. then C the function give TRUE if not FALSE. C======================================================================= SUBROUTINE TESTCO(NEWCON) INCLUDE 'All.f' LOGICAL NEWCON C NEWCON =.TRUE. EEXCI = EPH(NTOT(inp),inp) DO I = NCOM(inp)+1, NTOT(inp)-1 C Verify if it exist the same configuration in the many-body basis C First: Compare p-h Energie IF (DABS(EEXCI-EPH(I,inp)).LT.1D-5) THEN DO IH = 1, NHPMAX if (inp.eq.1) then IF (MBSTn(IH,1,I).NE.MBSTn(IH,1,NTOT(inp))) GOTO 10 IF (MBSTn(IH,2,I).NE.MBSTn(IH,2,NTOT(inp))) GOTO 10 else IF (MBSTp(IH,1,I).NE.MBSTp(IH,1,NTOT(inp))) GOTO 10 IF (MBSTp(IH,2,I).NE.MBSTp(IH,2,NTOT(inp))) GOTO 10 end if IF (IH.EQ.NHPMAX) GOTO 20 ENDDO ENDIF 10 CONTINUE ENDDO RETURN 20 NEWCON =.FALSE. C END C======================================================================= C This subroutine receive NHPMAX numbers and arrange them by C descending order C Results is store in NAFAR C======================================================================= SUBROUTINE ARRANG(I1,I2,I3,I4,I5) INCLUDE 'All.f' C NAFAR(1) = I1 NAFAR(2) = I2 NAFAR(3) = I3 NAFAR(4) = I4 NAFAR(5) = I5 DO I = 1, NHPMAX - 1 DO J = I + 1, NHPMAX IF (NAFAR(I).LT.NAFAR(J)) THEN NTMP = NAFAR(I) NAFAR(I) = NAFAR(J) NAFAR(J) = NTMP ENDIF ENDDO ENDDO C END C======================================================================= c C======================================================================= C c Calculate matrix element between the single states C======================================================================= DOUBLE PRECISION FUNCTION EMAD(N1,N2,N3,N4,IFPN) INCLUDE 'All.f' c emad=1!fhtda(inp) c return C c (IFPN.EQ.0) no pn pairing IF (IFPN.EQ.1) INP=1 !neutrons IF (IFPN.EQ.2) INP=2 !protons IF (N1.LT.0) THEN JZA = SS2K(-N1,1,INP) IPA = SSIP(-N1,1,INP) ELSE JZA = SS2K(N1,2,INP) IPA = SSIP(N1,2,INP) ENDIF IF (N2.LT.0) THEN JZB = SS2K(-N2,1,INP) IPBB = SSIP(-N2,1,INP) ELSE JZB = SS2K(N2,2,INP) IPBB = SSIP(N2,2,INP) ENDIF IF (N3.LT.0) THEN JZC = SS2K(-N3,1,INP) IPCC = SSIP(-N3,1,INP) ELSE JZC = SS2K(N3,2,INP) IPCC = SSIP(N3,2,INP) ENDIF IF (N4.LT.0) THEN JZD = SS2K(-N4,1,INP) IPD = SSIP(-N4,1,INP) ELSE JZD = SS2K(N4,2,INP) IPD = SSIP(N4,2,INP) ENDIF C ETMP = 0D0 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 DO IH = 1, NGHP DO IL = 1, NGL PHIPA = SPF(IH,IL,N1,INP) PHIPB = SPF(IH,IL,N2,INP) PHIPC = SPF(IH,IL,N3,INP) PHIPD = SPF(IH,IL,N4,INP) PHIMA = SMF(IH,IL,N1,INP) PHIMB = SMF(IH,IL,N2,INP) PHIMC = SMF(IH,IL,N3,INP) PHIMD = SMF(IH,IL,N4,INP) ETMP = ETMP+(PHIMA*PHIPB*PHIMC*PHIPD+PHIPA*PHIMB*PHIPC*PHIMD- & PHIMA*PHIPB*PHIPC*PHIMD-PHIPA*PHIMB*PHIMC*PHIPD)* & GH(IH)*GL(IL) ENDDO ENDDO ETMP = ETMP*FHTDA(INP)/BZP2/(2D0*PI) 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 ETMP = ETMP*WFU(N1,INP)*WFU(N2,INP)*WFU(N3,INP)*WFU(N4,INP) END IF EMAD = ETMP END C======================================================================= C houda c Calculate matrix element between the single states C======================================================================= DOUBLE PRECISION FUNCTION EMAQ(N1,N2,N3,N4) INCLUDE 'All.f' C IF (N1.LT.0) THEN JZA = SS2K(-N1,1,inp) IPA = SSIP(-N1,1,inp) ELSE JZA = SS2K(N1,2,inp) IPA = SSIP(N1,2,inp) ENDIF IF (N2.LT.0) THEN JZB = SS2K(-N2,1,inp) IPBB = SSIP(-N2,1,inp) ELSE JZB = SS2K(N2,2,inp) IPBB = SSIP(N2,2,inp) ENDIF IF (N3.LT.0) THEN JZC = SS2K(-N3,1,inp) IPCC = SSIP(-N3,1,inp) ELSE JZC = SS2K(N3,2,inp) IPCC = SSIP(N3,2,inp) ENDIF IF (N4.LT.0) THEN JZD = SS2K(-N4,1,inp) IPD = SSIP(-N4,1,inp) ELSE JZD = SS2K(N4,2,inp) IPD = SSIP(N4,2,inp) ENDIF C QTMP = 0D0 C =S_mu *-* C Parity IF (IPA.EQ.IPCC.AND.IPBB.EQ.IPD) THEN C Omega IF (JZA+JZB.EQ.JZC+JZD) THEN IF (IABS(JZA-JZC).LE.4) THEN QTMP = QTMP + QHF(N1,N3)*QHF(N2,N4) ENDIF ENDIF ENDIF C Parity IF (IPA.EQ.IPD.AND.IPBB.EQ.IPCC) THEN C Omega IF (JZA+JZB.EQ.JZC+JZD) THEN IF (IABS(JZA-JZD).LE.4) THEN QTMP = QTMP - QHF(N1,N4)*QHF(N2,N3) ENDIF ENDIF ENDIF QTMP = QTMP * Q00(INP) * 1D-0 EMAQ = QTMP END C******************************************************************* c Calculate matrix element between the single states C======================================================================= DOUBLE PRECISION FUNCTION EMAO(N1,N2,N3,N4) INCLUDE 'All.f' C IF (N1.LT.0) THEN JZA = SS2K(-N1,1,inp) IPA = SSIP(-N1,1,inp) ELSE JZA = SS2K(N1,2,inp) IPA = SSIP(N1,2,inp) ENDIF IF (N2.LT.0) THEN JZB = SS2K(-N2,1,inp) IPBB = SSIP(-N2,1,inp) ELSE JZB = SS2K(N2,2,inp) IPBB = SSIP(N2,2,inp) ENDIF IF (N3.LT.0) THEN JZC = SS2K(-N3,1,inp) IPCC = SSIP(-N3,1,inp) ELSE JZC = SS2K(N3,2,inp) IPCC = SSIP(N3,2,inp) ENDIF IF (N4.LT.0) THEN JZD = SS2K(-N4,1,inp) IPD = SSIP(-N4,1,inp) ELSE JZD = SS2K(N4,2,inp) IPD = SSIP(N4,2,inp) ENDIF C OTMP = 0D0 C =S_mu *-* C Parity IF (IPA.EQ.IPCC.AND.IPBB.EQ.IPD) THEN C Omega IF (JZA+JZB.EQ.JZC+JZD) THEN IF (IABS(JZA-JZC).LE.4) THEN OTMP = OTMP + OCHF(N1,N3)*OCHF(N2,N4) ENDIF ENDIF ENDIF C Parity IF (IPA.EQ.-IPD.AND.IPBB.EQ.-IPCC) THEN C Omega IF (JZA+JZB.EQ.JZC+JZD) THEN IF (IABS(JZA-JZD).LE.4) THEN OTMP = OTMP - OCHF(N1,N4)*OCHF(N2,N3) ENDIF ENDIF ENDIF OTMP = OTMP * OC00(INP) * 1D-0 EMAO = OTMP END C======================================================================= C To check C Output the main HTDA results on a file C======================================================================= SUBROUTINE OPHTDA INCLUDE 'All.f' C EPHT = EPHM(1) + EPHM(2) ECORER = EVALUE(1,1) + EVALUE(1,2) ECOR1 = EVALUE(1,1) - EPHM(1) ECOR2 = EVALUE(1,2) - EPHM(2) ECORT = ECOR1 + ECOR2 OPEN(UNIT=16,FILE='htdainfo') K2N = NUK2(1) IPN = NUIP(1) K2P = NUK2(2) IPP = NUIP(2) NURK2T = K2N + K2P NURIPT = IPN * IPP WRITE(16,*) NURK2T,PTY(NURIPT),K2N,PTY(IPN),K2P,PTY(IPP), & EPHM,EPHT,ECOR1,ECOR2,ECORT,ECORER CLOSE(16) C c 100 FORMAT(3(I3,2H,2,A,1H),6(1H,F7.3,1H,1H),1H,F7.3,5H,/, c & 6H,hline) END C======================================================================= C To check C Intrinsic state of the system C======================================================================= SUBROUTINE INTST(N,EAMP,EPHMP,EVECC,NFILE,PFILE) INCLUDE 'All.f' C C SUBROUTINE PARAMETERS C CHARACTER*(*) NFILE,PFILE DOUBLE PRECISION EAMP, EPHMP, EVECC(N) C C Local variable C INTEGER I, II, JJ IF (INP.EQ.1) THEN OPEN(UNIT=30,FORM='UNFORMATTED',FILE=NFILE) ELSEIF (INP.EQ.2) THEN OPEN(UNIT=30,FORM='UNFORMATTED',FILE=PFILE) ENDIF WRITE(30) N, EAMP, EPHMP WRITE(30) (EVECC(I),I=1,N) if (inp.eq.1) then WRITE(30) (((MBSTN(II,JJ,I),II=1,NHPMAX),JJ=1,2),I=1,N) else WRITE(30) (((MBSTN(II,JJ,I),II=1,NHPMAX),JJ=1,2),I=1,N) end if CLOSE(30) C END C======================================================================= C Count the number of pair of a given configuration C======================================================================= SUBROUTINE CNPAIR(I,IORD,NDPH,NORDER,NPHOLE,NPPAR,NPAIR) INCLUDE 'All.f' C IF (NORDER.EQ.1) THEN IORDH = IORD IORDP = IORD + NDPH ELSE IORDH = IORD + NDPH IORDP = IORD ENDIF NPHOLE = 0 NPPAR = 0 if (inp.eq.1) then C Count number of pair of hole configuration 1 IF (IORDH.GE.2) THEN IF (MOD(MBST(IORDH,1,I),2).NE.0) THEN IF (MBST(IORDH-1,1,I).EQ.MBST(IORDH,1,I)+1) THEN NPHOLE = NPHOLE + 1 IORDH = IORDH - 2 ELSE IORDH = IORDH - 1 ENDIF ELSE IORDH = IORDH - 1 ENDIF GOTO 1 ENDIF C Count number of pair of particles configuration 2 IF (IORDP.GE.2) THEN IF (MOD(MBST(IORDP,2,I),2).NE.0) THEN IF (MBST(IORDP-1,2,I).EQ.MBST(IORDP,2,I)+1) THEN NPPAR = NPPAR + 1 IORDP = IORDP - 2 ELSE IORDP = IORDP - 1 ENDIF ELSE IORDP = IORDP - 1 ENDIF GOTO 2 ENDIF end if if (inp.eq.2) then C Count number of pair of hole configuration 11 IF (IORDH.GE.2) THEN IF (MOD(MBST(IORDH,1,I),2).NE.0) THEN IF (MBST(IORDH-1,1,I).EQ.MBST(IORDH,1,I)+1) THEN NPHOLE = NPHOLE + 1 IORDH = IORDH - 2 ELSE IORDH = IORDH - 1 ENDIF ELSE IORDH = IORDH - 1 ENDIF GOTO 11 ENDIF C Count number of pair of particles configuration 22 IF (IORDP.GE.2) THEN IF (MOD(MBST(IORDP,2,I),2).NE.0) THEN IF (MBST(IORDP-1,2,I).EQ.MBST(IORDP,2,I)+1) THEN NPPAR = NPPAR + 1 IORDP = IORDP - 2 ELSE IORDP = IORDP - 1 ENDIF ELSE IORDP = IORDP - 1 ENDIF GOTO 22 ENDIF end if C NPAIR = MIN0(NPHOLE,NPPAR) C END C======================================================================= C To check C Calculate the correlated ground state C======================================================================= SUBROUTINE NFIRSTSTATES 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 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 C Diagonalise the matrix C 10 CALL ITIME(TIME) PRINT 100, TIME WRITE(20,100) TIME c CALL newDIAG c &(NVALMAX,NE,NTOT1(inp),-neigen,INP,INF,EVEC1,EVALUE) 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 REWIND(IFIL) c CALL NEWLAN(ZV,ZPREC,NTOT1(inp),NEIGEN,NLOOP,IFIL) c c************ reading eigenvectors from ifil ********************* REWIND(IFIL) DO I=1,NEIGEN READ(IFIL)(EVEC1(K,I),K=1,NTOT1(INP)) WRITE(60,*) (K,EVEC1(K,I),K=1,NTOT1(INP)) END DO c***************************************************************** c c rearragement due to time-reversal symmetry conservation c DO L=1,NEIGEN DO K=1,NTOT(INP) EVEC(K) =0D0 VAB(K,L)=0D0 END DO DO I=1,NTOT1(INP) K =PHIPLUS(I,1) KBAR=PHIPLUS(I,2) K1 =IABS(KBAR) 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 C C MAXPI : Index of the configuration of Maximal probability C PROBM : Maximal Probability C PROBL = 1D-3 CALL VMAXC(NTOT(inp),EVEC,MAXPI) CALL VMAXCS(NTOT(inp),EVEC,11,MAXPIS) EPHM(INP) = EPH(MAXPI,INP) DO K=1,NEIGEN EVALUE(K,INP)=ZV(K) ENDDO PRINT 101, EVALUE(1,INP) WRITE(20,101) EVALUE(1,INP) PROBM = EVEC(MAXPI)**2 PRINT 102, MAXPI, & ((MBST(I,J,MAXPI),I = 1, NHPMAX), J = 1, 2), & EPH(MAXPI,inp),PROBM WRITE(20,102) MAXPI, & ((MBST(I,J,MAXPI),I = 1, NHPMAX), J = 1, 2), & EPH(MAXPI,inp),PROBM WRITE(20,103) 20,PROBL, &(K-1,MAXPIS(K),((MBST(I,J,MAXPIS(K)),I = 1, NHPMAX), J = 1, 2), & EPH(MAXPIS(K),inp),EVEC(MAXPIS(K))**2,K=2,11) C Determine the contribution of each order of configuration DO I = 0, 10 NCONTR(I) = 0 CONTRI(I) = 0D0 NCONTP(I) = 0 CONTRP(I) = 0D0 ENDDO DO I = 1, NTOT(inp) IF (I.NE.MAXPI) THEN !Do not take the principal configuration IORD = IORDMA 1 IF (MBST(IORD,NORDER,I).NE.0) GOTO 2 IORD = IORD-1 IF (IORD.EQ.0) GOTO 2 GOTO 1 2 NCONTR(IORD) = NCONTR(IORD) + 1 CONTRI(IORD) = CONTRI(IORD) + EVEC(I)**2 CALL CNPAIR(I,IORD,NDPH,NORDER,NPHOLE,NPPAR,NPAIR) IF (NPAIR.GT.0) THEN NCONTP(IORD) = NCONTP(IORD) + 1 CONTRP(IORD) = CONTRP(IORD) + EVEC(I)**2 ENDIF ENDIF ENDDO C NNOTM = 0 CNOTM = 0D0 NNOTMP = 0 CNOTMP = 0D0 DO I = 0, 10 IF (NCONTR(I).NE.0) THEN NNOTM = NNOTM + NCONTR(I) CNOTM = CNOTM + CONTRI(I) NNOTMP = NNOTMP + NCONTP(I) CNOTMP = CNOTMP + CONTRP(I) PRINT 107,I,CONTRI(I),NCONTR(I),CONTRP(I),NCONTP(I) WRITE(20,107) I,CONTRI(I),NCONTR(I),CONTRP(I),NCONTP(I) ENDIF ENDDO PRINT 108,CNOTM,NNOTM,CNOTMP,NNOTMP WRITE(20,108) CNOTM,NNOTM,CNOTMP,NNOTMP C IF (MNHTDA.EQ.1.AND.MNHF.EQ.0) THEN C C Write out the intrinsic state to be used in a calculation C that take into account the rotation C CALL INTST(NTOT(inp),EVALUE(1,INP),EPHM(INP), & EVECBIS,'nstate0','pstate0') C C Write out htda infomation for ground state C (We do this task only when the calculation is based on C a previously calculated htda one) C Calculate Ep-h(n,p), Ecorr(n,p) C C IF (INP.EQ.2.AND.KOIP.EQ.1) CALL OPHTDA ENDIF 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),/, &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),0 P F12.3,2P F10.4,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======================================================================= C To check C Density in coordinate space for 0+ state C======================================================================= SUBROUTINE DENSIS(IBLOCK) INCLUDE 'All.f' C N = NVB(IBLOCK) NHBM1 = NHB(IBLOCK) - 1 DO I = 1, N IA = NHBM1 + I NSA = MS(IA) LA = ML(IA) XLA = DBLE(LA) NZA = MZ(IA) NRA = MR(IA) X41 = BP*XLA*DBLE(NSA)/BZ DO J = 1, I Z = A(I,J) IF (HTDA) THEN ZKP = AKP(I,J) ZKN = AKN(I,J) ZHF = AHF(I,J) ENDIF IF (DABS(Z).GE.1D-7.OR.(DABS(ZHF).GE.1D-7.AND.HTDA).OR. & (DABS(ZKN).GE.1D-7.AND.HTDA).OR. & (DABS(ZKP).GE.1D-7.AND.HTDA)) THEN IB = NHBM1 + J NSB = MS(IB) LB = ML(IB) XLB = DBLE(LB) NZB = MZ(IB) NRB = MR(IB) X32 = BZ2*DBLE(NZA+NZB+1)+BP2*DBLE(2*(NRA+NRB+LA+1)) IF (NSA.NE.NSB) THEN LC = (LA + LB - 1)/2 C SAMSB : (Sigma_A minus Sigma_B)*2 SAMSB = DBLE(NSA-NSB) DO IL = 1, NGL ETA = XL(IL) ETALC = 1D0 IF (LC.NE.0) ETALC = ETA**LC QL00 = QL(IL,NRA,LA)*QL(IL,NRB,LB) QL01 = QL(IL,NRA,LA)*QL1(IL,NRB,LB) QL10 = QL1(IL,NRA,LA)*QL(IL,NRB,LB) DO IH = 1, NGHP QH01 = QH(IH,NZA)*QH1(IH,NZB) QH10 = QH1(IH,NZA)*QH(IH,NZB) X10 = QH01*QL10 - QH10*QL01 X11 = QL00*(XLA*QH01 + XLB*QH10) X1 = X10*SAMSB - 2D0*X11 IF (DABS(Z).GE.1D-7) & DIVJ(IH,IL,INP) = DIVJ(IH,IL,INP) + Z*ETALC*X1 IF (HTDA.AND.DABS(ZHF).GE.1D-7) & DJHF(IH,IL,INP) = DJHF(IH,IL,INP) + ZHF*ETALC*X1 ENDDO ENDDO ELSE XLL = XLA*XLB FAC = 2D0 IF (IA.EQ.IB) FAC = 1D0 DO IL = 1, NGL ETA = XL(IL) ETALA= 1D0 IF (LA.NE.0) ETALA= ETA**LA QL00 = QL(IL,NRA,LA)*QL(IL,NRB,LB) QL01 = QL(IL,NRA,LA)*QL1(IL,NRB,LB) QL10 = QL1(IL,NRA,LA)*QL(IL,NRB,LB) QL11 = QL1(IL,NRA,LA)*QL1(IL,NRB,LB) DO IH = 1, NGHP XI = XH(IH) QH00 = QH(IH,NZA)*QH(IH,NZB) QH11 = QH1(IH,NZA)*QH1(IH,NZB) X1 = QL00*QH00*ETALA X21 = QH11*QL00*ETA*BZ2 X22 = (QL11 + XLL*QL00)*QH00*BP2 X2 = (X21+X22)*ETALA/ETA X31 = BZ2*XI*XI + BP2*ETA X3 = (X31-X32)*QH00*QL00*ETALA X4 = QH00*(QL01+QL10)*ETALA/ETA*X41 IF (DABS(Z).GE.1D-7) THEN RHO(IH,IL,INP) = RHO(IH,IL,INP) + Z*X1*FAC TAU(IH,IL,INP) = TAU(IH,IL,INP) + Z*X2*FAC DRHO(IH,IL,INP) = DRHO(IH,IL,INP) + Z*X3*FAC DIVJ(IH,IL,INP) = DIVJ(IH,IL,INP) + Z*X4*FAC ENDIF IF (HTDA.AND.DABS(ZKP).GE.1D-7) & ROKP(IH,IL,INP) = ROKP(IH,IL,INP) + ZKP*X1*FAC IF (HTDA.AND.DABS(ZKN).GE.1D-7) & ROKN(IH,IL,INP) = ROKN(IH,IL,INP) + ZKN*X1*FAC IF (HTDA.AND.DABS(ZHF).GE.1D-7) THEN RHF(IH,IL,INP) = RHF(IH,IL,INP) + ZHF*X1*FAC THF(IH,IL,INP) = THF(IH,IL,INP) + ZHF*X2*FAC DRHF(IH,IL,INP) = DRHF(IH,IL,INP) + ZHF*X3*FAC DJHF(IH,IL,INP) = DJHF(IH,IL,INP) + ZHF*X4*FAC ENDIF ENDDO ENDDO ENDIF ENDIF ENDDO ENDDO C END C======================================================================= C Diagonalization la matrice A NxN C Maximal dimension : NM C======================================================================= SUBROUTINE DIAGO(A,NM,N) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A(NM,*) DIMENSION D(N),E(N),Z(N,N) DATA EPS/9D-20/ DATA TOL/1D-36/ C Z(1,1) = 1D0 IF (N.LE.1) GOTO 10 DO I = 1, N DO J = 1, N Z(I,J) = A(I,J) ENDDO ENDDO DO II = 2, N L = N - II I = L + 2 F = Z(I,I-1) G = 0D0 IF (L.NE.0) THEN DO K = 1, L G = G + Z(I,K)**2 ENDDO ENDIF H = G + F*F IF (G.LE.TOL) THEN E(I) = F H = 0D0 ELSE L = L + 1 G =-DSIGN(DSQRT(H),F) E(I) = G H = H - F*G Z(I,I-1) = F - G F = 0D0 DO J = 1, L Z(J,I) = Z(I,J)/H G = 0D0 DO K = 1, J G = G + Z(J,K)*Z(I,K) ENDDO J1 = J + 1 IF (J1.LE.L) THEN DO K = J1, L G = G + Z(K,J)*Z(I,K) ENDDO ENDIF E(J) = G/H F = F + G*Z(J,I) ENDDO HH = F/(H+H) DO J = 1, L F = Z(I,J) G = E(J) - HH*F E(J) = G DO K = 1, J Z(J,K) = Z(J,K) - F*E(K) - G*Z(I,K) ENDDO ENDDO ENDIF D(I) = H ENDDO D(1) = 0D0 E(1) = 0D0 DO I = 1, N L = I - 1 IF (DABS(D(I)).GE.1D-36) THEN IF (L.NE.0) THEN DO J = 1, L G = 0D0 DO K = 1, L G = G + Z(I,K)*Z(K,J) ENDDO DO K = 1, L Z(K,J) = Z(K,J) - G*Z(K,I) ENDDO ENDDO ENDIF ENDIF D(I) = Z(I,I) Z(I,I) = 1D0 IF (L.NE.0) THEN DO J = 1, L Z(I,J)=0D0 Z(J,I)=0D0 ENDDO ENDIF ENDDO DO I = 2, N E(I-1) = E(I) ENDDO E(N) = 0D0 B = 0D0 F = 0D0 DO L = 1, N J = 0 H = EPS * (DABS(D(L))+DABS(E(L))) IF (B.LT.H) B = H DO M = L, N IF (DABS(E(M)).LE.B) GOTO 1 ENDDO 1 IF (M.NE.L) THEN 2 IF (J.EQ.30) THEN WRITE(20,*) 'Problem in Diagonalization a block' CLOSE(20) STOP 'Problem in Diagonalization a block' ENDIF J = J + 1 P = (D(L+1)-D(L))/(2D0*E(L)) R = DSQRT(P*P + 1D0) H = D(L) - E(L)/(P+DSIGN(R,P)) DO I = L, N D(I) = D(I) - H ENDDO F = F + H P = D(M) C = 1D0 S = 0D0 M1 = M - 1 ML = M1 + L DO II = L, M1 I = ML - II G = C*E(I) H = C*P IF (DABS(P).GE.DABS(E(I))) THEN C = E(I)/P R = DSQRT(C*C+1D0) E(I+1) = S*P*R S = C/R C = 1D0/R ELSE C = P/E(I) R = DSQRT(C*C+1D0) E(I+1) = S*E(I)*R S = 1D0/R C = C/R ENDIF P = C*D(I) - S*G D(I+1) = H + S*(C*G+S*D(I)) DO K = 1, N H = Z(K,I+1) Z(K,I+1) = S*Z(K,I) + C*H Z(K,I) = C*Z(K,I) - S*H ENDDO ENDDO E(L) = S*P D(L) = C*P IF (DABS(E(L)).GT.B) GOTO 2 ENDIF D(L) = D(L) + F ENDDO N1 = N - 1 DO I = 1, N1 K = I P = D(I) II = I + 1 DO J = II, N IF (D(J).LT.P) THEN K = J P = D(J) ENDIF ENDDO IF (K.NE.I) THEN D(K) = D(I) D(I) = P DO J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P ENDDO ENDIF ENDDO DO I = 1, N A(I,I) = D(I) ENDDO 10 CONTINUE C END C======================================================================= C To check C Density in the single states basis C This can be used for the self-consistant calculation C of the odd nuclei C======================================================================= SUBROUTINE NDENSS INCLUDE 'All.f' C C Initialize the density matrix in sp basis. C C rho_ij = C =Sum_m Sum_n Chi_m Chi_n C Phi_m and Phi_n have the same K/pi so the different between Phi_n C and Phi_m is 0+ => could differ from 0 only C when i and j belonging to the same block K/pi C So we can discribe rho by two different matrices corresponding to C K>0 and K<0. C open(unit=4,file='rho') DO I = 1, NOV DO J = I, NOV DENS(I,J) = 0D0 DENS(J,I) = 0D0 DENSN(I,J) = 0D0 DENSN(J,I) = 0D0 ENDDO ENDDO IF (INP.EQ.1) NF0=24 IF (INP.EQ.2) NF0=44 C print *,'dens=',dens(i,j),ndiff1 C If we don't have pairing by HTDA the distribution of rho C is a Fermi distribution DO II = 1, NUN(INP)/2 C DENS and DENSN is indexed using order before the rearrangement of C energies I = NBRE(II) DENS(I,I) = 1D0 DENSN(I,I) = 1D0 ENDDO C Contribution to non-diagonal terms of density matrix REWIND(NF0) READ(NF0) NDIFF1,NDIF0H,NDIF0P DO II = 1, NDIFF1 READ(NF0) IPHAS,ISWT,JSWT,I,J ITMP = NBRE(IABS(ISWT)) JTMP = NBRE(IABS(JSWT)) IF (ISWT.GT.0.AND.JSWT.GT.0) THEN C 2 states with K>0 IF (IPHAS.EQ.1) THEN DENS(ITMP,JTMP) = DENS(ITMP,JTMP) + EVEC(I)*EVEC(J) ELSE IF (IPHAS.EQ.-1) THEN DENS(ITMP,JTMP) = DENS(ITMP,JTMP) - EVEC(I)*EVEC(J) ELSE PRINT *,'Strange IPHAS is neither 1 nor -1' ENDIF ELSEIF (ISWT.LT.0.AND.JSWT.LT.0) THEN C 2 states with K<0 IF (IPHAS.EQ.1) THEN DENSN(ITMP,JTMP) = DENSN(ITMP,JTMP) + EVEC(I)*EVEC(J) ELSE IF (IPHAS.EQ.-1) THEN DENSN(ITMP,JTMP) = DENSN(ITMP,JTMP) - EVEC(I)*EVEC(J) ELSE PRINT *,'Strange IPHAS is neither 1 nor -1' ENDIF ELSE PRINT *,'Strange : one state with K>0 and the other with K<0' ENDIF ENDDO C Contribution to the diagonal terms of density matrix C Substraction by hole states DO II = 1,NDIF0H READ(NF0) ISWT,I ITMP = NBRE(IABS(ISWT)) IF (ISWT.GT.0) THEN C state with K>0 DENS(ITMP,ITMP) = DENS(ITMP,ITMP) - EVEC(I)*EVEC(I) C print *,'dens=',dens(itmp,itmp) ENDIF IF (ISWT.LT.0) THEN C state with K<0 DENSN(ITMP,ITMP) = DENSN(ITMP,ITMP) - EVEC(I)*EVEC(I) ENDIF ENDDO C Addition by particle states DO II = 1,NDIF0P READ(NF0) ISWT,I ITMP = NBRE(IABS(ISWT)) IF (ISWT.GT.0) THEN C state with K>0 DENS(ITMP,ITMP) = DENS(ITMP,ITMP) + EVEC(I)*EVEC(I) C print *,'dens=',dens(itmp,itmp) ENDIF IF (ISWT.LT.0) THEN C state with K<0 DENSN(ITMP,ITMP) = DENSN(ITMP,ITMP) + EVEC(I)*EVEC(I) C print *,'dens=',dens(itmp,itmp) ENDIF ENDDO C We have stored the different between Phi_m and Ph_n where n>m C So to take in to account the contribution of Phi_m and Ph_n where n0 or K<0 matrix C DMAND: Density maximal non diagonal C DMIND: Density maximal non diagonal NNDT = (ISUP-IINF+1)*(ISUP-IINF) EPS = 1D-7 NNZNDP = 0 NNZNDM = 0 DMANDP = 0D0 DMINDP = 1D0 DMANDM = 0D0 DMINDM = 1D0 C We will use A(NW,NW) to store windowed matrix element of rho C in canonic basis DO II = IINF,ISUP C IA = II-IINF+1 C I = NBRE(II) DO JJ = II+1, ISUP J = NBRE(JJ) C For K>0 IF (DENS(J,I).EQ.0) THEN DENS(J,I) = DENS(I,J) ELSE DENS(I,J) = DENS(J,I) + DENS(I,J) DENS(J,I) = DENS(I,J) ENDIF C JA = JJ-IINF+1 C A(IA,JA) = DENS(I,J) A(JA,IA) = DENS(I,J) C IF (DABS(DENS(I,J)).GT.EPS) NNZNDP = NNZNDP + 2 IF (DENS(I,J).GT.DMANDP) DMANDP = DENS(I,J) IF (DENS(I,J).LT.DMINDP) DMINDP = DENS(I,J) C For K<0 IF (DENSN(J,I).EQ.0) THEN DENSN(J,I) = DENSN(I,J) ELSE DENSN(I,J) = DENSN(J,I) + DENSN(I,J) DENSN(J,I) = DENSN(I,J) ENDIF C AHF(IA,JA) = DENSN(I,J) AHF(JA,IA) = DENSN(I,J) C IF (DABS(DENSN(I,J)).GT.EPS) NNZNDM = NNZNDM + 2 IF (DENSN(I,J).GT.DMANDM) DMANDM = DENSN(I,J) IF (DENSN(I,J).LT.DMINDM) DMINDM = DENSN(I,J) ENDDO C A(IA,IA) = DENS(I,I) AHF(IA,IA) = DENSN(I,I) C ENDDO C Write information to draw the density matrix in excitation window DO II = IINF, ISUP DO JJ = IINF, ISUP I = NBRE(II) J = NBRE(JJ) C For K>0 WRITE(30,*) 'K>0',DENS(I,J), II-IINF+1,JJ-IINF+1 C For K<0 WRITE(30,*) 'K<0',DENSN(I,J), II-IINF+1,JJ-IINF+1 ENDDO ENDDO CLOSE(30) C Write out the number of nucleons in each bloc before filling C approximation C IF (INP.EQ.1) THEN OPEN(UNIT=30,FILE='nblockfillingbefore') ELSEIF (INP.EQ.2) THEN OPEN(UNIT=30,FILE='pblockfillingbefore') ENDIF WRITE(30,*) PAP = 0D0 PAM = 0D0 DO ICB = 1, NOB PNP = 0D0 PNM = 0D0 IM = NHB(ICB) - 1 N = NVB(ICB) DO I = 1, N PNP = PNP + DENS(I+IM,I+IM) PNM = PNM + DENSN(I+IM,I+IM) ENDDO WRITE(30,*) ICB,K2B(ICB),PTY(IPB(ICB)),PNP,N, & -K2B(ICB),PTY(IPB(ICB)),PNM PAP = PAP + PNP PAM = PAM + PNM TRABP(ICB) = PNP TRABM(ICB) = PNM ENDDO WRITE(30,*) PAP,PAM CLOSE(30) C RESO = 0D0 RES = 0D0 RESON = 0D0 RESN = 0D0 DO I = 1, NUN(INP)/2+NPL(inp) RES = RES + DENS(NBRE(I),NBRE(I)) RESN = RESN + DENSN(NBRE(I),NBRE(I)) ENDDO DO I = 1, NUN(INP)/2 RESO = RESO + DENS(NBRE(I),NBRE(I)) RESON = RESON + DENSN(NBRE(I),NBRE(I)) ENDDO TRAKP = RES TRAKM = RESN CC PRINT 100, 'For K > 0',RESO,TRAKP CC PRINT 100, 'For K < 0',RESON,TRAKM WRITE(20,100) 'For K > 0',RESO,TRAKP WRITE(20,100) 'For K < 0',RESON,TRAKM C IF (INP.EQ.1) THEN OPEN(UNIT=30,FILE='HTDAndis') ELSEIF (INP.EQ.2) THEN OPEN(UNIT=30,FILE='HTDApdis') ENDIF C Canonical basis NTAIL = NHL(inp)+NPL(inp) CALL DIAGO(A,NW,NTAIL) CALL DIAGO(AHF,NW,NTAIL) C C K>0 SUIVI = 0D0 WRITE(30,*) (EHF(NUN(INP)/2)+EHF(NUN(INP)/2+1))/2,WFX,WFMU DO I = ISUP, IINF, -1 RHOII = DENS(NBRE(I),NBRE(I)) VI = DSQRT(RHOII) UI = DSQRT(1-RHOII) SUIVI = SUIVI + UI*VI WRITE(30,*) 'K>0 ',I,EHF(I),RHOII,VI,UI,UI*VI,IPB(MB(NBRE(I))) ENDDO WRITE(30,*) SUIVI C Canonical basis write(30,*) write(30,*) WRITE(30,*) '#Canonical basis' SUIVI = 0D0 SV2 = 0D0 DO I = 1, NTAIL c IF (A(I,I).LE.1.E-8) A(I,I)=0D0 RHOII = A(I,I) c PRINT*,INP,RHOII write(4,*)inp,rhoii VI = DSQRT(RHOII) UI = DSQRT(1-RHOII) SUIVI = SUIVI + UI*VI SV2 = SV2 + RHOII IR = ISUP - I + 1 WRITE(30,*) 'K>0 ',IR,EHF(IR),RHOII,VI,UI,UI*VI, & IPB(MB(NBRE(IR))) ENDDO WRITE(30,*) SUIVI SUUIVI(1,INP) = SUIVI SUV2(1,INP) = SV2 + IINF - 1 C write(30,*) write(30,*) WRITE(30,*) NNZNDP,NNDT,DMANDP,DMINDP C K<0 SUIVI = 0D0 DO I = ISUP, IINF, -1 RHOII = DENSN(NBRE(I),NBRE(I)) VI = DSQRT(RHOII) UI = DSQRT(1-RHOII) SUIVI = SUIVI + UI*VI WRITE(30,*) 'K<0 ',I,EHF(I),RHOII,VI,UI,UI*VI,IPB(MB(NBRE(I))) ENDDO WRITE(30,*) SUIVI C Canonical basis write(30,*) write(30,*) WRITE(30,*) '#Canonical basis' SUIVI = 0D0 SV2 = 0D0 DO I = 1, NTAIL c IF (AHF(I,I).LE.1E-8) AHF(I,I)=0D0 RHOII = AHF(I,I) * PRINT *,INP,RHOII c write(44,*)inp,rhoii VI = DSQRT(RHOII) UI = DSQRT(1-RHOII) SUIVI = SUIVI + UI*VI SV2 = SV2 + RHOII IR = ISUP - I + 1 WRITE(30,*) 'K<0 ',IR,EHF(IR),RHOII,VI,UI,UI*VI, & IPB(MB(NBRE(IR))) ENDDO WRITE(30,*) SUIVI SUUIVI(2,INP) = SUIVI SUV2(2,INP) = SV2 + IINF - 1 write(30,*) write(30,*) C WRITE(30,*) NNZNDM,NNDT,DMANDM,DMINDM CLOSE(30) C 100 FORMAT(6X,A,/, & 6X,24H#Sum up to Fermi level :,F21.4,/, & 6X,15H#Trace de Rho :,20X,F10.4) 101 FORMAT(14H#Fermi level :,F8.3,5X,2HX=,F5.2,5X,3HMu=,F5.2,/, &5X,6H#Index,5X,6HEnergy,7X,6Hrho_ii,11X,2HVi,11X,2HUi,8X,5HUi*Vi, &7H Parity) 102 FORMAT(A,I3,5(F13.4),I5) 103 FORMAT(43X,18H#Sum of U_i*V_i : ,F18.4) 107 FORMAT(47H#Bloc K/pi Number of particles Size of block, & 28H #K/pi Number of particles) 108 FORMAT(I3,I4,2H/2,A1,F21.5,I10,I8,2H/2,A1,F21.5) 109 FORMAT(21H#Number of particles:,F11.5,F42.5) 110 FORMAT(23H#Non-diagonal elements:,/,5X,11H#Non-zeros:,I5,1H/,I5, &5X, &5HMax =,F10.6,5X,5HMin =,F10.6) 111 FORMAT(A,F14.8,I5,I5) C END C======================================================================= SUBROUTINE Q20OVERLAP INCLUDE 'All.f' DIMENSION QWIN(NV,NV) DIMENSION PRHOP(NV,NV),PRHOM(NV,NV) C C Initialize the pseudo density matrix in sp basis. C C pseudo-rho_ij = C =Sum_m Sum_n Chi_m Chi'_n C Phi_m and Phi_n have the same K/pi so the different between Phi_n C and Phi_m is 0+ => could differ from 0 only C when i and j belonging to the same block K/pi C So we can describe pseudo-rho by two different matrices corresponding to C K>0 and K<0. C cc-->> IF (KOMAT.EQ.1) THEN CALL CALQHF ENDIF cc<<-- IMIN=NUN(INP)/2-NHL(inp)+1 IMAX=NUN(INP)/2+NPL(inp) DO I = IMIN,IMAX DO J = IMIN,IMAX QWIN(I,J) = 0D0 ENDDO ENDDO DO K = 1,NBW DO I=1,NSIB(K) N1 = ISB(K,I) DO J=1,NSIB(K) N2 = ISB(K,J) QWIN(ISW(N1,inp),ISW(N2,inp)) = 2D0 * QHF(N1,N2) ENDDO ENDDO ENDDO c cc-->> Check parity with respect to Time Reveral Operator eps1=1.d-4 eps2=5.d-5 do k=1,neigen ibreak=0 ieven=0 iodd=0 do i=1,ntot(inp) j=iabs(mbreversed(i,inp)) tsgn=dfloat(mbreversed(i,inp)/j) write(48,*) 'T |',i,'> = ',tsgn,' |',j,'>' write(48,*) ((mbst(m,n,i),m=2,1,-1),n=1,2),' --> ', # ((mbst(m,n,j),m=2,1,-1),n=1,2) xx=vab(i,k) yy=vab(j,k) if (i.eq.j) then if (dabs(xx).gt.eps1) then if (tsgn.gt.0.d0) ieven=ieven+1 if (tsgn.lt.0.d0) iodd=iodd+1 endif elseif (dabs(xx)+dabs(yy).gt.2.*eps1) then c if (k.eq.1) write(49,*) i,xx,j,yy if (dabs(tsgn*xx-yy).lt.eps2) then ieven=ieven+1 elseif (dabs(tsgn*xx+yy).lt.eps2) then iodd=iodd+1 else write(49,*) xx,'|',i,'> et ',yy,'|',j,'> ', # 'me les brise' ibreak=ibreak+1 endif endif enddo write(49,*) 'For ',k,'-th HTDA state : ',ieven,iodd,ibreak enddo cc<<-- DO K=2,NEIGEN DO I = IMIN,IMAX DO J = IMIN,IMAX PRHOP(I,J) = 0D0 PRHOM(I,J) = 0D0 ENDDO ENDDO C Contribution to non-diagonal terms of pseudo density matrix REWIND(NF0) READ(NF0) NDIFF1,NDIF0H,NDIF0P DO II = 1, NDIFF1 READ(NF0) IPHAS,ISWT,JSWT,I,J ITMP = (IABS(ISWT)) JTMP = (IABS(JSWT)) XX = DFLOAT(IPHAS)*EVEC(I)*VAB(J,K) YY = DFLOAT(IPHAS)*EVEC(J)*VAB(I,K) C if ((dabs(xx).ge.1.d-5).or.(dabs(yy).ge.1.d-5)) then C write(49,*) iphas,iswt,jswt,i,j C write(49,*) xx,yy C endif IF (ISWT.GT.0.AND.JSWT.GT.0) THEN PRHOP(ITMP,JTMP)=PRHOP(ITMP,JTMP) + XX PRHOP(JTMP,ITMP)=PRHOP(JTMP,ITMP) + YY ELSEIF (ISWT.LT.0.AND.JSWT.LT.0) THEN PRHOM(ITMP,JTMP)=PRHOM(ITMP,JTMP) + XX PRHOM(JTMP,ITMP)=PRHOM(JTMP,ITMP) + YY ENDIF ENDDO C Contribution to the diagonal terms of pseudo density matrix C Substraction of hole states/addition of particle states XHP = -1D0 DO II = 1,NDIF0H+NDIF0P READ(NF0) ISWT,I ITMP = (IABS(ISWT)) XX = XHP*EVEC(I)*VAB(I,K) cbktr IF (ISWT.GT.0) THEN PRHOP(ITMP,ITMP) = PRHOP(ITMP,ITMP) + XX ELSEIF (ISWT.LT.0) THEN PRHOM(ITMP,ITMP) = PRHOM(ITMP,ITMP) + XX ENDIF IF (II.EQ.NDIF0H) XHP=1D0 ENDDO c C Compute the Trace TRACEP=0D0 TRACEM=0D0 DO I=IMIN,IMAX DO J=IMIN,IMAX TRACEP=TRACEP+QWIN(I,J)*PRHOP(I,J) TRACEM=TRACEM+QWIN(I,J)*PRHOM(I,J) ENDDO ENDDO c if (tracep*tracem.lt.0) tracem=-tracem OVERLAP = TRACEP + TRACEM c Q10E(K,INP) = OVERLAP ENDDO C END C======================================================================= C Calculate for HTDA method C======================================================================= SUBROUTINE CALH00 INCLUDE 'All.f' C DOUBLE PRECISION SG(2),SQXL(NL),RENOHF(2) DATA SG/1D0,-1D0/ C NFERMI = NUN(INP)/2 C C Densities Initialization (for rho(Phi) where Phi is the HTDA vacuum) C DO I = 1, NGHP DO J = 1, NGL RHF(I,J,INP) = 0D0 THF(I,J,INP) = 0D0 DRHF(I,J,INP) = 0D0 DJHF(I,J,INP) = 0D0 ENDDO ENDDO C DO ICB = 1, NOB N = NVB(ICB) IM = NHB(ICB) - 1 C C --- Density in the assymptotic basis (the alpha basis) C C Density for the vacuum in HTDA calculation DO I = 1, N DO J = 1, N XHF = 0D0 DO K = 1, N IF (NARE(K+IM).LE.NFERMI) THEN XHFT = COEF(I,K,ICB,INP)*COEF(J,K,ICB,INP) XHF = XHF + XHFT + XHFT ENDIF ENDDO AHF(I,J) = XHF ENDDO ENDDO C C Density in coordinate space C---------------------- N = NVB(ICB) IM = NHB(ICB) - 1 DO I = 1, N IA = IM + I NSA = MS(IA) LA = ML(IA) XLA = DBLE(LA) NZA = MZ(IA) NRA = MR(IA) X41 = BP*XLA*DBLE(NSA)/BZ DO J = 1, I ZHF = AHF(I,J) IF (DABS(ZHF).GE.1D-7) THEN IB = IM + J NSB = MS(IB) LB = ML(IB) XLB = DBLE(LB) NZB = MZ(IB) NRB = MR(IB) X32 = BZ2*DBLE(NZA+NZB+1)+BP2*DBLE(2*(NRA+NRB+LA+1)) IF (NSA.NE.NSB) THEN LC = (LA + LB - 1)/2 C SAMSB : (Sigma_A minus Sigma_B)*2 SAMSB = DBLE(NSA-NSB) DO IL = 1, NGL ETA = XL(IL) ETALC = 1D0 IF (LC.NE.0) ETALC = ETA**LC QL00 = QL(IL,NRA,LA)*QL(IL,NRB,LB) QL01 = QL(IL,NRA,LA)*QL1(IL,NRB,LB) QL10 = QL1(IL,NRA,LA)*QL(IL,NRB,LB) DO IH = 1, NGHP QH01 = QH(IH,NZA)*QH1(IH,NZB) QH10 = QH1(IH,NZA)*QH(IH,NZB) X10 = QH01*QL10 - QH10*QL01 X11 = QL00*(XLA*QH01 + XLB*QH10) X1 = X10*SAMSB - 2D0*X11 DJHF(IH,IL,INP) = DJHF(IH,IL,INP) + ZHF*ETALC*X1 ENDDO ENDDO ELSE XLL = XLA*XLB FAC = 2D0 IF (IA.EQ.IB) FAC = 1D0 DO IL = 1, NGL ETA = XL(IL) ETALA= 1D0 IF (LA.NE.0) ETALA= ETA**LA QL00 = QL(IL,NRA,LA)*QL(IL,NRB,LB) QL01 = QL(IL,NRA,LA)*QL1(IL,NRB,LB) QL10 = QL1(IL,NRA,LA)*QL(IL,NRB,LB) QL11 = QL1(IL,NRA,LA)*QL1(IL,NRB,LB) DO IH = 1, NGHP XI = XH(IH) QH00 = QH(IH,NZA)*QH(IH,NZB) QH11 = QH1(IH,NZA)*QH1(IH,NZB) X1 = QL00*QH00*ETALA X21 = QH11*QL00*ETA*BZ2 X22 = (QL11 + XLL*QL00)*QH00*BP2 X2 = (X21+X22)*ETALA/ETA X31 = BZ2*XI*XI + BP2*ETA X3 = (X31-X32)*QH00*QL00*ETALA X4 = QH00*(QL01+QL10)*ETALA/ETA*X41 RHF(IH,IL,INP) = RHF(IH,IL,INP) + ZHF*X1*FAC THF(IH,IL,INP) = THF(IH,IL,INP) + ZHF*X2*FAC DRHF(IH,IL,INP) = DRHF(IH,IL,INP) + ZHF*X3*FAC DJHF(IH,IL,INP) = DJHF(IH,IL,INP) + ZHF*X4*FAC ENDDO ENDDO ENDIF ENDIF ENDDO ENDDO ENDDO C C Density Normalisation C DO IH = 1, NGHP XI2 = XH(IH)**2 DO IL = 1, NGL POWER = XI2 + XL(IL) X = BZP2PI*DEXP(-POWER) RHF(IH,IL,INP) = X*RHF(IH,IL,INP) THF(IH,IL,INP) = X*THF(IH,IL,INP) DRHF(IH,IL,INP) = 2D0*(THF(IH,IL,INP) + X*DRHF(IH,IL,INP)) DJHF(IH,IL,INP) = X*BPZ*DJHF(IH,IL,INP) ENDDO ENDDO C XHF = 0D0 DO I = 1, NGHP Z1 = GH(I) DO J = 1, NGL XHF = XHF + Z1*GL(J)*RHF(I,J,INP) ENDDO ENDDO C RENOHF(INP) = 5D-1*QP*XHF/BZP2 TMP = XHF/BZP2PI RENOHF(INP) = TMP + TMP C C Renormalisation C YHF = DBLE(NUN(INP))/RENOHF(INP) DO I = 1, NGHP DO J = 1, NGL RHF(I,J,INP) = YHF*RHF(I,J,INP) THF(I,J,INP) = YHF*THF(I,J,INP) DRHF(I,J,INP) = YHF*DRHF(I,J,INP) DJHF(I,J,INP) = YHF*DJHF(I,J,INP) ENDDO ENDDO IF (INP.EQ.2) THEN C C Coulomb for H00 C SG(1) = 1D0 SG(2) =-1D0 DO I = 1, NGL SQXL(I) = DSQRT(XL(I))/BP ENDDO DO I = 1, NGHP Z = XH(I)/BZ DO J = 1, NGL R = SQXL(J) SHF = 0D0 DO M = 1, 2 ZKK = SG(M)/BZ DO K = 1, NGHP GHK = GH(K) ZK = ZKK*XH(K) DO L = 1, NGL RL = SQXL(L) Y = (Z-ZK)**2 + (R+RL)**2 X = 1D0 - 4D0*RL*R/Y RJ = DSQRT(Y) E = 1D0 + X*(0.4630151+0.1077812*X) IF (X.GT.1D-6) THEN E = E - X*(0.2452727+0.0412496*X)*DLOG(X) ENDIF SHF = SHF + GHK*GL(L)*DRHF(K,L,INP)*E*RJ ENDDO ENDDO ENDDO VCHF(I,J) = EV*SHF/BZP2 ENDDO ENDDO C Calculate H00 COEFF =-75D-2*EV*(3D0/PI)**(1D0/3D0) H00 = 0D0 DO I = 1, NGHP Z = GH(I) DO J = 1, NGL Z1 = Z*GL(J) RHO1 = RHF(I,J,1) RHO2 = RHF(I,J,2) SRHO = RHO1+RHO2 TAU1 = THF(I,J,1) TAU2 = THF(I,J,2) STAU = TAU1+TAU2 DRHO1= DRHF(I,J,1) DRHO2= DRHF(I,J,2) SDRHO= DRHO1+DRHO2 DIVJ1= DJHF(I,J,1) DIVJ2= DJHF(I,J,2) SDIVJ= DIVJ1+DIVJ2 X = B1*SRHO**2 + B2*(RHO1*RHO1+RHO2*RHO2) X = X + B3*SRHO*STAU + B4*(RHO1*TAU1+RHO2*TAU2) X = X + (SRHO**XBET3)*(B7*SRHO**2+B8*(RHO1*RHO1+RHO2*RHO2)) X = X + HB0*STAU X = X + B5*SRHO*SDRHO+B6*(RHO1*DRHO1+RHO2*DRHO2) X = X + B9*(SRHO*SDIVJ+RHO1*DIVJ1+RHO2*DIVJ2) X = X + 5D-1*VCHF(I,J)*RHO2 + COEFF*(RHO2**(4D0/3D0)) H00 = H00 + Z1*X ENDDO ENDDO C X = 1D0/BZP2PI X = X + X C X: constant to take integral C 2D0 comme from the fact that we have dealed with only z>0 H00 = X * H00 C Calculate "les champs" H_HF(rho(Phi)) CALL CALCHS ENDIF C END C======================================================================= C Calculate "les champs" H_HF(rho(Phi)) C======================================================================= SUBROUTINE CALCHS INCLUDE 'All.f' C C BP4 = BP2*BP2 Y = EV*(3D0/PI)**(1D0/3D0) DO K = 1, 2 M = 3 - K C FK is used to add the coulomb field in the case of protons FK = DBLE(K-1) DO I = 1,NGHP XI = XH(I) Z3 = Q*XI*XI Z2 = 2D0*Z3 DO J = 1, NGL ETA = XL(J) RHOK= RHF(I,J,K) RHOM= RHF(I,J,M) RHOA= (RHOK+RHOM)**XBET3 TAUK= THF(I,J,K) TAUM= THF(I,J,M) DRHOK= DRHF(I,J,K) DRHOM= DRHF(I,J,M) DIVJK= DJHF(I,J,K) DIVJM= DJHF(I,J,M) X = 2D0*(B1*(RHOM+RHOK)+B2*RHOK)+B3*(TAUK+TAUM)+B4*TAUK X = X + 2D0*(B5*(DRHOK+DRHOM)+B6*DRHOK) X = X + B7*(2D0+XBET3)*RHOA*(RHOK+RHOM) X = X + B8*RHOA*(2D0*RHOK*(RHOK+RHOM) + & XBET3*(RHOK*RHOK+RHOM*RHOM))/(RHOK+RHOM) X = X + B9*(2D0*DIVJK+DIVJM) + FK*(VCHF(I,J)-Y*RHOK**(1D0/3D0)) CC X = X + XQ2*(Z2-ETA)/BP2 CC X = X + XQ4*(375D-3*ETA*ETA+Z3*Z3-3D0*ETA*Z3)/BP4 C CE(I,J,K) = X EM(I,J,K) = HB0 + B3*(RHOK+RHOM) + B4*RHOK SO(I,J,K) =-B9*(RHOM + RHOK + RHOK) ENDDO ENDDO ENDDO C END C======================================================================= C Calculate many-body matrix elements of residual interaction C======================================================================= FUNCTION RINEW(IA,IB) INCLUDE 'All.f' DIMENSION IEXOP(0:4) 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. CALL newMBCOMP(IA,IB,NPH,IEXOP) if (NPH.LT.0) then RINEW=ETMP RETURN endif IFPN=0 C C data written in different files for neutrons and protons C 25-diff1n,26-diff0hn,27-diff0p,125-diff1p,126-diff0hp,127-diff0pp c if (NPH.EQ.0) then ccc ccc Cas ou |A> et |B> sont égaux ccc iph=1 ! Somme_p et |B> diffèrent par une particule ccc if (NPH.EQ.1) then if (inp.eq.1) $ WRITE(25) IEXOP(0),ISW(IP1,inp),ISW(IT1,inp),IA,IB if (inp.eq.2) $ WRITE(45) IEXOP(0),ISW(IP1,inp),ISW(IT1,inp),IA,IB NDIFF1=NDIFF1+1 iph=1 ! Somme_p ... do while ( (MBST(iph,2,IA).NE.0) .AND. (iph.LE.NHPMAX) ) NP=MBST(iph,2,IA) ETMP=ETMP+EMAD(NP,IP1,NP,IT1,IFPN) iph=iph+1 enddo iph=1 ! ... - Somme_t do while ( (MBST(iph,1,IA).NE.0) .AND. (iph.LE.NHPMAX) ) NT=-MBST(iph,1,IA) ETMP=ETMP-EMAD(NT,IP1,NT,IT1,IFPN) iph=iph+1 enddo ccc ccc Cas ou |A> et |B> diffèrent par deux particules ccc else IP2=IEXOP(2) IT2=IEXOP(4) c print 12,IB,IEXOP(0),IA,IP1,IP2,IT1,IT2 ETMP=-EMAD(IP1,IP2,IT1,IT2,IFPN) endif RINEW = IEXOP(0)*ETMP c RETURN c END cccc SUBROUTINE MBCOMP(IA,IB,NPH,IEXOP) INCLUDE 'All.f' c DIMENSION IEXOP(0:4) DO I=0,4 IEXOP(I)=0 ENDDO c nhia=0 npia=0 do while (MBST(nhia+1,1,IA).NE.0) nhia=nhia+1 enddo do while (MBST(npia+1,2,IA).NE.0) npia=npia+1 enddo nhib=0 npib=0 do while (MBST(nhib+1,1,IB).NE.0) nhib=nhib+1 enddo do while (MBST(npib+1,2,IB).NE.0) npib=npib+1 enddo cc-->> On determine si IA et IB sont dans la petite fenetre ou non c itst=0 : grande fenetre c itst=1 : petite fenetre c itst=1 if((nhia.ne.0.and.mbst(1,1,ia).lt.2*(nhl(inp)-nhld(inp))+1).or. # (nhib.ne.0.and.mbst(1,1,ib).lt.2*(nhl(inp)-nhld(inp))+1) .or. # (npia.ne.0.and.mbst(npia,2,ia).gt.2*(npld(inp))) .or. # (npib.ne.0.and.mbst(npib,2,ib).gt.2*(npld(inp))))then itst=0 endif if (ia.eq.ib) Then NPH=0 RETURN endif cc<<-- if (nhia-npia.NE.nhib-npib) then NPH=-1 RETURN endif c c |IA> est la configuration d ordre le plus bas, pris comme vide c |IB> est la configuration d ordre le plus élevé c iposa=1 iposb=1 ipha=nhia iphb=nhib IEXOP(0)=1 ishift=0 do while ( (ipha.GE.1) .OR. (iphb.GE.1) ) if ( (iposa.GT.3) .OR. (iposb.GT.3) ) then NPH=-1 RETURN endif ista=MBST(ipha,1,IA) if (ipha.LT.1) ista=9000 istb=MBST(iphb,1,IB) if (iphb.LT.1) istb=9000 if (ista.EQ.istb) then ipha=ipha-1 iphb=iphb-1 elseif (ista.LT.istb) then IEXOP(iposb)=-ista iposb=iposb+1 ipha=ipha-1 if (MOD(ista+ishift+1,2).NE.0) IEXOP(0)=-IEXOP(0) else ! (istb.LT.ista) IEXOP(2+iposa)=-istb iposa=iposa+1 iphb=iphb-1 if (MOD(istb+ishift+1,2).NE.0) IEXOP(0)=-IEXOP(0) endif ishift=ishift-1 enddo c ipha=npia iphb=npib do while ( (ipha.GE.1) .OR. (iphb.GE.1) ) if ( (iposa.GT.3) .OR. (iposb.GT.3) ) then NPH=-1 RETURN endif ista=MBST(ipha,2,IA) if (ipha.LT.1) ista=9000 istb=MBST(iphb,2,IB) if (iphb.LT.1) istb=9000 if (ista.EQ.istb) then ipha=ipha-1 iphb=iphb-1 ishift=ishift+1 elseif (ista.LT.istb) then IEXOP(2+iposa)=+ista iposa=iposa+1 ipha=ipha-1 if (MOD(ishift,2).NE.0) IEXOP(0)=-IEXOP(0) else ! (istb.LT.ista) IEXOP(iposb)=+istb iposb=iposb+1 iphb=iphb-1 if (MOD(ishift,2).NE.0) IEXOP(0)=-IEXOP(0) endif enddo NPH=iposa-1 c if (NPH.EQ.2) THEN if (IEXOP(0).LT.0) THEN IEXOP(0)=1 else jph=IEXOP(3) IEXOP(3)=IEXOP(4) IEXOP(4)=jph endif endif c RETURN END cccc SUBROUTINE newMBCOMP(IA,IB,NPH,IEXOP) INCLUDE 'All.f' 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 c On determine si les deux etats sont dans la petite fenetre ou non c ITST=0 : grande fenetre c ITST=1 : petite fenetre c DIMENSION IEXOP(0:4*NHPMAX) dimension nab(2),nhpab(2,2) dimension nwick(-nhpmax:nhpmax,2),ncontr(2) cc-->> dimension mwick(-nhpmax:nhpmax,2) cc<<-- c NAB(1)=IA NAB(2)=IB NPH=0 ITST=1 c do iab=1,2 j=3-2*iab do ihp=1,2 i=NHPMAX do while (MBST(i,ihp,nab(iab)).eq.0.and.i.ge.1) nwick(i*j,ihp)=0 i=i-1 enddo nhpab(ihp,iab)=i do i=1,nhpab(ihp,iab) nwick(i*j,ihp)=MBST(nhpab(ihp,iab)-i+1,ihp,nab(iab)) enddo if (nhpab(ihp,iab).ne.0) then if (ihp.eq.1) then if (nwick(j,1).lt.2*(nhl(inp)-nhld(inp))+1) ITST=0 else if (mbst(1,2,nab(iab)).gt.2*npld(inp)) ITST=0 endif endif enddo if (NAB(1).eq.NAB(2)) RETURN ! états many-body identiques enddo c NPH=-1 if (nhpab(1,1)-nhpab(2,1).ne.nhpab(1,2)-nhpab(2,2)) RETURN ! differents nombres de nucléons if (iabs(nhpab(1,1)-nhpab(2,1)).gt.2) RETURN ! avec |n-m| > 2 c IEXOP(0)=nhpab(2,2)*(nhpab(1,1)+nhpab(1,2)) IEXOP(0)=0 IEXOP(2)=0 IEXOP(4)=0 cc-->> C rewind(54) C write(54) nwick cc<<-- MPH=0 do ihp=1,2 ncontr(ihp)=0 i=-nhpab(ihp,2) do while (i.le.-1) j=nhpab(ihp,1)-ncontr(ihp) do while (j.ge.1) if (nwick(j,ihp).eq.nwick(i,ihp)) then ncontr(ihp)=ncontr(ihp)+1 IEXOP(0)=IEXOP(0)-i+j do k=i,-nhpab(ihp,2)+ncontr(ihp),-1 nwick(k,ihp)=nwick(k-1,ihp) enddo nwick(-nhpab(ihp,2)+ncontr(ihp)-1,ihp)=0 do k=j,nhpab(ihp,1)-ncontr(ihp) nwick(k,ihp)=nwick(k+1,ihp) enddo nwick(nhpab(ihp,1)-ncontr(ihp)+1,ihp)=0 j=1 elseif (nwick(j,ihp).lt.nwick(i,ihp)) then j=1 endif j=j-1 enddo i=i+1 enddo MPH=MPH+nhpab(ihp,ihp)-ncontr(ihp) if (MPH.gt.2) RETURN ! difference greater than 2p-2h enddo c IEXOP(0)=IEXOP(0)+ncontr(1)+nhpab(2,2)*nhpab(2,1) IEXOP(0)=1-2*mod(IEXOP(0),2) NPH=MPH do ihp=1,2 j=3-2*ihp m=2*(ihp-1) do i=1,nhpab(ihp,1)-ncontr(ihp) m=m+1 IEXOP(m)=-j*nwick(i,ihp) enddo do i=-nhpab(3-ihp,2)+ncontr(3-ihp),-1 m=m+1 IEXOP(m)=j*nwick(i,3-ihp) enddo enddo C write(53,12) (IEXOP(I),I=0,4) 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 C======================================================================= C HTDA method to deal the pairing C======================================================================= SUBROUTINE HTDAME 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 test with one spectra c c if (inp.eq.1) CALL NEWEXCWIN1(DELE,DELE,DELED,DELED) 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 C Determine THE MANY-BODY BASIS {|npmh>} CALL NEWMBBASE C do i=1,ntot C do j=i,ntot C call newMBCOMP(i,j,kk,kk,iexop) C enddo C enddo C stop IF (KOMAT.EQ.0) THEN CALL HTDAMA ELSE READ(NF0) NDIFF1,NDIF0H,NDIF0P ENDIF C CALCULATE THE GROUND STATE + LOW LYING EXCITED STATES CALL NFIRSTSTATES cc-->> call calqhf xx=0d0 do l=1,nbw do i=1,nsib(l) nnh=isb(l,i) if (nnh.lt.0) then do j=1,nsib(l) nnp=isb(l,j) if (nnp.gt.0) then qelem=qhf(nnh,nnp)*2. energ=ehf(isw(nnp,inp))-ehf(isw(nnh,inp)) contrib=4*qelem*qelem*energ write(50,50) nnh,nnp,qelem,energ,contrib xx=xx+contrib endif enddo endif enddo enddo 50 format('<',i3,'| Q_20 |',i3,'> = ',f8.4,' ; E_ph = ',f7.3, # ' ==> ',f8.2) c IF (MOD(KORINT,2).EQ.1) THEN IF (NEIGEN.GT.1) THEN CALL Q20OVERLAP ENDIF c ENDIF cc<<-- C CALCULATE THE DENSITY IN SP. BASIS. c CALL ITIME(TIME) c PRINT 101, TIME c WRITE(20,101) TIME c CALL NDENSS 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 To check C Reconstruction of density in coordinate space C======================================================================= SUBROUTINE RECOND INCLUDE 'All.f' C NF = NUN(INP)/2 C C Densities Initialization C DO I = 1, NGHP DO J = 1, NGL RHO(I,J,INP) = 0D0 TAU(I,J,INP) = 0D0 DRHO(I,J,INP) = 0D0 DIVJ(I,J,INP) = 0D0 IF (HTDA) THEN ROKP(I,J,INP) = 0D0 ROKN(I,J,INP) = 0D0 ENDIF ENDDO ENDDO IF (HTDA) THEN DO I = 1, NGHP DO J = 1, NGL RHF(I,J,INP) = 0D0 THF(I,J,INP) = 0D0 DRHF(I,J,INP) = 0D0 DJHF(I,J,INP) = 0D0 ENDDO ENDDO ENDIF C DO ICB = 1, NOB N = NVB(ICB) IM = NHB(ICB) - 1 C C --- Density in the assymptotic basis (the alpha basis) C IF (HTDA) THEN C Density for the vacuum in HTDA calculation DO I = 1, N DO J = 1, N XHF = 0D0 DO K = 1, N IF (NARE(K+IM).LE.NF) THEN XHFT = COEF(I,K,ICB,INP)*COEF(J,K,ICB,INP) XHF = XHF + XHFT + XHFT ENDIF ENDDO AHF(I,J) = XHF ENDDO ENDDO ENDIF C IF (.NOT.HTDA) THEN DO I = 1, N DO J = 1, N X = 0D0 DO K = 1, N XT = DENS(K+IM,K+IM)*COEF(I,K,ICB,INP)*COEF(J,K,ICB,INP) X = X + XT + XT ENDDO A(I,J) = X ODEN(I+IM,J+IM,INP) = A(I,J) ENDDO ENDDO ELSE DO I = 1, N DO J = 1, N X1 = 0D0 X2 = 0D0 DO II = 1, N DO JJ = 1, N X1 = X1 + DENS(II+IM,JJ+IM)*COEF(I,II,ICB,INP)* &COEF(J,JJ,ICB,INP) X2 = X2 + DENSN(II+IM,JJ+IM)*COEF(I,II,ICB,INP)* &COEF(J,JJ,ICB,INP) ENDDO ENDDO X = X1 + X2 AKP(I,J) = X1 AKN(I,J) = X2 A(I,J) = X ODEN(I+IM,J+IM,INP) = A(I,J) ENDDO ENDDO ENDIF C C Density in coordinate space C---------------------- CALL DENSIS(ICB) c****** C CALL NEWDENSIS(ICB) C---------------------- ENDDO C C Write out the number of nucleons in each bloc after filling C approximation C IF (INP.EQ.1) THEN OPEN(UNIT=30,FILE='nblockfillingafter') ELSEIF (INP.EQ.2) THEN OPEN(UNIT=30,FILE='pblockfillingafter') ENDIF WRITE(30,*) PA = 0D0 DO ICB = 1, NOB PN = 0D0 IM = NHB(ICB) - 1 N = NVB(ICB) DO I = 1, N PN = PN + ODEN(I+IM,I+IM,INP) ENDDO WRITE(30,*) ICB,K2B(ICB),PTY(IPB(ICB)),PN,N PA = PA + PN ENDDO WRITE(30,*) PA CLOSE(30) C C Density Normalisation C C WRITE(15,*) NCL(INP) DO IH = 1, NGHP XI2 = XH(IH)**2 DO IL = 1, NGL POWER = XI2 + XL(IL) X = BZP2PI*DEXP(-POWER) RHO(IH,IL,INP) = X*RHO(IH,IL,INP) TAU(IH,IL,INP) = X*TAU(IH,IL,INP) DRHO(IH,IL,INP) = 2D0*(TAU(IH,IL,INP) + X*DRHO(IH,IL,INP)) DIVJ(IH,IL,INP) = X*BPZ*DIVJ(IH,IL,INP) IF (HTDA) THEN ROKP(IH,IL,INP) = X*ROKP(IH,IL,INP) ROKN(IH,IL,INP) = X*ROKN(IH,IL,INP) RHF(IH,IL,INP) = X*RHF(IH,IL,INP) THF(IH,IL,INP) = X*THF(IH,IL,INP) DRHF(IH,IL,INP) = 2D0*(THF(IH,IL,INP) + X*DRHF(IH,IL,INP)) DJHF(IH,IL,INP) = X*BPZ*DJHF(IH,IL,INP) ENDIF ENDDO ENDDO C**** C CALL NORMALISERHO C 100 FORMAT(46HBloc K/pi Number of particles Size of block) 101 FORMAT(I3,I4,2H/2,A1,F21.5,I10) 102 FORMAT(20HNumber of particles:,F11.5) END C======================================================================= C Resolution of the stationary Schroedinger equation C======================================================================= SUBROUTINE SCHRD INCLUDE 'All.f' C DIMENSION VEC(N,N) C C Find all the eigenvectors and eigenvalue of the one-body Hamiltonian C (Block diagonal) C C ICB : I Counter for Blocks WRITE(20,100) NCL(INP) C PRINT 100,NCL(INP) WRITE(20,101) C PRINT 101 DO ICB = 1, NOB N = NVB(ICB) C Hamiltonian matrix elements calculation for one block C---------------------- CALL CAHFME(ICB) C Diagonalization the block C------------------------ CALL DIAGON(N,ICB) C------------------------ DO I = 1, N J = NHB(ICB) + I - 1 EHF(J) = A(I,I) NBRE(J)= J ENDDO ENDDO WRITE(20,102) C**************************************************** c***************************************************************** C******************************************** C PRINT 102 C Sort single states according to increasing energies C Calculate parity and the projection of total momentum of each C single state DO I = 1, NOV-1 DO J = I+1, NOV IF (EHF(I).GT.EHF(J)) THEN X = EHF(I) EHF(I) = EHF(J) EHF(J) = X ITMP = NBRE(I) NBRE(I)= NBRE(J) NBRE(J)= ITMP ENDIF ENDDO ENDDO C IF (HTDA) THEN DO I = 1, NOV NARE(NBRE(I)) = I ENDDO ENDIF C Density distribution without pairing correction NF = NUN(INP)/2 DO I = 1, NOV DO J = 1, NOV DENS(I,J) = 0D0 ENDDO ENDDO C HF00(INP) = 0D0 C DO I = 1, NF K = NBRE(I) DENS(K,K) = 1D0 HF00(INP) = HF00(INP) + EHF(I) ENDDO C Have to double HF00 because we have taken into account only C state with K>0 HF00(INP) = HF00(INP) + HF00(INP) C Write out the levels IF (INP.EQ.1) THEN OPEN(UNIT = 30,FILE='nlevels') ELSE OPEN(UNIT = 30,FILE='plevels') ENDIF WRITE(30,*) DO I = NOV, 1, -1 K = NBRE(I) IB = MB(K) WRITE(30,*) I,EHF(I),K2B(IB),PTY(IPB(IB)) ENDDO CLOSE(30) C Density distribution with pairing correction C Pairing Calculation or Densities in the base of individual states IF (.NOT.HTDA.AND.MNBCS.GT.0) THEN CALL BCS ELSEIF (HTDA) THEN C Calculate CALL CALH00 ENDIF C IF (MNBCS.EQ.0) THEN XYX = 0D0 DO I = 1, NF XYX = XYX + EHF(I) + EHF(I) ENDDO ESOMIN = ESOMIN + XYX WRITE(20,*) 'Without BCS ESOMIN = ',ESOMIN,'XYX = ',XYX ENDIF C WRITE(20,103) C Reconstruction of density in coordinate space IF (HF) CALL RECOND C 100 FORMAT(/,15HResolution for ,A) 101 FORMAT(39H-Matrix calculation and diagonalization) 102 FORMAT(48H-Density distribution in the single states basis) 103 FORMAT(43H-Calculation of density in coordinate space) 104 FORMAT(4(F19.11)) 105 FORMAT(5HIndex,9X,6HEnergy,9X,4HK/Pi) c 106 FORMAT(I5,1H&,F15.3,1H&,I10,4H/2$^,A,10H$hline) C END C======================================================================= C Second HF iteration for HTDA energies C======================================================================= SUBROUTINE NSCHRD INCLUDE 'All.f' C C DO ICB = 1, NOB N = NVB(ICB) CALL CAHFME(ICB) CALL CALEHF(ICB) ENDDO C C Sort single states according to increasing energies C Calculate parity and the projection of total momentum of each C single state DO I = 1, NOV-1 DO J = I+1, NOV IF (EHF(I).GT.EHF(J)) THEN X = EHF(I) EHF(I) = EHF(J) EHF(J) = X ITMP = NBRE(I) NBRE(I)= NBRE(J) NBRE(J)= ITMP ENDIF ENDDO ENDDO C DO I = 1, NOV NARE(NBRE(I)) = I ENDDO C Density distribution without pairing correction NFERMI = NUN(INP)/2 DO I = 1, NOV DO J = 1, NOV DENS(I,J) = 0D0 ENDDO ENDDO C TMP = 0D0 C DO I = 1, NFERMI K = NBRE(I) DENS(K,K) = 1D0 TMP = TMP + EHF(I) ENDDO C Have to double HF00 because we have taken into account only C K>0 HF00(INP) = TMP + TMP C Write out the levels IF (INP.EQ.1) THEN OPEN(UNIT = 30,FILE='htdanlevels') ELSE OPEN(UNIT = 30,FILE='htdaplevels') ENDIF WRITE(30,*) DO I = NOV, 1, -1 K = NBRE(I) IB = MB(K) WRITE(30,*) I,EHF(I),K2B(IB),PTY(IPB(IB)) ENDDO CLOSE(30) IF (KOPN.EQ.0) CALL HTDAME IF (KOPN.EQ.1) CALL HTDAMEPN C CALCULATE THE DENSITY IN SP. BASIS. IF (KOPN.EQ.0) CALL NDENSS IF (KOPN.EQ.0) CALL RECOND IF (KOPN.EQ.1.AND.INP.EQ.2) THEN DO INP=1,2 CALL NDENSS CALL RECOND END DO END IF C 100 FORMAT(/,15HResolution for ,A) 101 FORMAT(39H-Matrix calculation and diagonalization) 102 FORMAT(48H-Density distribution in the single states basis) 103 FORMAT(43H-Calculation of density in coordinate space) 104 FORMAT(4(F19.11)) 105 FORMAT(5HIndex,9X,6HEnergy,9X,4HK/Pi) c 106 FORMAT(I5,1H&,F15.3,1H&,I10,4H/2$^,A,10H$\\\\ \\hline) C END C======================================================================= C To check C Root mean square radii - Quadrupolar and Hexadecapolar moments C Densities renormalisation C======================================================================= SUBROUTINE MOMENT INCLUDE 'All.f' DIMENSION BE2(2),VOL(4),RENORM(2),RENOHF(2) C DO K = 1, 2 X = 0D0 XR = 0D0 XQ = 0D0 XQ4 = 0D0 DO I = 1, NGHP Z = Q*XH(I)*XH(I) Z1 = GH(I) ZQ = 2D0*Z ZQ4 = 2D0*ZQ*ZQ DO J = 1, NGL Z2 = XL(J) Z3 = Z1*GL(J)*RHO(I,J,K) X = X + Z3 XR = XR + Z3*(Z2+Z) XQ = XQ + Z3*(ZQ-Z2) XQ4= XQ4+ Z3*(ZQ4+Z2*(-12D0*ZQ + 3D0*Z2)) ENDDO ENDDO RENORM(K) = 5D-1*QP*X/BZP2 RMS(K) = DSQRT(XR/(X*BP2)) RMS2 = RMS(K)*RMS(K) QM(K) = 5D-1*QP*XQ/(BZP2*BP2) BE2(K) = DSQRT(5D0*QP)*QM(K)/(DBLE(NUN(K))*RMS2*10D0) DEF(K) = QM(K)/(DBLE(NUN(K))*RMS2) Q40(K) = DSQRT(9D0*QP)*XQ4/(16D0*BZP2*BP2*BP2) VOL(K) = (2D0*DBLE(NUN(K))*RMS2 - & QM(K))*DSQRT(DBLE(NUN(K))*RMS2+QM(K)) ENDDO R = DSQRT(RMS(2)**2+64D-2) VOL(4) = DSQRT((DBLE(NUN(1))*RMS(1)*RMS(1)+ & DBLE(NUN(2))*RMS(2)*RMS(2))/DBLE(NUN(1)+NUN(2))) Q20 = QM(1) + QM(2) Q4M = Q40(1) + Q40(2) WRITE(20,100) RMS,R,VOL(4) C This line for look in the detail of proton=>for the test of r2 WRITE(20,*) RMS,R,VOL(4) C PRINT 100,RMS,R,VOL(4) WRITE(20,101) QM(1)/1D2,QM(2)/1D2,Q20/1D2 C PRINT 101,QM(1)/1D2,QM(2)/1D2,Q20/1D2 WRITE(20,103) Q40(1)/1D4,Q40(2)/1D4,Q4M/1D4 WRITE(20,104) BE2(1),BE2(2) C C Renormalisation C DO K = 1, 2 Y = DBLE(NUN(K))/RENORM(K) DO I = 1, NGHP DO J = 1, NGL RHO(I,J,K) = Y*RHO(I,J,K) TAU(I,J,K) = Y*TAU(I,J,K) DRHO(I,J,K) = Y*DRHO(I,J,K) DIVJ(I,J,K) = Y*DIVJ(I,J,K) ENDDO ENDDO ENDDO C 100 FORMAT(8H-Radii :,/, &5X,10HNeutrons =,F13.4,3H fm,9X,9HProtons =,F13.4,3H fm, &5X,10HCharge =,F13.4,3H fm,9X,9HMass =,F13.4,3H fm) 101 FORMAT(20H-Quadrupole moments:,/,5X,10HQneutrons=,F13.4,2H b, &10X,9HQprotons=,F13.4,2H b,/,40X,9HQmass =,F13.4,2H b) 102 FORMAT(24HThe true charged radius:,F13.5) 103 FORMAT(22H-Hexadecapole moments:,/,5X,10HHneutrons=,F13.4,2Hb2, &10X,9HHprotons=,F13.4,2Hb2,/,40X,9HHmass =,F13.4,2Hb2) 104 FORMAT(6H-Beta2:,/,5X,10HBneutrons=,F13.4, &10X,9HBprotons=,F13.4) C END C======================================================================= C To check C Direct Coulombian potential calculation C======================================================================= SUBROUTINE COULOM INCLUDE 'All.f' DIMENSION SG(2),SQXL(NL) C SG(1) = 1D0 SG(2) =-1D0 DO I = 1, NGL SQXL(I) = DSQRT(XL(I))/BP ENDDO C DO I = 1, NGHP Z = XH(I)/BZ DO J = 1, NGL R = SQXL(J) S = 0D0 DO M = 1, 2 ZKK = SG(M)/BZ DO K = 1, NGHP GHK = GH(K) ZK = ZKK*XH(K) DO L = 1, NGL RL = SQXL(L) Y = (Z-ZK)**2 + (R+RL)**2 X = 1D0 - 4D0*RL*R/Y RJ = DSQRT(Y) E = 1D0 + X*(0.4630151+0.1077812*X) IF (X.GT.1D-6) THEN E = E - X*(0.2452727+0.0412496*X)*DLOG(X) ENDIF S = S + GHK*GL(L)*DRHO(K,L,2)*E*RJ ENDDO ENDDO ENDDO VC(I,J) = EV*S/BZP2 ENDDO ENDDO C END C======================================================================= C To check C Energy density and total energy C======================================================================= SUBROUTINE ENERGY INCLUDE 'All.f' c c OPEN(UNIT=10,FILE='intensite',ACCESS='APPEND') C COEFF =-75D-2*EV*(3D0/PI)**(1D0/3D0) EKIN = 0D0 EVOL = 0D0 ESUR = 0D0 ESOB = 0D0 EDIR = 0D0 EXCH = 0D0 EREA = 0D0 T3 = (24D0*B7 + 12D0*B8)*2D0/3D0 DO I = 1, NGHP Z = GH(I) DO J = 1, NGL Z1 = Z*GL(J) RHO1 = RHO(I,J,1) RHO2 = RHO(I,J,2) SRHO = RHO1+RHO2 TAU1 = TAU(I,J,1) TAU2 = TAU(I,J,2) STAU = TAU1+TAU2 DRHO1= DRHO(I,J,1) DRHO2= DRHO(I,J,2) SDRHO= DRHO1+DRHO2 DIVJ1= DIVJ(I,J,1) DIVJ2= DIVJ(I,J,2) SDIVJ= DIVJ1+DIVJ2 X = B1*SRHO**2 + B2*(RHO1*RHO1+RHO2*RHO2) X = X + B3*SRHO*STAU + B4*(RHO1*TAU1+RHO2*TAU2) X = X + (SRHO**XBET3)*(B7*SRHO**2+B8*(RHO1*RHO1+RHO2*RHO2)) EVOL = EVOL + X*Z1 EKIN = EKIN + Z1*STAU ESUR = ESUR + Z1*(B5*SRHO*SDRHO+B6*(RHO1*DRHO1+RHO2*DRHO2)) ESOB = ESOB + Z1*B9*(SRHO*SDIVJ+RHO1*DIVJ1+RHO2*DIVJ2) EDIR = EDIR + Z1*5D-1*VC(I,J)*RHO2 EXCH = EXCH + Z1*(RHO2**(4D0/3D0)) EREA = EREA + Z1*RHO1*RHO2*SRHO*T3 ENDDO ENDDO C X = 2D0*PI/BZP2 C 2D0 comme from the fact that we have dealed with only z>0 EKIN = X*EKIN*HB0 EVOL = X*EVOL ESUR = X*ESUR ESOB = X*ESOB EDIR = X*EDIR EXCH = X*EXCH*COEFF ECOUL= EDIR + EXCH ETOT = EKIN + EVOL + ESUR + ESOB + ECOUL + EPAIR C HHF00= HF00(1) + HF00(2) EREA =-X*EREA/8D0 + EXCH/3D0 C WRITE(20,100) EKIN,EVOL,ESUR,ESOB,ECOUL,EDIR,EXCH IF (.NOT.HTDA.AND.MNBCS.GT.0) WRITE(20,101) EPAIR +EPAIR PRINT 102,ETOT,HF00(1),HF00(2),HHF00,EREA,KORINT WRITE(20,102) ETOT,HF00(1),HF00(2),HHF00,EREA,KORINT C IF (HTDA) THEN IF (KOPN.EQ.0) THEN EHTDA = H00 + EVALUE(1,1) + EVALUE(1,2) PRINT 103, H00,(EVALUE(1,K),K=1,2),EHTDA WRITE(20,103) H00,(EVALUE(1,K),K=1,2),EHTDA ELSE EHTDA = H00 +EVALUEPN(1) PRINT 106, H00, EVALUEPN(1),EHTDA WRITE(20,106) H00,EVALUEPN(1),EHTDA END IF ENDIF C 100 FORMAT(30X,7H EKIN =,F13.5,5X,6HEVOL =,F13.5,/, & 30X,7H ESURF=,F13.5,5X,6HESO =,F13.5,/, &4X,7H ECOUL=,F13.5,6X,7H EDIR =,F13.5,5X,6HECOEX=,F13.5) 101 FORMAT(30X,13H EPAIR (BCS)=,F31.5) 102 FORMAT(31X,43(1H-),/,30X,7H ETOT =,F37.5,5H(MeV),//, & 30X,7H ,18X,7H(n) :,F12.5,/, & 55X,7H(p) :,F12.5,/,55X,19(1H-),/, & 55X,7HTotal :,F12.5,5H(MeV),/, & 31X,6HEREAR=,F37.5,5H(MeV),/,'INTERAC=',I2) 103 FORMAT(/,31X,15H =,15X,F13.5,/, & 31X,15HElanczos(n) =,15X,F13.5,/, & 31X,15HElanczos(p) =,15X,F13.5,/,31X,43(1H-),/, & 31X,15HEHTDA =,15X,F13.5,5H(MeV)) 106 FORMAT(/,31X,15H =,15X,F13.5,/, & 31X,15HElanczos(pn) =,15X,F13.5,/,31X,43(1H-),/, & 31X,15HEHTDA =,15X,F13.5,5H(MeV)) 104 FORMAT(20H-Energy Correction :,/, &5X,22H (n):,F12.5,/,23X,4H(p):,F12.5,/, &23X,16(1H-),/,23X,F16.5,5H(MeV), &5X,8HE shift:,F12.5,5H(MeV)) 105 FORMAT(40HDelta E'=<2|(VHF1-dHF1)-(VHF2-dHF2)|2> =,F12.5,/, & 40HC1 =,F12.5,/, & 40HC2 =,F12.5,/, & 40HDelta E = Delta E' + C1 - C2 =,F12.5,/, & 40HCorrected HTDA energy =,F12.5,5H(MeV)) c END C======================================================================= C To check C One-Body Hamiltonian Calculation C======================================================================= SUBROUTINE FIELDS INCLUDE 'All.f' C C XIM = MIX IF (HF) THEN XIM = HFM ELSEIF (HTDA) THEN XIM = HTDAM ENDIF Q40T = Q40(1) + Q40(2) BP4 = BP2*BP2 Y = EV*(3D0/PI)**(1D0/3D0) YMU = 1D0-XIM XQ2 = 2D0*C2*(Q20-Q2) XQ4 = 2D0*C4*(Q40T-Q4) C TPENTE =-C2*(Q20-Q2) PENTE = TPENTE + TPENTE C IF (HTDA) THEN XQ2 = 2D0*C2*(QM(1)+QM(2)-Q2) XQ4 = 2D0*C4*(Q40(1)+Q40(2)-Q4) ENDIF DO K = 1, 2 M = 3-K C FK is used to add the coulomb field in the case of protons FK = DBLE(K-1) DO I = 1,NGHP XI = XH(I) Z3 = Q*XI*XI Z2 = 2D0*Z3 DO J = 1, NGL ETA = XL(J) RHOK= RHO(I,J,K) RHOM= RHO(I,J,M) RHOA= (RHOK+RHOM)**XBET3 TAUK= TAU(I,J,K) TAUM= TAU(I,J,M) DRHOK= DRHO(I,J,K) DRHOM= DRHO(I,J,M) DIVJK= DIVJ(I,J,K) DIVJM= DIVJ(I,J,M) X = 2D0*(B1*(RHOM+RHOK)+B2*RHOK)+B3*(TAUK+TAUM)+B4*TAUK X = X + 2D0*(B5*(DRHOK+DRHOM)+B6*DRHOK) X = X + B7*(2D0+XBET3)*RHOA*(RHOK+RHOM) X = X + B8*RHOA*(2D0*RHOK*(RHOK+RHOM) + & XBET3*(RHOK*RHOK+RHOM*RHOM))/(RHOK+RHOM) X = X + B9*(2D0*DIVJK+DIVJM) + FK*(VC(I,J)-Y*RHOK**(1D0/3D0)) X = X + XQ2*(Z2-ETA)/BP2 X = X + XQ4*(375D-3*ETA*ETA+Z3*Z3-3D0*ETA*Z3)/BP4 CE(I,J,K) = XIM*X + YMU*CE(I,J,K) X = HB0 + B3*(RHOK+RHOM) + B4*RHOK EM(I,J,K) = XIM*X + YMU*EM(I,J,K) X =-B9*(RHOM + RHOK + RHOK) SO(I,J,K) = XIM*X + YMU*SO(I,J,K) ENDDO ENDDO ENDDO C WRITE(20,100) PENTE PRINT 100,PENTE C 100 FORMAT('Pente ',F12.4) C END C======================================================================= C To check C Write out HF potentials C======================================================================= SUBROUTINE WRFIEL INCLUDE 'All.f' C OPEN(UNIT=21,FORM='UNFORMATTED',FILE='HFfields.out') REWIND (21) WRITE(21) NGHP,NGL WRITE(21) (((EM(I,J,K),I = 1, NGHP),J = 1, NGL),K = 1, 2) WRITE(21) (((CE(I,J,K),I = 1, NGHP),J = 1, NGL),K = 1, 2) WRITE(21) (((SO(I,J,K),I = 1, NGHP),J = 1, NGL),K = 1, 2) c CLOSE(21) C END C======================================================================= C To check C Write out HF potentials C======================================================================= SUBROUTINE OUTPUT INCLUDE 'All.f' dimension contrib(2),sumrule(2),triplecom(2) C OPEN(UNIT = 30,FILE='basicinfos') OPEN(UNIT=3,FILE='M2',ACCESS='APPEND') OPEN(UNIT=4,FILE='r2',ACCESS='APPEND') C OPEN(UNIT=4,FILE='Q10',ACCESS='APPEND') IF (MNHTDA.GT.0) THEN KN = NUK2(1) KP = NUK2(2) KT = KN + KP IPN = NUIP(1) IPP = NUIP(2) IPT = IPN * IPP WRITE(30,*) ETOTAL = EHTDA Q20N = QM(1) Q20P = QM(2) Q20T = Q20N + Q20P Q40N = Q40(1) Q40P = Q40(2) Q40T = Q40N + Q40P RN = RMS(1) RP = RMS(2) SUVPN = SUUIVI(1,1) SUVPP = SUUIVI(1,2) SUVMN = SUUIVI(2,1) SUVMP = SUUIVI(2,2) SV2PN = SUV2(1,1) SV2PP = SUV2(1,2) SV2MN = SUV2(2,1) SV2MP = SUV2(2,2) WRITE(30,*) KT,PTY(IPT),ETOTAL,Q20N,Q20P,Q20T,Q40N,Q40P,Q40T, & KN,PTY(IPN),KP,PTY(IPP),RN,RP,SUVPN,SUVPP,SUVMN,SUVMP C********* c do inp=1,2 c sumrule(inp)=0 c do k=2,neigen c WRITE(*,*) k-1,Q10E(k,inp),EVALUE(k,inp), c # H00+EVALUE(k,inp)+evalue(1,3-inp) c write(3,*)k-1,DELE,Q10E(k,inp),EVALUE(k,inp), c # H00+EVALUE(k,inp)+evalue(1,3-inp) c enddo c triplecom(inp)=8d0*hb0*(2d0*dfloat(nun(inp))*RMS(inp)**2+QM(inp)) c enddo c do k=2,neigen c do inp=1,2 c contrib(inp)=2*q10e(k,inp)**2*(evalue(k,inp)-evalue(1,inp)) c sumrule(inp)=sumrule(inp)+contrib(inp) c enddo c write(64,12) k,contrib,(evalue(k,inp)-evalue(1,inp),inp=1,2) c enddo c write(64,*) c write(64,13) sumrule,triplecom,0,0 c write(64,14) (sumrule(i)/triplecom(i),i=1,2), c # (sumrule(i)/triplecom(i),i=1,2) 12 format(2x,i3,2x,2(f9.3,2x),2x,2(f7.3,2x)) 13 format('total',2x,2(f9.0,2x),/,'theo.',2x,2(f9.0,2x) # ,' (vide corrélé)',/,'theo.',2x,2(f9.0,2x) # ,' (vide HF)') 14 format(' % ',2x,2(2p f9.2,2x),' (vide corrélé)',/,' % ', # 2x,2(f9.2,2x),' (vide HF)') C***************** WRITE(30,*) WRITE(30,*) SV2PN,SV2PP,SV2MN,SV2MP WRITE(30,*) WRITE(30,*) KT,PTY(IPT),KN,PTY(IPN),KP,PTY(IPP), & H00,(EVALUE(1,K),K=1,2),EHTDA write(30,*) 'sumhtda=', evalue(1,1)+evalue(1,2) ENDIF CLOSE(30) C For the accuracy, store informations in binary file OPEN(UNIT = 30,FORM='UNFORMATTED',FILE='bbasicinfos') IF (MNHTDA.GT.0) THEN WRITE(30) ETOTAL,Q20N,Q20P,Q20T,Q40N,Q40P,Q40T,RN,RP, & SUVPN,SUVPP,SUVMN,SUVMP,SV2PN,SV2PP,SV2MN,SV2MP ENDIF CLOSE(30) C CALL WRFIEL C c 100 FORMAT(23HEssential informations:,/,6H,hline,/, c &58H,Calcul Energy$$(n)$$(p)$$(t), c &40H,$$(n)$$(p)$$(t)) c 101 FORMAT(1H$,I2,3H/2^,A1,2H$&,F9.3,3(1H&,F12.2),3(1H&,F12.1),2H, c &/,1H$,I2,3H,2^,A1,12H,nuotimes ,I2,3H,2^,A1,4H_pi,2H$&, c &2(1H&,F8.5),4(1H&,F7.4),2H) c 102 FORMAT(30X,11H&& r(n)&,4X,5Hr(p)&, c &8HSuv+(n)&,8HSuv+(p)&,8HSuv-(n)&,7HSuv-(p),2H\\\\) c 103 FORMAT(4(1H&,F7.4),/, c &8H&Sv2+(n),8H&Sv2+(p),8H&Sv2-(n),8H&Sv2-(p)) c 104 FORMAT(6H\\hline,/,33HEtat &$K^\\pi_n\\otimes K^\\pi_p$,6X,1H&, c &48H $H_{00}$&$E_\\text{lan.}(n)$&$E_\\text{lan.}(p)$&, c &43H$E_\\text{HTDA}$& $E^*$&$\\Delta E$&${E^*}$\\\\) c 105 FORMAT(1H$,I2,3H/2^,A1,2H$&,1H$,I2,3H/2^,A1,12H_\\nu\\otimes , c &I2,3H/2^,A1,4H_\\pi,2H$&,F9.3,2(1H&,F18.3),1H&,F15.3,1H&,F6.3,1H&, c &F10.3,1H&,F8.3,2H\\\\) c 106 FORMAT(6H\\hline,/,33HEtat &$K^\\pi_n\\otimes K^\\pi_p$,6X,1H&, c &48H $H_{00}$&$E_\\text{lan.}(n)$&$E_\\text{lan.}(p)$&, c &42H$E_\\text{HTDA}$& $E^*$&$E^*_\\text{exp.}$\\\\) c 107 FORMAT(1H$,I2,3H/2^,A1,2H$&,1H$,I2,3H/2^,A1,12H_\\nu\\otimes , c &I2,3H/2^,A1,4H_\\pi,2H$&,F9.3,2(1H&,F18.3),1H&,F15.3,1H&,6X,1H&, c &17X,2H\\\\) END C============================================================================== C This subroutine determine the convergence in function of 3 final C iteration and a predefined tolerance C============================================================================== SUBROUTINE CON3LI(RBBL,RBL,RL,DR,RCONV) LOGICAL RCONV DOUBLE PRECISION RBBL,RBL,RL,DR C IF (DABS(RBBL-RBL).LE.DR.AND.DABS(RBL-RL).LE.DR) THEN RCONV = .TRUE. ELSE RCONV = .FALSE. ENDIF RBBL = RBL RBL = RL C END C============================================================================== C This soubroutine use to check if one want to stop the running C program in an unusual way. C - If the file 'STOP' does not exist => do nothing (don't touch C NMNHF,NMNHTD) C - If 'STOP' existed => C + if cannot read NMNHF, NMNHTD Stop immediately (MNHF=-1) C + if NMNHF>=ITHF and NMNHTD>=ITHTDA then change MNHF, MNHTDA C + if the precedent is not verified then do nothing C============================================================================== SUBROUTINE SSTOP(MNHF,MNHTDA,ITHF,ITHTDA) C MNHFO = MNHF MNHTDO = MNHTDA OPEN(UNIT=19,FILE='STOP',STATUS='OLD',ERR=20) READ(19,*,ERR=10,END=10) MNHFT,MNHTDT IF (MNHFT.GE.ITHF.AND.MNHTDT.GE.ITHTDA) THEN MNHF = MNHFT MNHTDA = MNHTDT ENDIF GOTO 11 10 MNHF =-1 11 CLOSE (19) OPEN(UNIT=19,FILE='REPLYforSTOP') WRITE(19,100) ITHF+ITHTDA IF (MNHF.EQ.-1) THEN WRITE(19,101) ELSE WRITE(19,102) MNHF, MNHTDA ENDIF CLOSE(19) 20 CONTINUE C 100 FORMAT(14HWe are now at ,I3,11H iteration.) 101 FORMAT(51HThe program will stop at the end of this iteration.) 102 FORMAT(32HNew values for MNHF and MNHTDA :,I4,I10) END C============================================================================== C This soubroutine do the main calculation and is used for the C simplex method C============================================================================== SUBROUTINE CALCUL(BZRO,QZRO,CHI2,ENEW,Q0NEW,Q40CA) INCLUDE 'All.f' DOUBLE PRECISION OEM(NH/2,NL,2),OCE(NH/2,NL,2),OSO(NH/2,NL,2) C C OPEN(UNIT=15,FILE='density') C C HTDAconverg is used to store differents physics quantities in C a HTDA calculation, so the file is used only when we do HTDA C method IF (MNHTDA.GT.0) THEN OPEN(UNIT=91,FILE='HTDAconverg') C Different temporary files for HTDA calculation C HTDA states that differ by one nucleon C File in interaction with LASO package OPEN(UNIT=28,FORM='UNFORMATTED',FILE='iovec') ENDIF C HFconverg, like HTDAconverg, is used when we do HF and/or HF+BCS C method IF (MNHF.GT.0) THEN OPEN(UNIT=90,FILE='HFconverg') ENDIF OPEN(UNIT=80,FILE='convergence') OPEN(UNIT=81,FILE='noscilation') OPEN(UNIT=82,FILE='poscilation') C------------------------------------------------------------------ C This part is added to treat the case one want a HF calculation at C plusieurs mesh point in (b,q) C------------------------------------------------------------------ DOMESH = .FALSE. NMESH = 0 NMESHO = 0 IF (ICTRL.EQ.1) THEN C If we do simplex then B0 and Q is taken from simplex algorithm C and we need to recontruct the basis and the polynomials B0 = BZRO Q = QZRO CALL BASE CALL POLYHL GOTO 10 ENDIF C If the file "mesh" exist, and contain good information C we will used the mesh in B0,Q from this file. OPEN(UNIT=17,FILE='mesh',STATUS='OLD',ERR=10) READ(17,*,END = 10,ERR = 10) B0MIN,B0MAX,NB0,QMIN,QMAX,NQ WRITE(17,*) IF (NQ.LE.1) THEN WRITE(20,*) 'Number of meshpoint for q is negative' CLOSE (20) STOP 'Number of meshpoint for q is negative' ENDIF IF (NB0.LE.1) THEN WRITE(20,*) 'Number of meshpoint for b is negative' CLOSE (20) STOP 'Number of meshpoint for b is negative' ENDIF DOMESH = .TRUE. NMESH = NB0*NQ ENEMIN = 0D0 NMESHO = NMESH DQ = 0D0 IF (NQ.GT.1) THEN DQ = (QMAX-QMIN)/(NQ-1) ENDIF DB0 = 0D0 IF (NB0.GT.1) THEN DB0 = (B0MAX-B0MIN)/(NB0-1) ENDIF PRINT *,'Do calculation at ',NMESH,' meshpoints in (b,q)' WRITE(20,*) 'Do calculation at ',NMESH,' meshpoints in (b,q)' 10 IF (NMESH.GT.0.AND.ICTRL.EQ.0) THEN NMESHT = NMESHO - NMESH + 1 Q = QMIN + DQ*DBLE(MOD(NMESHT-1,NQ)) B0 = B0MIN + DB0*DBLE((NMESHT-1)/NQ) CALL BASE CALL POLYHL ENDIF PRINT *,'Beta_0 = ',B0,' q = ',Q WRITE(20,*) 'Beta_0 = ',B0,' q = ',Q C WRITE (20,100) IF (MOD(NUN(1),2).NE.0.OR.MOD(NUN(2),2).NE.0) & PRINT *,'This code is used only for pair-pair nuclids' C Initialise logical variable for execution of the program C EO2 = 2D0 EO1 = 1D0 Q202 = 2D10 Q201 = 1D10 Q402 = 2D10 Q401 = 1D10 C COE1 = .FALSE. COE2 = .FALSE. COQ201 = .FALSE. COQ202 = .FALSE. COQ401 = .FALSE. COQ402 = .FALSE. C COHF = .FALSE. COHTDA = .FALSE. C HF = .TRUE. HTDA = .FALSE. EHTDA = 0D0 IF (MNHF.EQ.0) THEN HF = .FALSE. HTDA = .TRUE. ENDIF ITHTDA = 0 ITHF = 0 C---------------- CALL INIPOT C**** NOUVEAU CALCUL DES POLYNOMES C CALL NEWPOLYHL C***** NOUVEAU CALCUL DE LA DENSITÉ C CALL NEWDENSIS(IBLOCK) C---------------- 90 IF (HF) ITHF = ITHF + 1 IF (HTDA) ITHTDA = ITHTDA + 1 C ITERR : Real Iteration ITERR = ITHF + ITHTDA C PRINT 103, ITERR WRITE(20,103) ITERR EPAIR = 0D0 ESOMIN = 0D0 C C Resolution of the stationary Schroedinger equation C WRITE(20,104) C PRINT 104 C DO INP = 1, 2 C Store the old values of CE, EM, SO IF (HTDA) THEN DO I = 1, NGHP DO J = 1, NGL OCE(I,J,INP) = CE(I,J,INP) OEM(I,J,INP) = EM(I,J,INP) OSO(I,J,INP) = SO(I,J,INP) ENDDO ENDDO ENDIF C---------------- CALL SCHRD C---------------- ENDDO C Second HF iteration for HTDA energies IF (HTDA) THEN C Restore values from stored potentiels DO INP = 1, 2 DO I = 1, NGHP DO J = 1, NGL CE(I,J,INP) = OCE(I,J,INP) EM(I,J,INP) = OEM(I,J,INP) SO(I,J,INP) = OSO(I,J,INP) ENDDO ENDDO OPEN(UNIT=25,FORM='UNFORMATTED',FILE='diff1n') OPEN(UNIT=26,FORM='UNFORMATTED',FILE='diff0hn') OPEN(UNIT=27,FORM='UNFORMATTED',FILE='diff0pn') OPEN(UNIT=45,FORM='UNFORMATTED',FILE='diff1p') OPEN(UNIT=46,FORM='UNFORMATTED',FILE='diff0hp') OPEN(UNIT=47,FORM='UNFORMATTED',FILE='diff0pp') OPEN(UNIT=24,FORM='UNFORMATTED',FILE='diffn') OPEN(UNIT=44,FORM='UNFORMATTED',FILE='diffp') c c REWIND(INF) c REWIND(24) c REWIND(44) CALL SSTOP(MNHF,MNHTDA,ITHF,ITHTDA) CALL NSCHRD c CLOSE(25,STATUS='DELETE') CLOSE(26,STATUS='DELETE') CLOSE(27,STATUS='DELETE') CLOSE(45,STATUS='DELETE') CLOSE(46,STATUS='DELETE') CLOSE(47,STATUS='DELETE') CLOSE(24) CLOSE(44) CLOSE(INF) ENDDO ENDIF C C Root mean square radii - Quadrupole and Hexadecapole moments C Densities renormalisation WRITE(20,105) C PRINT 105 C---------------- CALL MOMENT C---------------- C Direct Coulombian potential calculation WRITE(20,106) C PRINT 106 C---------------- CALL COULOM C---------------- C Energy density and total energy WRITE(20,107) CC PRINT 107 C---------------- CALL ENERGY C---------------- C One-Body Hamiltonian Calculation WRITE(20,108) CC PRINT 108 C------------ ---- CALL FIELDS C---------------- C RCHARG = DSQRT(RMS(2)**2+64D-2) IF (.NOT.HTDA) THEN ECONV = ETOT Q20T = Q20 Q40T = Q40(1) + Q40(2) CALL CON3LI(EO2,EO1,ECONV,EPSE,COE1) CALL CON3LI(Q202,Q201,Q20T,EPQ20,COQ201) CALL CON3LI(Q402,Q401,Q40T,EPQ40,COQ401) COHF = COE1.AND.COQ201.AND.COQ401 WRITE(80,114) ITERR,RCHARG,Q20,ETOT WRITE(90,116) ITHF,ETOT,Q20,RCHARG ELSE ECONV = EHTDA ECOR = EVALUE(1,1) + EVALUE(1,2) Q20T = QM(1) + QM(2) Q40T = Q40(1) + Q40(2) CALL CON3LI(EO2,EO1,ECONV,EPSE,COE2) CALL CON3LI(Q202,Q201,Q20T,EPQ20,COQ202) CALL CON3LI(Q402,Q401,Q40T,EPQ40,COQ402) COHTDA = COE2.AND.COQ202.AND.COQ402 WRITE(80,115) ITERR,RCHARG,Q20,EHTDA WRITE(91,118) ITHTDA,(EVALUE(1,K),K=1,2),ECOR,H00,EHTDA,ETOT, & Q20T,RCHARG,QM,Q40,RMS ENDIF WRITE(20,109) ITERR,Q20T,EO1 PRINT 109, ITERR,Q20T,EO1 C Here we will write out the fields in order to be able to use it C for another calculation in the case the program stop accidently. CALL WRFIEL C C C In the case of long calculations we give a possibility to stop the C running programme from exterier. C To do that, one need to put a file named 'STOP' in the C directory that the programme is running. CALL SSTOP(MNHF,MNHTDA,ITHF,ITHTDA) IF (MNHF.EQ.-1) GOTO 91 IF (.NOT.HTDA) THEN IF (.NOT.COHF) THEN C Don't converged before HTDA IF (ITHF.LT.MNHF) THEN ! Try again GOTO 90 ELSE ! Try HTDA HF = .FALSE. IF (MNHTDA.GT.0) THEN HTDA = .TRUE. GOTO 90 ENDIF ENDIF ELSE C converged before HTDA HF = .FALSE. IF (MNHTDA.GT.0) THEN HTDA = .TRUE. GOTO 90 ENDIF ENDIF ELSE IF (.NOT.COHTDA) THEN C Don't converged after HTDA IF (ITHTDA.LT.MNHTDA) GOTO 90 ! Try again ENDIF ENDIF C---------------- 91 CALL OUTPUT C---------------- C IF (NMESH.GT.0.AND.ICTRL.EQ.0) THEN IF (ENEMIN.GT.ETOT) THEN ENEMIN = ETOT B0MINE = B0 QMINE = Q Q2MINE = Q20/1D2 ITMINE = ITHF+ITHTDA ENDIF WRITE(17,112) B0,Q,Q2/1D2,Q20/1D2,ETOT,ITHF+ITHTDA IF (NMESH.NE.1) THEN NMESH = NMESH - 1 GOTO 10 ENDIF ENDIF C Write (b,q) corresponding to the minimal energy C only in the case we have more than one point in the mesh IF (NMESHO.GT.0) THEN IF (NMESHO.GT.1) & WRITE(17,113) B0MINE,QMINE,ENEMIN,Q2/1D2,Q2MINE,ITMINE NMESH = NMESHO CLOSE(17) ENDIF C Close files opened in HTDA calculation IF (MNHTDA.GT.0) THEN CLOSE(91) close(25) close(26) close(27) close(29) CLOSE(28,STATUS='DELETE') CLOSE(35,STATUS='DELETE') ENDIF IF (MNHF.GT.0) THEN C Close files opened in HF or HF+BCS calculation CLOSE(90) ENDIF CLOSE(80) CLOSE(81) CLOSE(82) C ENEW = ETOT Q0NEW = Q20 CHI2 = ENEW + 15000D0 Q40CA = Q40(1) + Q40(2) C 100 FORMAT(//,26HExecution of the program :,/,26(1H=),/) 101 FORMAT(26HThe vacuum is the nuclid :,/, &A,1H-,I3,2H :,5X,A,8X,4HN = ,I3,5X,4HZ = ,I3,10X, & 4HA = ,I3) 102 FORMAT(9HHF Part :,/,9(1H-)) 103 FORMAT(19H*Iteration Number :,I3) 104 FORMAT(/,40X,39HSolve the one-body Hamiltonian equation, &/,40X,39(1H-)) 105 FORMAT(/,17X,24HRoot mean square radii -, &38H Quadrupolar and Hexadecapolar moments,/,17X,62(1H-)) 106 FORMAT(/,62X,17HCoulomb Potential,/,62X,17(1H-)) 107 FORMAT(/,67X,12HTotal Energy,/,67X,12(1H-)) 108 FORMAT(/,30X,49HEffective Mass, central and spin-orbit potentials, &/,30X,49(1H-)) 109 FORMAT(/,79(1H-),/,18HEnd of iteration :,I3,5X,5HQ20 =,F11.6, & 4H fm2,5X,14HTotal energy =,F12.5,4H MeV,/,79(1H-)) 111 FORMAT(A) 112 FORMAT(3(F10.4),4H => ,F10.4,F17.5,I10) 113 FORMAT(23HCase of minimal energy:,/,3Hb =,F8.4,4X,3Hq =,F10.4, &7X,3HE =,F17.5,/,4Hmu =,F10.4,4H => ,4HQ20=,F10.4,I5,7H iters.) 114 FORMAT(18HEnd of iteration :,I3,/, &9HRcharge =,F10.2,3H fm,5X,5HQ20 =,F11.2,3H fm, &5X,6HETOT =,F12.5,4H MeV,/,79(1H-)) 115 FORMAT(18HEnd of iteration :,I3,/, &9HRcharge =,F10.2,3H fm,5X,5HQ20 =,F11.2,3H fm, &5X,6HEHTDA=,F12.5,4H MeV,/,79(1H-)) 116 FORMAT(I3,F11.3,F9.2,F7.3) 117 FORMAT(I3,3(F9.4),3(F11.3),F9.2,F7.3) 118 FORMAT(I3,14(F11.4)) C END C============================================================================== C SUBROUTINES for SIMPLEX calculations C============================================================================== C This is used by simplex for enter the data and return the results C PARAM : input data C Chi2 : output C ICTRL : = 1 b,q from simplex method C = 0 b,q from file data or meshsimplex C============================================================================== SUBROUTINE FCN(PARAM,CHI2) INCLUDE 'All.f' DIMENSION PARAM(10) C BZRO = PARAM(1) QZRO = PARAM(2) CALL CALCUL(BZRO,QZRO,CHI2,ENEW,Q0NEW,Q40CA) IF (ICTRL.EQ.1) THEN WRITE (18,100) BZRO,QZRO,Q0NEW/1D2,Q40CA/1D4,ENEW,ITHF+ITHTDA ENDIF C 100 FORMAT(2Hb=,F7.4,1X,2Hq=,F7.4,1X,4HQ20=,F8.2,1Hb,1X, &4HQ40=,F10.4,2Hb2,1X,2HE=,F11.4,I4,7H iters.) END C============================================================================== C Initialisation for simplex method C============================================================================== SUBROUTINE INITS INCLUDE 'All.f' CHARACTER*10 PNAME(10) CHARACTER*78 TITLE C C--------------------------------- C Print parameters starting values C--------------------------------- FMIN = 1D20 OPEN(UNIT=18,FILE='simplexres') OPEN(UNIT=30,FILE = 'fixedata/simplexdat',STATUS='OLD',ERR=900) GOTO 901 900 PRINT *,'Error opening file "simplexdat" in directory "fixedata"' PRINT *,'File is not existed' PRINT *,'' WRITE(20,*) & 'Error opening file "simplexdat" in directory "fixedata"' WRITE(20,*) 'File is not existed' WRITE(20,*) CLOSE (20) STOP 901 READ(30,105) TITLE,TITLE WRITE (18,100) TITLE READ(30,105) TITLE READ(30,*) NSP WRITE(18,101) NSP N = NSP + 1 READ(30,105) TITLE DO I = 1, NSP READ(30,106) PNAME(I) READ(30,*) SIMP(1,I),STEP(I),BINF(I),BSUP(I),GMAXER(I) ENDDO READ(30,105) TITLE READ(30,*) MAXTER READ(30,105) TITLE READ(30,*) GMAXER(NSP+1) CLOSE(30) C DO I = 1, NSP WRITE (18,102) PNAME(I),SIMP(1,I),STEP(I),BINF(I),BSUP(I), & GMAXER(I) ENDDO WRITE (18,103) GMAXER(NSP+1) C 100 FORMAT(/,80(1H*),/,1H*,19X, &39HS I M P L E X M I N I M I S A T I O N,20X,1H*,/,1H*, &78X,1H*,/,1H*,A78,1H*,/,80(1H*),//) 101 FORMAT(22HNumber of parameters :,I3,//, &38H Parameter Value Pas , &39H Inf.B. Sup.B. Tolerance,/) 102 FORMAT(A10,5(3X,F10.5)) 103 FORMAT(/,28HTolerance on the Function : ,F12.5,///, & 21X,37HC O N T R O L E O F S I M P L E X/) 105 FORMAT(A78) 106 FORMAT(A10) END C============================================================================== C Gives LOW/HIGH in each parameter C============================================================================== SUBROUTINE ORDER(LOW,HIGH) INCLUDE 'All.f' INTEGER HIGH DIMENSION HIGH(10),LOW(10) C N = NSP+1 DO J = 1, N DO I = 1, N IF (SIMP(I,J).LT.SIMP(LOW(J),J)) LOW(J) = I IF (SIMP(I,J).GT.SIMP(HIGH(J),J)) HIGH(J) = I ENDDO ENDDO C END C============================================================================== C Controle the bornes of variation of parameters C============================================================================== SUBROUTINE BKTRL(AA) INCLUDE 'All.f' DIMENSION AA(10) C DO I = 1, N IF (AA(I).GT.BSUP(I)) AA(I) = BSUP(I) IF (AA(I).LT.BINF(I)) AA(I) = BINF(I) ENDDO C END C============================================================================== C VXNEXT in place of the worst vertex C============================================================================== SUBROUTINE NWVRTX(NITER,VXNEXT,HIGH) INCLUDE 'All.f' INTEGER HIGH DIMENSION VXNEXT(10),HIGH(10) C IUNIT = 50 DO I = 1, N SIMP(HIGH(N),I) = VXNEXT(I) ENDDO IF (VXNEXT(N).LT.FMIN/5D0) THEN FMIN = VXNEXT(N) WRITE (IUNIT,100) NITER,FMIN C WRITE (IUNIT,110) (NAME(I),VXNEXT(I),I=1,M) ENDIF C 100 FORMAT(1H ,'ITERATION:',I4,' CHI2 = ',G12.5) C110 FORMAT(1H ,36X,A10,' = ',G12.5) END C============================================================================= C Report results of simplex C============================================================================= SUBROUTINE REPORT(NITER,VLMEAN,ERROR) INCLUDE 'All.f' DIMENSION VLMEAN(10),ERROR(10) C WRITE (18,100) NITER DO J = 1, N WRITE (18,101) (SIMP(J,I), I = 1, N) ENDDO WRITE (18,102) WRITE (18,101) (VLMEAN(I), I = 1, N) WRITE (18,103) WRITE (18,101) (ERROR(I), I = 1, N) C 100 FORMAT(/,26X,27HS I M P L E X R E S U L T,///, &32HMinimisation has converged after, I4,11H iterations,//, &22HThe final simplex is :) 101 FORMAT(5F15.8) 102 FORMAT(/,13HThe mean is :) 103 FORMAT(/,35HThe estimated fractional error is :) END c============================================================================== PROGRAM AxialCalculation INCLUDE 'All.f' DIMENSION PARAM(10) LOGICAL DONE INTEGER HIGH DIMENSION HIGH(10),LOW(10),VXNEXT(10),CENTER(10),VLMEAN(10) DIMENSION PS(10),QS(10),ERROR(10) DATA ALFA,BETA,GAMMA/1D0,5D-1,2D0/ C OPEN(UNIT=20,FILE='result') C CALL MEMCOM C INITIALITION OF CONSTANTS : Pi,4Pi,eV,2**x,n!,options messages CALL INICTS C READ INPUT DATA CALL REDATA C---------------- IF (ICTRL.EQ.0) THEN ! No simplex C In this case, read B0 and Q from "data" file or calculate them C from the "mesh" file PARAM(1) = B0 PARAM(2) = Q CALL FCN(PARAM,CHI2) ELSEIF (ICTRL.EQ.1) THEN ! Simplex ROOT2 = DSQRT(2D0) CALL INITS C Compute the first vertex DO I = 1, NSP PARAM(I) = SIMP(1,I) ENDDO CALL FCN(PARAM,CHI2) N = NSP+1 SIMP(1,N) = CHI2 C Compute offset of the vertexes of the starting simplex DO I = 1, NSP SQRN = DSQRT(DBLE(N)) PS(I) = STEP(I)*(SQRN+DBLE(NSP)-1D0)/(DBLE(NSP)*ROOT2) QS(I) = STEP(I)*(SQRN-1D0)/(DBLE(NSP)*ROOT2) ENDDO C All vertexes of the starting simplex DO I = 2, N DO J = 1, NSP SIMP(I,J) = SIMP(1,J) + QS(J) ENDDO SIMP(I,I-1) = SIMP(1,I-1) + PS(I-1) DO K = 1, NSP PARAM(K) = SIMP(I,K) ENDDO CALL FCN(PARAM,CHI2) SIMP(I,N) = CHI2 ENDDO C Preset the vertex (Chi2 value) DO I = 1, N LOW(I) = 1 HIGH(I) = 1 ENDDO C-------------------------- CALL ORDER(LOW,HIGH) C-------------------------- C STARTING ITERATION NITER = 0 999 DONE = .TRUE. NITER = NITER + 1 DO I = 1, N CENTER(I) = 0D0 ENDDO C Compute centroid excluding the worst DO I = 1, N IF (I.NE.HIGH(N)) THEN DO J = 1, NSP CENTER(J) = CENTER(J) + SIMP(I,J) ENDDO ENDIF ENDDO C First attempt to reflect the worst (in VXNEXT) DO I = 1, N CENTER(I) = CENTER(I)/NSP VXNEXT(I) = (1D0+ALFA)*CENTER(I) - ALFA*SIMP(HIGH(N),I) ENDDO CALL BKTRL(VXNEXT) DO K = 1, NSP PARAM(K) = VXNEXT(K) ENDDO CALL FCN(PARAM,CHI2) VXNEXT(N) = CHI2 C C MAIN TEST LOOP C C Better than the best? IF (VXNEXT(N).LE.SIMP(LOW(N),N) ) THEN C OK Accepted CALL NWVRTX(NITER,VXNEXT,HIGH) C VXNEXT expanded DO I = 1, NSP VXNEXT(I) = GAMMA*SIMP(HIGH(N),I)+(1D0-GAMMA)*CENTER(I) ENDDO CALL BKTRL(VXNEXT) DO K = 1, NSP PARAM(K) = VXNEXT(K) ENDDO CALL FCN(PARAM,CHI2) VXNEXT(N) = CHI2 C Still better IF (VXNEXT(N).LE.SIMP(LOW(N),N)) THEN C Expansion accepted CALL NWVRTX(NITER,VXNEXT,HIGH) ENDIF C If not better than the best ELSE IF (VXNEXT(N).LE.SIMP(HIGH(N),N) ) THEN C Better than the worst CALL NWVRTX(NITER,VXNEXT,HIGH) C Worst than worst ELSE C Then contract DO I = 1, NSP VXNEXT(I) = BETA*SIMP(HIGH(N),I)+(1D0-BETA)*CENTER(I) ENDDO CALL BKTRL(VXNEXT) DO K = 1, NSP PARAM(K) = VXNEXT(K) ENDDO CALL FCN(PARAM,CHI2) VXNEXT(N) = CHI2 IF (VXNEXT(N).LE.SIMP(HIGH(N),N) ) THEN C Contraction accepted CALL NWVRTX(NITER,VXNEXT,HIGH) C If still bad ELSE C Shrink all bad vertexes DO I = 1, N DO J = 1, NSP SIMP(I,J) = (SIMP(I,J)+SIMP(LOW(N),J))*BETA ENDDO DO K = 1, NSP PARAM(K) = SIMP(I,K) ENDDO CALL FCN(PARAM,CHI2) SIMP(I,N) = CHI2 ENDDO ENDIF ENDIF ENDIF CALL ORDER(LOW,HIGH) C Check for convergence DO J = 1, N ERROR(J) = (SIMP(HIGH(J),J)-SIMP(LOW(J),J))/SIMP(HIGH(J),J) IF (DONE) THEN IF (ERROR(J).GT.GMAXER(J)) THEN DONE = .FALSE. ENDIF ENDIF ENDDO C Test for end of iteration IF ((.NOT.DONE).AND.(NITER.LT.MAXTER)) GOTO 999 C Average each parameter DO I = 1, N VLMEAN(I) = 0D0 DO J = 1, N VLMEAN(I) = VLMEAN(I) + SIMP(J,I) ENDDO VLMEAN(I) = VLMEAN(I) / N ENDDO C Print simplex minimisation results CALL REPORT(NITER,VLMEAN,ERROR) C Print user's results CALL FCN(VLMEAN,CHI2) CLOSE (18) ELSE WRITE(20,*)'Error in simplex option' CLOSE (20) STOP 'Error in simplex option' ENDIF C CLOSE(20) C END C============================================================================== C C Lanczos diagonalisation method C C The folowing two subroutine are used like to "communicate" with C Lanczos others subroutines and functions. C C============================================================================== C C======================================================================= C This is the matrix multiply for the matrix A C======================================================================= SUBROUTINE OP(N,M,P,QQ) INCLUDE 'All.f' DOUBLE PRECISION P(N,M),QQ(N,M),VABL(NE) DIMENSION J(NE) COMMON/NFIL/NFILE C DO II = 1, N DO JJ = 1, M QQ(II,JJ) = 0D0 ENDDO ENDDO REWIND (NFILE) DO LINE = 1, N C Read the matrix elements on line LINE READ(NFILE) KMAX, (VABL(IK),J(IK),IK = 1, KMAX) IF (J(1).EQ.LINE) THEN DIAG = 5D-1 ELSE DIAG = 1D0 ENDIF DO IK = 1, KMAX DO IM = 1, M QQ(LINE,M) = QQ(LINE,M) + VABL(IK)*P(J(IK),M)*DIAG QQ(J(IK),M) = QQ(J(IK),M) + VABL(IK)*P(LINE,M)*DIAG ENDDO DIAG = 1D0 ENDDO ENDDO C END C======================================================================= C This subroutine stores and recalls vectors. C======================================================================= SUBROUTINE IOVECT(N,M,Q,J,K) INTEGER N,M,J,K,I,L DOUBLE PRECISION Q DIMENSION Q(N,M) C IF (J.EQ.M) REWIND(28) IF (K.EQ.0) WRITE(28) ((Q(I,L),I = 1, N),L = 1, M) IF (K.EQ.1) READ(28) ((Q(I,L),I = 1, N),L = 1, M) C END C****************************************************************************** C****************************************************************************** C LE CALCUL RPA C****************************************************************************** C****************************************************************************** C======================================================================= C Determine the many-body basis for HTDA calculation C======================================================================= SUBROUTINE NEWMBBASE INCLUDE 'All.f' dimension ict(4),jct(4),kct(4) C INTEGER TIME(3) C NTOT(inp) = 0 CALL ITIME(TIME) WRITE(20,100) TIME PRINT 100,TIME C DO I = 1, NHPMAX NALL(I,INP) = 0 NSS(INP)=0 NFOUR(INP)=0 ENDDO C Find all the configuration of order 0, 1 and 2 CALL CONFIG(NUK2(INP),NUIP(INP),0,0,ECHTDA) if (kopair.ne.1) then CALL CONFIG(NUK2(INP),NUIP(INP),1,1,ECHTDA+400.) end if CALL CONFIG(NUK2(INP),NUIP(INP),2,2,ECHTDA) CALL CONFIG(NUK2(INP),NUIP(INP),4,4,ECHTDA) c CALL CONFIG(NUK2(INP),NUIP(INP),6,6,ECHTDA) c IF (INP.EQ.1) THEN DO I=1,NHPMAX DO J=1,2 DO K=1,NE MBST(I,J,K)=MBSTN(I,J,K) END DO END DO END DO ELSE DO I=1,NHPMAX DO J=1,2 DO K=1,NE MBST(I,J,K)=MBSTP(I,J,K) END DO END DO END DO END IF c C If one want to stop taking the configuration of order more than 2 C choose the command below GOTO 3 C Find configuration whose Eph is minimal CALL VMIN(NTOT(inp),EPH,IMINPH) C Put this configuration in the common /MPHC/ DO I = 1, NHPMAX IF (MBST(I,1,IMINPH).NE.0) THEN IHC(I) = MBST(I,1,IMINPH) ELSE GOTO 1 ENDIF ENDDO c 1 NHC = I - 1 DO I = 1, NHPMAX IF (MBST(I,2,IMINPH).NE.0) THEN IPC(I) = MBST(I,2,IMINPH) ELSE GOTO 2 ENDIF ENDDO 2 NPC = I - 1 C Add bases based on this configuration NCOM(inp) = NTOT(inp) NCOM3 = NCOM(inp) DO IO = 3, 4 C Decide type of configurations that we have to add to the existed one IF (NHC.EQ.1) THEN IF (NPC.EQ.1) THEN C Ordre 1=> add ordre 3 IF (IO.EQ.3) & CALL CONFIG(NUK2(INP),NUIP(INP),3,3,ECHTDA) ENDIF ELSEIF (NHC.EQ.2) THEN IF (NPC.EQ.2) THEN C Ordre 2=> add ordre 3, 4 IF (IO.EQ.3) & CALL CONFIG(NUK2(INP),NUIP(INP),3,3,ECHTDA) ENDIF ENDIF IF (NCOM(inp).NE.NTOT(inp)) THEN WRITE(20,102) IO+NHC-MIN(NHC,NPC),IO+NPC-MIN(NHC,NPC), & NTOT(inp)-NCOM(inp) PRINT 102, IO+NHC-MIN(NHC,NPC),IO+NPC-MIN(NHC,NPC), & NTOT(inp)-NCOM(inp) ENDIF NCOM(inp) = NTOT(inp) ENDDO C label used if one need to stop to take the configuration of C order more than 2 3 CONTINUE NCOM(inp) = NCOM3 PRINT 101, NTOT(inp) WRITE(20,101) NTOT(inp) 100 FORMAT(4X,42H*Determine the many-body basis {|npmh>} at,19X, &I2.2,1H:,I2.2,1X,I2.2,1Hs) 101 FORMAT(6X,38(1H-),/,6X,27HSize of the many-body basis,I11) 102 FORMAT(6X,25HNumber of configuration |,I1,1Hh,I1,3Hp>:,I7) RETURN 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 CONFIG(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 NCONFS = 0 NSS(inp) = 0 NPS = 2*NPL(inp) NHS = 2*NHL(inp) cc-->> NPSD = 2*NPLD(inp) NHSD = 2*NHLD(inp) cc-->> 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 |0p-0h> C======================================================================= IF (NHOLE.EQ.0) THEN IF (NPAR.EQ.0) THEN IF (KCONF2.EQ.0) THEN IF (IPCONF.EQ.1) THEN NALL(1,inp) = NALL(1,inp) + 1 NSS(inp) = NSS(inp) + 1 NTOT(inp) = NTOT(inp) + 1 IF (NTOT(inp).GT.NE) GOTO 999 DO I = 1, NHPMAX DO J = 1, 2 MBSTP(I,J,NTOT(inp)) = 0 MBSTN(I,J,NTOT(inp)) = 0 ENDDO ENDDO EPH(NTOT(inp),inp) = 0D0 cc-->> mbreversed(ntot(inp),inp)=ntot(inp) cc<<-- ENDIF ENDIF GOTO 5 ENDIF ENDIF 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=ntot(inp)+1 if (mod(ict1,2).eq.0) mbmax=ntot(inp) 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(2,inp) = NALL(2,inp) + 1 EEXCI = SSE(ICT2,2,inp) - SSE(ICT1,1,inp) IF (EEXCI.LE.CUTOFF+EREF) THEN NSS(inp) = NSS(inp) + 1 NTOT(inp) = NTOT(inp) + 1 IF (NTOT(inp).GT.NE) GOTO 999 DO I = 1, NHPMAX DO J = 1, 2 if (inp.eq.1) then MBSTN(I,J,NTOT(inp)) = 0 else MBSTP(i,j,ntot(inp)) = 0 end if ENDDO ENDDO if (inp.eq.1) then MBSTn(1,1,NTOT(inp)) = ICT1 MBSTn(1,2,NTOT(inp)) = ICT2 else mbstp(1,1,ntot(inp)) = ict1 mbstp(1,2,ntot(inp)) = ict2 end if EPH(NTOT(inp),inp) = 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 mbreversed(ntot(inp),inp)=isig*mbrev mbreversed(mbrev,inp)=isig*ntot(inp) mbrev=mbmax c write(60,*)ict1,ict2,ntot(inp),mbrev, c & mbreversed(ntot(inp),inp),mbreversed(mbrev,inp) endif else if (mbstp(1,2,mbrev).eq.jct2) then mbreversed(ntot(inp),inp)=isig*mbrev mbreversed(mbrev,inp)=isig*ntot(inp) 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 5 ENDIF !1p ENDIF !1h 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=ntot(inp)+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(3,inp) = NALL(3,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 NTOT(inp) = NTOT(inp) + 1 IF (NTOT(inp).GT.NE) GOTO 999 DO I = 1, NHPMAX DO J = 1, 2 if (inp.eq.1) then MBSTn(I,J,NTOT(inp)) = 0 else mbstp(i,j,ntot(inp)) = 0 end if ENDDO ENDDO if (inp.eq.1) then MBSTN(2,1,NTOT(INP)) = ICT1 MBSTN(1,1,NTOT(INP)) = ICT2 MBSTN(2,2,NTOT(INP)) = ICT3 MBSTN(1,2,NTOT(INP)) = ICT4 ELSE MBSTP(2,1,NTOT(INP)) = ICT1 MBSTP(1,1,NTOT(INP)) = ICT2 MBSTP(2,2,NTOT(INP)) = ICT3 MBSTP(1,2,NTOT(INP)) = ICT4 END IF EPH(NTOT(inp),inp) = EEXCI c 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 mbreversed(ntot(inp),inp)=ntot(inp) ! 1paire_h 1paire_p else mbreversed(ntot(inp),inp)=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 mbreversed(ntot(inp),inp)=isig*mbrev mbreversed(mbrev,inp)=isig*ntot(inp) 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 mbreversed(ntot(inp),inp)=isig*mbrev mbreversed(mbrev,inp)=isig*ntot(inp) mbrev=mbmax endif mbrev=mbrev+1 enddo end if endif cc<<-- ENDIF !eexci 2 CONTINUE ENDIF !spin ENDIF !parity ENDDO !ict4 ENDDO !ict3 cc-->> if (inp.eq.1) then if ((mod(ict1,2).eq.1).and.(ict2.eq.ict1+1)) then do mb1=mbmin,ntot(inp) if (mbreversed(mb1,inp).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.ntot(inp)) if ((mbstn(1,2,mbrev).eq.jct3).and. # (mbstn(2,2,mbrev).eq.jct4)) then mbreversed(mb1,inp)=isig*mbrev mbreversed(mbrev,inp)=isig*mb1 mbrev=ntot(inp) 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,ntot(inp) if (mbreversed(mb1,inp).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.ntot(inp)) if ((mbstp(1,2,mbrev).eq.jct3).and. # (mbstp(2,2,mbrev).eq.jct4)) then mbreversed(mb1,inp)=isig*mbrev mbreversed(mbrev,inp)=isig*mb1 mbrev=ntot(inp) endif mbrev=mbrev+1 enddo endif enddo endif end if c write(60,*)ntot(inp),ict1,ict2,ict3,ict4,mbreversed(mb1,inp) cc<<-- ENDDO !ict2 cc-->> if (mod(ict1,2).eq.1) then mbmax=ntot(inp) endif cc<<-- ENDDO ! ict1 GOTO 5 ENDIF !2h ENDIF !2p 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 NTOT(inp) = NTOT(inp) + 1 IF (NTOT(inp).GT.NE) GOTO 999 DO I = 1, NHPMAX DO J = 1, 2 if (inp.eq.1) then MBSTn(I,J,NTOT(inp)) = 0 else mbstp(i,j,ntot(inp)) = 0 end if ENDDO ENDDO if (inp.eq.1) then MBSTN(2,1,NTOT(INP)) = ICT1 MBSTN(1,1,NTOT(INP)) = ICT2 MBSTN(2,2,NTOT(INP)) = ICT3 MBSTN(1,2,NTOT(INP)) = ICT4 mbstn(3,1,ntot(inp)) = ict5 mbstn(4,1,ntot(inp)) = ict6 mbstn(3,2,ntot(inp)) = ict7 mbstn(4,2,ntot(inp)) = ict8 ELSE MBSTP(2,1,NTOT(INP)) = ICT1 MBSTP(1,1,NTOT(INP)) = ICT2 MBSTP(2,2,NTOT(INP)) = ICT3 MBSTP(1,2,NTOT(INP)) = ICT4 mbstp(3,1,ntot(inp)) = ict5 mbstp(4,1,ntot(inp)) = ict6 mbstp(3,2,ntot(inp)) = ict7 mbstp(4,2,ntot(inp)) = ict8 END IF C EPH(NTOT(inp),inp) = EEXCI c mbreversed(ntot(inp),inp)=ntot(inp) 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 5 ENDIF !2h ENDIF !2p C======================================================================= C |3p-3h> C !!!! uniquement dans la petite fenetre !!!! C======================================================================= IF (NHOLE.EQ.3) THEN IF (NPAR.EQ.3) THEN IF (NHC.EQ.1) THEN C C Based on |1p-1h> (I1=>J1) C KKNOW = SS2K(J1,2,inp) - SS2K(I1,1,inp) IPKNOW = SSIP(J1,2,inp) * SSIP(I1,1,inp) EKNOW = SSE(J1,2,inp) - SSE(I1,1,inp) cc-->> DO ICT1 = NHS-NHSD+1, NHS - 1 IF (ICT1.EQ.I1) GOTO 10 DO ICT2 = ICT1 + 1, NHS IF (ICT2.EQ.I1) GOTO 11 KT12 = SS2K(ICT1,1,inp) + SS2K(ICT2,1,inp) DO ICT3 = 1, NPSD - 1 IF (ICT3.EQ.J1) GOTO 12 DO ICT4 = ICT3 + 1, NPSD cc-->> IF (ICT4.EQ.J1) GOTO 13 KT34 = SS2K(ICT3,2,inp) + SS2K(ICT4,2,inp) IF (KT34-KT12+KKNOW.EQ.KCONF2) THEN IP12 = SSIP(ICT1,1,inp) * SSIP(ICT2,1,inp) IP34 = SSIP(ICT3,2,inp) * SSIP(ICT4,2,inp) IF (IP12*IP34*IPKNOW.EQ.IPCONF) THEN EEXCI = EKNOW + & SSE(ICT3,2,inp) + SSE(ICT4,2,inp) & - SSE(ICT1,1,inp) - SSE(ICT2,1,inp) IF (EEXCI.LE.CUTOFF+EREF) THEN NCONFS = NCONFS + 1 NSS(inp) = NSS(inp) + 1 NTOT(inp) = NTOT(inp) + 1 IF (NTOT(inp).GT.NE) GOTO 999 DO I = 1, NHPMAX DO J = 1, 2 if (inp.eq.1) then MBSTn(I,J,NTOT(inp)) = 0 else mbstp(i,j,ntot(inp)) = 0 end if ENDDO ENDDO CALL ARRANG(ICT1,ICT2,I1,0,0) if (inp.eq.1) then MBSTn(3,1,NTOT(inp)) = NAFAR(3) MBSTn(2,1,NTOT(inp)) = NAFAR(2) MBSTn(1,1,NTOT(inp)) = NAFAR(1) else mbstp(3,1,ntot(inp)) = nafar(3) mbstp(2,1,ntot(inp)) = nafar(2) mbstp(1,1,ntot(inp)) = nafar(1) end if CALL ARRANG(ICT3,ICT4,J1,0,0) if (inp.eq.1) then MBSTn(3,2,NTOT(inp)) = NAFAR(3) MBSTn(2,2,NTOT(inp)) = NAFAR(2) MBSTn(1,2,NTOT(inp)) = NAFAR(1) else mbstp(3,2,ntot(inp)) = nafar(3) mbstp(2,2,ntot(inp)) = nafar(2) mbstp(1,1,ntot(inp)) = nafar(1) end if EPH(NTOT(inp),inp) = EEXCI CALL TESTCO(NEWCON) IF (.NOT.NEWCON) THEN NSS(inp) = NSS(inp) - 1 NTOT(inp) = NTOT(inp) - 1 ENDIF ENDIF ENDIF ENDIF 13 CONTINUE ENDDO 12 CONTINUE ENDDO 11 CONTINUE ENDDO 10 CONTINUE ENDDO GOTO 5 ELSEIF (NHC.EQ.2) THEN C C Based on |2p-2h> (I1,I2=>J1,J2) C KKNOW = SS2K(J1,2,inp) + SS2K(J2,2,inp) & - SS2K(I1,1,inp) - SS2K(I2,1,inp) IPKNOW = SSIP(J1,2,inp) * SSIP(J2,2,inp) & * SSIP(I1,1,inp) * SSIP(I2,1,inp) EKNOW = SSE(J1,2,inp) + SSE(J2,2,inp) & - SSE(I1,1,inp) - SSE(I2,1,inp) cc-->> DO ICT1 = NHS-NHSD+1, NHS IF (ICT1.EQ.I1.OR.ICT1.EQ.I2) GOTO 14 KT1 = SS2K(ICT1,1,inp) DO ICT2 = 1, NPSD cc-->> IF (ICT2.EQ.J1.OR.ICT2.EQ.J2) GOTO 15 KT2 = SS2K(ICT2,2,inp) IF (KT2-KT1+KKNOW.EQ.KCONF2) THEN IP1 = SSIP(ICT1,1,inp) IP2 = SSIP(ICT2,2,inp) IF (IP1*IP2*IPKNOW.EQ.IPCONF) THEN EEXCI = EKNOW + & SSE(ICT2,2,inp) - SSE(ICT1,1,inp) IF (EEXCI.LE.CUTOFF+EREF) THEN NCONFS = NCONFS + 1 NSS(inp) = NSS(inp) + 1 NTOT(inp) = NTOT(inp) + 1 IF (NTOT(inp).GT.NE) GOTO 999 DO I = 1, NHPMAX DO J = 1, 2 if (inp.eq.1) then MBSTn(I,J,NTOT(inp)) = 0 else mbstp(i,j,ntot(inp)) = 0 end if ENDDO ENDDO CALL ARRANG(ICT1,I1,I2,0,0) if (inp.eq.1) then MBSTn(3,1,NTOT(inp)) = NAFAR(3) MBSTn(2,1,NTOT(inp)) = NAFAR(2) MBSTn(1,1,NTOT(inp)) = NAFAR(1) else mbstp(3,1,ntot(inp)) = nafar(3) mbstp(2,1,ntot(inp)) = nafar(2) mbstp(1,1,ntot(inp)) = nafar(1) end if CALL ARRANG(ICT2,J1,J2,0,0) if (inp.eq.1) then MBSTn(3,2,NTOT(inp)) = NAFAR(3) MBSTn(2,2,NTOT(inp)) = NAFAR(2) MBSTn(1,2,NTOT(inp)) = NAFAR(1) else mbstp(3,2,ntot(inp)) = nafar(3) mbstp(2,2,ntot(inp)) = nafar(2) mbstp(1,1,ntot(inp)) = nafar(1) end if EPH(NTOT(inp),inp) = EEXCI CALL TESTCO(NEWCON) IF (.NOT.NEWCON) THEN NSS(inp) = NSS(inp) - 1 NTOT(inp) = NTOT(inp) - 1 ENDIF ENDIF ENDIF ENDIF 15 CONTINUE ENDDO 14 CONTINUE ENDDO GOTO 5 ENDIF ENDIF ENDIF 5 IF (NHOLE+NPAR.LT.6) THEN WRITE(20,100) & NHOLE,NPAR,NSS(inp),NALL(MIN0(NPAR,NHOLE)+1,inp) C***** c WRITE(40,*) NHOLE,NPAR,NSS C***** PRINT 100, NHOLE,NPAR,NSS(inp),NALL(MIN0(NPAR,NHOLE)+1,inp) WRITE(20,100) NHOLE,NPAR,NSS(inp),NCONFS c WRITE(20,*) ' (based on [', c & (IHC(I), I = 1, NHC),'=>',(IPC(I), I = 1, NPC),'])' c PRINT 100, NHOLE,NPAR,NSS(inp),NCONFS c PRINT *,' (based on [', c & (IHC(I), I = 1, NHC),'=>',(IPC(I), I = 1, NPC),'])' ENDIF if (nhole+npar.eq.8) then write(20,100) nhole,npar, nfour(inp),nall(4,inp) print 100, nhole,npar,nfour(inp),nall(4,inp) end if c RETURN 999 WRITE(20,*) 'Number of configuration is too big',NTOT(inp),NE CLOSE(20) STOP 'Number of configuration is too big' C 100 FORMAT(6X,25HNumber of configuration |,I1,1Hh,I1,3Hp>:,I7,1H/,I8) C END C C======================================================================= C Determine the matrix for HTDA calculation+tr correction C======================================================================= SUBROUTINE HTDAMA 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 elementHTDAMA 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 NT1=1 NPLUS=0 DO I=1,NTOT(INP) IBAR=MBREVERSED(I,INP) IF (I.LE.IABS(IBAR)) THEN ! PHIPLUS(NT1,1)=I ! state |Psi> IF(IBAR.EQ.I) THEN PHIPLUS(NT1,2)=0 !state + NPLUS=NPLUS+1 ELSE PHIPLUS(NT1,2)=IBAR ! state - END IF NT1=NT1+1 END IF END DO NTOT1(INP)=NT1-1 print*, 'ntot1=', ntot1(inp) C CALL ITIME(TIME) WRITE(20,100) TIME PRINT 100,TIME C C NDIFF1 = 0 NDIF0H = 0 NDIF0P = 0 cc<<-- C C EIDEN: Identity Matrice * EIDEN c DO I = 1, NTOT1(inp) 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, NTOT1(inp) 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=RINEW(K,L) ! c IF(L1.NE.0) E2=RINEW(K,L1) ! c IF(K1.NE.0) E3=RINEW(K1,L) !<\bar psi|V|psi> c IF(K1.NE.0.AND.L1.NE.0) E4=RINEW(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 IF (I.EQ.J) THEN c DIAGE = EMTMP + EPH(K,INP) + 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 c ENDDO c IF (INP.EQ.1) INF=29 c IF (INP.EQ.2) INF=39 c WRITE(INF) KMAX, (VABL(J),ICOL(J), J = 1, KMAX) c ENDDO 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, NTOT1(inp) C C |i+>=|k>+\bar|k> C K =PHIPLUS(I,1) KBAR=PHIPLUS(I,2) K1 =IABS(KBAR) C DIAGEI(I)=EPH(K,INP)+RINEW(K,K) c KMAX = 0 cc DO J = I+1, NTOT1(inp) 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=RINEW(K,L) ! IF(L1.NE.0) E2=RINEW(K,L1) ! IF(K1.NE.0) E3=RINEW(K1,L) !<\bar psi|V|psi> IF(K1.NE.0.AND.L1.NE.0) E4=RINEW(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,NTOT1(INP)) CLOSE(33) CLOSE(35) C PRINT 101, EPHTDA WRITE(20,101) EPHTDA PRINT 102, NDIFF1,NDIF0H,NDIF0P WRITE(20,102) NDIFF1,NDIF0H,NDIF0P REWIND(25) REWIND(26) REWIND(27) REWIND(45) REWIND(46) REWIND(47) c if (inp.eq.1) then WRITE(24) NDIFF1,NDIF0H,NDIF0P DO I=1,NDIFF1 READ (25) N1,N2,N3,N4,N5 WRITE(24) N1,N2,N3,N4,N5 ENDDO DO I=1,NDIF0H READ (26) N1,N2 WRITE(24) N1,N2 ENDDO DO I=1,NDIF0P READ (27) N1,N2 WRITE(24) N1,N2 ENDDO end if if (inp.eq.2) then WRITE(44) NDIFF1,NDIF0H,NDIF0P DO I=1,NDIFF1 READ (45) N1,N2,N3,N4,N5 WRITE(44) N1,N2,N3,N4,N5 ENDDO DO I=1,NDIF0H READ (46) N1,N2 WRITE(44) N1,N2 ENDDO DO I=1,NDIF0P READ (47) N1,N2 WRITE(44) N1,N2 ENDDO end if 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