C####################################################################### C########## ################## C########## PROGRAM GAMPN ################## C########## ################## C########## VERSION: JULY, 1991 ################## C########## ################## C########## RECEIVED AT ANU FROM PAUL SEMMES ################## C########## 4-MAR-1993 A.E. STUCHBERY ################## C########## ################## C########## PC version AES May 2007 ################## C###### I have replaced time with elapsed time ############## C########## ################## C########## ################## C####################################################################### C*** ORIGINAL PROGRAM DESCRIBED IN NUCL PHYS A307 (1978) 189 **** C AND PHYSICA SCRIPTA 8 (1973) 17 * C * C PLEASE REPORT ERRORS, INCONSISTENCIES OR PROBLEMS TO: * C * C INGEMAR RAGNARSSON PAUL SEMMES * C DEPT OF MATH PHYSICS OR PHYSICS DEPT, BOX 5051 * C LTH, BOX 118 TENESSEE TECH. UNIV. * C S-22100 LUND, SWEDEN COOKEVILLE, TN 38505 USA * C * C * C*********************************************************************** C C VERSION: JULY, 1991 C C DISTRIBUTED AT THE OAK RIDGE THEORY WORKSHOP, AUG 1991, FOR USE C WITH OTHER MODIFIED-OSCILLATOR CODES ASYRMO (DIAGONALIZES THE C PARTICLE + TRIAXIAL ROTOR HAMILTONIAN) AND PROBAMO (CALCULATES C M1/E2 MATRIX ELEMENTS FROM THE ENERGY EIGENVECTORS OF ASYRMO) C C C GAMPN CALCULATES OF SINGLE-PARTICLE ENERGIES AND VARIOUS MATRIX C ELEMENTS IN A MODIFIED OSCILLATOR (NILSSON) POTENTIAL. C PROGRAM MAINLY DESIGNED TO PROVIDE INPUT FOR SUBSEQUENT ODD-A C AND ODD-ODD PARTICLE-ROTOR CALCULATIONS C C SINGLE-PARTICLE TRANSFORMATION: /N l j omega> basis ---> C /N l lambda sigma> basis WRITTEN ON FILE2 (NEEDED FOR ODD-ODD C CALCULATIONS WITH A PROTON-NEUTRON RESIDUAL INTERACTION) C C DATA FOR VMI CALCULATIONS (IN ASYRMO/ASYRPN) WRITTEN ON FILE16 C DATA FOR PART-ROTOR PROGRAMS (ASYRPN AND PROBA) WRITTEN ON FILE17 C C ALSO POSSIBLE TO CALCULATE ENERGIES IN A CRANKED POTENTIAL C (OMROT .NE. 0) C C*********************************************************************** C C DESCRIPTION OF THE INPUT DATA (FREE FORMAT UNLESS SPECIFIED): C C CARD 1. FILE2, FILE16, FILE17 C OUTPUT FILES TO BE USED FOR SUBSEQUENT APPLICATIONS C (RECALL THAT ON A VAX FREE FORMAT CHARACTER C DATA MUST BE ENCLOSED WITH APOSTROPHES) C FILE2 CONTAINS THE S.P. WAVE FUNCTIONS IN THE C |N,L,LAMBDA,SIGMA> BASIS NEEDED FOR Vpn CALCULATIONS C FILE16 CONTAINS THE S.P. WAVEFUNCTIONS NEEDED FOR THE VMI C OPTION IN ASYR C FILE17 CONTAINS THE USUAL OUTPUT NEEDED FOR ASYR, IE, C SINGLE PARTICLE MATRIX ELEMENTS LIKE , , etc C C 2. ISTRCH, ICORR, IREC C ISTRCH=0: CALCULATION IN SPHERICAL COORDINATES C (DEF PARAM DELTA2, GAMMAD, DELTA4, DELTA6) C ISTRCH=1: CALCULATION IN STRETCHED COORDINATES (STANDARD) C (DEF PARAM EPS2, GAMMA, EPS4, EPS6) C (NOTE THAT THE SHAPES GENERATED ARE DIFFE- C RENT IN THE TWO CASES; FURTHERMORE THE C L*S AND L*L TERMS ARE SOMEWHAT DIFFERENT C WITH L DEF IN STRETCHED OR SPHER COORDINATES) C C ICORR=0: NORMAL CALCULATION OF J+/J-/JZ/J+J_, ETC, IE, C APPROXIMATE WITH J+/J-/JZ (STRETCHED) C ICORR=1: INCLUDE CORRECTION TERMS FOR J+/J-/JZ FROM THE C STRETCHED BASIS (SEE APPENDIX, NUCL PHYS A307,189) C AND CONSTRUCT CORRECTED J+J-, JZ2, ETC FROM THESE C C IREC = 0: USUAL TREATMENT OF JX2, JY2, JZ2 IN ASYR, IE, AS C ONE-BODY OPERATORS WITH UU' + VV' PAIRING FACTORS C IREC = 1: "RECOIL" TERMS JX2, JY2, JZ2 ARE TREATED AS TWO- C BODY OPERATORS IN ASYRMO => MORE j+, j_, jz C MATRIX ELEMENTS MUST BE CALCULATED HERE (FOR THE C SUMS OVER INTERMEDIATE STATES IN ASYRMO) C C ** RECOMMENDATIONS: FOR NORMAL DEFORMATIONS, THE DIFFERENCES BETWEEN C *** THE STRETCHED AND PHYSICAL j+, j_, jz OPERATORS C *** ARE NEGLIGIBLE ==> ICORR=0 IS STANDARD. C *** SIMILARLY, THE SIMPLER ONE-BODY TREATMENT OF C *** THE RECOIL TERM SEEMS QUITE ADEQUATE ==> IREC=0 C *** IS STANDARD. (NOTE THAT IF IREC=1 IS SELECTED C *** HERE, IREC=0 CAN STILL BE CHOSEN IN ASYRMO.) C C 3. NKAMY, NNEUPR C NKAMY .LE. 0: SAME KAPPA AND MY FOR ALL N-SHELLS C NKAMY .GT. 1: DIFFERENT KAPPA AND MY FOR DIFFERENT N-SHELLS C (SEE CARD 3 BELOW) C FOR NKAMY .GT. 1:(NKAMY-1) .GE. MAX(NPROT,NNEUTR); (CARD 9) C (NKAMY-1) .LE. 10 C NNEUPR=1: CALCULATION FOR PROTONS C NNEUPR=-1: CALCULATION FOR NEUTRONS C NNEUPR=2: CALCULATION FOR PROTONS AND NEUTRONS (FOR ODD-ODD) C C 4. KAPPAP(I), MYP(I), KAPPAN(I), MYN(I) C OMITTED IF NKAMY .LE. 1; OTHERWISE NKAMY CARDS NEEDED C KAPPAP AND MYP (KAPPAN AND MYN) VALUES OF L.S AND L**2 C COUPLING PARAMETERS FOR PROTONS (NEUTRONS) FOR THE C DIFFERENT SHELLS, N=0,1,2,3,..., (NKAMY-1) C C 4A. NKAMYL C NUMBER OF N-SHELLS WITH KAPPA AND MY L-DEPENDENT C C CARDS 3B-3C REPEATED NKAMYL TIMES (OMIT FOR NKAMYL=0): C 4B. NSHELL C N QUANTUM NUMBER FOR SHELLS WITH L-DEPENDENCE C 4C. KAPPAP(N,L), MYP(N,L), KAPPAN(N,L), MYN(N,L) C L=0,2,.. OR L=1,3,.. (N=NSHELL) C C *** FOR "STANDARD" KAPPA'S & MU 'S, SEE NUCL PHYS A436 (1985) 14 *** C C 5. KAPPAP, MYP, KAPPAN, MYN C THIS CARD MAY NEVER BE OMITTED BUT IS DUMMY FOR NKAMY .GT. 1 C FOR NKAMY .LE. 1 USED AS CARDS 3 BUT REFERRING TO ALL SHELLS C C 6. NUU, IPKT, NOYES, ITRANS C NUU: NUMBER OF OSCILLATOR SHELLS COUPLED, NUU=4 STAND VALUE C (REDUNDANT FOR EPS4=EPS6=OMROT=0 AND ISTRCH=1 BECAUSE C THE N-SHELLS ARE PURE IN THIS CASE) C IPKT: NUMBER OF DEFORMATIONS, CF. CARD 9 BELOW C NOYES=1 IN PART-ROT CALCULATIONS C (FOR NOYES=0 ONLY SINGLE PARTICLE ENERGIES BUT NO C WAVE-FUNCTIONS OR MATRIX ELEMENTS ARE CALCULATED) C ITRANS: CONTROLS OPTION OF TRANFORMING THE LEVELP(J) AND C LEVELN(K) PROTON AND NEUTRON ORBITALS TO THE C /N,L,LAM,SIGMA> BASIS (NEEDED IF SUBSEQUENT Vpn C CALCULATIONS ARE DESIRED, OUTPUT TO FILE2) C ITRANS=0: NO SUCH TRANSFOMATION C ITRANS=1: WAVEFUNCTIONS ARE TRANSFORMED, BUT NOT PRINTED C ITRANS=2: W.F.'S ARE TRANSFORMED AND PRINTED OUT C C 7. EMIN, EMAX C FOR ENERGIES BETWEEN EMIN AND EMAX (IN OSCILLATOR UNITS) THE C WAVE-FUNCTIONS ARE PRINTED C C 8a. IPARP, NORBP, (LEVELP(J),J=1,NORBP) FORMAT(A1,I2,15I3) C IPARP: PARITY CONSIDERED FOR PROTONS, "+" OR "-" C NORBP: NUMBER OF PROTON ORBITALS CONSIDERED IN THE C PARTICLE-ROTOR CALCULATION, NORBP .LE. 15 C LEVELP(J): "NUMBERS" ON THE ORBITALS INCLUDED C (CALCULATED FROM THE "BOTTOM" AND FOR THE PARITY "IPARP") C C 8b. IPARN,NORBN,(LEVELN(J),J=1,NORBN) C ANALOGOUS TO 7a BUT FOR NEUTRONS C C 9. Z,A C Z: PROTON NUMBER; A: MASS NUMBER C FOR Z AND/OR N=A-Z PARTICLES (Z AND A EVEN), THE MICROSCOPIC C QUADRUPOLE MOMENTS Q0 AND Q2 ARE CALCULATED (WITH A SHARP C FERMI SURFACE). AN ASYMMETRY PARAMETER GAMMA CORRESPONDING C TO THE MATTER DISTRIBUTION IS THEN DEDUCED C C 10. IPKT CARDS: C EPS2, GAMMA, EPS4, EPS6, OMROT, NPROT, NNEUTR, NSHELP, NSHELN C (ISTRCH = 1) C OR C DELTA2,GAMMAD,DELTA4,DELTA6,OMROT,NPROT,NNEUTR,NSHELP,NSHELN C (ISTRCH = 0) C FOUR DEFORMATION PARAMETERS C (EPS6 (DELTA6) IS NOT DEFINED TO HAVE THE CORRECT C TRANSFORMATION PROPERTIES FOR GAMMA .NE. 0 AND SHOULD C ONLY BE USED FOR SMALL GAMMA, SAY GAMMA .LE. 15) C OMROT: ROTATIONAL FREQUENCY (=0 IN PART-ROTOR CALC) C NPROT: NUMBER OF OSCILLATOR SHELLS INCLUDED FOR PROTONS C NNEUTR: NUMBER OF OSCILLATOR SHELLS INCLUDED FOR NEUTRONS C MAX NUMBER OF SHELLS IN PRESENT VERSION: 10 C NSHELP,NSHELN:USED IN THE TRANSFORMATION TO /N,L,LAM,SIG> C BASIS IF THE N-SHELLS ARE PURE. NSHELP IS THEN THE C N-SHELL FOR THE LEVELP(J) PROTON ORBITALS, NSHELN C FOR THE LEVELN(K) NEUTRON ORBITALS C C*********************************************************************** C COMMON /COUPL/ NUU COMMON /LU/ LU COMMON /PI/ PI COMMON /CMN1/ KAPPA,OMOM,MY,L2F(18) COMMON /CMN2/ EPS1,EPS2,EPS3,EPS4,EPS5,EPS6,EPS22 COMMON /CMN7/ EPSP,GAMMA COMMON /QSTR/ Q COMMON/E/ EMIN,EMAX COMMON/LEVEL/ NU,LEVEL(15) COMMON /NUCL/ Z,A COMMON /PAR/ IPAR COMMON/OMROT/ OMROT COMMON/EQFAC/ EFAC,QFAC,ISO COMMON/ITRANS/ ITRANS,NSHELL DIMENSION LEVELP(15),LEVELN(15) COMMON/BAS/ ISTRCH,ICORR,IREC COMMON/XKAMY/ XKAPPA(13,7),XMY(13,7),AVKAMY(13) INTEGER Z,A REAL KAPPA,MY,KAPPAP,MYP,KAPPAN,MYN REAL XKAP(13),XMYP(13),XKAN(13),XMYN(13), $ XKAPL(13,7),XMYPL(13,7),XKANL(13,7),XMYNL(13,7), $ AVKMYP(13),AVKMYN(13) CHARACTER IPAR*4,IPARP*4,IPARN*4,FILE2*24,FILE16*24,FILE17*24 integer count,len,status,iostat character(len=132) inputfile count = command_argument_count () if(count>0)then call get_command_argument (1, inputfile, len, iostat) open(unit=5,file=inputfile,status='old',form='formatted', 1 iostat=iostat) endif open(unit=6,file='GAMPN.out',status='unknown',form='formatted') LU=6 PI=3.1415926536 READ(5,*)FILE2,FILE16,FILE17 C C OPEN FILES (APPROPRIATE FOR VAX) C OPEN(2,FILE=FILE2, FORM='FORMATTED', STATUS='unknown') OPEN(16,FILE=FILE16,FORM='UNFORMATTED',STATUS='unknown') OPEN(17,FILE=FILE17,FORM='UNFORMATTED',STATUS='unknown') ! OPEN(2,FILE=FILE2, FORM='FORMATTED', STATUS='NEW') ! OPEN(16,FILE=FILE16,FORM='UNFORMATTED',STATUS='NEW') ! OPEN(17,FILE=FILE17,FORM='UNFORMATTED',STATUS='NEW') C C CALCULATIN OF FACTORIALS, GAMMA FUNCTION AND C AVERAGE WITHIN SHELLS OF L*L C CALL HELP C EPS1=0.0 EPS3=0.0 EPS5=0.0 READ(5,*) ISTRCH,ICORR,IREC READ(5,*) NKAMY,NNEUPR IF (NKAMY .LE. 1) GO TO 51 write(6,16) 16 FORMAT(///,' KAPPAP MYP KAPPAN MYN') DO 52 IN=1,NKAMY READ(5,*)XKAP(IN),XMYP(IN),XKAN(IN),XMYN(IN) NSHELL=IN-1 write(6,15) NSHELL,XKAP(IN),XMYP(IN),XKAN(IN),XMYN(IN) 15 FORMAT(I5,4F10.4) NL = (IN+1)/2 DO 57 IL = 1,NL XKAPL(IN,IL) = XKAP(IN) XKANL(IN,IL) = XKAN(IN) XMYPL(IN,IL) = XMYP(IN) XMYNL(IN,IL) = XMYN(IN) AVKMYP(IN) = XMYP(IN)*XKAP(IN) AVKMYN(IN) = XMYN(IN)*XKAN(IN) 57 CONTINUE 52 CONTINUE READ(5,*) NKAMYL IF (NKAMYL .EQ. 0) GO TO 61 write(6,18) ' SHELLS WITH KAPPA L-DEPENDENT', $' N L KAPPAP MYP KAPPAN MYN' 18 FORMAT(/,A,/,TR19,A,/) DO 63 I = 1,NKAMYL READ(5,*) NSHELL IN = NSHELL+1 NL = NSHELL/2 + 1 SUMMYP = 0. SUMMYN = 0. SUMORB = 0. DO 62 IL = 1,NL READ(5,*)XKAPL(IN,IL),XMYPL(IN,IL),XKANL(IN,IL),XMYNL(IN,IL) IF (MOD(NSHELL,2) .EQ. 0) THEN L = 2*(IL-1) ELSE L = 2*IL-1 ENDIF SUMMYP = SUMMYP + L*(L+1)*(2*L+1)*XMYPL(IN,IL)*XKAPL(IN,IL) SUMMYN = SUMMYN + L*(L+1)*(2*L+1)*XMYNL(IN,IL)*XKANL(IN,IL) SUMORB = SUMORB + L*(L+1)*(2*L+1) write(6,17) NSHELL,L, $ XKAPL(IN,IL),XMYPL(IN,IL),XKANL(IN,IL),XMYNL(IN,IL) 17 FORMAT(15X,2I6,4F10.4) 62 CONTINUE AVKMYP(IN) = SUMMYP/SUMORB AVKMYN(IN) = SUMMYN/SUMORB 63 CONTINUE 61 CONTINUE 51 CONTINUE READ(5,*)KAPPAP,MYP,KAPPAN,MYN READ(5,*)NUU,IPKT,NOYES,ITRANS READ(5,*)EMIN,EMAX READ(5,1952) IPARP,NORBP,(LEVELP(J),J=1,NORBP) READ(5,1952) IPARN,NORBN,(LEVELN(J),J=1,NORBN) 1952 FORMAT(A1,I2,15I3) READ(5,*)Z,A N=A-Z DO 1 I=1,IPKT ! IRK=MSLEFT(XX) READ(5,*)EPSP,GAMMA,EPS4,EPS6,OMROT,NPROT,NNEUTR,NSHELP,NSHELN IF(NPROT.LT.3.OR.NNEUTR.LT.3) GO TO 111 IF(NPROT.GT.12.OR.NNEUTR.GT.12) GO TO 111 GAMMAR=GAMMA*PI/180.0 EPS2=EPSP*COS(GAMMAR) EPS22=EPSP*SIN(GAMMAR)/SQRT(0.75) C C CALCULATION OF VOLUME CONSERVATION FACTOR, OMOM C CALL OMO(OMOM) C C*************SINGLE PARTICLE ENERGY LEVELS CALCULATION IF (I .NE. 1) write(6,12) 12 FORMAT(1H1) IF (NNEUPR .LT. 0) GO TO 71 IF (ISTRCH .EQ. 1) THEN write(6,3) ELSEIF (ISTRCH .EQ. 0) THEN write(6,13) ENDIF 3 FORMAT(//,' KAPPA MY EPS GAMMA EPS4 ', *' EPS6 W0/W00 NMAX COUPL OMROT EFAC QFAC') 13 FORMAT(//,' KAPPA MY DELTA GAMMAD DELTA4 ', *'DELTA6 W0/W00 NMAX COUPL OMROT EFAC QFAC') ISO=1 IPAR = IPARP KAPPA=KAPPAP MY=MYP NPROT1 = NPROT + 1 DO 53 IN=1,NPROT1 NL = (IN-1)/2+1 DO 58 IL=1,NL IF (NKAMY .GE. 2) THEN XKAPPA(IN,IL) = XKAPL(IN,IL) XMY(IN,IL) = XMYPL(IN,IL) AVKAMY(IN) = AVKMYP(IN) ELSE XKAPPA(IN,IL) = KAPPAP XMY(IN,IL) = MYP AVKAMY(IN) = MYP*KAPPAP ENDIF 58 CONTINUE 53 CONTINUE C********FACTORS FOR TRANSFORMATION FROM OSC UNITS TO PHYSICAL UNITS EFAC=41.0/FLOAT(A)**(1./3.)*(1.-FLOAT(N-Z)/A/3.)*OMOM QFAC=41.5/EFAC NU = NORBP DO 56 IP=1,NU LEVEL(IP) = LEVELP(IP) 56 CONTINUE NSHELL=NSHELP IF (NKAMY .LE. 0) write(6,100) KAPPA,MY, *EPSP,GAMMA,EPS4,EPS6, OMOM,NPROT,NUU,OMROT,EFAC,QFAC 100 FORMAT(2X,F7.4,2F8.3,F7.1,F8.3,F9.3,F11.5,2I6,F12.5,2F10.5,/) IF (NKAMY .GT. 1) write(6,101) ' SEE TABLE ', *EPSP,GAMMA,EPS4,EPS6, OMOM,NPROT,NUU,OMROT,EFAC,QFAC 101 FORMAT(2X,A,F8.3,F7.1,F8.3,F9.3,F11.5,2I6,F12.5,2F10.5,/) IF (ISTRCH .EQ. 0) write(6,201) IF (ISTRCH .EQ. 1) write(6,202) 201 FORMAT(/,20X,' **** PROTONS ****',20X, & ' **** SPHERICAL COORDINATES ****') 202 FORMAT(/,20X,' **** PROTONS ****',20X, & ' **** STRETCHED COORDINATES ****') C CALCULATION OF ENERGY LEVELS, MATRIX ELEMENTS ETC + OUTPUT C CALL ENERLJ(OMOM,NPROT,NOYES) C ! JRK=MSLEFT(XX) TID=(IRK-JRK)*0.001 write(6,10) TID 10 FORMAT(1X,16H TOTAL TIME:,F6.1,4H SEC) IF (NNEUPR .EQ. 1) GO TO 1 71 CONTINUE C CALCULATIONS FOR NEUTRONS C ! IRK=MSLEFT(XX) ISO=2 IPAR = IPARN KAPPA = KAPPAN MY = MYN NNEUT1 = NNEUTR + 1 DO 55 IN = 1,NNEUT1 NL = (IN-1)/2+1 DO 59 IL=1,NL IF (NKAMY .GT. 2) THEN XKAPPA(IN,IL) = XKANL(IN,IL) XMY(IN,IL) = XMYNL(IN,IL) AVKAMY(IN) = AVKMYN(IN) ELSE XKAPPA(IN,IL) = KAPPAN XMY(IN,IL) = MYN AVKAMY(IN) = MYN*KAPPAN ENDIF 59 CONTINUE 55 CONTINUE NU = NORBN DO 54 IN=1,NU LEVEL(IN) = LEVELN(IN) 54 CONTINUE NSHELL=NSHELN C********FACTORS FOR TRANSFORMATION FROM OSC UNITS TO PHYSICAL UNITS EFAC=41.0/FLOAT(A)**(1./3.)*(1.+FLOAT(N-Z)/A/3.)*OMOM QFAC=41.5/EFAC write(6,12) IF (ISTRCH .EQ. 1) THEN write(6,3) ELSEIF (ISTRCH .EQ. 0) THEN write(6,13) ENDIF IF (NKAMY .LE. 0) write(6,100) KAPPA,MY, *EPSP,GAMMA,EPS4,EPS6, OMOM,NNEUTR,NUU,OMROT,EFAC,QFAC IF (NKAMY .GT. 1) write(6,101) ' SEE TABLE ', *EPSP,GAMMA,EPS4,EPS6, OMOM,NNEUTR,NUU,OMROT,EFAC,QFAC IF (ISTRCH .EQ. 0) write(6,203) IF (ISTRCH .EQ. 1) write(6,204) 203 FORMAT(/,20X,' **** NEUTRONS ****',19X, & ' **** SPHERICAL COORDINATES ****') 204 FORMAT(/,20X,' **** NEUTRONS ****',19X, & ' **** STRETCHED COORDINATES ****') CALL ENERLJ(OMOM,NNEUTR,NOYES) C c JRK=MSLEFT(XX) TID=(IRK-JRK)*0.001 write(6,10) TID 1 CONTINUE GO TO 112 111 CONTINUE 113 FORMAT(1X,30HINPUT CARD INCORRECTLY PUNCHED) write(6,113) STOP 112 CONTINUE 1946 FORMAT( ) STOP END C C####################################################################### C SUBROUTINE QNUMLJ(BB,NP,L,JJ,JOMEGA,SN) C C SET UP OF BASIS FUNCTIONS (ONE FOR EACH CALL): C ON OUTPUT: NP=N, L=L, JJ=2*J, JOMEGA=2*OMEGA C INPUT: BB=NUMBERING OF WAVE-FUNCTIONS C SN=MAX NUMBER OF SHELLS (FOR GIVEN PARITY) C INTEGER B,SN,ORDN,F,D,BB F=0 B=BB OMEGA=SN+0.5 IF((SN/2)*2.NE.SN) OMEGA=SN-0.5 ORDN=SN-OMEGA+1.6 NP=SN I1=ORDN I2=B-1 IF(I1.GT.I2) GO TO 2 D=I1 1 CONTINUE F=F+ORDN OMEGAM=-NP+0.5 IF((SN/2)*2.NE.SN) OMEGAM=-NP-0.5 IF(ABS(OMEGA-OMEGAM).GT.0.1) GO TO 10 NP=NP-2 OMEGA=NP+0.5 IF((NP/2)*2.NE.NP) OMEGA=NP-0.5 GO TO 11 10 CONTINUE OMEGA=OMEGA-2.0 11 CONTINUE ORDN=NP-ABS(OMEGA)+1.6 D=D+ORDN IF(D.LE.I2) GO TO 1 2 CONTINUE B=B-F RJ=FLOAT(NP)+0.5 IF((B/2)*2.EQ.B) GO TO 3 S=-0.5 GO TO 4 3 CONTINUE S=0.5 4 CONTINUE RJ=RJ+1.0-FLOAT(B) L=RJ+S+0.1 JJ=2.*RJ+0.1 JOMEGA=2.*OMEGA+0.1*SIGN(1.0,OMEGA) RETURN END C C####################################################################### C SUBROUTINE QTME C CALCULATION OF QUADRUPOLE TRANSITION MATRIX ELEMENTS WITHIN THOSE C ORBITALS CONSIDERED IN PART-ROT CALCULATION C***** QUADRUPOLE OPERATORS DEFINED AS: C Q20=SQRT(16*PI/5)*R2*Y20 (=2Z2-X2-Y2) C Q22=SQRT(4*PI/5)*R2*(Y22+Y2-2) (=SQRT(3/2)*(X2-Y2)) C Q20, Q22, R2 CALCULATED IN SPHERICAL BASIS (PHYSICAL VALUES) C Q20S, Q22S AND RHO2 SAME MATRIX ELEMENTS IN STRETCHED BASIS C (ISTRCH=1) C RJZ: MATRIX ELEMENTS OF JZ C KV1=N, KV2=L, KV3=2*J, KV4=2*OMEGA, KATT=PARITY ('+' OR '-') C COMMON/CMN2/ EPS1,EPS2,EPS3,EPS4,EPS5,EPS6,EPS22 COMMON/DIM/ IDIM,KATT CHARACTER KATT*4 COMMON/KVANT/ KV1(161),KV2(161),KV3(161),KV4(161) COMMON/EIG/ EIGVEC(161,161),VALU(161) COMMON/BAS/ ISTRCH,ICORR,IREC COMMON/FACTOR/ GAC,HAC COMMON/QUADR/ Q20S(15,15),Q22S(15,15),Q20(15,15),Q22(15,15), *RHO2(15,15),R2(15,15),RJZ(15,15) COMMON/LEVEL/ NU,LEVEL(15) INTEGER RNP,RLP,RJP,OMEGAP,RN,RL,RJ,OMEGA DATA PI/3.1415926536/ FAC=(1.+EPS2/3.+EPS22/2.)*(1.+EPS2/3.-EPS22/2.)*(1.-2.*EPS2/3.) IF (NU .LE. 15) GO TO 4 write(6,10) 10 FORMAT(1H1,'QTME STOP') STOP 4 CONTINUE DO 3 MYP=1,NU DO 3 MY =1,NU Q20S(MYP,MY)=0.0 Q22S(MYP,MY)=0.0 Q20(MYP,MY)=0.0 Q22(MYP,MY)=0.0 RHO2(MYP,MY)=0.0 R2(MYP,MY)=0.0 RJZ(MYP,MY)=0.0 3 CONTINUE DO 2 IP=1,IDIM RNP=KV1(IP) RLP=KV2(IP) RJP=KV3(IP) OMEGAP=KV4(IP) DO 2 I=1,IDIM RN=KV1(I) RL=KV2(I) RJ=KV3(I) OMEGA=KV4(I) C***************** IS THIS COMBINATION OF BASIS VECTORS SIGNIFICANT? *** C EIGMAX = 1.E-8 C DO 25 MMY=1,NU C DO 24 NNY=MMY,NU C MY = LEVEL(MMY) C NY = LEVEL(NNY) C EIGABS = ABS(EIGVEC(IP,MY)*EIGVEC(I,NY)) C EIGMAX = MAX(EIGABS,EIGMAX) C 24 CONTINUE C 25 CONTINUE C IF (EIGMAX .LT. 1.E-5) GO TO 2 C C IF ((ABS(EPS4)+ABS(EPS6)) .LT. 0.001 .AND. RNP .NE. RN .AND. C * (ISTRCH .EQ. 1)) GO TO 2 IF(IABS(RNP-RN).GT.2.OR.IABS(RJP-RJ).GT.4.OR.IABS(RLP-RL).GT.2. $OR.IABS(OMEGAP-OMEGA).GT.4) GO TO 6 CGC=0.0 CGC1=0.0 CGC2=0.0 RME=0.0 HLZ=0.0 HJZ=0.0 C CALCULATION OF RADIAL MATRIX ELEMENTS, R2 C RME=RMATEL(2*RNP,2*RLP,2*RN,2*RL,2) C IF(OMEGAP.EQ.OMEGA) CGC=CLEBI(RJ,4,RJP,-1,0,-1)* $CLEBI(RJ,4,RJP,OMEGA,0,OMEGA) IF(OMEGAP.EQ.OMEGA+4) CGC1=CLEBI(RJ,4,RJP,-1,0,-1)* $CLEBI(RJ,4,RJP,OMEGA,4,OMEGAP) IF(OMEGAP.EQ.OMEGA-4) CGC2=CLEBI(RJ,4,RJP,-1,0,-1)* $CLEBI(RJ,4,RJP,OMEGA,-4,OMEGAP) IF(RNP.EQ.RN.AND.RLP.EQ.RL.AND.RJP.EQ.RJ.AND.OMEGAP.EQ.OMEGA) $HJZ=FLOAT(OMEGA)/2. IF(RNP.EQ.RN.AND.RLP.EQ.RL.AND.OMEGAP.EQ.OMEGA) HLZ= $(OMEGA-1.)/2.*CLEBI(2*RL,1,RJP,OMEGA-1,1,OMEGA)* $CLEBI(2*RL,1,RJ,OMEGA-1,1,OMEGA)+ $(OMEGA+1.)/2.*CLEBI(2*RL,1,RJP,OMEGA+1,-1,OMEGA)* $CLEBI(2*RL,1,RJ,OMEGA+1,-1,OMEGA) DO 5 MYP=1,NU DO 5 MY=1,NU NYP=LEVEL(MYP) NY=LEVEL(MY) Q20S(MYP,MY)=Q20S(MYP,MY)+EIGVEC(IP,NYP)*EIGVEC(I,NY)* $2.0*RME*SQRT((RJ+1.)/(RJP+1.))*CGC Q22S(MYP,MY)=Q22S(MYP,MY)+EIGVEC(IP,NYP)*EIGVEC(I,NY)* $RME*SQRT((RJ+1.)/(RJP+1.))*(CGC1+CGC2) IF(RLP.EQ.RL.AND.RJP.EQ.RJ.AND.OMEGAP.EQ.OMEGA) $RHO2(MYP,MY)=RHO2(MYP,MY)+EIGVEC(IP,NYP)*EIGVEC(I,NY)*RME RJZ(MYP,MY)=RJZ(MYP,MY)+EIGVEC(IP,NYP)*EIGVEC(I,NY)* $(HJZ+GAC*HLZ-HAC/SQRT(1.5)*(CGC1-CGC2)*SQRT((RJ+1.)/(RJP+1.))*RME) 1 CONTINUE 5 CONTINUE 6 CONTINUE 2 CONTINUE IF (ISTRCH .EQ. 0) GO TO 21 DO 7 MYP=1,NU DO 7 MY=1,NU Q20(MYP,MY)=((1.+EPS2/3.-EPS22*EPS22/6.)*Q20S(MYP,MY)+ $EPS22/SQRT(6.)*(1.-2.*EPS2/3.)*Q22S(MYP,MY)+ $(2.*EPS2/3.*(1.+EPS2/3.)-EPS22*EPS22/6.)*RHO2(MYP,MY))/FAC Q22(MYP,MY)=((1.+EPS2/3.)*Q22S(MYP,MY)+ $EPS22/SQRT(24.)*Q20S(MYP,MY)- $EPS22/SQRT(6.)*RHO2(MYP,MY))/FAC*(1.-2.*EPS2/3.) R2(MYP,MY)=(((1.+EPS2/3.)*(1.-EPS2/3.)-1./3.*(EPS22/2.)**2)* $RHO2(MYP,MY)+((1.+EPS2/3.)*EPS2/3.-1./3.*(EPS22/2.)**2)* $Q20S(MYP,MY)-EPS22/2.*(1.-2.*EPS2/3.)*SQRT(2./3.)*Q22S(MYP,MY))/ $FAC 7 CONTINUE GO TO 22 21 CONTINUE DO 23 MYP=1,NU DO 23 MY=1,NU Q20(MYP,MY) = Q20S(MYP,MY) Q22(MYP,MY) = Q22S(MYP,MY) R2(MYP,MY) = RHO2(MYP,MY) 23 CONTINUE 22 CONTINUE WRITE(17) NU,(LEVEL(J),J=1,NU),IDIM,(VALU(J),J=1,IDIM), $((Q20(MYP,MY),MY=1,NU),MYP=1,NU),((Q22(MYP,MY),MY=1,NU),MYP=1,NU) RETURN END C C####################################################################### C SUBROUTINE Q20Q22 C C SAME AS QTME BUT ONLY DIAG MATR ELEM FOR ALL LOW-ENERGY ORBITALS C***** QUADRUPOLE MOMENTS IN UNITS OF 'HBAR/MW0' C Q20=SQRT(16*PI/5)*R2*Y20 (=2Z2-X2-Y2) C Q22=SQRT(4*PI/5)*R2*(Y22+Y2-2) (=SQRT(3/2)*(X2-Y2)) C COMMON /CMN2/ EPS1,EPS2,EPS3,EPS4,EPS5,EPS6,EPS22 COMMON /DIM/ IDIM,KATT CHARACTER KATT*4 COMMON /KVANT/ KV1(161),KV2(161),KV3(161),KV4(161) COMMON/EIG/ EIGVEC(161,161),VALU(161) COMMON /QMOM/ Q20S(70),Q22S(70),Q20(70),Q22(70),RHO2(70) $,R2(70),RJZ(70) COMMON /NVEC/ NVEC,NQVEC COMMON/FACTOR/ GAC,HAC COMMON/BAS/ ISTRCH,ICORR,IREC COMMON/OMROT/ OMROT INTEGER RNP,RLP,RJP,OMEGAP,RN,RL,RJ,OMEGA DATA PI/3.1415926536/ FAC=(1.+EPS2/3.+EPS22/2.)*(1.+EPS2/3.-EPS22/2.)*(1.-2.*EPS2/3.) DO 3 NY=1,NQVEC Q20S(NY)=0.0 Q22S(NY)=0.0 Q20(NY)=0.0 Q22(NY)=0.0 RHO2(NY)=0.0 R2(NY)=0.0 RJZ(NY)=0.0 3 CONTINUE DO 2 IP=1,IDIM RNP=KV1(IP) RLP=KV2(IP) RJP=KV3(IP) OMEGAP=KV4(IP) DO 2 I=1,IDIM RN=KV1(I) RL=KV2(I) RJ=KV3(I) OMEGA=KV4(I) C***************** IS THIS COMBINATION OF BASIS VECTORS SIGNIFICANT? *** C EIGMAX = 1.E-8 C DO 25 NY=1,NVEC C EIGABS = ABS(EIGVEC(IP,NY)*EIGVEC(I,NY)) C EIGMAX = MAX(EIGABS,EIGMAX) C 25 CONTINUE C IF (EIGMAX .LT. 1.E-5) GO TO 2 C C ABCOUP = ABS(EPS4) + ABS(EPS6) + ABS(OMROT) C IF ((ABCOUP) .LT. 0.001 .AND. RNP .NE. RN .AND. ISTRCH .EQ. 1) C *GO TO 2 IF(IABS(RNP-RN).GT.2.OR.IABS(RJP-RJ).GT.4.OR.IABS(RLP-RL).GT.2. $OR.IABS(OMEGAP-OMEGA).GT.4) GO TO 6 CGC=0.0 CGC1=0.0 CGC2=0.0 RME=0.0 HLZ=0.0 HJZ=0.0 RME=RMATEL(2*RNP,2*RLP,2*RN,2*RL,2) IF(OMEGAP.EQ.OMEGA) CGC=CLEBI(RJ,4,RJP,-1,0,-1)* $CLEBI(RJ,4,RJP,OMEGA,0,OMEGA) IF(OMEGAP.EQ.OMEGA+4) CGC1=CLEBI(RJ,4,RJP,-1,0,-1)* $CLEBI(RJ,4,RJP,OMEGA,4,OMEGAP) IF(OMEGAP.EQ.OMEGA-4) CGC2=CLEBI(RJ,4,RJP,-1,0,-1)* $CLEBI(RJ,4,RJP,OMEGA,-4,OMEGAP) IF(RNP.EQ.RN.AND.RLP.EQ.RL.AND.RJP.EQ.RJ.AND.OMEGAP.EQ.OMEGA) $HJZ=FLOAT(OMEGA)/2. IF(RNP.EQ.RN.AND.RLP.EQ.RL.AND.OMEGAP.EQ.OMEGA) HLZ= $(OMEGA-1.)/2.*CLEBI(2*RL,1,RJP,OMEGA-1,1,OMEGA)* $CLEBI(2*RL,1,RJ,OMEGA-1,1,OMEGA)+ $(OMEGA+1.)/2.*CLEBI(2*RL,1,RJP,OMEGA+1,-1,OMEGA)* $CLEBI(2*RL,1,RJ,OMEGA+1,-1,OMEGA) DO 5 NY=1,NQVEC Q20S(NY)=Q20S(NY)+EIGVEC(IP,NY)*EIGVEC(I,NY)* $2.0*RME*SQRT((RJ+1.)/(RJP+1.))*CGC Q22S(NY)=Q22S(NY)+EIGVEC(IP,NY)*EIGVEC(I,NY)* $RME*SQRT((RJ+1.)/(RJP+1.))*(CGC1+CGC2) IF(RLP.EQ.RL.AND.RJP.EQ.RJ.AND.OMEGAP.EQ.OMEGA) $RHO2(NY)=RHO2(NY)+EIGVEC(IP,NY)*EIGVEC(I,NY)*RME RJZ(NY)=RJZ(NY)+EIGVEC(IP,NY)*EIGVEC(I,NY)* $(HJZ+GAC*HLZ-HAC/SQRT(1.5)*(CGC1-CGC2)*SQRT((RJ+1.)/(RJP+1.))*RME &*(RNP-RN)/2.) 5 CONTINUE 6 CONTINUE 2 CONTINUE IF (ISTRCH .EQ. 0) GO TO 21 DO 7 NY=1,NQVEC Q20(NY)=((1.+EPS2/3.-EPS22*EPS22/6.)*Q20S(NY)+ $EPS22/SQRT(6.)*(1.-2.*EPS2/3.)*Q22S(NY)+ $(2.*EPS2/3.*(1.+EPS2/3.)-EPS22*EPS22/6.)*RHO2(NY))/FAC Q22(NY)=((1.+EPS2/3.)*Q22S(NY)+ $EPS22/SQRT(24.)*Q20S(NY)- $EPS22/SQRT(6.)*RHO2(NY))/FAC*(1.-2.*EPS2/3.) R2(NY)=(((1.+EPS2/3.)*(1.-EPS2/3.)-1./3.*(EPS22/2.)**2)* $RHO2(NY)+((1.+EPS2/3.)*EPS2/3.-1./3.*(EPS22/2.)**2)*Q20S(NY)- $EPS22/2.*(1.-2.*EPS2/3.)*SQRT(2./3.)*Q22S(NY))/FAC 7 CONTINUE GO TO 22 21 CONTINUE DO 23 NY=1,NQVEC Q20(NY) = Q20S(NY) Q22(NY) = Q22S(NY) R2(NY) = RHO2(NY) 23 CONTINUE 22 CONTINUE RETURN END C C####################################################################### C SUBROUTINE DECOUP C C***** CALCULATION OF GENERALIZED DECOUPLING FACTORS ETC. C IN THE CASE OF STRETCHED COORDINATES (ISTRCH=1) THE DECOUPLING C FACTORS ARE APPROXIMATED BY THEIR "STRETCHED VALUES" C (CF. APPENDIX, NUCL PHYS 307 (1978) 189 C IF ICORR=1, // MATRIX ELEMENTS ARE CORRECTED AS C DESCRIBED IN THE APPENDIX, NPA 307,189 (1978) C AND CORRECTED D, E, AND JZ2 ARE CONTRUCTED FROM C J+, J-, JZ BY SUMMING OVER INTERMEDIATE STATES C C GENERALISED DECOUPLING FACTORS: APLB, C, D, E C APLB = A+B = , C = , D = , E = J-**2> C NORM: NORMALIZATION INTEGRAL C JZ = , JZ2 = C SNOLL = , SPLUS = , SMINUS = C R2SZ = <(R**2)*SZ>,.... R2SMIN = <(R**2)*S->,... ETC C D2 = -(SQRT(10)/2)(S*C(2))**1 - TENSOR OF RANK 1 C *** SPHERICAL COMPONENTS CALC FOR D2, S+,S-,J+,J- FOR S AND J ** C C(2M) = SQRT(4*PI/5)*Y(2M) WHERE Y(2M) SPHERICAL HARMONIC C Q20 ETC. TRANSFERRED FROM QTME FOR PRINTING C C IF IREC=1, THE 'RECOIL TERMS' (J+J-, J+J+, JZ2, ETC) ARE TREATED AS C TWO BODY TERMS IN ASYR, WHICH REQUIRES THAT THE MATRIX ELEMENTS OF C J+, J- AND JZ BE CALCULATED FOR A COMPLETE SET OF INTERMEDIATE STATES. C COMMON/KVANT/ KV1(161),KV2(161),KV3(161),KV4(161) COMMON/EIG/ EIGVEC(161,161),VALU(161) COMMON/DIM/ IDIM,KATT CHARACTER KATT*4 COMMON/QUADR/ Q20S(15,15),Q22S(15,15),Q20(15,15),Q22(15,15), *RHO2(15,15),R2(15,15),RJZ(15,15) COMMON/MATREL/ APLB(15,15),C(15,15),D(15,15),E(15,15), *NORM(15,15),JZ(15,15),JZ2(15,15), *SNOLL(15,15),SPLUS(15,15),SMINUS(15,15), *R2SZ(15,15),R2SPLU(15,15),R2SMIN(15,15), *R4SZ(15,15),R4SPLU(15,15),R4SMIN(15,15), *R2JZ(15,15),R2JPLU(15,15),R2JMIN(15,15), *R4JZ(15,15),R4JPLU(15,15),R4JMIN(15,15), *R2DZ(15,15),R2DPLU(15,15),R2DMIN(15,15), *R4DZ(15,15),R4DPLU(15,15),R4DMIN(15,15) REAL JZ,JZ2,NORM C C**** EQUIVA IS INTRODUCED TO REFER TO ALL MATRIX ELEMENTS C**** SIMULTANEOUSLY: INITIALIZE TO ZERO, MAKE SYMMETRIC AND OUTPUT REAL EQUIVA(15,15,28) EQUIVALENCE (APLB,EQUIVA) C C*** JZR, JPLR, AND JMNR ARE NEEDED TO CALCULATE MATRIX ELEMENTS OF C*** THE RECOIL TERM IN ASYR IF IT IS TREATED AS A TWO BODY TERM. C*** SZR, ...,FZR ARE NEEDED ALSO IF THE CORRECTION DUE TO STRETCHING C*** IS DESIRED C C**** FMNR, FPLR, FZR AS DEFINED IN APP., NUCL PHYS A307, 189 C**** NOTE THAT FPLR IS PROPORTIONAL TO Y2,-1 C**** FMNR Y2,+1 C**** FZR Y2,+2 - Y2,-2 C DIMENSION JPLR(161,15),JMNR(161,15),JZR(161,15), & SPLR(161,15),SMNR(161,15),SZR(161,15), & FPLR(161,15),FMNR(161,15),FZR(161,15) REAL JZR,JPLR,JMNR REAL JPLRS C CQS** Q2TEST = Q20S, SZ=SNOLL, SPLU=-SPLUS/SQRT(2), SMIN=SMINUS/SQRT(2) CQS** INTRODUCED AS TEST CASES CQS REAL Q2TEST(15,15), CQS *SZ(15,15),SPLU(15,15),SMIN(15,15) INTEGER RNP,RLP,RJP,OMEGAP,RN,RL,RJ,OMEGA COMMON/LEVEL/ NU,LEVEL(15) COMMON/BAS/ ISTRCH,ICORR,IREC COMMON/CMN7/ EPSP,GAMMA DATA PI/3.1415926536/ ICORRS=ICORR*ISTRCH IF (ICORRS .EQ. 1) write(6,101) 101 FORMAT(1H ,20(1H*),' NOTE: ICORR=1 => J+, J-, JZ, J+J-, J+J-,' & ' JZ2, CORRECTIONS DUE TO STRETCHED BASIS APPLIED ',20(1H*),///) write(6,100) 100 FORMAT(1H ,46(1H*),' SINGLE PARTICLE MATRIX ELEMENTS ',47(1H*), *///' ENY EMY A+B C D E JZ ', *' S0 S+ S- R2DZ R2D+ R2D- ', CC *' JZ2 S+ S- R2DZ R2D+ R2D- ', *' Q20S Q22S Q20 Q22 RHO2 R2',/) C 100 FORMAT(1H ,42(1H*),' SINGLE PARTICLE MATRIX ELEMENTS ',42(1H*), C *///' ENY EMY A+B C D E ', C *' ', C *' ',/, C *' ', C *' ', C *' Q20S Q22S Q20 Q22 RHO2 R2',/) IF (NU .LE. 15) GO TO 10 write(6,80) NU 80 FORMAT(1H1,'NYIC=',I3) STOP 10 CONTINUE DO 1001 I=1,NU DO 1001 J=1,NU DO 1001 K=1,28 EQUIVA(I,J,K) = 0. 1001 CONTINUE IF (IREC .EQ. 1 .OR. ICORRS .EQ. 1) THEN DO 1002 I=1,IDIM DO 1002 J=1,NU JZR (I,J) = 0.0 JPLR(I,J) = 0.0 JMNR(I,J) = 0.0 SZR (I,J) = 0.0 SPLR(I,J) = 0.0 SMNR(I,J) = 0.0 FZR (I,J) = 0.0 FPLR(I,J) = 0.0 FMNR(I,J) = 0.0 1002 CONTINUE ENDIF CQS DO 1000 NN=1,NU CQS DO 1000 MM=1,NU CQS SZ(NN,MM)=0.0 CQS SMIN(NN,MM)=0.0 CQS SPLU(NN,MM)=0.0 CQS Q2TEST(NN,MM)=0.0 C1000 CONTINUE C C***** COMPUTE MATRIX ELEMENTS C***** AND C C***** FIRST LOOP OVER BASIS FUNCTIONS (INITIAL STATE) C DO 11 IP = 1,IDIM RNP = KV1(IP) RLP = KV2(IP) RJP = KV3(IP) OMEGAP = KV4(IP) C C***** SECOND LOOP OVER BASIS FUNCTIONS (FINAL STATE) C DO 12 I=1,IDIM RN = KV1(I) RL = KV2(I) RJ = KV3(I) OMEGA = KV4(I) C***************** IS THIS COMBINATION OF BASIS VECTORS SIGNIFICANT? *** C EIGMAX = 1.E-8 C DO 21 MMY=1,NU C DO 22 NNY=MMY,NU C MY = LEVEL(MMY) C NY = LEVEL(NNY) C EIGABS = ABS(EIGVEC(IP,MY)*EIGVEC(I,NY)) C EIGMAX = MAX(EIGABS,EIGMAX) C 22 CONTINUE C 21 CONTINUE C IF (EIGMAX .LT. 1.E-5) GO TO 12 C CQ IF (ABS(RL-RLP) .GT. 2) GO TO 12 CC IF (ABS(RJ-RJP) .GT. 2) GO TO 12 C C*********************************************RADIAL MATRIX ELEMENT***** IF (RN .NE. RNP .OR. RL .NE. RLP) THEN RAD0 = 0. ELSE RAD0 = 1. ENDIF RAD2 = RMATEL(2*RNP,2*RLP,2*RN,2*RL,2) CALL MRHOT(RNP,RLP,0,RN,RL,0,4,RAD4) C C**** REDUCED MATRIX ELEMENT FOR TENSOR OPERATOR D *** SEE EDMONDS FOR * C *** DEFINITIONS ***** REDMD = REDMD1(2*RL,1,RJ,2*RLP,1,RJP) C C**** REDUCED MATRIX ELEMENT OF S USED AS A TEST CASE *************** CS CS IF (RL .EQ. RLP) REDMS = 0.5*SQRT(3.)*(-1.)**((RJP-OMEGAP)/2)* CS *(CLEBI(2*RL,1,RJ,0,1,1)*CLEBI(2*RLP,1,RJP,0,1,1)- CS *CLEBI(2*RL,1,RJ,2,-1,1)*CLEBI(2*RLP,1,RJP,2,-1,1))/ CS *CLEBI(RJ,RJP,2,1,-1,0) C IF (OMEGAP .NE. OMEGA) GO TO 41 C C**** MATRIX ELEM OF SCALAR OPERATORS AND ZERO COMP OF TENSOR OPERATORS C X = (RJ*(RJ+2.) - OMEGA*OMEGA)/4. XA = OMEGA/2. XB = FLOAT(OMEGA*OMEGA)/4. XD = (-1.)**((RJP-OMEGAP)/2)/SQRT(3.)* *CLEBI(RJ,RJP,2,OMEGA,-OMEGAP,0) CG = CLEBI(2*RL,1,RJP,OMEGAP-1, 1,OMEGAP)* * CLEBI(2*RL,1,RJ ,OMEGA -1, 1,OMEGA )- * CLEBI(2*RL,1,RJP,OMEGAP+1,-1,OMEGAP)* * CLEBI(2*RL,1,RJ ,OMEGA +1,-1,OMEGA ) CQ MJ=OMEGA CQ XQ=2.*SQRT((2*RLP+1.)/(2*RL+1.))* CQ *(CLEBI(2*RL,1,RJ,MJ-1,1,MJ)*CLEBI(2*RLP,1,RJP,MJ-1,1,MJ)* CQ *CLEBI(4,2*RLP,2*RL,0,MJ-1,MJ-1) + CQ *CLEBI(2*RL,1,RJ,MJ+1,-1,MJ)*CLEBI(2*RLP,1,RJP,MJ+1,-1,MJ)* CQ *CLEBI(4,2*RLP,2*RL,0,MJ+1,MJ+1))* CQ *CLEBI(4,2*RLP,2*RL,0,0,0) C******************************* START LOOPS OF NILSSON ORBITALS ******* C DO 31 MMY=1,NU DO 32 NNY=MMY,NU MY = LEVEL(MMY) NY = LEVEL(NNY) EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY) IF (RNP .NE. RN .OR. RLP .NE. RL .OR. RJP .NE. RJ) GO TO 51 NORM(NNY,MMY) = NORM(NNY,MMY) + EIGV12 D (NNY,MMY) = D (NNY,MMY) + EIGV12*X JZ2(NNY,MMY) = JZ2(NNY,MMY) + EIGV12*XB 51 CONTINUE IF (RLP .NE. RL) GO TO 52 SNOLL(NNY,MMY) = SNOLL(NNY,MMY) + EIGV12*CG*0.5*RAD0 R2SZ (NNY,MMY) = R2SZ (NNY,MMY) + EIGV12*CG*0.5*RAD2 R4SZ (NNY,MMY) = R4SZ (NNY,MMY) + EIGV12*CG*0.5*RAD4 CS SZ(NNY,MMY) = SZ(NNY,MMY) + EIGV12*XD*REDMS*RAD0 52 CONTINUE IF (RJP .NE. RJ) GO TO 53 JZ (NNY,MMY) = JZ (NNY,MMY) + EIGV12*XA*RAD0 R2JZ(NNY,MMY) = R2JZ(NNY,MMY) + EIGV12*XA*RAD2 R4JZ(NNY,MMY) = R4JZ(NNY,MMY) + EIGV12*XA*RAD4 53 CONTINUE R2DZ(NNY,MMY) = R2DZ(NNY,MMY) + EIGV12*XD*REDMD*RAD2 R4DZ(NNY,MMY) = R4DZ(NNY,MMY) + EIGV12*XD*REDMD*RAD4 CQ Q2TEST(NNY,MMY) = Q2TEST(NNY,MMY) + EIGV12*XQ*RAD2 32 CONTINUE 31 CONTINUE C IF (IREC .EQ. 1 .OR. ICORRS .EQ. 1) THEN DO 131 MMY=1,NU MY = LEVEL(MMY) IF (EIGVEC(IP,MY) .EQ. 0.0) GO TO 132 caes DO 132 NY=1,IDIM DO NY=1,IDIM ! AES EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY) IF (RLP .EQ. RL) & SZR(NY,MMY) = SZR(NY,MMY) + EIGV12*CG*0.5*RAD0 IF (RJP .EQ. RJ) & JZR(NY,MMY) = JZR(NY,MMY) + EIGV12*XA*RAD0 ENDDO ! AES 132 CONTINUE 131 CONTINUE ENDIF C C********************************** END LOOPS OVER NILSSON ORBITALS **** 41 CONTINUE IF (OMEGAP .NE. -OMEGA+2) GO TO 42 C C**** (+1) COMPONENT OF TENSOR OPERATORS ****************************** C X = (-1.)**((RJ-OMEGA)/2)*0.5*SQRT((RJ+OMEGA)*(RJ-OMEGA+2.)) XD = (-1.)/SQRT(3.)*CLEBI(RJ,RJP,2,OMEGA,OMEGAP,2) CGX=CLEBI(2*RL,1,RJP,OMEGAP-1,1,OMEGAP)* *CLEBI(2*RL,1,RJ,-(OMEGA-1),-1,-OMEGA)*(-1.)**((RJ-OMEGA)/2) C IF (IABS(RNP-RN) .EQ. 2 .AND. IABS(RLP-RL) .LE. 2) THEN CGF = (-1)**((RJP-OMEGAP)/2)*FLOAT(RNP-RN)/SQRT(1.5)* &SQRT((RJP+1.0)/(RJ+1.0))*CLEBI(RJP,4,RJ,-OMEGAP,2,OMEGA)* &CLEBI(RJP,4,RJ,-1,0,-1) ELSE CGF = 0.0 ENDIF C****************************** START LOOPS OVER NILSSON ORBITALS ****** C DO 33 MMY=1,NU DO 34 NNY=MMY,NU MY = LEVEL(MMY) NY = LEVEL(NNY) EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY) IF (RLP .NE. RL) GO TO 62 SPLUS(NNY,MMY) = SPLUS(NNY,MMY) + EIGV12*CGX*RAD0 R2SPLU(NNY,MMY) = R2SPLU(NNY,MMY) + EIGV12*CGX*RAD2 R4SPLU(NNY,MMY) = R4SPLU(NNY,MMY) + EIGV12*CGX*RAD4 CS SPLU(NNY,MMY) = SPLU(NNY,MMY) + EIGV12*XD*REDMS*RAD0 62 CONTINUE IF (RJP .NE. RJ) GO TO 63 APLB(NNY,MMY) = APLB(NNY,MMY) + EIGV12*X*RAD0 R2JPLU(NNY,MMY) = R2JPLU(NNY,MMY) + EIGV12*X*RAD2 R4JPLU(NNY,MMY) = R4JPLU(NNY,MMY) + EIGV12*X*RAD4 63 CONTINUE R2DPLU(NNY,MMY) = R2DPLU(NNY,MMY) + EIGV12*XD*REDMD*RAD2 R4DPLU(NNY,MMY) = R4DPLU(NNY,MMY) + EIGV12*XD*REDMD*RAD4 34 CONTINUE 33 CONTINUE C IF (IREC .EQ. 1 .OR. ICORRS .EQ. 1) THEN DO 133 MMY=1,NU MY = LEVEL(MMY) IF (EIGVEC(IP,MY) .EQ. 0.0) GO TO 134 DO NY=1,IDIM EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY) IF (RLP .EQ. RL) & SPLR(NY,MMY) = SPLR(NY,MMY) + EIGV12*CGX*RAD0 IF (RJP .EQ. RJ) & JPLR(NY,MMY) = JPLR(NY,MMY) + EIGV12*X*RAD0 FMNR(NY,MMY) = FMNR(NY,MMY) + EIGV12*CGF*RAD2 ENDDO 134 CONTINUE 133 CONTINUE ENDIF C C***************************** END LOOPS OVER NILSSON ORBITALS ********* 42 CONTINUE IF (OMEGAP .NE. -OMEGA-2) GO TO 43 C C***** (-1) COMPONENT OF VECTOR OPERATORS *********************** C X = (-1.)**((RJ-OMEGA)/2)*0.5*SQRT((RJ-OMEGA)*(RJ+OMEGA+2.)) XD = (-1.)/SQRT(3.)*CLEBI(RJ,RJP,2,OMEGA,OMEGAP,-2) CGX=CLEBI(2*RL,1,RJP,OMEGAP+1,-1,OMEGAP)* *CLEBI(2*RL,1,RJ,-(OMEGA+1),1,-OMEGA)*(-1.)**((RJ-OMEGA)/2) IF (IABS(RNP-RN) .EQ. 2 .AND. IABS(RLP-RL) .LE. 2) THEN CGF = (-1)**((RJP-OMEGAP)/2)*FLOAT(RNP-RN)/SQRT(1.5)* &SQRT((RJP+1.0)/(RJ+1.0))*CLEBI(RJP,4,RJ,-OMEGAP,-2,OMEGA)* &CLEBI(RJP,4,RJ,-1,0,-1) ELSE CGF = 0.0 ENDIF C******************************* START LOOPS OVER NILSSON ORBITALS ***** C DO 35 MMY=1,NU DO 36 NNY=MMY,NU MY = LEVEL(MMY) NY = LEVEL(NNY) EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY) IF (RLP .NE. RL) GO TO 72 SMINUS(NNY,MMY) = SMINUS(NNY,MMY) + EIGV12*CGX*RAD0 R2SMIN(NNY,MMY) = R2SMIN(NNY,MMY) + EIGV12*CGX*RAD2 R4SMIN(NNY,MMY) = R4SMIN(NNY,MMY) + EIGV12*CGX*RAD4 CS SMIN(NNY,MMY) = SMIN(NNY,MMY) + EIGV12*XD*REDMS*RAD0 72 CONTINUE IF (RJP .NE. RJ) GO TO 73 C (NNY,MMY) = C (NNY,MMY) + EIGV12*X*RAD0 R2JMIN(NNY,MMY) = R2JMIN(NNY,MMY) + EIGV12*X*RAD2 R4JMIN(NNY,MMY) = R4JMIN(NNY,MMY) + EIGV12*X*RAD4 73 CONTINUE R2DMIN(NNY,MMY) = R2DMIN(NNY,MMY) + EIGV12*XD*REDMD*RAD2 R4DMIN(NNY,MMY) = R4DMIN(NNY,MMY) + EIGV12*XD*REDMD*RAD4 36 CONTINUE 35 CONTINUE C IF (IREC .EQ. 1 .OR. ICORRS .EQ. 1) THEN DO 135 MMY=1,NU MY = LEVEL(MMY) IF (EIGVEC(IP,MY) .EQ. 0.0) GO TO 136 DO NY=1,IDIM EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY) IF (RLP .EQ. RL) & SMNR(NY,MMY) = SMNR(NY,MMY) + EIGV12*CGX*RAD0 IF (RJP .EQ. RJ) & JMNR(NY,MMY) = JMNR(NY,MMY) + EIGV12*X*RAD0 FPLR(NY,MMY) = FPLR (NY,MMY) + EIGV12*CGF*RAD2 ENDDO 136 CONTINUE 135 CONTINUE ENDIF C C***************************** END LOOPS OVER NILSSON ORBITALS ********* 43 CONTINUE IF (OMEGAP .NE. OMEGA-4 .OR. ICORRS .EQ. 1) GO TO 44 C C***** -2 AND +2 COMPONENTS OF TENSOR OPERATORS ****************** C IF (RJ .NE. RJP .OR. RN .NE. RNP .OR. RL .NE. RLP) GO TO 44 IF ((RJ+OMEGA-2) .LT. 0) GO TO 44 X = SQRT((RJ+OMEGA)*(RJ-OMEGA+2.)*(RJ+OMEGA-2)*(RJ-OMEGA+4.))/4. C******************************** START LOOPS OVER NILSSON ORBITALS **** C DO 38 MMY=1,NU DO 37 NNY=MMY,NU MY = LEVEL(MMY) NY = LEVEL(NNY) EIGV21 = EIGVEC(IP,NY)*EIGVEC(I,MY) EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY) E(NNY,MMY) = E(NNY,MMY) + (EIGV12+EIGV21)*X 37 CONTINUE 38 CONTINUE C C********************************** END LOOPS OVER NILSSON ORBITALS **** 44 CONTINUE C C***** CORRECTION TO JZ DUE TO STRETCHING OF BASIS, NONAXIALITY C IF (ICORRS .NE. 1 .OR. IABS(RNP-RN) .NE. 2) GO TO 45 MU = OMEGA - OMEGAP IF (IABS(MU) .NE. 4 .OR. IABS(RLP-RL) .GT. 2) GO TO 45 CGF = FLOAT(RNP-RN)/SQRT(6.)*SQRT((RJP+1.0)/(RJ+1.0))* * CLEBI(RJP,4,RJ,-1,0,-1)*CLEBI(RJP,4,RJ,OMEGAP,MU,OMEGA) IF (MU .EQ. 4) CGF = -CGF C******************************** START LOOPS OVER NILSSON ORBITALS **** C DO 139 MMY=1,NU MY = LEVEL(MMY) IF (EIGVEC(IP,MY) .EQ. 0.0) GO TO 140 DO NY=1,IDIM EIGV12 = EIGVEC(IP,MY)*EIGVEC(I,NY) FZR(NY,MMY) = FZR(NY,MMY) + EIGV12*CGF*RAD2 ENDDO 140 CONTINUE 139 CONTINUE C C********************************** END LOOPS OVER NILSSON ORBITALS **** 45 CONTINUE C 12 CONTINUE 11 CONTINUE C C***************************** END LOOPS OVER BASIS FUNCTIONS ********* C**** CORRECT J+/J-/JZ DUE TO STRETCHED BASIS IF ICORR=1 C IF (ICORRS .EQ. 1) THEN GAMMAR = GAMMA*PI/180.0 OMX = 1.-2./3.*EPSP*COS(GAMMAR+2.*PI/3.) OMY = 1.-2./3.*EPSP*COS(GAMMAR-2.*PI/3.) OMZ = 1.-2./3.*EPSP*COS(GAMMAR) F1 = SQRT(OMZ/OMY) F2 = SQRT(OMZ/OMX) C FAC1 = (F1 + 1./F1 + 1./F2 + F2)/4. FAC2 = (F1 + 1./F1 - 1./F2 - F2)/4. FAC3 = (F1 - 1./F1 + 1./F2 - F2)/4. FAC4 = (F1 - 1./F1 - 1./F2 + F2)/4. FAC5 = (F2/F1 + F1/F2)/2. FAC6 = (F2/F1 - F1/F2)/2. C DO 125 MMY=1,NU DO 125 NY=1,IDIM JPLRS = JPLR(NY,MMY) JPLR(NY,MMY) = FAC1*JPLRS + (1.-FAC1)*SPLR(NY,MMY) & +FAC2*JMNR(NY,MMY) - FAC2 *SMNR(NY,MMY) & +FAC3*FPLR(NY,MMY) + FAC4 *FMNR(NY,MMY) JMNR(NY,MMY) = FAC1*JMNR(NY,MMY) + (1.-FAC1)*SMNR(NY,MMY) & +FAC2*JPLRS - FAC2 *SPLR(NY,MMY) & +FAC3*FMNR(NY,MMY) + FAC4 *FPLR(NY,MMY) JZR (NY,MMY) = FAC5*JZR(NY,MMY) + (1.-FAC5)* SZR(NY,MMY) & +FAC6*FZR(NY,MMY) 125 CONTINUE DO 126 MMY=1,NU DO 126 NNY=1,NU D (NNY,MMY) = 0.0 E (NNY,MMY) = 0.0 JZ2(NNY,MMY) = 0.0 126 CONTINUE DO 129 MMY=1,NU MY = LEVEL(MMY) DO 128 NNY=MMY,NU NY = LEVEL(NNY) APLB(NNY,MMY) = JPLR(NY,MMY) C (NNY,MMY) = JMNR(NY,MMY) JZ (NNY,MMY) = JZR (NY,MMY) DO 127 IK=1,IDIM D (NNY,MMY) = D(NNY,MMY) + (JPLR(IK,MMY)*JPLR(IK,NNY) & + JMNR(IK,MMY)*JMNR(IK,NNY))/2. E (NNY,MMY) = E(NNY,MMY) + JPLR(IK,MMY)*JMNR(IK,NNY) & + JMNR(IK,MMY)*JPLR(IK,NNY) JZ2(NNY,MMY) = JZ2(NNY,MMY) + JZR(IK,MMY)*JZR(IK,NNY) 127 CONTINUE 128 CONTINUE 129 CONTINUE ENDIF C************** START LOOPS OVER NILSSON ORBITALS FOR PRINTING ******** C DO 1 MMY=1,NU MY = LEVEL(MMY) DO 1 NNY=MMY,NU AVR2S= $(ABS(R2SZ(NNY,MMY))+ABS(R2SPLU(NNY,MMY))+ABS(R2SMIN(NNY,MMY)))/ $(ABS(SNOLL(NNY,MMY))+ABS(SPLUS(NNY,MMY))+ABS(SMINUS(NNY,MMY))) AVR4S= $(ABS(R4SZ(NNY,MMY))+ABS(R4SPLU(NNY,MMY))+ABS(R4SMIN(NNY,MMY)))/ $(ABS(SNOLL(NNY,MMY))+ABS(SPLUS(NNY,MMY))+ABS(SMINUS(NNY,MMY))) NY = LEVEL(NNY) IF(MY.EQ.NY) write(6,200) VALU(NY),KATT,APLB(NNY,MMY), $C(NNY,MMY),D(NNY,MMY),E(NNY,MMY),JZ(NNY,MMY), CC $JZ2(NNY,MMY),SPLUS(NNY,MMY),SMINUS(NNY,MMY),AVR2S,AVR4S, $SNOLL(NNY,MMY),SPLUS(NNY,MMY),SMINUS(NNY,MMY),AVR2S,AVR4S, $R2DZ(NNY,MMY),R2DPLU(NNY,MMY),R2DMIN(NNY,MMY), $Q20S(NNY,MMY),Q22S(NNY,MMY),Q20(NNY,MMY),Q22(NNY,MMY), $RHO2(NNY,MMY),R2(NNY,MMY) IF(NY.GT.MY) write(6,300) VALU(NY),KATT,VALU(MY),KATT, $APLB(NNY,MMY),C(NNY,MMY),D(NNY,MMY),E(NNY,MMY),JZ(NNY,MMY), CC $JZ2(NNY,MMY),SPLUS(NNY,MMY),SMINUS(NNY,MMY),AVR2S,AVR4S, $SNOLL(NNY,MMY),SPLUS(NNY,MMY),SMINUS(NNY,MMY),AVR2S,AVR4S, $R2DZ(NNY,MMY),R2DPLU(NNY,MMY),R2DMIN(NNY,MMY), $Q20S(NNY,MMY),Q22S(NNY,MMY),Q20(NNY,MMY),Q22(NNY,MMY), $RHO2(NNY,MMY),R2(NNY,MMY) CC PRINT 423,APLB(NNY,MMY), , CC $C(NNY,MMY),D(NNY,MMY),E(NNY,MMY),JZ(NNY,MMY),JZ2(NNY,MMY) 423 FORMAT(' ',10F9.4) 200 FORMAT(/,1X,F5.2,A1,' MY=NY',2F7.3,2F7.2,15F6.2) 300 FORMAT(1X,2(F5.2,A1),2F7.3,2F7.2,15F6.2) DO 3 K=1,28 EQUIVA(MMY,NNY,K) = EQUIVA(NNY,MMY,K) 3 CONTINUE 1 CONTINUE C C************* END LOOPS OVER NILSSON ORBITALS FOR PRINTING *********** C WRITE(17) KATT,IREC,NU,(LEVEL(J),J=1,NU), $(((EQUIVA(MYP,MY,I),MY=1,NU),MYP=1,NU),I=1,28) C IF (ISTRCH .EQ. 1) GO TO 9 WRITE(17) IDIM,(KV1(I),KV3(I),KV4(I),I=1,IDIM) DO 2 MMY=1,NU MY=LEVEL(MMY) WRITE(17) (EIGVEC(I,MY),I=1,IDIM) 2 CONTINUE 9 CONTINUE C IF (IREC .EQ. 1) WRITE(17) ((JPLR(NY,MMY),MMY=1,NU),NY=1,IDIM), & ((JMNR(NY,MMY),MMY=1,NU),NY=1,IDIM), & ((JZR(NY,MMY),MMY=1,NU),NY=1,IDIM) RETURN END C C####################################################################### C FUNCTION REDMD1(L,S,J,LP,SP,JP) C C**** REDUCED MATRIX ELEMENT FOR TENSOR OPERATOR D1 OF RANK ONE FORMED * C**** FROM THE COUPLING OF C2 AND S; D1 = -SQRT(10)/2*(C2*S)1 * C**** DEFINITIONS FROM EDMOMDS: ANGULAR MOMENTUM IN QUANT MECHANICS * C C**** L S,J,LP,SP,JP HAVE TWICE THEIR PHYSICAL VALUES, I.E. S=SP=1 INTEGER L,S,J,LP,SP,JP REAL NINEJ IF ((LP+L) .GE. 4 .AND. ABS(LP-L) .LE. 4 * .AND.(JP+J) .GE. 2 .AND. ABS(JP-J) .LE. 2) GO TO 1 REDMD1 = 0. RETURN 1 CONTINUE C C**** SIX-J SYMBOLS FROM CLOSED FORMULAS IN EDMONDS, (6.3.3) AND (6.3.4) C IF (JP .EQ. (LP-1)) THEN XT = (J+LP+5.)*(LP+3.-J)/4. XN = LP*(LP+1.)*12. SIXJA = (-1.)**((J+LP+3)/2)*SQRT(XT/XN) ELSEIF (JP .EQ. (LP+1)) THEN XT = (J+3.-LP)*(J+LP-1.)/4. XN = (LP+1.)*(LP+2.)*12. SIXJA = (-1.)**((J+LP+3)/2)*SQRT(XT/XN) ELSE STOP 6 ENDIF C IF (J .EQ. (L-1)) THEN XT = (LP+L+6.)*(L+4.-LP)/4. XN = L*(L+1.)*20. SIXJB = (-1.)**((L+LP)/2)*SQRT(XT/XN) ELSEIF (J .EQ. (L+1)) THEN XT = (LP+4.-L)*(LP+L-2.)/4. XN = (L+1.)*(L+2.)*20. SIXJB = (-1.)**((L+LP)/2)*SQRT(XT/XN) ELSE STOP 66 ENDIF C SIXJC = SQRT(3.)/6 C C******************** NINE-J SYMBOL CALCULATED FROM EDMONDS (6.4.17) NINEJ = -SIXJA*SIXJB/SIXJC/3. C C*** REDUCED MATRIX ELEMENTS OF C2 =SQRT(4*PI/5)*Y2 AND S C REDMC2 = SQRT(L+1.)*CLEBI(L,4,LP,0,0,0) REDMS1 = SQRT(1.5) C C************************ REDUCED MATRIX ELEMENT OF PRODUCT OPERATOR D1, C ********************* SEE EDMOMDS (7.1.5) ***** C = -SQRT(10.)/2. REDMD1 = C*REDMC2*REDMS1*SQRT((JP+1.)*(J+1.)*3.)*NINEJ C RETURN END C C####################################################################### C SUBROUTINE ENERLJ(OMOM,SNI,NOYES) C C***** CALCULATION OF SINGLE PARTICLE ENERGIES, WAVE FUNCTIONS C AND MATRIX ELEMENTS C ****** MAX 10 N-SHELLS ******************************************* C CV COMMON/LEVEL/NU,LEVEL(15) CV COMMON /LU/ LU COMMON /COUPL/ NUU COMMON /CMN2/ EPS1,EPS2,EPS3,EPS4,EPS5,EPS6,EPS22 COMMON/CMN7/ EPSP,GAMMA COMMON /QSTR/ Q COMMON/CMN1/ KAPPA,DUMMY,MY,L2F(18) COMMON/E/ EMIN,EMAX COMMON/KVANT/ KV1(161),KV2(161),KV3(161),KV4(161) C COMMON/EIG/ EIGVEC(203,70),VALU(203) COMMON/EIG/ EIGVEC(161,161),VALU(161) COMMON/DIM/ IDIM,KATT COMMON /NUCL/ Z,A COMMON /QMOM/ Q20S(70),Q22S(70),Q20(70),Q22(70),RHO2(70) $,R2(70),RJZ(70) COMMON/EQFAC/ EFAC,QFAC,ISO COMMON /PAR/ IPAR COMMON /NVEC/ NVEC,NQVEC COMMON/OMROT/ OMROT COMMON/FACTOR/ FAC,HAC COMMON/XKAMY/ XKAPPA(13,7),XMY(13,7),AVKAMY(13) COMMON/ITRANS/ ITRANS,NSHELL COMMON/BAS/ ISTRCH,ICORR,IREC ! DIMENSION AM(13041) DIMENSION AM(0:25425) DIMENSION E(572),EIGVAL(70) DIMENSION KVA(161),LEVNUM(280) DIMENSION EB(280) INTEGER SNI,SN,Q,RN,RL,RJ,RNP,RLP,RJP,OMEGA,OMEGAP DIMENSION DIM(2),TIME(2),TID(2) INTEGER DIM,DIME INTEGER Z,A,Z2 INTEGER QQ REAL KAPPA,MY COMMON/Q2022/ QQ20(280),QQ22(280) $,RR2(280),RRJZ(280) CHARACTER IPLUS*4,IMINUS*4,KATT*4,IPAR*4,KVANT(2)*4,KVB(280)*4 CHARACTER KV(728)*4 DATA MATDIM/161/ DATA IPLUS,IMINUS/'+','-'/ C DATA DELTA/0.0316/ DATA DELTA/0.000316/ DATA PI/3.1415926536/ SN=SNI IF (SNI .LE. 10) GO TO 200 199 FORMAT(1X,' N IS GREATER THAN MAX VALUE 10 ... ENERLJ ') write(6,199) STOP 200 CONTINUE 100 FORMAT(/4X,'OMROT=',F5.3) GAMMAR=GAMMA*PI/180.0 OMX=1.-2./3.*EPSP*COS(GAMMAR+2.*PI/3.) OMY=1.-2./3.*EPSP*COS(GAMMAR-2.*PI/3.) OMZ=1.-2./3.*EPSP*COS(GAMMAR) FAC=0.5*(SQRT(OMY/OMX)+SQRT(OMX/OMY))-1. HAC=0.5*(SQRT(OMY/OMX)-SQRT(OMX/OMY)) OMROT=OMROT/OMOM Q=0 QQ=0 IR=1 IF(ABS(OMROT).GT.0.0001) IR=2 C C LOOP IN SIGNATURES (OMROT .GT. 0) C DO 201 ISYM=1,IR IS=ISYM IB=0 C C LOOP IN PARITIES C DO 201 NUM=1,2 7132 CONTINUE MATNUM=IABS(NUM-3*IB) SN=SNI+1-MATNUM IF((SN/2)*2.EQ.SN) GO TO 7133 IDIM=(SN+1)*(2*SN+7)*(SN+3)/24 KATT=IMINUS GO TO 7135 7133 CONTINUE IDIM=(SN+2)*(2*SN+3)*(SN+4)/24 KATT=IPLUS 7134 CONTINUE 7135 CONTINUE IF(NUM.EQ.1.AND.KATT.EQ.IPAR.AND.MATNUM.EQ.1) THEN IB=1 GO TO 7132 ENDIF IF(NOYES.EQ.0) GO TO 101 IF (NUM .EQ. 2 .OR. IS .EQ. 2) write(6,397) 397 FORMAT(//) 101 CONTINUE DIM(MATNUM)=IDIM KVANT(MATNUM)=KATT KLAFS=IDIM*(IDIM+1)/2 DO 2 KK=1,KLAFS AM(KK)=0.0 2 CONTINUE DO 1 I=1,IDIM CALL QNUMLJ(I,RN,RL,RJ,OMEGA,SN) KV1(I)=RN KV2(I)=RL KV3(I)=RJ KV4(I)=OMEGA*(3-2*IS) 1 CONTINUE ! IRK=MSLEFT(XX) KK=0 C C LOOP OVER BASIS FUNCTIONS, CREATION OF HAM MATR ELEMENTS C DO 210 I=1,IDIM RNP=KV1(I) RLP=KV2(I) RJP=KV3(I) OMEGAP=KV4(I) C C SECOND LOOP OVER BASIS FUNCTIONS (J .GE. I) C DO 210 J=I,IDIM RN=KV1(J) RL=KV2(J) RJ=KV3(J) OMEGA=KV4(J) KK=KK+1 IF(RNP.GT.RN+NUU.OR.OMEGAP.GT.OMEGA+8.OR. $ OMEGAP.LT.OMEGA-8) GO TO 303 E20=0.0 E22=0.0 E40=0.0 E42=0.0 E44=0.0 E60=0.0 ELZ=0.0 END=0.0 IF(OMEGAP.EQ.OMEGA+8.OR.OMEGAP.EQ.OMEGA-8) GO TO 305 IF(OMEGAP.EQ.OMEGA+4.OR.OMEGAP.EQ.OMEGA-4) GO TO 306 IF(RNP.EQ.RN) E20=-2.0*EPS2/3.0*CLEBI(RJ,4,RJP,-1,0,-1)* $ CLEBI(RJ,4,RJP,OMEGA,0,OMEGAP) IF (ISTRCH .EQ. 1) GO TO 901 IF(RNP.EQ.RN+2) E20=-2.0*EPS2/3.0*CLEBI(RJ,4,RJP,-1,0,-1)* $ CLEBI(RJ,4,RJP,OMEGA,0,OMEGAP) 901 CONTINUE C E40=EPS4*(COS(1.5*GAMMAR)**2+0.375*SIN(1.5*GAMMAR)**2)* C**** EPS4 FOR GAMMA .NE. 0 OR 60 AS DEFINED E.G. IN NUCL PHYS C A436 (1985) 14 ###### E40, E42 AND E44 #### C SEE ALSO FUNCTION FVCON ########### C C**** ALSO NOTE TYPO IN THAT PAPER : V4 TERM OMITS THE FACTOR C 2*SQRT(9/4*PI), REQUIRED TO GIVE 2*EPS4*P4 FOR GAMMA=0 C AS IS CONVENTIONAL, E.G. NUCL. PHYS. A129 (1969) 445 C E40=EPS4*((5.*COS(GAMMAR)**2+1.)/6.)* $ CLEBI(RJ,8,RJP,-1,0,-1)*CLEBI(RJ,8,RJP,OMEGA,0,OMEGAP) E60=EPS6*CLEBI(RJ,12,RJP,-1,0,-1)* $ CLEBI(RJ,12,RJP,OMEGA,0,OMEGAP) IF(RNP.EQ.RN.AND.RLP.EQ.RL) ELZ=-OMROT*FAC* $ ((OMEGA-1.)/2.*CLEBI(2*RL,1,RJP,OMEGA-1,1,OMEGA)* $ CLEBI(2*RL,1,RJ,OMEGA-1,1,OMEGA)+ $ (OMEGA+1.)/2.*CLEBI(2*RL,1,RJP,OMEGA+1,-1,OMEGA)* $ CLEBI(2*RL,1,RJ,OMEGA+1,-1,OMEGA)) GO TO 304 306 CONTINUE IF(RNP.EQ.RN) E22=EPS22/SQRT(6.0)*CLEBI(RJ,4,RJP,-1,0,-1)* $ (CLEBI(RJ,4,RJP,OMEGA,4,OMEGAP)+ $ CLEBI(RJ,4,RJP,OMEGA,-4,OMEGAP)) IF (ISTRCH .EQ. 1) GO TO 902 IF(RNP.EQ.RN+2) E22=EPS22/SQRT(6.0)*CLEBI(RJ,4,RJP,-1,0,-1)* $ (CLEBI(RJ,4,RJP,OMEGA,4,OMEGAP)+ $ CLEBI(RJ,4,RJP,OMEGA,-4,OMEGAP)) 902 CONTINUE IF(RNP.EQ.RN+2) $ END= OMROT/SQRT(1.5)*HAC*CLEBI(RJ,4,RJP,-1,0,-1)* $ (CLEBI(RJ,4,RJP,OMEGA,4,OMEGAP)- $ CLEBI(RJ,4,RJP,OMEGA,-4,OMEGAP)) C E42=-EPS4*SQRT(10.)/8.*SIN(1.5*GAMMAR)**2* E42=-EPS4*SQRT(30.)/12.*SIN(2.*GAMMAR)* $ CLEBI(RJ,8,RJP,-1,0,-1)* $ (CLEBI(RJ,8,RJP,OMEGA,4,OMEGAP)+ $ CLEBI(RJ,8,RJP,OMEGA,-4,OMEGAP)) GO TO 304 305 CONTINUE C E44=EPS4*SQRT(70.)/16.*SIN(1.5*GAMMAR)**2* E44=EPS4*SQRT(70.)/12.*SIN(GAMMAR)**2* $ CLEBI(RJ,8,RJP,-1,0,-1)* $ (CLEBI(RJ,8,RJP,OMEGA,8,OMEGAP)+ $ CLEBI(RJ,8,RJP,OMEGA,-8,OMEGAP)) 304 CONTINUE CALL MRHOT(RNP,RLP,0,RN,RL,0,2,RME) ARJ=FLOAT(RJ) ARJP = FLOAT(RJP) ARL = FLOAT(RL) AM(KK) = RME*SQRT((ARJ+1.)/(ARJP+1.))* * (E20+E22+E40+E42+E44+E60+END)+ELZ IF(I.EQ.J) AM(KK)=AM(KK)+FLOAT(RN)+1.5-XKAPPA(RN+1,RL/2+1)/ * OMOM*(ARJ*(ARJ+2.)/4.-ARL*(ARL+1.)-0.75)- * XKAPPA(RN+1,RL/2+1)/OMOM*XMY(RN+1,RL/2+1)*ARL*(ARL+1.)+ * AVKAMY(RN+1)/OMOM*FLOAT(L2F(RN+1))- * OMROT*FLOAT(OMEGA)/2. C C AM IS SYMMETRICAL MATRIX TO BE DIAGONALIZED C (ONLY STORED FOR (J .GE. I) C AM(KK)=-AM(KK) 303 CONTINUE 210 CONTINUE C C c JRK=MSLEFT(XX) TIME(MATNUM)=(IRK-JRK)*0.001 IDI=IDIM NQVEC=MIN0(IDI,70) IF (IREC .EQ. 1) THEN NVEC=IDIM ELSE NVEC=NQVEC ENDIF IF(NOYES.EQ.0) GO TO 11 C***** PRINT AV BASFUNKTIONER DIME=DIM(MATNUM) IMIN=(DIME+4)/5 KVR1=3-2*ISYM write(6,9) KVR1 9 FORMAT(4X,'SYMMETRY R1=',I2,/) write(6,7) KVANT(MATNUM) 7 FORMAT(/,' BASIS FUNCTIONS /N,L,J,OMEGA>, PARITY ',A1,//) DO 113 K=1,IMIN IMAX = K + 4*IMIN IF (IMAX .GT. DIME) IMAX = IMAX - IMIN write(6,16)(J,KV1(J),KV2(J),KV3(J),KV4(J),J=K,IMAX,IMIN) 16 FORMAT(1X,5(I3,1X,1H/,2I2,I3,2H/2,I3,2H/2,1H>,3X)) 113 CONTINUE 11 CONTINUE IF(NOYES.EQ.0) NVEC=0 ! IRK=MSLEFT(XX) C C DIAGONALIZATION C CALL BIGMAX(AM,VALU,IDIM,IDI,NVEC,EIGVEC,MATDIM) C c JRK=MSLEFT(XX) TID(MATNUM)=(IRK-JRK)*0.001 DO 405 J=1,IDI VALU(J)=-VALU(J) 405 CONTINUE IF(NOYES.EQ.0) GO TO 12 C***** PRINT AV VAVE-FUNCTIONS write(6,398) KATT,EMIN,EMAX 398 FORMAT(//,' SINGLE PARTICLE (PARITY',A1,') WAVE FUNCTIONS ', * 'IN THE ENERGY INTERVAL',F4.1,3H)/11X)) 600 FORMAT(/1X,2HE=,F6.4,2X,5(7(F9.6,1H/,I3,1H>)/11X)) 403 CONTINUE 404 CONTINUE CV CV PRINT OF SINGLE-PARTICLE WAVEFUNCTIONS FOR VMI CV IF (KATT .NE. IPAR) GOTO 801 WRITE(16) IDIM DO 800 II=1,NU IIII=LEVEL(II) IF (IPAR .EQ. IPLUS) WRITE(16) (KV1(III)/2+1,KV2(III), $KV3(III),KV4(III),EIGVEC(III,IIII),III=1,IDIM) IF (IPAR .EQ. IMINUS) WRITE(16) ((KV1(III)+1)/2,KV2(III), $KV3(III),KV4(III),EIGVEC(III,IIII),III=1,IDIM) 800 CONTINUE 801 CONTINUE CV CV C*** TRANSFORM TO /N,L,LAMBDA,SIGMA> BASIS (NECESSARY FOR Vpn CALCN) CV ITEST=NUU*(ABS(EPS4)+ABS(EPS6)) IF (ITEST .EQ. 0) THEN NMAX=NSHELL IDIM1=((NSHELL+1)*(NSHELL+2))/2 ELSE IDIM1=IDIM NMAX=SN ENDIF IF (KATT .EQ. IPAR .AND. ITRANS .NE. 0) & CALL TRANSF(NMAX,IDIM1) C C***** CALCULATION OF QUADRUPOLE MATRIX ELEMENTS C SUBROUTINES Q20Q22 AND QTME C ! IRKA=MSLEFT(XX) CALL Q20Q22 ! IRK = MSLEFT(XX) IF(KATT .EQ. IPAR .AND. OMROT .LT. 0.0001) THEN CALL QTME ENDIF c JRKA=MSLEFT(XX) SEC = (IRKA-IRK)*0.001 SEC1 = (IRK-JRKA)*0.001 write(6,55)SEC,SEC1 55 FORMAT(/,1X,8H Q20Q22:,F6.2,13H SEC QTME:,F6.2,4H SEC) CT WRITE(1,*)'NQVEC,QQ',NQVEC,QQ DO 204 I=1,NQVEC IQ=I+QQ EB(IQ)=VALU(I) KVB(IQ)=KATT QQ20(IQ)=Q20(I) QQ22(IQ)=Q22(I) RR2(IQ)=R2(I) RRJZ(IQ)=RJZ(I) 204 CONTINUE QQ=QQ+NQVEC IF(NUM.NE.2.) GO TO 406 IF(IS.NE.IR) GO TO 406 CT WRITE(1,*)'NQVEC,QQ,NUM,IS,IR,Z',NQVEC,QQ,NUM,IS,IR,Z CT WRITE(1,*)'EB,KVB',(EB(I),I=1,QQ),(KVB(I),I=1,QQ) C C SORTING ACCORDING TO SINGLE-PARTICLE ENERGIES C CALL SORT4(EB,KVB,QQ20,QQ22,RR2,RRJZ,QQ) LEVPOS=0 LEVNEG=0 DO 309 I=1,QQ IF (KVB(I).EQ.IPLUS) THEN LEVPOS=LEVPOS+1 LEVNUM(I)=LEVPOS ELSE LEVNEG=LEVNEG+1 LEVNUM(I)=LEVNEG ENDIF 309 CONTINUE C CT WRITE(1,*)'EB,KVB',(EB(I),I=1,QQ),(KVB(I),I=1,QQ) CT WRITE(1,*)'A',A N=A-Z C********QUADRUPOLE MOMENTS OF NUCLEUS Q0=0.0 Q2=0.0 IF (ISO .EQ. 1) THEN Z2=Z/2 ELSEIF (ISO .EQ. 2) THEN Z2=N/2 ENDIF IZ=Z2*IS DO 205 I=1,IZ Q0=Q0+QQ20(I)*2./FLOAT(IS) Q2=Q2+QQ22(I)*2./FLOAT(IS) 205 CONTINUE GAM=180./PI*ATAN(-SQRT(2.)*Q2/Q0) IF (ISO .EQ. 1) THEN write(6,35) Q0,Q2,GAM,Z,A ELSEIF (ISO .EQ. 2) THEN write(6,35) Q0,Q2,GAM,N,A ENDIF 35 FORMAT(/1X,3HQ0=,F8.3,' HBAR/MW0',5X,3HQ2=,F8.3,' HBAR/MW0', * 5X,6HGAMMA=,F5.1,' DEG FOR',I3,' PARTICLES; A=',I3,/) 406 CONTINUE IF(KATT.NE.IPAR) GO TO 202 C C***** CALC OF GENERALIZED DECOUPLING FACTORS, SUBROUTINE DECOUP ! IRKA = MSLEFT(XX) IF(KATT.EQ.IPAR.AND.OMROT.LT.0.0001) THEN write(6,399) CALL DECOUP ENDIF ! JRKA=MSLEFT(XX) SEC=(IRKA-JRKA)*0.001 write(6,65) SEC 65 FORMAT(/,1X,'DECOUP:',F6.2,' SEC') 399 FORMAT(1H1) 19 CONTINUE 202 CONTINUE 12 CONTINUE DO 203 I=1,IDI IQ=I+Q E(IQ)=VALU(I) KV(IQ)=KATT 203 CONTINUE Q=Q+IDI 201 CONTINUE IF(NOYES.EQ.0) GO TO 17 write(6,8) 8 FORMAT(///,1X,2(' # ENERGY +/-(#) ', *' ',4X),/) IMIN = (QQ+1)/2 IF (IMIN .GT. 40) IMIN = 40 DO 111 K=1,IMIN IMAX = K + IMIN IF (IMAX .GT. QQ) IMAX = IMIN write(6,14) (J,EB(J),KVB(J),LEVNUM(J),QQ20(J),QQ22(J),RR2(J), * RRJZ(J),J = K,IMAX,IMIN) C IF (MOD(K,5) .EQ. 0) PRINT 24 24 FORMAT(/) 14 FORMAT(1X,2(I3,F9.4,2X,A1,'(#',I2,')',4F9.4,3X)) 111 CONTINUE WRITE(17) EPSP,GAMMA,EPS4,EPS6,OMOM,QQ, $(EB(I),I=1,QQ),(KVB(I),I=1,QQ), $(QQ20(I),I=1,QQ),(QQ22(I),I=1,QQ),(RR2(I),I=1,QQ),(RRJZ(I),I=1,QQ) 17 CONTINUE C C***** PRINT OF S P ENERGIES AND QUANTUM NUMBERS - NOYES = 0 IF(NOYES.NE.0) GO TO 114 CALL SORT(E,KV,Q) IF(NOYES.NE.0) write(6,399) write(6,6) 6 FORMAT(1X,' LEVEL,ENERGY,PARITY',/) IMIN=MIN0(44,Q) DO 112 K=1,IMIN IMAX=K+176 IMAX=MIN0(IMAX,Q)/44*44+K IF(IMAX.GT.Q) IMAX=IMAX-44 write(6,15) (J,E(J),KV(J),J=K,IMAX,44) 15 FORMAT(1X,5(I7,F9.4,4X,A1,3X)) 112 CONTINUE 114 CONTINUE C***** PRINT AV TIDER. write(6,40) 40 FORMAT(//) DO 10 I=1,2 write(6,30) KVANT(I),DIM(I),TIME(I),TID(I) 30 FORMAT(4X,9H PARITY: ,A1,15H DIMENSION:,I4,12H MATRIX: $ ,F6.1,4H SEC,21H DIAGONALISATION:,F6.1,4H SEC) 10 CONTINUE RETURN END C C*********************************************************************** C SUBROUTINE QNUMLL(NUM,NP,LP,LAMP,ISIGP,NMAX) C C SETS UP /N,L,LAMBDA,SIGMA> BASIS FUNCTIONS C (ONE FOR EACH CALL) BUT WITH THE RESTRICTION C OMEGA = LAMBDA + SIGMA =..., -3/2, +1/2, +5/2, ... C ON OUTPUT: NP=N, LP=L, LAMP=LAMBDA, ISIGP=2*SIGMA C INPUT: NUM=NUMBERING ON BASIS FUNCTION C NMAX=MAXIMUM N-SHELL (PARITY=(-1)**NMAX) C NUMP=NUM NTRY=NMAX 1 CONTINUE NTEST=(NTRY+1)*(NTRY+2)/2 IF(NUMP .LE. NTEST) GO TO 2 NUMP=NUMP-NTEST NTRY=NTRY-2 GO TO 1 2 CONTINUE NP=NTRY C LTRY=NP 3 CONTINUE LTEST=2*LTRY +1 IF(NUMP .LE. LTEST) GO TO 4 NUMP=NUMP-LTEST LTRY=LTRY-2 GO TO 3 4 CONTINUE LP=LTRY LAMP=LP+1-NUMP ISIGP=-1 IF((LAMP/2)*2 .EQ. LAMP) ISIGP=1 RETURN END C C*********************************************************************** C SUBROUTINE TRANSF(NMAX,IDIM1) C C TRANSFORMS S.P. STATES FROM THE /N,L,J,OMEGA> BASIS C TO THE /N,L,LAMBDA,SIGMA> BASIS (CALLED FROM ENERLJ) C COMMON/KVANT/ KV1(161),KV2(161),KV3(161),KV4(161) COMMON/EIG/ EIGVEC(161,161),VALU(161) COMMON/DIM/ IDIM,KATT COMMON/LEVEL/ NU,LEVEL(15) COMMON/ITRANS/ ITRANS,NSHELL DIMENSION KVL1(161),KVL2(161),KVL3(161),KVL4(161),KDOM(161) DIMENSION T1(161,15),AV(15),AM(161) INTEGER OMEGA,OMEGAP CHARACTER KATT*4 DATA DELTA/0.00316/ DO 1 I=1,IDIM1 CALL QNUMLL(I,N,L,LAM,ISIG,NMAX) KVL1(I)=N KVL2(I)=L KVL3(I)=LAM KVL4(I)=ISIG 1 CONTINUE DO 3 I=1,IDIM1 DO 3 JJ=1,NU T1(I,JJ)=0.0 3 CONTINUE DO 150 IJ=1,NU JJ=LEVEL(IJ) AV(IJ)=VALU(JJ) DO 140 I=1,IDIM IF (ABS(EIGVEC(I,JJ)) .LT. 0.00001) GO TO 139 N = KV1(I) L = KV2(I) J = KV3(I) OMEGA = KV4(I) DO 130 K=1,IDIM1 NP = KVL1(K) LP = KVL2(K) LAM = KVL3(K) ISIG = KVL4(K) OMEGAP = 2*LAM + ISIG IF (NP.NE.N .OR. LP.NE.L .OR. OMEGAP.NE.OMEGA) & GO TO 129 CG1 = CLEBI(2*L,1,J,2*LAM,ISIG,OMEGA) T1(K,IJ) = T1(K,IJ) + CG1*EIGVEC(I,JJ) 129 CONTINUE 130 CONTINUE 139 CONTINUE 140 CONTINUE 149 CONTINUE 150 CONTINUE IF (ITRANS .EQ. 1) GO TO 406 C***** PRINT AV BASFUNKTIONER write(6,398) 398 FORMAT(//5X,'SELECTED SINGLE-PARTICLE STATES EXPANDED IN ', & '/N,L,LAMBDA,SIGMA> BASIS ',//) IMIN=(IDIM1+4)/5 write(6,7) 7 FORMAT(/,' BASIS FUNCTIONS /N,L,LAMBDA,SIGMA> ',//) DO 113 K=1,IMIN IMAX = K + 4*IMIN IF (IMAX .GT. IDIM1) IMAX = IMAX - IMIN write(6,16)(J,KVL1(J),KVL2(J),KVL3(J),KVL4(J),J=K,IMAX, & IMIN) 16 FORMAT(1X,5(I3,1X,1H/,4I3,2H/2,1H>,3X)) 113 CONTINUE C***** PRINT OUT WAVE-FUNCTIONS K=0 DO 404 J=1,NU KK=0 DO 402 K=1,IDIM1 IF(ABS(T1(K,J)) .LT. DELTA) GO TO 401 KK=KK+1 AM(KK)=T1(K,J) KDOM(KK)=K 401 CONTINUE 402 CONTINUE KKMAX=KK write(6,600) AV(J), (AM(KK),KDOM(KK),KK=1,KKMAX) 600 FORMAT(/2X,'E=',F6.4,1X,5(7(F9.6,1H/,I3,1H>)/11X)) 404 CONTINUE 406 CONTINUE C*** WRITE TO FILE2 WRITE(2,*) IDIM1 DO 410 J=1,IDIM1 WRITE(2,*) KVL1(J),KVL2(J),KVL3(J),KVL4(J) 410 CONTINUE WRITE(2,*) NU WRITE(2,*) (LEVEL(I),I=1,NU) DO 420 K=1,NU WRITE(2,*) (T1(J,K),J=1,IDIM1) 420 CONTINUE RETURN END C C####################################################################### C FUNCTION RMATEL(N,L,NP,LP,LAMBDA) C C**** N, L, NP, LP ARE INTEGERS WITH TWICE THEIR PHYSICAL VALUE. C RADIAL MATRIX ELEMENTS OF R**LAMBDA (R DIMENSIONLESS) C (LAMBDA = 2,3; CLOSED FORMULAS USED) C N/2 ,NP/2 - NUMBER OF OSCILLATOR QUANTA C L/2, LP/2 - ORBITAL ANGULAR MOMENTUM C RMATEL=0. IF (LAMBDA.EQ.2) GO TO 200 IF (LAMBDA.EQ.3) GO TO 300 write(6,100)LAMBDA 100 FORMAT(/,'RADIAL MATRIX ELEMENTS OF R',2H**,I2, *' CAN BE SUPPLIED BY MRHOT') STOP 200 IF (NP-N) 210,220,230 210 IF (NP.NE.N-4) GO TO 9 IF (LP-L) 211,212,213 211 IF (LP.NE.L-4) GO TO 9 RMATEL=SQRT(FLOAT((N+L+2)*(N+L-2)))/4. GO TO 9 212 RMATEL=SQRT(FLOAT((N-L)*(N+L+2)))/4. GO TO 9 213 IF (LP.NE.L+4) GO TO 9 RMATEL=SQRT(FLOAT((N-L)*(N-L-4)))/4. GO TO 9 220 IF (LP-L) 221,222,223 221 IF (LP.NE.L-4) GO TO 9 RMATEL=SQRT(FLOAT((N-L+4)*(N+L+2)))/2. GO TO 9 222 RMATEL=FLOAT(N+3)/2. GO TO 9 223 IF (LP.NE.L +4) GO TO 9 RMATEL=SQRT(FLOAT((N-L)*(N+L+6)))/2. GO TO 9 230 IF (NP.NE.N+4) GO TO 9 IF (LP-L) 231,232,233 231 IF (LP.NE.L-4) GO TO 9 RMATEL=SQRT(FLOAT((N-L+4)*(N-L+8)))/4. GO TO 9 232 RMATEL=SQRT(FLOAT((N-L+4)*(N+L+6)))/4. GO TO 9 233 IF (LP.NE.L+4) GO TO 9 RMATEL=SQRT(FLOAT((N+L+6)*(N+L+10)))/4. GO TO 9 300 ND=NP-N LD=LP-L NA=IABS(ND) LA=IABS(LD) IF (NA.EQ.2) GO TO 301 IF (NA.EQ.6) GO TO 303 GO TO 9 301 IF (LA.EQ.2) GO TO 310 IF (LA.EQ.6) GO TO 320 GO TO 9 310 IF (LD*ND) 311,9,312 311 RMATEL=FLOAT(3*N+L+10+ND)*SQRT(FLOAT(N-L+2+ND))/8. GO TO 9 312 RMATEL= FLOAT(3*N-L+8+ND)*SQRT(FLOAT(N+L+4+ND))/8. GO TO 9 320 IF (LD*ND) 321,9,322 321 RMATEL= SQRT(FLOAT((N+L+4-ND)*(N-L+2*ND)*(N-L+4+2*ND)))*3./8. GO TO 9 322 RMATEL=SQRT(FLOAT((N+L+2+2*ND)*(N+L+6+2*ND)*(N-L+2-ND)))*3./8. GO TO 9 303 IF (LA.EQ.2) GO TO 330 IF (LA.EQ.6) GO TO 340 GO TO 9 330 ND=ND/3 IF (LD*ND) 331,9,332 331 RMATEL= SQRT(FLOAT((N+L+4+ND)*(N-L+2*ND)*(N-L+4+2*ND)))/8. GO TO 9 332 RMATEL=SQRT(FLOAT((N+L+2+2*ND)*(N+L+6+2*ND)*(N-L+2+ND)))/8. GO TO 9 340 ND=ND/3 IF (LD*ND) 341,9,342 341 RMATEL=SQRT(FLOAT((N-L+4*ND)*(N-L+4+4*ND)*(N-L+2+ND)))/8. GO TO 9 342 RMATEL= SQRT(FLOAT((N+L+2+4*ND)*(N+L+6+4*ND)*(N+L+4+ND)))/8. 9 RETURN END C C####################################################################### C SUBROUTINE SORT(E,KV,NI) C C SORTING OF E(1:NI), KV(1:NI) ACCORDING TO INCREASING VALUES IN E C DIMENSION KV(1),E(1) INTEGER D,I,J REAL Y character KV*4,Y1*4 N=NI D=2**INT(1.4427*ALOG(FLOAT(N-1)))-1 341 IF (D .LE. 0) RETURN I=1 342 J=I IPD=I+D Y=E(IPD) Y1=KV(IPD) 343 IF (Y.LT. E(J)) GO TO 344 345 JPD=J+D E(JPD)=Y KV(JPD)=Y1 I=I+1 IF (I+D .LE. N) GO TO 342 D=(D-1)/2 GO TO 341 344 JPD=J+D E(JPD)=E(J) KV(JPD)=KV(J) J=J-D IF (J .GE. 1) GO TO 343 GO TO 345 END C C####################################################################### C SUBROUTINE SORT4(E,KV1,KV2,KV3,KV4,KV5,NI) C C SORTING OF E, KV1,...., KV5 ACCORDING TO INCREASING VALUES IN E C NI: DIMENSION OF E, KV1,..., KV5 C DIMENSION E(1),KV1(1),KV2(1),KV3(1),KV4(1),KV5(1) REAL KV2,KV3,KV4,KV5 CHARACTER Y1*4,KV1*4 CT WRITE(1,*)'UTSKRIFT I SORT4,NI',NI N=NI M=2**INT(1.4427*ALOG(FLOAT(N-1)))-1 1 IF(M.LE.0) RETURN I=1 2 J=I K=I+M Y=E(K) Y1=KV1(K) Y2=KV2(K) Y4=KV4(K) Y5=KV5(K) Y3=KV3(K) 3 IF(Y.LT.E(J)) GO TO 5 4 L=J+M E(L)=Y KV1(L)=Y1 KV2(L)=Y2 KV4(L)=Y4 KV5(L)=Y5 KV3(L)=Y3 I=I+1 IF(I+M.LE.N) GO TO 2 M=(M-1)/2 GO TO 1 5 L=J+M E(L)=E(J) KV1(L)=KV1(J) KV2(L)=KV2(J) KV4(L)=KV4(J) KV5(L)=KV5(J) KV3(L)=KV3(J) J=J-M IF(J.GE.1) GO TO 3 GO TO 4 END C C####################################################################### C FUNCTION BVCON(FVCON,XA,XB,YA,YB,NX,NY,X,WX,Y,WY,IX,IY) C C DOUBLE INTEGRATION SUBROUTINE FOR VOLUME CONSERVATION CONDITION. C FVCON IS THE EXTERNAL FUNCTION TO BE INTEGRATED. C X,WX ARE THE ARRAYS OF NORMALIZED POINTS AND WEIGHTS FOR X-VARI C Y,WY ARE THE ARRAYS OF NORMALIZED POINTS AND WEIGHTS FOR Y-VARI C DIMENSION X(1),WX(1),Y(1),WY(1) C CALCULATION OF THE DENORMALIZING X AND Y SLOPES HX=(XB-XA)/(2.*IX) HY=(YB-YA)/(2.*IY) BVCON=0.0 DO 10 L=1,IY C CALCULATION OF THE DENORMALIZING Y INTERCEPT BY=YA+(2.*L-1.)*HY DO 10 K=1,IX C CALCULATION OF THE DENORMALIZING X INTERCEPT BX=XA+(2.*K-1.)*HX DO 10 J=1,NY C CALCULATION OF THE DENORMALIZED POINT Y A=HY*Y(J)+BY DO 10 I=1,NX C CALCULATION OF THE DENORMALIZED POINT X AND SUMMATION OF THE C APPROXIMATING SUM 10 BVCON=BVCON+WY(J)*WX(I)*FVCON(HX*X(I)+BX,A) C ADJUSTMENT OF THE CALCULATED SUM. BVCON=HX*HY*BVCON RETURN END C C####################################################################### C FUNCTION CLEBI(I1,I2,I3,N1,N2,N3) C CALCULATION OF CLEBSCH-GORDAN COEFF C COMMON /PTR/ FCT(57) INTEGER Z,ZMIN,ZMAX J1=I1 J2=I2 J =I3 N=57 M1=N1 M2=N2 M=-N3 CC=0.0 JSUM=J1+J2+J JM1 =J1-IABS(M1) JM2 =J2-IABS(M2) JM3 =J -IABS(M ) IF((MOD(JSUM,2).NE.0).OR.(MOD(JM1,2).NE.0).OR.(MOD(JM2,2).NE.0) 1.OR.(MOD(JM3,2).NE.0).OR.(JM1.LT.0).OR.(JM2.LT.0).OR.(JM3.LT.0)) 2GO TO 1 IF((M1+M2+M.NE.0).OR.(J.GT.J1+J2).OR.(J.LT.IABS(J1-J2))) GO TO 1 ZMIN=0 IF(J-J2+M1.LT.0) ZMIN=-J+J2-M1 IF(J-J1-M2+ZMIN.LT.0) ZMIN=-J+J1+M2 ZMAX=J1+J2-J IF(J2+M2-ZMAX.LT.0) ZMAX=J2+M2 IF(J1-M1-ZMAX.LT.0) ZMAX=J1-M1 JA=(J1+M1)/2+1 JB=JA-M1 JC=(J2+M2)/2+1 JD=JC-M2 JE=(J +M )/2+1 JF=JE-M JG=(J1+J2-J)/2+1 JH=JA+JB-JG JI=JC+JD-JG JJ=JE+JF+JG-1 IF(JJ.GT.N) GO TO 5 IA=ZMIN/2 IB=JG-IA+1 IC=JB-IA+1 ID=JC-IA+1 IE=JA-JG+IA IF=JD-JG+IA FASE=1.0 IF(MOD(IA,2).EQ.0) FASE=-FASE Z =ZMIN ARII=FLOAT(J+1) ARI=(FCT(JA)+FCT(JB)+FCT(JC)+FCT(JD)+FCT(JE)+FCT(JF)+FCT(JG) 1 +FCT(JH)+FCT(JI)-FCT(JJ)+ALOG(ARII))/2.0 2 IA=IA+1 IB=IB-1 IC=IC-1 ID=ID-1 IE=IE+1 IF=IF+1 FASE=-FASE ARP=-FCT(IA)-FCT(IB)-FCT(IC)-FCT(ID)-FCT(IE)-FCT(IF) ARS=ARI+ARP IF(ARS.GT.80.0) GO TO 5 IF(ARS.LT.-80.0) GO TO 10 CC=CC+FASE*EXP(ARS) 10 CONTINUE Z=Z+2 IF(Z.LE.ZMAX) GO TO 2 1 CLEBI=CC RETURN C 3 WRITE(6,4) C 4 FORMAT(26H0 ERROR IN ARGUMENT OF VCC) C GO TO 7 5 WRITE(61,6) STOP 6 FORMAT(29H0 ERROR - FACTORIALS EXCEEDED) C 7 WRITE(61,8)I1,I2,I3,N1,N2,N3 C 8 FORMAT(10I10) C RETURN END C C####################################################################### C FUNCTION FVCON(XX,YY) C C FUNCTION TO BE INTEGRATED IN VOLUME CALC (STRETCHED COORD) C COMMON /CMN2/ EPS1,EPS2,EPS3,EPS4,EPS5,EPS6,EPS22 COMMON/CMN7/ EPSP,GAMMA X=XX Y=YY PI=3.1415926536 GAMMAR=GAMMA*PI/180.0 C ANM=1.0-EPS2/3.0*(3.*X*X-1.)+0.5*EPS22*(1.-X*X)*COS(2.*Y)+ C $EPS4/4.0*(COS(1.5*GAMMAR)**2+0.375*SIN(1.5*GAMMAR)**2)* C $(35.*X*X*X*X-30.*X*X+3.)-5.*EPS4/8.*SIN(1.5*GAMMAR)**2*(7.*X*X-1.) C $*(1.-X*X)*COS(2.*Y)+35.*EPS4/32.*SIN(1.5*GAMMAR)**2*(1.-X*X)* C $(1.-X*X)*COS(4.*Y) ANM=1.0-EPS2/3.0*(3.*X*X-1.)+0.5*EPS22*(1.-X*X)*COS(2.*Y)+ $EPS4/4.0*((5.*COS(GAMMAR)**2+1.)/6.)* $(35.*X*X*X*X-30.*X*X+3.) $-5.*EPS4/12.*SQRT(3.)*SIN(2.*GAMMAR)*(7.*X*X-1.) $*(1.-X*X)*COS(2.*Y)+35.*EPS4/24.*SIN(GAMMAR)**2*(1.-X*X)* $(1.-X*X)*COS(4.*Y) FVCON=1./(4.*PI)/SQRT((1.+EPS2/3.+EPS22/2.)*(1.+EPS2/3.-EPS22/2.)* $(1.-2.*EPS2/3.))/(ANM**(1.5)) RETURN END C C####################################################################### C SUBROUTINE HELP COMMON /LU/ LU COMMON /PTR/ F(57) COMMON /GM/ GM(150) COMMON /CMN1/ KAPPA,OMOM,MY,L2F(18) LU = 6 C*****COMPUTATION OF F(N+1)=LN(N^) F(1) = 0. DO 101 I=1,56 AII = FLOAT(I) F(I+1) = F(I) + ALOG(AII) 101 CONTINUE C*****COMPUTATION OF (L2)AV = L2F(N+1) L2F(1) = 0 DO 106 I = 1,17 L2F(I+1) = L2F(I) + I + 1 106 CONTINUE C*****COMPUTATION OF GAMMA-FUNCTION IN INTEGER AND HALF-INTEGER POINTS. GM(52)=1.0 GM(51)=1.77245385 DO 34 I=1,28 C KOLLA OM 28 R EN TILLR CKLIGT H\G GR NS^ GM(2*I+52)=I*GM(2*I+50) 34 GM(2*I+51)=(I-0.5)*GM(2*I+49) GM(49)=-2.0*GM(51) DO 35 I=1,20 IR=-2*I+49 IIR=IR+2 35 GM(IR)=GM(IIR)/(-I-0.5) RETURN END C C####################################################################### C SUBROUTINE LGAUSS(X,W,N) C C SUBROUTINE LGAUSS COMPUTES THE ZEROES OF THE LEGENDRE POLYNOMIAL C AND THEIR ASSOCIATED WEIGHTS FOR A GAUSSIAN QUADRATURE C DIMENSION X(1),W(1) IF(N-1) 1,2,3 C REQUEST FOR A ZERO POUNT FORMULA IS MEANINGLESS SO RETURN 1 write(6,11) 11 FORMAT(//,'REQUEST FOR ZERO POINT FORMULA IN LGAUSS') STOP 2 X(1)=0.0 W(1)=2.0 RETURN C FOR A ONE POINT FORMULA SEND BACK RESULTS WITHOUT COMPUTING 3 R=N G=-1. C THE INITIAL GUESS FOR THE SMALLEST ROOT OF P(N) IS TAKEN AS -1. DO 147 I=1,N TEST=-2. IC=N+1-I C WHENEVER WE FIND A ROOT OF THE POLYNOMIAL, ITS NEGATIVE IS ALSO A C THE INDEX IC TELLS WHERE TO STORE THE OTHER ROOT IF(IC.LT.I) GO TO 150 4 S=G T=1. U=1. V=0. C EVALUATION OF THE N-TH LEGENDRE POLYNOMIAL AND ITS FIRST DERIVATIV C WHERE U=DS/DX C V=DT/DX C DP=DP/DX DO 50 K=2,N A=K P=((2.0*A-1.0)*S*G-(A-1.0)*T)/A DP=((2.0*A-1.0)*(S+G*U)-(A-1.0)*V)/A V=U U=DP T=S 50 S=P IF (ABS(TEST+G) .LT. 0.00000005) GO TO 100 IF (ABS((TEST-G)/(TEST+G)).LT.0.0000005) GO TO 100 SUM=0. IF(I.EQ.1) GO TO 52 C THE FOLLOWING COMPUTES THE REDUCED LEGENDRE POLYNOMIAL AND ITS DER DO 51 K=2,I 51 SUM=SUM+1./(G-X(K-1)) 52 TEST=G G=G-P/(DP-P*SUM) GO TO 4 100 X(IC)=-G X(I)=G W(I)=2./(R*T*DP) W(IC)=W(I) 147 G=G-R*T/((R+2.)*G*DP+R*V-2.*R*T*SUM) 150 RETURN END C C####################################################################### C SUBROUTINE MRHOT(NNP,LLP,LAMP,NN,LL,LAM,T,ME) C C RADIAL MATRIX ELEMENT, ME, OF R**T (R DIMENSIONLESS) C N,NP - NUMBER OF OSCILLATOR QUANTA C L,LP - ORBITAL ANGULAR MOMENTUM C LAM, LAMP - PROJECTION OF L AND LP (DUMMY IN THIS CASE) C COMMON /LU/ LU REAL ME 1 FORMAT(12I6) INTEGER T,TA,E,C,D,SLUT COMMON /GM/ GM(150) 300 FORMAT(7H DANGER) 301 FORMAT(31H TERMINATION SOMETHING IS WRONG) 302 FORMAT(33H IS NOT THIS A WONDERFUL COMPILER) NANP=NNP LALP=LLP LAAMP=LAMP NAN=NN LAL=LL LAAM=LAM TA=T IF(LAAM-LAAMP) 231,201,231 201 CONTINUE NP=(NANP-LALP)/2+1 N=(NAN-LAL)/2+1 IF(NP-N)202,203,203 202 NPNP=NP NP=N N=NPNP L=LALP LALP=LAL LAL=L 203 CONTINUE C=(LAL-LALP-TA)/2 D=(LALP-LAL-TA)/2 E=D+NP-N-1 IF(2*C.EQ.(LAL-LALP-TA)) GO TO 401 ME=0.0 DO 3 K=1,N K1=2*(N-K+1)+50 K2=LAL-LALP-TA+2*(K-1)+50 K3=LALP-LAL-TA+2*(K+NP-N-1)+50 K4=2*(NP-N+K)+50 K5=2*(N-K+1)+LAL+LALP+1+TA+50 K6=LAL-LALP-TA+50 K7=LALP-LAL+50-TA 3 ME=ME+GM(2*N+50)*GM(2*NP+50)/GM(K1)*GM(K2)/GM(2*K+50)*GM(K3) */GM(K4)*GM(K5)/GM(K6)/GM(K7) K1=2*N+2*LAL+51 K2=2*NP+2*LALP+51 ME=ME/SQRT(GM(2*N+50)*GM(K1)*GM(2*NP+50)*GM(K2))*(-1)**(NP-N) GO TO 230 401 CONTINUE IF(C.LE.0.AND.D.LE.0) GO TO 204 GO TO 208 204 CONTINUE IF(E.GT.-1) GO TO 206 GO TO 209 206 SLUT=0 GO TO 220 209 CONTINUE IF(C.LE.E+1) GO TO 207 GO TO 221 207 SLUT=-E GO TO 220 221 SLUT=1-C GO TO 210 208 CONTINUE IF(C.LE.0.AND.D.GT.0) GO TO 211 GO TO 212 211 SLUT=1-C GO TO 213 212 CONTINUE IF(C.GT.0.AND.D.LE.0) GO TO 214 GO TO 215 214 CONTINUE IF(E.GT.-1) GO TO 216 GO TO 217 216 SLUT=0 GO TO 218 217 SLUT=-E 218 CONTINUE GO TO 219 215 CONTINUE WRITE(LU,300) WRITE(LU,301) WRITE(LU,302) STOP 219 CONTINUE 210 CONTINUE 220 CONTINUE 213 CONTINUE IF(SLUT.EQ.0) GO TO 231 IF(SLUT.LT.0) GO TO 222 GO TO 223 222 WRITE(LU,300) WRITE(LU,301) WRITE(LU,302) STOP 223 ME=0.0 IF(SLUT.GT.N) SLUT=N DO 234 K=1,SLUT KK1=2*(N-K+1)+50 KK2=2*(N-K+1)+LAL+LALP+1+TA+50 KK3=2*(NP-N+K)+50 KK4=2-2*C+50 KK5=2*(2-K-C)+50 KK6=2*(C+K-1)+50 KK7=2*C+50 IF(C.LE.0) GO TO 224 GO TO 225 224 GLI=(-1)**(K+1)*GM(KK4)/GM(KK5) GO TO 226 225 GLI=GM(KK6)/GM(KK7) 226 CONTINUE KK8=2-2*D+50 KK9=2*(2-D-K-NP+N)+50 KK10=2*(D+K+NP-N-1)+50 KK11=2*D+50 IF(D.LE.0.) GO TO 227 GO TO 228 227 GLA=GM(KK8)/GM(KK9)*(-1)**(1+K+NP-N) GO TO 229 228 GLA=GM(KK10)/GM(KK11) 229 CONTINUE ME=ME+GM(2*N+50)/GM(2*K+50)*GM(2*NP+50)/GM(KK1)*GM(KK2)/GM(KK3)* *GLI*GLA 234 CONTINUE KK12=2*(N+LAL)+51 KK13=2*(NP+LALP)+51 ME=ME*(-1)**(N-NP)/SQRT(GM(2*N+50)*GM(2*NP+50)*GM(KK12)*GM(KK13)) GO TO 230 231 ME=0.0 230 CONTINUE RETURN END C C####################################################################### C SUBROUTINE OMO(OMOM) C C VOLUME CONSERVATION CALCULATION C COMMON /CMN2/ EPS1,EPS2,EPS3,EPS4,EPS5,EPS6,EPS22 DIMENSION Z(9),WZ(9) COMMON/BAS/ ISTRCH,ICORR,IREC EXTERNAL FVCON C C VOLUME CONSERVATION FACTOR CALCULATED WITH DELTA4, DELTA6 C = 0 IN THE CASE OF SPHERICAL COORDINATES C IF (ISTRCH .EQ. 0) GO TO 6 IF(ABS(EPS4).LT.0.0001) GO TO 4 PI=3.1415926536 CALL LGAUSS(Z,WZ,9) NX=1 NY=1 A=BVCON(FVCON,-1.0,0.0,0.0,PI/2.,9,9,Z,WZ,Z,WZ,NX,NY) OMOM=(8.*A)**(1.0/3.0) RETURN 4 OMOM=((1.+EPS2/3.+EPS22/2.)*(1.+EPS2/3.-EPS22/2.)*(1.-2.*EPS2/3.)) $**(-1.0/3.0) RETURN 6 CONTINUE FAC = 2.*EPS2/3. OMOM=((1.+FAC+EPS22)*(1.+FAC-EPS22)*(1.-2.*FAC))**(-1./6.) RETURN END C C####################################################################### C SUBROUTINE BIGMAX(A,VALU,NN,NEIG,NVEC,VECTR,MDIM) C C MATRIX DIAGONALISATION USING HOUSEHOULDER'S METHOD C C A: MATRIX TO BE DIAGONALIZED - TRIANGULAR PART STORED IN C ONE DIMENSION; (NN+1)*NN/2 C VALU: EIGENVALUES ORDERED WITH LARGEST ONE FIRST C NN: DIMENSION OF MATRIX; "NN * NN" C NEIG: NUMBER OF EIGENVALUES CALCULATED C NVEC: NUMBER OF EIGENVECTORS CALCULATED C VECTR: EIGENVECTORS C MDIM: FIRST DIMENSION AS DECLARED FOR VECTR C ! DIMENSION A(1),VALU(1),VALL(225),UPPERD(225),DIAG(225),V(225), ! $SQGAMP(225),INTER(1),TMULT(1),T(225,3),SQGAMR(225),UPPER2(225) ! DIMENSION VECTR(MDIM,1),IFANA(225) ! EQUIVALENCE (TMULT,INTER,SQGAMP,VALL),(T,SQGAMR),(T(1,2),UPPER2) ! EQUIVALENCE (TMULTP,IMULTP),(COMPL1,ICMPL1) DIMENSION A(0:25425),VALU(1), 1 VALL(225),UPPERD(225),DIAG(225),V(0:2250), 1 SQGAMP(225),INTER(225),TMULT(225),T(225,3) 1 ,SQGAMR(0:255),UPPER2(0:450) DIMENSION VECTR(MDIM,1),IFANA(225) EQUIVALENCE (TMULT,INTER,SQGAMP,VALL) equivalence (T,SQGAMR(1)),(T(1,2),UPPER2(1)) EQUIVALENCE (TMULTP,IMULTP),(COMPL1,ICMPL1) v(0)=0. sqgamr(0)=0. upper2(0)=0. A(0)=0. ! NZ=0 N=NN IF(N.LE.1) GO TO 41 NP1=N+1 NM1=N-1 NM2=N-2 NT2P1=N*2+1 IX=0 DO 10 I=1,NM2 SIGMA2=0. IP1=I+1 DO 1 J=IP1,N IJ=IX+J 1 SIGMA2=SIGMA2+A(IJ)**2 SIGMA=SQRT(SIGMA2) II=IX+I DIAG(I)=A(II) IIP1=IX+I+1 UPPERD(I)=-SIGN(SIGMA,A(IIP1)) UPPER2(I)=SIGMA2 IF(ABS(SIGMA).GT.ABS(A(IIP1)))GO TO 2 UPPERD(I)=A(IIP1) A(IIP1)=0. GO TO 10 2 A(IIP1)=SQRT(1.+ABS(A(IIP1))/SIGMA) SQTGAM=-SIGN(SIGMA*A(IIP1),UPPERD(I)) IP2=I+2 DO 3 J=IP2,N IJ=IX+J 3 A(IJ)=A(IJ)/SQTGAM JK1=I*(2*N-I-1)/2 JX=JK1 IIX=JK1 DO 5 J=IP1,N SQGAMP(J)=0. JK=JK1+J DO 4 K=IP1,J IK=IX+K SQGAMP(J)=SQGAMP(J)+A(JK)*A(IK) 4 JK=JK+N-K IF(J.EQ.N)GO TO 6 CALL LOOP1(J+2,NP1,SQGAMP(J),A(JX),A(IX)) 5 JX=JX+N-J 6 DELGAM=0. DO 7 J=IP1,N IJ=IX+J 7 DELGAM=DELGAM+A(IJ)*SQGAMP(J) DGO2=.5*DELGAM DO 8 J=IP1,N IJ=IX+J 8 SQGAMR(J)=SQGAMP(J)-DGO2*A(IJ) DO 9 II=IP1,N III=IX+II CALL LOOP2(A(IIX),A(IX),SQGAMR(NZ),SQGAMR(II),A(III),II+1,NP1) 9 IIX=IIX+N-II 10 IX=IX+N-I M=N*(N+1)/2 UPPERD(NM1)=A(M-1) UPPER2(NM1)=UPPERD(NM1)**2 DIAG(NM1)=A(M-2) DIAG(N)=A(M) ENORM=AMAX1(ABS(DIAG(1))+ABS(UPPERD(1)),ABS(DIAG(N))+ABS(UPPERD(NM 11))) DO 11 I=2,NM1 ENRTMP=ABS(DIAG(I))+ABS(UPPERD(I))+ABS(UPPERD(I-1)) 11 IF(ENRTMP.GT.ENORM)ENORM=ENRTMP DO 12 I=1,NEIG VALU(I)=ENORM 12 VALL(I)=-ENORM DO 24 I=1,NEIG 13 ROOT=.5*(VALU(I)+VALL(I)) IF(ROOT.EQ.VALL(I).OR.ROOT.EQ.VALU(I))GO TO 24 NAGREE=0 PM2=0. PM1=1. DO 21 J=1,N IF(PM2.NE.0.)GO TO 15 14 PM1=SIGN(1.,PM1) GO TO 17 15 IF(PM1.NE.0.)GO TO 17 16 P=-SIGN(1.,PM2) PM2=0. IF(UPPER2(J-1))18,14,18 17 P=DIAG(J)-ROOT-UPPER2(J-1)*PM2/PM1 PM2=1. 18 IF(P)21,19,20 19 PM2=PM1 IF(PM2)21,20,20 20 NAGREE=NAGREE+1 21 PM1=P DO 23 J=I,NEIG IF(J.LE.NAGREE)GO TO 22 IF(VALU(J).LE.ROOT)GO TO 13 VALU(J)=ROOT GO TO 23 22 VALL(J)=ROOT 23 CONTINUE GO TO 13 24 CONTINUE IF(NVEC.EQ.0)GO TO 39 EPSLON=ENORM*1.E-6 ICMPL1=NOT(1) DO 38 I=1,NVEC DO 25 J=1,N V(J)=1. T(J,2)=DIAG(J)-VALU(I) IF(J.EQ.N)GO TO 26 T(J,3)=UPPERD(J) 25 T(J+1,1)=UPPERD(J) 26 T(N,3)=0. DO 29 J=1,N IF(ABS(T(J,2)).LT.1.E-15)T(J,2)=EPSLON T(J,1)=T(J,2) T(J,2)=T(J,3) T(J,3)=0. IF(J.EQ.N)GO TO 30 INTER(J)=0 IFANA(J)=0 JP1=J+1 IF(ABS(T(JP1,1)).LE.ABS(T(J,1)))GO TO 28 CO INTER(J)=1 IFANA(J)=1 DO 27 K=1,3 TEMP=T(J,K) T(J,K)=T(JP1,K) 27 T(JP1,K)=TEMP 28 TMULTP=T(JP1,1)/T(J,1) CO TMULT(J)=IOR(INTER(J),IAND(IMULTP,ICMPL1)) TMULT(J)=TMULTP T(JP1,2)=T(JP1,2)-TMULTP*T(J,2) 29 T(JP1,3)=T(JP1,3)-TMULTP*T(J,3) 30 ITER=1 31 DO 32 J=1,N L=N+1-J V(L)=(V(L)-T(L,2)*V(L+1)-T(L,3)*V(L+2))/T(L,1) 32 CONTINUE VNORM=0. DO 33 L=1,N 33 VNORM=VNORM+V(L)**2 VNORM=SQRT(VNORM) DO 34 J=1,N 34 V(J)=V(J)/VNORM IF(ITER.EQ.2)GO TO 36 ITER=2 DO 35 L=2,N LM1=L-1 CO IF(IAND(INTER(LM1),1).EQ.0.)GO TO 35 IF(IFANA(LM1).EQ.0)GO TO 35 VTEMP=V(LM1) V(LM1)=V(L) V(L)=VTEMP 35 V(L)=V(L)-TMULT(LM1)*V(LM1) GO TO 31 36 IF(VNORM.EQ.0.)V(I)=1. IIX=(N*N-N-6)/2 DO 37 KK=1,NM2 IIP1=N-KK UTV=0. CALL LOOP3(UTV,A(IIX),V(NZ),IIP1+1,NP1) CALL LOOP4(A(IIX),V(NZ),NP1,IIP1+1,UTV) 37 IIX=IIX+IIP1-N-2 DO 100 J=1,N 100 VECTR(J,I)=V(J) 38 CONTINUE GO TO 39 41 VECTR(1,1) = 1.0 VALU(1) = A(1) 39 RETURN END C C####################################################################### C SUBROUTINES FOR BIGMAX: LOOP1, LOOP2, LOOP3 C SUBROUTINE LOOP1(JP2,NP1,SGAMPJ,AJX,AIX) DIMENSION AJX(1), AIX(1) DO 1 L=JP2,NP1 1 SGAMPJ=SGAMPJ+AJX(L)*AIX(L) RETURN END SUBROUTINE LOOP2(AIIX,AIX,S,SI,AIII,IP1,NP1) DIMENSION AIIX(1),AIX(1),S(1) DO 2 JJ=IP1,NP1 2 AIIX(JJ)=AIIX(JJ)-AIII*S(JJ)-SI*AIX(JJ) RETURN END SUBROUTINE LOOP3(UTV,AIIX,V,IIP2,NP1) DIMENSION AIIX(1), V(1) DO 3 J=IIP2,NP1 3 UTV=UTV+AIIX(J)*V(J) RETURN END SUBROUTINE LOOP4(AIIX,V,NP1,IIP2,UTV) DIMENSION AIIX(1),V(1) DO 4 K=IIP2,NP1 4 V(K)=V(K)-AIIX(K)*UTV RETURN END C C####################################################################### C APPROPRIATE FOR NORD COMPUTER C CC FUNCTION MSLEFT(XX) C DOUBLE INTEGER IUSE,TUSED C IUSE=TUSED(GOLEG) C MSLEFT=-IUSE*20 C RETURN C END C C*********************************************************************** C APPROPRIATE FOR VAX (FROM TORD) C ! FUNCTION MSLEFT(XX) ! ! LOGICAL STATUS ! DATA I0,IPID,JPI_CPU/0,0,1031/ ! IF(IPID.EQ.0)THEN ! STATUS = LIB$GETJPI(I0,IPID,,ISEX,,) ! IF(.NOT.STATUS)THEN ! WRITE(*,*)' something wrong with lib$getjpi!' ! ENDIF ! ENDIF ! STATUS = LIB$GETJPI(JPI_CPU,IPID,,ISEX,,) ! IF(.NOT.STATUS)THEN ! WRITE(*,*)' something wrong with lib$getjpi!' ! ENDIF ! MSLEFT=-ISEX*10 ! RETURN ! END C C*********************************************************************** C APPROPRIATE FOR VAX (FROM TORD) C c FUNCTION SECOND(XX) c LOGICAL STATUS c DATA I0,IPID,JPI_CPU/0,0,1031/ c IF(IPID.EQ.0)THEN c STATUS = LIB$GETJPI(I0,IPID,,ISEX,,) c IF(.NOT.STATUS)THEN c WRITE(*,*)' something wrong with lib$getjpi!' c ENDIF c ENDIF c STATUS = LIB$GETJPI(JPI_CPU,IPID,,ISEX,,) c IF(.NOT.STATUS)THEN c WRITE(*,*)' something wrong with lib$getjpi!' c ENDIF c SECOND=FLOAT(ISEX)/100. c RETURN c END cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc