C INCLUDE 'CALPGM.INC' C C Copyright (C) 1989, California Institute of Technology C All rights reserved. U. S. Government Sponsorship under C NASA Contract NAS7-918 is acknowledged. C C Herbert M. Pickett, 8 Mar 1989 C HMP, GETQN modified for bug in case of oblate sym top C HMP, GETQN modified for b type parity in symmetric top C HMP, INTENS modified to add ITYP = 6 C HMP, SYMKSQ modified for -KSQ C HMP, DPMAKE modified for no base C HMP, add SPECOP C HMP, new KROLL code C HMP, fix code for oblate off-diagonal elements C HMP, new code for I(TOT) representation C HMP, new flag to signal no diagonalization or proj. sort (IDIAG) C HMP, new code to use PCARD in SETOPT, new kmin in SETBLK C HMP, new option for fourier coefficient C HMP, new option for IDIAG=2 using FIXBLK C HMP, divided definition of dipole ITYP=5 by 2 C HMP, fix SETBLK for negative weights C HMP, new code when number of quanta .GT. 6 C HMP, new call for SETOPT C HMP, spins can now change with 'vibrational' state C HMP, ESYMWT now recognizes 'l' selection rules C HMP, modify call to HAMX for hamiltonian dump C HMP, negative ESYMWT (default) now diables ESYMWT feature C HMP, added type 9 and 10 dipoles to INTENS C HMP, correction to bug for I(TOT) weights in SETOPT C HMP, new options for IAX in SETOPT C HMP, constained derivatives calculated within HAMX C HMP, IDPAR is now REAL*8 C HMP, NVIB now can be <99 C HMP, N now can be <360 C HMP, fixed error in dipole matrix element phase for vib. dipoles C HMP, 'l' doubling changed to K-l basis, REMOVED IDIAG=2 option C HMP, E states have positive K C HMP, SETBLK now a function to return size of vibrational field C HMP, SETBLK passed PAR and rho now precomputed C HMP, fix phase error for sine operators in DIRCOS C HMP, fix error in SETINT, INTENS interaction C HMP, add option for SzNz interaction term C HMP, check for reasonable ESYMWT C HMP, new Fourier parameter id format C HMP, add new IDIAG=2 and IDIAG=3 option C HMP, fixed bug in SETINT for complex dipoles C HMP, remove code to recognize old Fourier parameter id format C SUBROUTINE HAMX (IBLK,NSIZD,IDPAR,PAR,EGY,T,DEDP,IFDUMP) INTEGER IBLK,NSIZD,IFDUMP REAL*8 IDPAR(0:*),PAR(*),EGY(*),T(NSIZD,*),DEDP(NSIZD,*) C C .. PACKAGE FOR 9 INTERACTING VIBRATIONAL STATES WITH MULTI-SPIN C C calculate energies and derivatives for block IBLK C C on entry: C IBLK= block number C NSIZE= block size and dimension C IDPAR= list of parameter identifiers ( element 0 is length ) C PAR = list of parameter values C NPARD= number of parameters with positive IDPAR C on return: C EGY = energy list C T = eigenvector matrix C DEDP = derivative of energy with respect to parameter C INTEGER NDPAR,NDMX,MAXS,MXSPIN,NDSPIN,I0,I1 INTEGER*4 IDTMP,L10 PARAMETER (NDPAR=750,NDMX=450,MAXS=64,MXSPIN=5,NDSPIN=MXSPIN*9) PARAMETER (I0=0,I1=1,L10=10) INTEGER NSIZE,IXCOM(6),JXCOM(6),ISCOM(0:MXSPIN+1),IS1,NTENS, & JSCOM(0:MXSPIN+1),LLCOM(MXSPIN+1),LSCOM(MXSPIN+1), & OLDPAR,GETLL,GETQQ,NOLD,NSBLK,IFF,IIWT,N,NPAR,KL, & IXX,JXX,IBASE,JBASE,KBGNI,KBGNJ,IV,JV,JSPN,NI,NJ, & I,IJD,ND,IPARX,ITY,NCOS,IZ,JZ,IPAR,NSQJ,ISPN, & IKQ,L,LD,KD,INS,SI1,SI2,SIZ,IERR,BLKSYM,NBLK,KDIAG INTEGER*2 SBKPTR(MAXS+1),KMIN(MAXS),IVS(MAXS),ITAU(MAXS) REAL*8 SQNN,SQJ(10),RMATRX,ELE,DN,ZPAR,PBAR,ZERO(2), & DNI,ZFAC(NDPAR),SPFAC(0:NDSPIN),WK(NDMX),DBAR(NDMX) INTEGER*4 IPTYP(NDPAR) INTEGER*2 NPARX,NJQ(NDPAR),IP(NDPAR),IPSYM(NDPAR), & IDX(NDMX),JDX(NDMX),IPDER(0:NDPAR),IVSYM INTEGER KAPPA,NSPIN,NTSPIN,IISPN(NDSPIN),ITPTR,ITPTRM,IDIAG,NVIB LOGICAL FIRST,PARSKP,ROLL COMMON/HSAVE/SQNN,SQJ,ZFAC,IPTYP,IPSYM,IP,NPARX,NJQ,IPDER COMMON /SQN/ KAPPA,NSPIN,NTSPIN,IISPN,ITPTR,ITPTRM,IDIAG,NVIB COMMON /SPFACV/ SPFAC,NOLD SAVE ZERO SAVE /HSAVE/,/SQN/,/SPFACV/ DATA ZERO/0.,0./ NSIZE=NSIZD C get F and sub-block structure NSBLK=GETQQ(MAXS,IBLK,IFF,IIWT,SBKPTR,KMIN,IVS) C zero hamiltonian matrix CALL DCLR(NSIZE,NSIZE,T,I1) NPAR=IDPAR(0) C zero derivative matrix N=IPDER(0) CALL DCLR(NSIZE,N,DEDP,I1) C set up hamiltonian FIRST=.TRUE. KDIAG=0 C loop on wang blocks 240 OLDPAR=0 DO 720 IXX=1,NSBLK IBASE=SBKPTR(IXX) N=SBKPTR(IXX+1)-IBASE KBGNI=KMIN(IXX) ISPN=IVS(IXX) CALL GETQS(ISPN,IFF,N,KBGNI,IXCOM,ISCOM,IV) IS1=ISCOM(0)+1 NI=IXCOM(3) ITAU(IXX)=KBGNI+KBGNI-IAND(NI+IXCOM(2),1) IF (ISCOM(1).NE.NOLD) THEN C set up powers of N(N+1) NOLD=ISCOM(1) DN=NI SQNN=DN*(NI+1) DO 200 I=2,10 200 SQJ(I)=SQJ(I-1)*SQNN DNI=NI+NI+1 SPFAC(0)=SQRT(SQNN*DNI) ENDIF DO 730 JXX=1,IXX C set up sub-block quanta IJD=IXX-JXX JBASE=SBKPTR(JXX) N=SBKPTR(JXX+1)-JBASE KBGNJ=KMIN(JXX) JSPN=IVS(JXX) CALL GETQS(JSPN,IFF,N,KBGNJ,JXCOM,JSCOM,JV) NJ=JXCOM(3) ND=NI-NJ IF(ND.GT.3) GO TO 730 IF(ND.NE.0) ND=-20 IDTMP=IV+JV+99*MIN(IV,JV) IVSYM=BLKSYM(IXCOM,JXCOM)+4*IDTMP C loop over parameters NBLK=0 DO 740 IPARX=1,NPARX IF(IPSYM(IPARX).NE.IVSYM) GO TO 740 IDTMP=IPTYP(IPARX) IF(IDTMP.GE.0) THEN IF(OLDPAR.NE.0) THEN C add diagonal or sub-diagonal to Hamiltonian C using composite parameter value PBAR DO 300 I=1,NCOS IZ=IDX(I) JZ=JDX(I) 300 T(IZ,JZ)=T(IZ,JZ)+PBAR*WK(I) OLDPAR=0 ENDIF C first occurance of direction cosine operator NCOS=-1 ELSE IF(NCOS.EQ.0) GO TO 740 C last cosine operator was zero IDTMP=-IDTMP ENDIF IPAR=IP(IPARX) NSQJ=ABS(NJQ(IPARX)) IF(NJQ(IPARX).LT.ND) THEN C sum parameters which are different only by N*(N+1) IF(PARSKP) GO TO 740 C jump if mother operator was zero ITY=MOD(IDTMP,L10) ELSE C new type of operator PARSKP=.TRUE. CALL IDPARS(IDTMP,IKQ,ITY,L,LD,KD,INS, & SI1,SI2,SIZ,KL) NTENS=SI1 C .. test tensor compatability IF(SIZ.NE.0) THEN CALL SZNZ(SIZ,NI,NJ,IXCOM,JXCOM,ISCOM,JSCOM) IF(SIZ.EQ.0) GO TO 750 ENDIF IF(GETLL(LD,KD,L,SI1,SI2,LLCOM,LSCOM,ISCOM, & JSCOM).NE.0) GO TO 750 C .. get direction cosines and K dependent part of operator IF(NCOS.LT.0) THEN NCOS=NDMX IF(ITY.GT.160) NCOS=-NCOS CALL DIRCOS(IXCOM,JXCOM,LD,KD,NCOS,WK, & IDX,JDX,IJD,KL) IF(NCOS.EQ.0) GO TO 750 IF(ITY.NE.0) THEN C .. special operator CALL SPECOP(ITY,NBLK,NSQJ,IKQ,KBGNI, & KBGNJ,NCOS,WK,IDX,JDX,IXCOM,PAR(IPAR)) CALL SPECFC(ITY,IV,JV,KD,KBGNI, & KBGNJ,NCOS,WK,IDX,JDX) ENDIF IF(SIZ.NE.0) CALL SZNZOP(NI,NJ,KBGNI, & KBGNJ,ISCOM,JSCOM,NCOS,WK,IDX,JDX) IF(IKQ.NE.0) THEN C ..find K*K dependence (anti-commutated) CALL SYMKSQ(IKQ,KBGNI, & KBGNJ,NCOS,WK,IDX,JDX) ENDIF IF(NCOS.EQ.0) GO TO 750 DO 220 I=1,NCOS IDX(I)=IDX(I)+IBASE 220 JDX(I)=JDX(I)+JBASE IF(FIRST) THEN PBAR=0. ELSE C .. begin by setting up reduced expectation values CALL DPMAKE(NSIZD,DBAR,T,NCOS,WK,IDX,JDX) ENDIF ENDIF C .. find K independent parts ZPAR=RMATRX(LD,L,IXCOM,JXCOM) C .. spherical tensor transformations CALL TENSOR(ZPAR,NTENS,ISCOM,JSCOM,LLCOM,LSCOM) IF(SIZ.LT.0) THEN CALL SZNZ(SIZ,NI,NJ,IXCOM,JXCOM,ISCOM,JSCOM) ENDIF C .. find N*(N+1) and N.S dependence CALL SYMNSQ(NSQJ,INS,ISCOM,JSCOM,IISPN(IS1),ZPAR) PARSKP=ABS(ZPAR).LT.1.E-30 IF(PARSKP) GO TO 740 ENDIF ELE=ZPAR*SQJ(NSQJ)*ZFAC(IPARX) IF(FIRST) THEN C .. sum into composite parameter PBAR=PBAR+ELE*PAR(IPAR) OLDPAR=IPAR ELSE C .. sum into derivative matrix I=IPDER(IPAR) IF(I.LT.0) THEN I=-I ELE=ELE*PAR(IPAR)/PAR(I) I=IPDER(I) ENDIF CALL DAXPY(NSIZE,ELE,DBAR,I1,DEDP(1,I),I1) ENDIF 750 IF(SIZ.LT.0) & CALL SZNZ(SIZ,NI,NJ,IXCOM,JXCOM,ISCOM,JSCOM) 740 CONTINUE 730 CONTINUE 720 CONTINUE IF (FIRST) THEN IF(OLDPAR.NE.0) THEN C add diagonal or sub-diagonal to Hamiltonian C using composite parameter value PBAR DO 800 I=1,NCOS IZ=IDX(I) JZ=JDX(I) 800 T(IZ,JZ)=T(IZ,JZ)+PBAR*WK(I) ENDIF ROLL=.FALSE. IF(IDIAG.LT.0.OR.NSIZE.LE.1.OR.IFDUMP.NE.0) THEN ND=NSIZE+1 CALL DCOPY(NSIZE,T,ND,EGY,I1) IF(IFDUMP.NE.0) RETURN DO 840 I=1,NSIZE CALL DCOPY(NSIZE,ZERO,I0,T(1,I),I1) 840 T(I,I)=1 ELSE IF(IDIAG.EQ.0) THEN CALL KROLL(NSIZD,T,NSBLK,SBKPTR,KMIN,KAPPA,ROLL) ELSE IF(IDIAG.EQ.2) THEN ND=NSIZE+1 CALL DCOPY(NSIZE,T,ND,EGY,I1) CALL ORDHAM(NSIZD,NSIZE,T,EGY,NSBLK,SBKPTR,JDX,KAPPA) ENDIF C diagonalize CALL HDIAG (NSIZD,NSIZE,T,EGY,WK,IERR) IF (IERR.NE.0) THEN WRITE (*,890) IERR,IBLK,NSIZE 890 FORMAT (' DIAGONALIZATION FAILURE',3I5) STOP ENDIF IF(IDIAG.EQ.0) THEN C energy sort of subblock CALL ORDBLK(NSIZD,NSIZE,T,EGY,NSBLK,SBKPTR,WK,IDX,-KAPPA) ELSE IF(IDIAG.EQ.1) THEN C projection sort of full block DO 841 I=1,NSIZE+1 841 JDX(I)=I CALL ORDBLK(NSIZD,NSIZE,T,EGY,NSIZE,JDX,WK,IDX,-KAPPA) ELSE IF(IDIAG.EQ.2) THEN CALL ORDBLK(NSIZD,NSIZE,T,EGY,NSBLK,SBKPTR,WK,IDX,-KAPPA) I=0 CALL ORDHAM(NSIZD,NSIZE,T,EGY,NSBLK,SBKPTR,JDX,I) ELSE IF(IDIAG.EQ.3) THEN C projection sort by vibration JXX=1 JDX(1)=1 DO 842 I=2,NSBLK IF(IVS(I-1)/4.NE.IVS(I)/4) THEN JXX=JXX+1 JDX(JXX)=SBKPTR(I) ENDIF 842 CONTINUE JDX(JXX+1)=SBKPTR(NSBLK+1) CALL ORDBLK(NSIZD,NSIZE,T,EGY,JXX,JDX,WK,IDX,-KAPPA) DO 843 I=1,NSBLK IXX=SBKPTR(I) WK(IXX)=ITAU(I) JZ=SBKPTR(I+1)-2 DO 844 IZ=IXX,JZ 844 WK(IZ+1)=WK(IZ)+4 843 CONTINUE I=-1 CALL ORDHAM(NSIZD,NSIZE,T,WK,JXX,JDX,IDX,I) I=0 CALL ORDHAM(NSIZD,NSIZE,T,EGY,JXX,JDX,IDX,I) ENDIF ENDIF C .. go back and get expectation values FIRST=.FALSE. GO TO 240 ELSE IF(ROLL) THEN CALL DCOPY(NSIZE,ZERO,I0,EGY,I1) I=IPDER(0) DO 850 IPAR=NPAR,1,-1 IF(IPDER(IPAR).GE.0) THEN ELE=PAR(IPAR) CALL DAXPY(NSIZE,ELE,DEDP(1,I),I1,EGY,I1) I=I-1 ENDIF 850 CONTINUE ENDIF RETURN END SUBROUTINE DCLR(N1,N2,VEC,IX) C formal parameters: INTEGER N1,N2,IX REAL*8 VEC(0:*) C clear a N1*N2 block INTEGER*4 N1L,NSQ,IV INTEGER N,NMX,I0 REAL*8 ZERO(2) PARAMETER (NMX=32766,I0=0) SAVE ZERO DATA ZERO/0.,0./ IV=0 N1L=N1 NSQ=N1L*N2 10 IF(NSQ.GT.NMX) THEN CALL DCOPY(NMX,ZERO,I0,VEC(IV),IX) IV=IV+NMX NSQ=NSQ-NMX GO TO 10 ENDIF N=NSQ CALL DCOPY(N,ZERO,I0,VEC(IV),IX) RETURN END SUBROUTINE SPECOP(ITY,NBLK,NSQJ,IKQ,KSI,KSJ,NCOS, & WK,IX,JX,IXCOM,PAR) C formal parameters: INTEGER ITY,NBLK,NSQJ,IKQ,KSI,KSJ,NCOS,IXCOM(*) INTEGER*2 IX(*),JX(*) REAL*8 WK(*),PAR C subroutine to make Euler transform operators C on entry: C ITY= type of operator C NBLK= flag for new block (=0 initially) C NSQJ=power of N*N+N +1 C IKQ= power of K*K C KSI,KSJ= starting values of K C NCOS= length of vector WK C WK= vector of a sub-diagonal C IX,JX= relative indices of elements in WK C on return: C WK= modified vector INTEGER I,J,KI,KJ,IST,KINC PARAMETER (KINC=2) REAL*8 ADEN(0:4),BDEN(0:4),AKISQ,AKJSQ,XI,XJ,DI,DJ,SQN SAVE ADEN,BDEN IF(NBLK.EQ.0) THEN DO 10 I=0,4 ADEN(I)=0 10 BDEN(I)=1 NBLK=1 ENDIF J=IAND(ITY,15)-2 IF(J.LT.0) RETURN IST=J/2 IF(IST+IST.NE.J) THEN IF(NSQJ.EQ.1.AND.IKQ.EQ.1) THEN ADEN(IST)=ADEN(IST)+PAR IKQ=0 ELSE IF(NSQJ.EQ.2.AND.IKQ.EQ.0) THEN ADEN(IST)=ADEN(IST)-PAR SQN=IXCOM(3) BDEN(IST)=BDEN(IST)+PAR*(SQN*SQN+SQN) ELSE IF(NSQJ.EQ.1.AND.IKQ.EQ.0) THEN IF(J.EQ.1) ITY=0 ENDIF NCOS=0 ELSE IF(IST.LE.4) THEN SQN=IXCOM(3) SQN=SQN*SQN+SQN DO 20 I=1,NCOS KI=KSI+IX(I)*KINC AKISQ=KI*KI DI=ADEN(IST)*AKISQ+BDEN(IST) AKISQ=AKISQ/DI KJ=KSJ+JX(I)*KINC AKJSQ=KJ*KJ DJ=ADEN(IST)*AKJSQ+BDEN(IST) AKJSQ=AKJSQ/DJ XI=0.5 XJ=0.5 C K*K part DO 22 J=1,IKQ XI=AKISQ*XI 22 XJ=AKJSQ*XJ AKISQ=SQN/DI-AKISQ AKJSQ=SQN/DJ-AKJSQ C N*N+N-K*K part DO 23 J=2,NSQJ XI=AKISQ*XI 23 XJ=AKJSQ*XJ C extra part IF(IST.GT.0) THEN XI=XI/DI XJ=XJ/DJ ENDIF WK(I)=WK(I)*(XI+XJ) 20 CONTINUE NSQJ=1 IKQ=0 ENDIF RETURN END SUBROUTINE SPECFC(ITY,IV,JV,KDEL,KSI,KSJ,NCOS,WK,IX,JX) C formal parameters: INTEGER ITY,IV,JV,KDEL,KSI,KSJ,NCOS INTEGER*2 IX(*),JX(*) REAL*8 WK(*) C subroutine to make Fourier coefficient operators C on entry: C ITY= type of operator C IV= v', JV=v" C KDEL= change in K C KSI,KSJ= starting values of K C NCOS= length of vector WK C WK= vector of a sub-diagonal C IX,JX= relative indices of elements in WK C on return: C WK= modified vector INTEGER KINC,I,J,KI,KJ,N PARAMETER (KINC=2) REAL*8 PFAC(0:98),DI,DJ,ANG SAVE PFAC IF(ITY.LT.16) THEN IF(ITY.GT.0) RETURN IF(ITY.LT.0) THEN DO 10 I=0,98 10 PFAC(I)=0. RETURN ENDIF ENDIF IF(ITY.EQ.0) THEN C multiply by pi/3 IF(IV.EQ.JV) PFAC(IV)=WK(1)*1.047197551D0 RETURN ENDIF J=ITY/16 IF(J.EQ.0) RETURN N=J IF(N.GT.10) N=N-10 DI=N*PFAC(IV) DJ=N*PFAC(JV) DO 40 I=1,NCOS KI=KSI+IX(I)*KINC KJ=KSJ+JX(I)*KINC IF(ABS(KI-KJ).NE.KDEL) KJ=-KJ ANG=DI*KI+DJ*KJ IF(J.LE.10) THEN WK(I)=WK(I)*COS( ANG ) ELSE WK(I)=WK(I)*SIN( ANG ) ENDIF 40 CONTINUE RETURN END SUBROUTINE SZNZ(SIZ,NI,NJ,IXCOM,JXCOM,ISCOM,JSCOM) C formal parameters: INTEGER SIZ,NI,NJ,IXCOM(*),JXCOM(*),ISCOM(0:*),JSCOM(0:*) IF(SIZ.EQ.1) THEN IF(NI.LE.0.AND.NJ.LE.0) SIZ=0 RETURN ELSE IF(SIZ.LT.0) THEN IF(SIZ.LE.-4) THEN JSCOM(1)=NJ+NJ JXCOM(3)=NJ ELSE ISCOM(1)=NI+NI IXCOM(3)=NI ENDIF ELSE IF(SIZ.GE.4) THEN IF(SIZ.EQ.4) THEN IF(NJ.LE.1) THEN SIZ=0 RETURN ENDIF JSCOM(1)=JSCOM(1)-2 JXCOM(3)=JXCOM(3)-1 ELSE IF(NJ.LE.0) THEN SIZ=0 RETURN ENDIF JSCOM(1)=JSCOM(1)+2 JXCOM(3)=JXCOM(3)+1 ENDIF ELSE IF(SIZ.EQ.2) THEN IF(NI.LE.1) THEN SIZ=0 RETURN ENDIF ISCOM(1)=ISCOM(1)-2 IXCOM(3)=IXCOM(3)-1 ELSE IF(NI.LE.0) THEN SIZ=0 RETURN ENDIF ISCOM(1)=ISCOM(1)+2 IXCOM(3)=IXCOM(3)+1 ENDIF SIZ=-SIZ RETURN END SUBROUTINE SZNZOP(NI,NJ,KSI,KSJ,ISCOM,JSCOM,NCOS,WK,IX,JX) C formal parameters: INTEGER NI,NJ,KSI,KSJ,ISCOM(0:*),JSCOM(0:*),NCOS INTEGER*2 IX(*),JX(*) REAL*8 WK(*) C subroutine to add SzNz C on entry: C ISCOM,JSCOM= spins C NCOS= length of vector WK C WK= vector of a sub-diagonal C IX,JX= relative indices of elements in WK C on return: C WK= modified vector INTEGER KINC,NNI,NNJ,NNI0,NNJ0,I,IFLG,KK,KI,KJ,LL, & LSCOM(2),LLCOM(2),KSCOM(0:2) PARAMETER (KINC=2) REAL*8 C3JJ,VALI,VALJ,TWK DATA LSCOM,LLCOM/2,0,2,0/ SAVE LSCOM,LLCOM IF(NCOS.LE.0) RETURN IFLG=0 NNI=NI+NI IF(NNI.NE.ISCOM(1)) IFLG=1 NNJ=NJ+NJ IF(NNJ.NE.JSCOM(1)) IFLG=-1 LL=2 IF(IFLG.LE.0) THEN NNJ0=JSCOM(1) VALJ=NNJ0+1 IF(NNJ0.NE.NNJ) VALJ=SQRT(VALJ*(NNJ+1)) IF(IAND(NNJ0+KSJ+KSJ,3).NE.0) VALJ=-VALJ KSCOM(0)=JSCOM(0) KSCOM(1)=NNJ KSCOM(2)=JSCOM(2) CALL TENSOR(VALJ,1,JSCOM,KSCOM,LLCOM,LSCOM) ENDIF IF(IFLG.GE.0) THEN NNI0=ISCOM(1) VALI=NNI0+1 IF(NNI0.NE.NNI) VALI=SQRT(VALI*(NNI+1)) IF(IAND(NNI+KSI+KSI,3).NE.0) VALI=-VALI KSCOM(0)=ISCOM(0) KSCOM(1)=NNI KSCOM(2)=ISCOM(2) CALL TENSOR(VALI,1,KSCOM,ISCOM,LLCOM,LSCOM) ENDIF DO 40 I=1,NCOS TWK=WK(I) WK(I)=0. IF(IFLG.LE.0) THEN KJ=KSJ+JX(I)*KINC KK=KJ+KJ WK(I)=VALJ*TWK*KJ*C3JJ(NNJ0,LL,NNJ,-KK,0,KK) ENDIF IF(IFLG.GE.0) THEN KI=KSI+IX(I)*KINC KK=KI+KI WK(I)=WK(I)+VALI*TWK*KI*C3JJ(NNI,LL,NNI0,-KK,0,KK) ENDIF 40 CONTINUE RETURN END FUNCTION BLKSYM(IXCOM,JXCOM) INTEGER BLKSYM,IXCOM(6),JXCOM(6) C get block symmetry from state symmetry INTEGER SYM(0:15),KSYM SAVE SYM DATA SYM/0,1,2,3, & 1,0,3,2, & 2,3,0,1, & 3,2,1,0/ KSYM=IXCOM(2) KSYM=KSYM+KSYM KSYM=KSYM+KSYM+JXCOM(2) BLKSYM=SYM(KSYM) RETURN END SUBROUTINE ORDHAM(NDM,NN,T,EGY,NBLK,ISBLK,JDX,IORD) INTEGER NDM,NN,NBLK,IORD INTEGER*2 ISBLK(*),JDX(*) REAL*8 T(NDM,*),EGY(*) C subroutine to order eigenvalues within sub-block like diagonal of C Hamiltonian C INTEGER I1,IBLK,I,IS,IX,II,IIS,IIX,IQ PARAMETER (I1=1) REAL*8 TMP IF(IORD.NE.0) THEN C set up jdx to match egy IF(IORD.GT.0) THEN TMP=-1 CALL DSCAL(NN,TMP,EGY,I1) ENDIF DO 5 I=1,NN 5 JDX(I)=I DO 10 IBLK=1,NBLK IS=ISBLK(IBLK) IX=ISBLK(IBLK+1)-2 IIX=IX+1 DO 11 I=IS,IX IQ=I TMP=EGY(I) IIS=I+1 DO 12 II=IIS,IIX IF(EGY(II).LT.TMP) THEN TMP=EGY(II) IQ=II ENDIF 12 CONTINUE IF(IQ.NE.I) THEN EGY(IQ)=EGY(I) EGY(I)= TMP II=JDX(IQ) JDX(IQ)=JDX(I) JDX(I)=II ENDIF 11 CONTINUE 10 CONTINUE ELSE C swap to follow energy IIX=NN-1 DO 20 I=1,NN IF(JDX(I).NE.I) THEN IIS=I+1 DO 21 II=IIS,IIX IF(JDX(II).EQ.I) GO TO 22 21 CONTINUE II=NN 22 JDX(II)=JDX(I) TMP=EGY(II) EGY(II)=EGY(I) EGY( I)=TMP CALL DSWAP(NN,T(1,I),I1,T(1,II),I1) ENDIF 20 CONTINUE ENDIF RETURN END SUBROUTINE KROLL(NSIZD,T,NSBLK,SBKPTR,KMIN,KAPPA,ROLL) INTEGER NSIZD,NSBLK,KAPPA INTEGER*2 SBKPTR(*),KMIN(*) REAL*8 T(NSIZD,*) LOGICAL ROLL C subroutine to make sure that diagonal elements of the Hamiltonian C are monotonically increasing (decreasing for oblate) C after K=KTHRSH INTEGER I,IXX,IBGN,IEND,K,IFLG,N,KTHRSH,KINC,I1,I0 PARAMETER (KTHRSH=5,KINC=2,I0=0,I1=1) REAL*8 TMP,TLAST,VALL,ZERO(2) SAVE ZERO DATA ZERO/0.,0./ ROLL=.FALSE. DO 1 IXX=1,NSBLK I=KTHRSH-KMIN(IXX) IBGN=SBKPTR(IXX) IF(I.GE.KINC) IBGN=IBGN+I/KINC IEND=SBKPTR(IXX+1)-1 IF(IBGN.GE.IEND) GO TO 1 TMP=T(IBGN,IBGN) IBGN=IBGN+1 DO 10 I=IBGN,IEND IFLG=I IF(KAPPA.GT.0) IFLG=-IFLG TLAST=TMP TMP=T(I,I) IF(TMP.GT.TLAST) IFLG=-IFLG IF(IFLG.GT.0) GO TO 12 10 CONTINUE 12 IF(IFLG.LE.0) GO TO 1 K=KMIN(IXX)+KINC*(IFLG-SBKPTR(IXX)) WRITE(*,*) ' roll-over at K = ',K ROLL=.TRUE. K=NSIZD VALL=1.E+15 IF(KAPPA.GT.0) VALL=-VALL TMP=TLAST DO 11 I=IFLG,IEND N=I-1 CALL DCOPY(N,ZERO,I0,T(I,1),K) N=K-I CALL DCOPY(N,ZERO,I0,T(I+1,I),I1) TMP=TMP+VALL 11 T(I,I)=TMP 1 CONTINUE RETURN END SUBROUTINE GETQS(IM,IFF,NSIZ,KBGN,IXCOM,ISCOM,IV) C formal parameters: INTEGER IM,IFF,NSIZ,KBGN,IXCOM(6),ISCOM(0:*),IV C subroutine to get quanta corresponding to sub-block C on entry: C IM= position in sub-block mold C IFF = F quantum number *2 C NSIZ= size of sub-block C KBGN= first K value C on return: C IXCOM = list of rotational containing C 1: DIMENSION OF SUB-BLOCK C 2: SYMMETRY CODE FOR BLOCK (0=A,1=Bx,2=By,3=Bz) C 3: 'N' QUANTUM NUMBER C 4: BEGINNING K C 5: L VALUE C 6: VIBRATIONAL QUANTA C ISCOM = list of angular momenta C IV = vibrational quantum number if L=0,1 C IV = vibrational quantum number-1 if L=-1 INTEGER MXSPIN,NDSPIN,NX PARAMETER (MXSPIN=5,NDSPIN=MXSPIN*9,NX=170) INTEGER*2 NVSM,JJB(0:NX,MXSPIN),MINF(0:NX),IIVPT(0:98),IISPT(0:98) : ,LVQN(0:98) INTEGER I,ISPX,LV INTEGER KAPPA,NSPIN,NTSPIN,IISPN(NDSPIN),ITPTR,ITPTRM,IDIAG,NVIB COMMON /SQN/ KAPPA,NSPIN,NTSPIN,IISPN,ITPTR,ITPTRM,IDIAG,NVIB COMMON /MOLDS/ NVSM,JJB,MINF,IIVPT,IISPT,LVQN SAVE /SQN/,/MOLDS/ C find list of angular momenta ISPX=IM/NVSM DO 10 I=1,NSPIN ISCOM(I)=JJB(ISPX,I) IF(I.NE.ITPTR) ISCOM(I)=ISCOM(I)+IFF 10 CONTINUE ISCOM(NSPIN+1)=IFF IXCOM(1)=NSIZ C find symmetry IXCOM(2)=IAND(IM,3) C find N IXCOM(3)=ISCOM(1)/2 IXCOM(4)=KBGN IV=(IM-NVSM*ISPX)/4 IXCOM(6)=IV LV=LVQN(IV) IF(LV.LT.0) THEN LV=-LV IV=IV-1 ENDIF IXCOM(5)=LV ISCOM(0)=IISPT(IV) IM=ISPX-IIVPT(IV) RETURN END SUBROUTINE IDPARS(IDPAR,KSQ,ITP,L,LD,KDEL,INS,SI1,SI2,SIZ,KLOFF) C formal parameters: INTEGER*4 IDPAR INTEGER KSQ,ITP,L,LD,KDEL,INS,SI1,SI2,SIZ,KLOFF C subroutine to parse parameter identifiers to subfields C on entry: C IDPAR= 10*( parameter ID/1000) + ITP C on return: C KSQ= power of K*K C ITP= parameter subtype C L= tensor order C LD= order of direction cosine ( L <= LD ) C KDEL= change in K C INS=power of N.S or fourier coefficient C SI1,SI2=spin identifiers INTEGER*4 L10K PARAMETER (L10K=10000) INTEGER IDHIGH,ITY,II,ISYM,IDLOW,IX SIZ=0 KLOFF=0 ITP=0 IDHIGH=IDPAR/L10K IF(IDHIGH.EQ.0) THEN IDLOW=IDPAR SI2=0 SI1=0 INS=0 ELSE IDLOW=IDPAR-L10K*IDHIGH SI2=IDHIGH/100 II=IDHIGH-100*SI2 SI1=II/10 INS=II-10*SI1 IF(SI2.GE.10) THEN II=SI2/10 SI2=SI2-10*II ITP=II*16 ENDIF ENDIF ITY=IDLOW/100 II =IDLOW-ITY*100 KSQ=II/10 IX=II-10*KSQ IF(IX.GE.5) THEN IX=IX-5 KLOFF=1 ENDIF LD=0 KDEL=0 IF(ITY.LE.3) THEN IF(ITY.NE.0) LD=2 ELSE IF(ITY.LT.20) THEN KDEL=ITY+ITY-6 LD=KDEL ELSE IF(ITY.LT.80) THEN ISYM=ITY/20 LD=ITY-ISYM*20+1 KDEL=LD-IAND(ISYM/2+LD,1) ELSE IF(ITY.LT.90) THEN C K energies coded as 8n with K=10*n+KSQ KSQ=-1-(ITY-80)*10-KSQ ELSE ITY=ITY-90 ITP=ITP+ITY+2 KDEL=ITY/2 KDEL=KDEL+KDEL IF(KDEL.NE.ITY) KDEL=0 ENDIF IF(SI1.EQ.0) THEN L=0 ELSE IF(SI2.EQ.0) THEN L=1 ELSE IF(LD.LT.2) THEN L=0 ELSE L=2 ENDIF IF(INS.GE.5) THEN INS=INS-5 SIZ=IX+1 ENDIF C .. try to use raising operator IF(LD.EQ.KDEL.AND.LD.GT.L) LD=L RETURN END FUNCTION GETLL(LD,KD,L,SI1,SI2,LLCOM,LSCOM,ISCOM,JSCOM) C formal parameters: INTEGER GETLL,LD,KD,L,SI1,SI2,LLCOM(*),LSCOM(*), : ISCOM(0:*),JSCOM(0:*) C function to find sequence of L values for coupling of spins C on entry: C LD = order of direction cosine C L = order of N dependent tensor C SI1, SI2 = position of spin values C ISCOM,JSCOM = list of angular momenta *2 C on return: C LLCOM = list of tensor orders for succesive couplings *2 C LSCOM = list of tensor orders for spins *2 C SI1, SI2 = position of spin values with possible offsets C return true if delta QN not consistent with L INTEGER LLJ,LLX,LLS,LLIT,I,NCHK INTEGER MXSPIN,NDSPIN PARAMETER (MXSPIN=5,NDSPIN=MXSPIN*9) INTEGER KAPPA,NSPIN,NTSPIN,IISPN(NDSPIN),ITPTR,ITPTRM,IDIAG,NVIB COMMON /SQN/ KAPPA,NSPIN,NTSPIN,IISPN,ITPTR,ITPTRM,IDIAG,NVIB SAVE /SQN/ GETLL=1 C check for scalar if delta N is 0 I=LD IF(I.LT.KD) I=KD IF(I.EQ.0) THEN IF(ISCOM(1).NE.JSCOM(1)) RETURN ELSE IF(I.GT.L) THEN IF(I+I.GT.ISCOM(1)+JSCOM(1)) RETURN ENDIF NCHK=NSPIN+1 DO 2 I=1,NCHK LSCOM(I)=0 IF(SI1.EQ.I) LSCOM(I)=2 IF(SI2.EQ.I) LSCOM(I)=LSCOM(I)+2 2 CONTINUE LLJ=L+L LLIT=0 IF(SI1.GE.ITPTRM) THEN C .. do ITOT tensor types , SI2 le SI1 LLIT=LSCOM(ITPTRM)+LSCOM(ITPTR) IF(LLIT.EQ.4) LLIT=LLJ ENDIF C generate list of tensor orders, checking quantum numbers along the C way DO 10 I=1,NCHK IF(I.EQ.ITPTR) LLJ=LLIT LLCOM(I)=LLJ IF(LLJ.LT.ABS(ISCOM(I)-JSCOM(I)) ) RETURN IF(LLJ.GT. ISCOM(I)+JSCOM(I) ) RETURN IF(I.NE.ITPTR) THEN LLX=LLJ LLS=LSCOM(I) IF(SI2+I.EQ.0) LLS=0 ELSE LLS=LLX ENDIF LLJ=ABS(LLS-LLJ) 10 CONTINUE GETLL=0 RETURN END FUNCTION RMATRX(LD,L,IXCOM,JXCOM) INTEGER LD,L,IXCOM(*),JXCOM(*),NI,NJ,NDIF,LT,LMIN,LMAX,LBGN REAL*8 RMATRX,DN,DNSQ,DLSQ,DEL,TMP C .. find N corrections to reduced matrix elements IF(LD.NE.L) THEN NI=IXCOM(3) NJ=JXCOM(3) DN=NI+NJ+1 LMAX=LD LMIN=L IF(LMAX.LT.LMIN) THEN LMIN=LMAX LMAX=L ENDIF NDIF=NI-NJ IF(NDIF.NE.0) DEL=NDIF*NDIF DNSQ=DN*DN TMP=1. LBGN=2 IF(LMIN.GT.1) LBGN=LMIN+1 DO 10 LT=LBGN,LD DLSQ=LT*LT TMP=TMP*0.25*(DNSQ-DLSQ) IF(NDIF.NE.0) TMP=TMP-TMP*DEL/DLSQ 10 CONTINUE IF(LMIN.EQ.0) TMP=TMP*DN*NI*(NI+1) RMATRX=SQRT(TMP) IF(L.EQ.0) RMATRX=RMATRX/DN ELSE RMATRX=1. ENDIF RETURN END SUBROUTINE SYMNSQ (INQ,INS,ISCOM,JSCOM,ISPN,Z) C formal parameters: INTEGER INQ,INS,ISCOM(0:*),JSCOM(0:*),ISPN REAL*8 Z C subroutine to find powers of N*(N+1) and N.S C on entry: C INQ= power of N*(N+1) + 1 C INS= power of N.S C ISCOM,JSCOM= list of angular momenta *2 C ISPN = S *2 C Z= matrix multipling factor C on return: C Z= matrix factor REAL*8 X,DOTNS,ZR INTEGER I LOGICAL MODF MODF=ISCOM(1).NE.JSCOM(1).AND.INQ.GT.1 ZR=Z IF(MODF) THEN C CORRECT ELE FOR COMMUTATION X=ISCOM(1)*(ISCOM(1)+2) X=JSCOM(1)*(JSCOM(1)+2)/X DO 10 I=2,INQ 10 ZR=ZR*X ENDIF IF(INS.GT.0) THEN DOTNS=0.125*( ISCOM(2)*(ISCOM(2)+2)-ISPN*(ISPN+2) & -ISCOM(1)*(ISCOM(1)+2) ) X =0.125*( JSCOM(2)*(JSCOM(2)+2)-ISPN*(ISPN+2) & -JSCOM(1)*(JSCOM(1)+2) ) DO 20 I=1,INS Z=Z*DOTNS 20 ZR=ZR*X MODF=.TRUE. ENDIF IF(MODF) Z=0.5*(Z+ZR) RETURN END SUBROUTINE SYMKSQ (IKQ,KSI,KSJ,N,WK,IX,JX) C formal parameters: INTEGER IKQ,KSI,KSJ,N INTEGER*2 IX(*),JX(*) REAL*8 WK(*) C subroutine to multiply right and left by a power of K*K C on entry: C IKQ= power of K*K C KSI,KSJ= starting values of K C N= length of vector WK C WK= vector of a sub-diagonal C IX,JX= relative indices of elements in WK C note: if IKQ < 0 then operator is K*K only when K=-1-IKQ C on return: C WK= modified vector INTEGER I,J,KI,KJ,KCMP,KINC,NT PARAMETER (KINC=2) REAL*8 AKISQ,AKJSQ,XI,XJ KCMP=-1-IKQ NT=N DO 20 I=1,NT KI=KSI+IX(I)*KINC AKISQ=KI*KI KJ=KSJ+JX(I)*KINC IF(KI.NE.KJ) THEN AKJSQ=KJ*KJ XI=0.5 XJ=0.5 DO 10 J=1,IKQ XI=AKISQ*XI 10 XJ=AKJSQ*XJ WK(I)=WK(I)*(XI+XJ) ELSE IF(IKQ.GT.0) THEN DO 15 J=1,IKQ 15 WK(I)=WK(I)*AKISQ ELSE IF(KI.EQ.KCMP) THEN IX(1)=IX(I) JX(1)=JX(I) N=1 RETURN ENDIF 20 CONTINUE IF(IKQ.LT.0) N=0 RETURN END SUBROUTINE DPMAKE(NSIZD,DP,T,N,WK,IDX,JDX) C find derivative contribution from sub-diagonal INTEGER NSIZD,N INTEGER*2 IDX(*),JDX(*) REAL*8 DP(*),T(*),WK(*) INTEGER I0,I1 PARAMETER (I0=0,I1=1) INTEGER NSIZE,I,IZ,JZ,K REAL*8 ELE,ZERO(2) SAVE ZERO DATA ZERO/0.,0./ NSIZE=NSIZD CALL DCOPY(NSIZE,ZERO,I0,DP,I1) DO 701 I=1,N IZ=IDX(I) JZ=JDX(I) ELE=WK(I) IF(IZ.NE.JZ) THEN ELE=ELE+ELE DO 700 K=1,NSIZE DP(K)=DP(K)+T(IZ)*T(JZ)*ELE IZ=IZ+NSIZE 700 JZ=JZ+NSIZE ELSE DO 710 K=1,NSIZE DP(K)=DP(K)+T(IZ)*T(IZ)*ELE 710 IZ=IZ+NSIZE ENDIF 701 CONTINUE RETURN END SUBROUTINE TENSOR(Z,NS,ISCOM,JSCOM,LLCOM,LSCOM) C formal parameters: INTEGER NS,ISCOM(0:*),JSCOM(0:*),LLCOM(*),LSCOM(*) REAL*8 Z C subroutine to find spherical tensor spin coupling coefficients C on entry: C Z= initial multiplier C NS= number of spins C ISCOM,JSCOM= list of angular momenta *2 C LLCOM= list of tensor orders *2 C LSCOM = list of tensor orders for spins *2 C on return: C Z= modified multiplier REAL*8 ZSQ,C6JJ,C9JJ,ONE INTEGER MXSPIN,NDSPIN PARAMETER (MXSPIN=5,NDSPIN=MXSPIN*9,ONE=1.) INTEGER SSF,SSI,ISGN,LLN,LLS,NNI,NNF,I,LLJ,JJI,JJF,N,ORGI,ORGF INTEGER KAPPA,NSPIN,NTSPIN,IISPN(NDSPIN),ITPTR,ITPTRM,IDIAG,NVIB COMMON /SQN/ KAPPA,NSPIN,NTSPIN,IISPN,ITPTR,ITPTRM,IDIAG,NVIB SAVE /SQN/ IF(NS.LE.0) RETURN IF(NS.EQ.ITPTRM) THEN N=ITPTR ELSE N=NS ENDIF ZSQ=1 ISGN=0 LLN=LLCOM(1) NNI=ISCOM(1) NNF=JSCOM(1) ORGI=ISCOM(0) ORGF=JSCOM(0) DO 300 I=1,N LLJ=LLCOM(I+1) JJI=ISCOM(I+1) JJF=JSCOM(I+1) SSF=IISPN(I+ORGF) SSI=IISPN(I+ORGI) IF(I.EQ.ITPTRM) THEN LLN=LSCOM(ITPTR) NNI=IISPN(ITPTR+ORGI) NNF=IISPN(ITPTR+ORGF) ELSE IF(I.EQ.ITPTR) THEN SSI=ISCOM(ITPTRM) SSF=JSCOM(ITPTRM) ENDIF IF(LLN.GT.0) THEN IF(LLJ.EQ.0) THEN C ....T * U ISGN=ISGN+JJI+SSI+NNF Z=Z*C6JJ(NNI,LLN,NNF,SSF,JJI,SSI) ELSE IF(LSCOM(I).NE.0) THEN C ....T X U ZSQ=ZSQ*(LLJ+ONE)*(JJI+ONE)*(JJF+ONE) IF(I.EQ.ITPTRM) THEN ZSQ=ZSQ*1.5 ELSE ISGN=ISGN+2 ZSQ=ZSQ*2.5 ENDIF LLS=LSCOM(I) Z=Z*C9JJ(NNF,SSF,JJF,LLN,LLS,LLJ,NNI,SSI,JJI) ELSE C ....U = 1 ISGN=ISGN+JJF+NNI+SSI+LLJ ZSQ=ZSQ*(JJI+ONE)*(JJF+ONE) Z=Z*C6JJ(NNI,LLJ,NNF,JJF,SSI,JJI) ENDIF ELSE IF(LLJ.GT.0) THEN C ....T = 1 ISGN=ISGN+JJI+NNI+SSF+LLJ ZSQ=ZSQ*(JJI+ONE)*(JJF+ONE) Z=Z*C6JJ(SSI,LLJ,SSF,JJF,NNI,JJI) ENDIF NNI=JJI NNF=JJF LLN=LLJ 300 CONTINUE ISGN=ISGN/2 IF(IAND(ISGN,1).NE.0) Z=-Z Z=Z*SQRT(ZSQ) RETURN END SUBROUTINE SETINT(IFDIAG,NSAV,IDIP) C formal parameters: INTEGER IFDIAG,NSAV INTEGER*4 IDIP(0:*),ITMP,L10K PARAMETER (L10K=10000) C subroutine to initialize intensity data INTEGER NDIP,I,II,IV1,IV2,MXSPIN,NDSPIN,NX,ICASE,ISYM,LV1,LV2 PARAMETER (MXSPIN=5,NDSPIN=MXSPIN*9,NX=170) INTEGER KAPPA,NSPIN,NTSPIN,IISPN(NDSPIN),ITPTR,ITPTRM,IDIAG,NVIB INTEGER*2 NVSM,JJB(0:NX,MXSPIN),MINF(0:NX),IIVPT(0:98),IISPT(0:98) : ,LVQN(0:98) COMMON /MOLDS/ NVSM,JJB,MINF,IIVPT,IISPT,LVQN COMMON /SQN/ KAPPA,NSPIN,NTSPIN,IISPN,ITPTR,ITPTRM,IDIAG,NVIB SAVE /SQN/,/MOLDS/ IFDIAG=IDIAG+1 NSAV=1 NDIP=IDIP(0) DO 1 I=1,NDIP ITMP=IDIP(I)/10 ISYM=IDIP(I)-ITMP*10 IF(NVIB.GT.9) THEN ICASE=ITMP/L10K II=ITMP-ICASE*L10K IV2=II/100 IV1=II-100*IV2 ELSE ICASE=ITMP/100 II=ITMP-ICASE*100 IV2=II/10 IV1=II-10*IV2 ENDIF C MAKE SURE VIBRATIONS ARE ORDERED IF(IV2.GT.IV1) THEN II=IV1 IV1=IV2 IV2=II ENDIF LV1=LVQN(IV1) IF(LV1.LT.0) THEN LV1=-LV1 IV1=IV1-1 ENDIF LV2=LVQN(IV2) IF(LV2.LT.0) THEN LV2=-LV2 IV2=IV2-1 ENDIF C .. CONVERT SYMMETRY IF(ISYM.GT.3) ISYM=0 IF(ISYM.NE.0.AND.KAPPA.LT.0) ISYM=4-ISYM IF(ICASE/10.EQ.8) THEN IF(LV1.EQ.0.OR.LV1.NE.LV2) ISYM=3-ISYM ENDIF IF(LV2.EQ.2.AND.LV1.NE.1) ISYM=4 IF(LV1.EQ.2.AND.LV2.NE.1) ISYM=4 IDIP(I)=(ICASE*L10K+100*IV2+IV1)*5+ISYM IF(LV1.NE.LV2) IDIP(I)=-IDIP(I) 1 CONTINUE RETURN END C FUNCTION INTENS (IBLK,ISIZ,JBLK,JSIZ,IDIP,DIP,S,NSIZD) C formal parameters: INTEGER INTENS,IBLK,JBLK,ISIZ,JSIZ,NSIZD INTEGER*4 IDIP(0:*) REAL*8 DIP(*),S(NSIZD,*) C function to calculate intensities C on entry: C IBLK,JBLK= block numbers C ISIZ,JSIZ= block sizes C IDIP= dipole id codes , IDIP(0)= number of dipoles C DIP= dipole values C on return: C S= transition dipole matrix in original basis C INTENS returns true is S is not 0 INTEGER*4 LTMP,L50K,IDV,IDVP PARAMETER (L50K=50000) INTEGER NDMX,MAXS,I1,MXSPIN,NDSPIN,ITY PARAMETER (NDMX=450,MAXS=64,MXSPIN=5,NDSPIN=MXSPIN*9,I1=1) INTEGER IXCOM(6),JXCOM(6),ISCOM(0:MXSPIN+1),JSCOM(0:MXSPIN+1), & LLCOM(MXSPIN+1),LSCOM(MXSPIN+1),SI2,KDV(0:3), & GETLL,GETQQ,SI1,NOLD,NBLKI,NBLKJ,IFF,JFF,IWT,JWT,N,KL, & NDIP,IXX,JXX,KBGNI,KBGNJ,IXTMP,IVV,JVV,NNI,NNJ,BLKSYM, & NX,I,ICASE,KD,L,LD,NCOS,K,IX,JX,ISYM,IFUP,IXDIP,I0, & IBASE,JBASE PARAMETER (I0=0) INTEGER*2 IBKPTR(MAXS+1),IKMIN(MAXS),IVS(MAXS),JBKPTR(MAXS+1), & JKMIN(MAXS),JVS(MAXS),IDX(NDMX),JDX(NDMX) REAL*8 SPFAC(0:NDSPIN),DRCS(NDMX),DD,DNI,RMATRX INTEGER KAPPA,NSPIN,NTSPIN,IISPN(NDSPIN),ITPTR,ITPTRM,IDIAG,NVIB COMMON /SQN/ KAPPA,NSPIN,NTSPIN,IISPN,ITPTR,ITPTRM,IDIAG,NVIB COMMON /SPFACV/ SPFAC,NOLD SAVE /SQN/,/SPFACV/ SAVE KDV DATA KDV/0,1,1,0/ INTENS=0 C get quantum information NBLKI=GETQQ(MAXS,IBLK,IFF,IWT,IBKPTR,IKMIN,IVS) NBLKJ=GETQQ(MAXS,JBLK,JFF,JWT,JBKPTR,JKMIN,JVS) IF(IWT.NE.JWT) RETURN IF(ABS(IFF-JFF).GT.2) RETURN IXX=IFF+JFF IF(IXX.LT.2) RETURN IF(IAND(IXX,1).NE.0) RETURN C clear dipole matrix CALL DCLR(ISIZ,JSIZ,S,I1) NDIP=IDIP(0) C loop over sub-blocks DO 390 IXX=1,NBLKI IBASE=IBKPTR(IXX) N=IBKPTR(IXX+1)-IBASE KBGNI=IKMIN(IXX) IXTMP=IVS(IXX) CALL GETQS(IXTMP,IFF,N,KBGNI,IXCOM,ISCOM,IVV) NNI=ISCOM(1) DNI=NNI+1 DO 380 JXX=1,NBLKJ JBASE=JBKPTR(JXX) N=JBKPTR(JXX+1)-JBASE KBGNJ=JKMIN(JXX) IXTMP=JVS(JXX) CALL GETQS(IXTMP,JFF,N,KBGNJ,JXCOM,JSCOM,JVV) NNJ=JSCOM(1) NX=MAX(NNI,NNJ)/2 ISYM=BLKSYM(IXCOM,JXCOM) C .. ISYM is the symmetry of the transition LTMP=IVV+JVV+ 99*MIN(IVV,JVV) IDV=ISYM+5*LTMP C.. set dipole for transition DO 200 I=1,NDIP IXDIP=IDIP(I) KL=0 IF(IXDIP.LT.0) THEN IXDIP=-IXDIP KL=1 ENDIF C.. isolate ICASE*10+SI1 ITY=IXDIP/L50K C.. eliminate trivial SIN(0) Fourier term IF(ITY.EQ.80) GO TO 200 IDVP=IXDIP-ITY*L50K IF(IDVP.NE.IDV) GO TO 200 C .. decode dipole id ICASE=ITY/10 SI1=ITY-10*ICASE IF(ICASE.EQ.7.OR.ICASE.EQ.8) SI1=0 SI2=0 KD=KDV(ISYM) L=1 IF(ICASE.EQ.1.OR.ICASE.EQ.2) THEN C .. anisotropies LD=2 IF(SI1.NE.0) L=2 IF(ICASE.EQ.1) THEN IF(ISYM.EQ.3) KD=2 ELSE IF(ISYM.EQ.0) KD=2 ENDIF ELSE IF(ISYM.NE.0) THEN C .. simple electric dipole LD=1 IF(ICASE.EQ.6) KD=KD+2 C .. neg SI2 signals vector cross product SI2=-SI1 ELSE C .. isotropic magnetic dipole LD=0 IF(SI1.NE.0) L=0 ENDIF C .. find matrix elements IF(GETLL(LD,KD,L,SI1,SI2,LLCOM,LSCOM,ISCOM,JSCOM) & .NE.0) GO TO 200 C check for real or imaginary dipoles IFUP=-1 DD=DIP(I) IF(ISYM.EQ.0) THEN C .. correct for magnetic dipole DD=DD*0.009274 IFUP=1 ENDIF NCOS=NDMX IF(ICASE.EQ.8) NCOS=-NCOS CALL DIRCOS(IXCOM,JXCOM,LD,KD,NCOS,DRCS,IDX,JDX,IFUP,KL) IF(NCOS.LE.0) GO TO 200 C .. correct for special cases IF(ICASE.EQ.2) THEN IF(SI1.EQ.0) DD=-DD ELSE IF(ICASE.EQ.3.OR.ICASE.EQ.9) THEN IF(NNI.EQ.NNJ) THEN K = NX*(NX+1) DD = DD*K IF(ICASE.EQ.9) DD = DD*K ELSE K = NX*NX DD = DD*K IF(ICASE.EQ.9) DD = DD*(K+1) ENDIF ELSE IF(ICASE.EQ.4) THEN CALL SYMKSQ(I1,KBGNI,KBGNJ,NCOS,DRCS,IDX,JDX) ELSE IF(ICASE.EQ.5.OR.ICASE.EQ.10) THEN IF(NNI.EQ.NNJ) GO TO 200 IF(IVV.GT.JVV) DD = -DD IF(NNI.LT.NNJ) DD = -DD DD = DD*NX IF(ICASE.EQ.10) DD = DD*(NX*NX) ELSE IF(ICASE.EQ.7.OR.ICASE.EQ.8) THEN ITY=(ITY-70)*16 CALL SPECFC(ITY,IVV,JVV,KD,KBGNI,KBGNJ,NCOS,DRCS, & IDX,JDX) ENDIF C .. correct for reduced matrix of N IF(LD.NE.L) THEN DD=DD*RMATRX(LD,L,IXCOM,JXCOM) ENDIF IF(SI1.GT.0) DD=DD*SPFAC(SI1) C .. couple dipoles through the spins CALL TENSOR(DD,NSPIN,ISCOM,JSCOM,LLCOM,LSCOM) DO 300 K=1,NCOS IX=IDX(K)+IBASE JX=JDX(K)+JBASE 300 S(IX,JX)=S(IX,JX)+DD*DRCS(K) INTENS=IFF+1 200 CONTINUE 380 CONTINUE 390 CONTINUE RETURN END FUNCTION SETOPT (LUIN,IQFMT,ITD,NAMFIL) C formal parameters: CHARACTER*80 NAMFIL INTEGER SETOPT,LUIN,IQFMT,ITD C this subroutine reads in quantum quantities C on entry: C LUIN is the fortran unit number to read from C on return: C IQFMT= quantum number format for catalog C ITD= 2 for linear, 3 for non-linear C NAMFIL contains a file name on which to find parameter names C SETOPT = number of cards read, -1 for end-of-file INTEGER*4 VSYM,ISPNX INTEGER KINC,MAXV,MAXQ,MXSPIN,NDSPIN,NX,NOPTS,J, & PCARD,NCARDS,MAXVS PARAMETER (KINC=2,MAXV=99,MAXQ=500,MXSPIN=5,NDSPIN=MXSPIN*9) PARAMETER (NOPTS=11,MAXVS=MAXV*4,NX=170) INTEGER*2 IIWT(MAXVS) REAL*8 RVEC(NOPTS) INTEGER*2 IWT(MAXQ),BLKPTR(MAXQ+1),MOLDVS(MAXQ) INTEGER NQNN,NQN,IFVIB,ESYMWT(0:98),LOPT,KNNMIN,KNNMAX,IAX, & IWTPL,IWTMN,NSSTAT(0:98),SI2,I,IVIB,II,NT,NSSMAX,K, & IREV(0:11),IOFF,NVIBM,KLSYM INTEGER NBKPJ,KNMIN(0:3,0:98),KNMAX(0:98),IXZ CHARACTER CTYP*1 CHARACTER*80 CARD INTEGER*2 NVSM,JJB(0:NX,MXSPIN),MINF(0:NX),IIVPT(0:98),IISPT(0:98) : ,LVQN(0:98) INTEGER KAPPA,NSPIN,NTSPIN,IISPN(NDSPIN),ITPTR,ITPTRM,IDIAG,NVIB COMMON /MOLDS/ NVSM,JJB,MINF,IIVPT,IISPT,LVQN COMMON /SQN/ KAPPA,NSPIN,NTSPIN,IISPN,ITPTR,ITPTRM,IDIAG,NVIB COMMON /MOLD/NBKPJ,KNMIN,KNMAX,IXZ,IWT,BLKPTR,MOLDVS COMMON /SQFMT/NQNN,NQN,IFVIB,KLSYM,ESYMWT SAVE /SQN/,/MOLD/,/MOLDS/,/SQFMT/ DATA IREV/2,0,2,0, 2,0,2,0, 3,3,-1,-1/ CALL CJJINI RVEC(1)=0 RVEC(2)=1 RVEC(3)=0 RVEC(4)=99 RVEC(5)=0 RVEC(6)=1 RVEC(7)=1 RVEC(8)=1 RVEC(9)=0 RVEC(10)=-1 RVEC(11)=0 NCARDS=0 ITD=2 VSYM=-1 C read option cards 1 IF(VSYM.LT.0) THEN READ(LUIN,'(A)',END=999) CARD NCARDS=NCARDS+1 RVEC(9)=0 I=NOPTS I=PCARD(CARD,RVEC,I) ISPNX= RVEC(1) LOPT= RVEC(2) KNNMIN=RVEC(3) KNNMAX=RVEC(4) IAX=RVEC(6) IWTPL= RVEC(7) IWTMN= RVEC(8) VSYM= RVEC(9) IVIB=ABS(LOPT) IF(KNNMAX.LT.0) KNNMAX=99 IF(KNNMIN.GT.KNNMAX) KNNMAX=KNNMIN IF(KNNMIN.NE.KNNMAX) ITD=3 IF(NCARDS.EQ.1) THEN CTYP=CARD(1:1) IF(CTYP.LT.'A'.AND.CTYP.GT.'9') CTYP=CARD(2:2) NAMFIL='sping.nam' IF(CTYP.GE.'A'.AND.CTYP.LE.'z') NAMFIL(5:5)=CTYP IXZ= RVEC(5) IDIAG= RVEC(11) IF(IDIAG.LT.0) IDIAG=-1 NVIB=MIN(IVIB,MAXV) NVIB=MAX(NVIB,1) NVIBM=NVIB-1 IVIB=0 NVSM=4*NVIB ITPTR=-IAX ITPTRM=MXSPIN+1 IF(LOPT.GE.0) THEN KAPPA=-1 ELSE KAPPA=1 ENDIF NQNN=1 IF(ISPNX.GE.0) NQNN=3 ELSE C if IVIB out of range read another card IF(IVIB.GE.NVIB) GO TO 1 ENDIF IAX=ABS(IAX) IF(KAPPA.GT.0) IAX=IAX+2-2*MOD(IAX+2,3) ISPNX=ABS(ISPNX) C set up minimum K values KNMAX(IVIB)=KNNMAX KNMIN(0,IVIB)=0 KNMIN(1,IVIB)=1 KNMIN(2,IVIB)=1 KNMIN(3,IVIB)=2 ESYMWT(IVIB)=RVEC(10) LVQN(IVIB)=MOD(ESYMWT(IVIB)/100,3) IF(LVQN(IVIB).EQ.2) THEN KNMIN(3,IVIB)=0 ELSE IF(LVQN(IVIB).EQ.1) THEN KNMIN(0,IVIB)=2 ELSE IF(KNNMAX.EQ.KNNMIN) THEN ESYMWT(IVIB)=-1 ENDIF DO 3 I=0,3 4 IF(KNMIN(I,IVIB).LT.KNNMIN) THEN KNMIN(I,IVIB)=KNMIN(I,IVIB)+KINC GO TO 4 ENDIF 3 CONTINUE C fill in defaults IF(NCARDS.EQ.1) THEN IVIB=-NVIB NT=0 CALL SETSP(ISPNX,IVIB,NT,NSSTAT) NSSMAX=NSSTAT(0) DO 10 I=1,NVIBM NSSTAT(I)=NSSTAT(0) KNMAX(I)=KNNMAX KNMIN(0,I)=KNMIN(0,0) KNMIN(1,I)=KNMIN(1,0) KNMIN(2,I)=KNMIN(2,0) KNMIN(3,I)=KNMIN(3,0) ESYMWT(I)=ESYMWT(0) LVQN(I)=LVQN(0) 10 CONTINUE ELSE CALL SETSP(ISPNX,IVIB,NT,NSSTAT) NSSMAX=MAX(NSSMAX,NSSTAT(IVIB)) ENDIF IF(KNNMAX.NE.0.AND.NQNN.EQ.1) NQNN=2 C set up weights CALL SETWT(IVIB,IIWT,IAX,IWTPL,IWTMN,VSYM) GO TO 1 ENDIF C end of card reading II=0 IF(LVQN(0).GT.0) II=1 C make sure +/- l values come in pairs DO 12 I=1,NVIBM IF(LVQN(I).GT.0) THEN IF(LVQN(I)+LVQN(I-1).EQ.3) THEN LVQN(I)=-LVQN(I) II=II-1 ELSE II=II+1 ENDIF ENDIF 12 CONTINUE IF(II.NE.0) THEN WRITE(*,*) 'L DOUBLETS SHOULD BE IN PAIRS' STOP ENDIF NT=NVSM-1 IAX=4*MOD(IAX+2,3) NBKPJ=0 NSSMAX=NSSMAX-1 DO 13 I=0,NSSMAX C set up spin weights and block code DO 15 II=0,NT IVIB=II/4 IF(I.LT.NSSTAT(IVIB)) THEN NBKPJ=NBKPJ+1 IF(NBKPJ.LE.MAXQ) THEN K=I+IIVPT(IVIB) MOLDVS(NBKPJ)=II+NVSM*K IOFF=1 IF(ITPTR.GT.0) THEN J=JJB(K,ITPTR) IF(IAND(J,3).EQ.2) IOFF=IREV(IAND(II,3)+IAX) ENDIF IWT(NBKPJ)=IIWT(II+IOFF) ENDIF ENDIF 15 CONTINUE 13 CONTINUE IF(NBKPJ.GT.MAXQ) THEN WRITE(*,*)'NUMBER OF SUB-BLOCKS PER F IS TOO BIG ',NBKPJ,MAXQ STOP ENDIF C set up format IF(NQNN.EQ.3) ITD=3 NQN=NQNN IFVIB=NVIBM IF(IFVIB.EQ.0) THEN IQFMT=NQN*101 ELSE NQN=NQN+1 IQFMT=NQN*101+1000 ENDIF IF(ITPTR.GT.0) IQFMT=IQFMT+2000 SI2=0 DO 2 I=2,NSPIN SI2=SI2+IAND(JJB(0,I)+MINF(0),1) SI2=SI2+SI2 2 CONTINUE SI2=SI2+IAND(MINF(0),1) IF(SI2.GT.9) SI2=IAND(SI2,7) NQN=NQN+NSPIN IF(NQN.LE.6) THEN IQFMT=IQFMT+SI2*10+NSPIN ELSE IQFMT=IQFMT+IAND(SI2,1)*10+4002 ENDIF SETOPT=NCARDS RETURN 999 SETOPT=-1 RETURN END FUNCTION SETBLK(LU,IDPAR,PAR,NBKPF,NEGY) C formal parameters: INTEGER SETBLK,LU,NBKPF,NEGY REAL*8 IDPAR(0:*),PAR(*) C this subroutine initializes quantum quantities C on entry: C LU= logical unit for listing C IDPAR= parameter id, IDPAR(0)= number of parameters C NBKPF= maximum value of F C NEGY = max dimensioned energy in main program C on return: C NBKPF= blocks per F C NEGY= possibly downsized NEGY to account for local storage C INTEGER NDPAR PARAMETER (NDPAR=750) INTEGER*4 IPTYP(NDPAR) INTEGER*2 NPARX,NJQ(NDPAR),IP(NDPAR),IPSYM(NDPAR), & IPDER(0:NDPAR),IVSYM INTEGER MAXQ,KINC,MXSPIN,NDSPIN,NX PARAMETER (MAXQ=500,KINC=2,MXSPIN=5,NDSPIN=MXSPIN*9,NX=170) REAL*8 SQNN,SQJ(10) INTEGER NPAR,I,ITY,IKQ,L,LD,KD,INS,SI1,SI2,SIZ,I0,IFF0, & IV2X,IV1X,NTSTAT,IFF,II,IM,JJ,JJT,BLKSYM,IV,K, & KK,IX,LV,NVIBM,LX,N,KNNMIN,GETLL,KDIAG,KL,NI,NJ, & NEVEN,NODD,NSIZE,IPARX,IXZT,ISYM INTEGER IXCOM(6),JXCOM(6),ISCOM(0:MXSPIN+1),JSCOM(0:MXSPIN+1), & LLCOM(MXSPIN+1),LSCOM(MXSPIN+1) INTEGER NBKPJ,KNMIN(0:3,0:98),KNMAX(0:98),IXZ INTEGER*2 IWT(MAXQ),BLKPTR(MAXQ+1),MOLDVS(MAXQ) INTEGER*4 IDTMP INTEGER NQNN,NQN,IFVIB,KLSYM,ESYMWT(0:98) REAL*8 AI(MXSPIN),ZFAC(NDPAR) INTEGER*2 NVSM,JJB(0:NX,MXSPIN),MINF(0:NX),IIVPT(0:98),IISPT(0:98) : ,LVQN(0:98) INTEGER KAPPA,NSPIN,NTSPIN,IISPN(NDSPIN),ITPTR,ITPTRM,IDIAG,NVIB COMMON/HSAVE/SQNN,SQJ,ZFAC,IPTYP,IPSYM,IP,NPARX,NJQ,IPDER COMMON /MOLDS/ NVSM,JJB,MINF,IIVPT,IISPT,LVQN COMMON /SQN/ KAPPA,NSPIN,NTSPIN,IISPN,ITPTR,ITPTRM,IDIAG,NVIB COMMON /MOLD/NBKPJ,KNMIN,KNMAX,IXZ,IWT,BLKPTR,MOLDVS COMMON /SQFMT/NQNN,NQN,IFVIB,KLSYM,ESYMWT SAVE /HSAVE/,/SQN/,/MOLD/,/MOLDS/,/SQFMT/ I0=0 NVIBM=NVIB-1 IF(NBKPF.GT.0) THEN DO 1 I=0,NVIBM II=IIVPT(I) C JJB(II,1) is the smallest N relative to F K=NBKPF-JJB(II,1)/2 IF(KNMAX(I).GT.K) KNMAX(I)=K 1 CONTINUE IFF0=NBKPF+NBKPF ELSE IFF0=200 ENDIF NPAR=IDPAR(0) CALL PASORT(NEGY,IDPAR,PAR,LVQN,KLSYM,IISPT,ESYMWT) C find interactions KDIAG=IDIAG IDIAG=-1 NTSTAT=NBKPJ DO 7 II=1,NTSTAT BLKPTR(II)=II-1 IM=MOLDVS(II) IFF=IFF0 CALL GETQS(IM,IFF,I0,I0,IXCOM,ISCOM,IV1X) IF(IAND(ISCOM(1),1).NE.0) THEN C special for mixed spin multiplicity (make N integer) IM=MOLDVS(II) IFF=IFF-1 CALL GETQS(IM,IFF,I0,I0,IXCOM,ISCOM,IV1X) ENDIF ISYM=IXCOM(2) NI=IXCOM(3) IF(KNMIN(3-ISYM,IV1X).GT.KNMAX(IV1X).AND. & KNMIN(ISYM,IV1X).GT.KNMAX(IV1X)) IWT(II)=0 IF(IWT(II).EQ.0) THEN BLKPTR(II)=NTSTAT GO TO 7 ENDIF DO 9 JJ=1,II IF(IWT(II).NE.IWT(JJ) ) GO TO 9 JJT=MOLDVS(JJ) CALL GETQS(JJT,IFF,I0,I0,JXCOM,JSCOM,IV2X) NJ=JXCOM(3) IDTMP=IV1X+IV2X+99*MIN(IV1X,IV2X) IVSYM=BLKSYM(IXCOM,JXCOM)+4*IDTMP IXZT=IXZ DO 10 K=1,NSPIN C check for spin multiplicity IF(IAND(ISCOM(K)+JSCOM(K),1).NE.0) GO TO 9 C check for connection restrictions IF(IAND(IXZT,1).NE.0) THEN IF( ISCOM(K).NE.JSCOM(K) ) GO TO 9 ENDIF IXZT=IXZT/2 10 CONTINUE IPARX=0 SIZ=0 C loop over parameters until connection found 11 IPARX=IPARX+1 IF(IPARX.LE.NPARX) THEN IF(IPSYM(IPARX).NE.IVSYM) GO TO 11 IDTMP=ABS(IPTYP(IPARX)) CALL IDPARS(IDTMP,IKQ,ITY,L,LD,KD,INS,SI1,SI2,SIZ,KL) IF(IAND(ITY,15).NE.0) GO TO 11 IF(SIZ.GT.1) CALL SZNZ(SIZ,NI,NJ,IXCOM,JXCOM, & ISCOM,JSCOM) K=GETLL(LD,KD,L,SI1,SI2,LLCOM,LSCOM,ISCOM,JSCOM) IF(SIZ.LT.0) CALL SZNZ(SIZ,NI,NJ,IXCOM,JXCOM, & ISCOM,JSCOM) IF(K.NE.0) GO TO 11 IF(KD.GT.0) IDIAG=KDIAG C II and JJ are coupled IF(BLKPTR(JJ).EQ.BLKPTR(II) ) THEN IF(IDIAG.LT.0) GO TO 11 ELSE IF(BLKPTR(II).GT.BLKPTR(JJ)) THEN C rename ptr II to ptr JJ KK=BLKPTR(II) IX=JJ ELSE C rename ptr JJ to ptr II KK=BLKPTR(JJ) IX=II ENDIF IDIAG=KDIAG DO 15 K=1,II IF(BLKPTR(K).EQ.KK) BLKPTR(K)=BLKPTR(IX) 15 CONTINUE ENDIF ENDIF 9 CONTINUE 7 CONTINUE C block connections now found IF(IDIAG.EQ.-1) WRITE(LU,*) 'NO DIAGONALIZATION NEEDED' IF(IDIAG.EQ. 0) WRITE(LU,*) 'ENERGY SORT OF WANG SUB-BLOCKS' IF(IDIAG.EQ. 1) WRITE(LU,*) 'EIGENVECTOR SORT OF STATES' IF(IDIAG.EQ. 2) WRITE(LU,*) 'EIGENVECTOR SORT OF K DOUBLETS' IF(ITPTR.GT. 0) WRITE(LU,*) & 'I(TOT) FOR LAST 2 SPINS IS LAST QUANTUM BEFORE F' IF(KAPPA.LT.0) THEN WRITE(LU,*) 'PROLATE ROTOR' ELSE WRITE(LU,*) 'OBLATE ROTOR' ENDIF WRITE(LU,*) ' V KMIN KMAX ESYMWT SPINS ' DO 18 IV=0,NVIBM II=IISPT(IV) DO 19 I=1,NSPIN 19 AI(I)=0.5*IISPN(I+II) KNNMIN=MIN(KNMIN(0,IV),KNMIN(1,IV),KNMIN(2,IV),KNMIN(3,IV)) WRITE(LU,100) IV,KNNMIN,KNMAX(IV),ESYMWT(IV),(AI(I),I=1,NSPIN) ESYMWT(IV)=MOD(ESYMWT(IV),100) 18 CONTINUE 100 FORMAT(4I5,9F5.1) C bubble sort block number LX=NTSTAT C DO WHILE( LX.GE.2) 13 IF( LX.GE.2 ) THEN LV=0 K=1 DO 20 I=2,LX IF(BLKPTR(K).GT.BLKPTR(I) ) THEN LV=K JJ =BLKPTR(I) BLKPTR(I)=BLKPTR(K) BLKPTR(K)=JJ JJ =MOLDVS(I) MOLDVS(I)=MOLDVS(K) MOLDVS(K)=JJ JJ =IWT(I) IWT(I)=IWT(K) IWT(K)=JJ ENDIF 20 K=I LX=LV GO TO 13 ENDIF C END DO WRITE(LU,101) 101 FORMAT(' BLOCK - WT - SYM - V - TSP - N - other quanta ', & '(rel. to F=0 )' ) C convert block label to pointer I=-1 N=0 NSIZE=0 NEVEN=0 NODD=0 DO 25 JJ=1,NTSTAT IF(BLKPTR(JJ).NE.I) THEN N=N+1 I=BLKPTR(JJ) BLKPTR(N)=JJ NSIZE=MAX(NSIZE,NEVEN,NODD) NEVEN=0 NODD=0 ENDIF JJT=MOLDVS(JJ) CALL GETQS(JJT,I0,I0,I0,IXCOM,ISCOM,IV) DO 30 II=1,NSPIN 30 AI(II)=0.5*ISCOM(II) ISYM=IXCOM(2) IV=IXCOM(6) IF(IWT(JJ).NE.0) THEN NEVEN=NEVEN+(KNMAX(IV)-KNMIN(ISYM,IV))/KINC+1 NODD =NODD +(KNMAX(IV)-KNMIN(3-ISYM,IV))/KINC+1 ENDIF WRITE(LU,102) N,IWT(JJ),ISYM,IV,JJT,(AI(II),II=1,NSPIN) 25 CONTINUE NSIZE=MAX(NSIZE,NEVEN,NODD) WRITE(LU,*) 'Maximum Dimension for Hamiltonian = ',NSIZE & ,', Limit in SPINV = ',NEGY BLKPTR(N+1)=NTSTAT+1 IF(I.EQ.NTSTAT) N=N-1 NBKPF=N NBKPJ=NBKPF IF(NSIZE.LT.NEGY) NEGY=NSIZE SETBLK=100 IF(NVIB.GT.9) SETBLK=10000 RETURN 102 FORMAT(1X,5I5,9F6.1) END SUBROUTINE PASORT(NSIZE,IDPAR,PAR,LVQN,KLSYM,IISPT,ESYMWT) C formal parameters: INTEGER NSIZE,KLSYM,ESYMWT(0:*) REAL*8 IDPAR(0:*),PAR(*) INTEGER*2 LVQN(0:*),IISPT(0:*) C C .. PACKAGE FOR 9 INTERACTING VIBRATIONAL STATES WITH MULTI-SPIN C C initialize spin factors and arrange operators C C on entry: C NSIZE= block size and dimension C IDPAR= list of parameter identifiers ( element 0 is length ) C INTEGER NDPAR,NDMX,MXSPIN,NDSPIN PARAMETER (NDPAR=750,NDMX=450,MXSPIN=5,NDSPIN=MXSPIN*9) INTEGER NOLD,I,IPAR,ITY,ITYI,IKQ,L,LD,KD,KSY,INS,IV12,KL, & SI1,SI2,NPARQ,ISWAP,K,KK,ITYP,IKQP,LP,LDP,ISYM,J, & KDP,INSP,SI1P,SI2P,IV12Q,NJQT,NNVIB,IFAC,IV1,IV2, & LV1,LV2,IV12D,SIZ,SIZP,KLP,IDMY(5),II,IJ INTEGER*4 IDTMP,IPTYP(NDPAR),IDVAL REAL*8 SQNN,SQJ(10),SPFAC2(0:NDSPIN),ZVAL,Z, & ZFAC(NDPAR),SPFAC(0:NDSPIN),VFAC,VDIV,ZV(2) INTEGER*2 NPARX,NJQ(NDPAR),IP(NDPAR),IPSYM(NDPAR), & IPDER(0:NDPAR) COMMON/HSAVE/SQNN,SQJ,ZFAC,IPTYP,IPSYM,IP,NPARX,NJQ,IPDER INTEGER KAPPA,NSPIN,NTSPIN,IISPN(NDSPIN),ITPTR,ITPTRM,IDIAG,NVIB COMMON /SQN/ KAPPA,NSPIN,NTSPIN,IISPN,ITPTR,ITPTRM,IDIAG,NVIB COMMON /SPFACV/ SPFAC,NOLD SAVE /HSAVE/,/SPFACV/,/SQN/ C KLSYM=0 C find reduced spin matrix elements DO 20 I=1,NTSPIN SPFAC(I)=0. SPFAC2(I)=0. IF(IISPN(I).GT.0) THEN Z=IISPN(I)+1 SPFAC(I)=0.5*SQRT(IISPN(I)*(Z+1)*Z) ENDIF IF(IISPN(I).GE.2) SPFAC2(I)=0.5*SQRT((Z+2)/(Z-2))/IISPN(I) 20 CONTINUE IF(NSIZE.GT.NDMX) NSIZE=NDMX C initialize SPECFC and DIRCOS NOLD=-1 CALL SPECFC(NOLD,I,I,I,I,I,I,ZV,IP,IP) I=0 CALL DIRCOS(IDMY,IDMY,I,I,I,ZFAC,IP,IP,I,I) SQJ(1)=1 C .. code parameter id for power of N*(N+1) and equivalence C .. NJQ = POWER + 1 (first occurrance of power) C .. IPTYP= ITY + 1 (first occurrance of cosine) C .. ITY = sequence number for hybrid operators C .. negate NJQ and IPTYP for successive occurrance C .. IP points to parameter IPAR=IDPAR(0) NNVIB=NVIB*100 IFAC=10 IF(NVIB.GT.9) IFAC=100 VFAC=IFAC*IFAC VDIV=1./VFAC KK=0 DO 50 I=1,IPAR IF(IDPAR(I).GE.0) THEN KK=KK+1 IPDER(I)=KK ITYI=I ELSE IPDER(I)=-ITYI ENDIF 50 CONTINUE IPDER(0)=KK NPARX=0 ITYI=0 IV12D=0 C while(IPAR.GT.0) 112 IF(IPAR.GT.0) THEN ITYI=0 C repeat until parameter subtypes exhausted (break= go to 110) 111 Z=ABS(IDPAR(IPAR))+0.5 IDVAL=VDIV*Z IV12Q=Z-IDVAL*VFAC IV2=IV12Q/IFAC IV1=IV12Q-IV2*IFAC IF(IV2.GT.IV1) THEN I=IV1 IV1=IV2 IV2=I ENDIF IV12Q=IV1+100*IV2 IF(IV12Q.EQ.909.AND.NVIB.LE.9) IV12Q=9999 IF(IV12Q.EQ.9999) THEN IV12=IV12D IV1=IV12/100 IV2=IV1 ELSE IF(IV1.GE.NVIB) GO TO 110 IV12=IV12Q ENDIF IF(IDVAL.EQ.9100) THEN I=0 ZV(1)=PAR(IPAR) CALL SPECFC(I,IV1,IV2,I,I,I,I,ZV,IP,IP) GO TO 110 ENDIF CALL IDPARI(IDVAL,ITYI,NJQT,ISYM,ZVAL) IF(ITYI.EQ.0) GO TO 110 CALL IDPARS(IDVAL,IKQ,ITY,L,LD,KD,INS,SI1,SI2,SIZ,KL) C setup spin factors I=IISPT(IV1) J=IISPT(IV2) IF(SI1.GT.0) THEN II=SI1+I IJ=SI1+J IF(IISPN(II).EQ.IISPN(IJ)) THEN IF(IISPN(II).EQ.0) GO TO 110 ZVAL=ZVAL*SPFAC(II) IF(SI2.EQ.SI1) THEN IF(IISPN(II).LT.2) GO TO 110 ZVAL=ZVAL*SPFAC2(II) ELSE IF(SI2.NE.0) THEN II=SI2+I IJ=SI2+J IF(IISPN(II).EQ.IISPN(IJ)) THEN IF(IISPN(II).EQ.0) GO TO 110 ZVAL=ZVAL*SPFAC(II) ELSE IF(IAND(IISPN(II)-IISPN(IJ),1).NE.0) THEN GO TO 110 ENDIF ENDIF ELSE IF(IAND(IISPN(II)-IISPN(IJ),1).NE.0) THEN GO TO 110 ENDIF ENDIF IF(SIZ.NE.0) THEN C .. setup for SzNz IF(IISPN(I+1).EQ.IISPN(J+1)) THEN IF(IISPN(I+1).EQ.0) GO TO 110 ZVAL=ZVAL*SPFAC(I+1) ELSE IF(IAND(IISPN(I+1)-IISPN(J+1),1).NE.0) THEN GO TO 110 ENDIF ZVAL=ZVAL*0.5 ENDIF LV1=LVQN(IV1) LV2=LVQN(IV2) C .. fixup ESYMWT if A state has wrong K dependence IF(IV1.EQ.IV2.AND.ESYMWT(IV1).GE.0) THEN IF(LV1.EQ.0.AND.MOD(KD,3).NE.0) ESYMWT(IV1)=-1 ENDIF C .. setup for l-doubling IF(LV1.LT.0) THEN LV1=-LV1 IV12=IV12-1 ENDIF IF(LV2.LT.0) THEN LV2=-LV2 IV12=IV12-100 ENDIF IF(LV1.NE.LV2) THEN C .. add flag for KL symmetry KL=1 IDVAL=IDVAL+5 KLSYM=1 ENDIF IF(LV1.EQ.2.AND.LV2.NE.1) GO TO 110 IF(LV2.EQ.2.AND.LV1.NE.1) GO TO 110 IF(LV1.NE.0.AND.LV1.EQ.LV2) THEN C .. change symmetry for Lz operator I=MAX(KD,LD) IF(ITY.GT.160) I=I+1 IF(IAND(I,1).NE.0) ISYM=3-ISYM ENDIF NPARQ=NPARX NPARX=NPARX+1 IF(NPARX.GT.NDPAR) THEN WRITE(*,*) ' HAMX: TOO MANY PARAMETERS ! ',NPARX,NDPAR STOP ENDIF IF(IV12.GT.9898) GO TO 110 IDTMP=IV12 KSY=ISYM+4*IDTMP ISWAP=0 IF(IAND(ITY,1).EQ.0) THEN ISWAP=NPARX DO 130 KK=NPARQ,1,-1 IF(IPSYM(KK).NE.KSY) GO TO 130 IDTMP=ABS(IPTYP(KK)) CALL IDPARS(IDTMP,IKQP,ITYP,LP,LDP,KDP,INSP,SI1P, & SI2P,SIZP,KLP) C .. check for same K dependence IF(IKQ.NE.IKQP) GO TO 130 IF(LD .NE. LDP) GO TO 130 IF(KD .NE. KDP) GO TO 130 IF(ITY.NE.ITYP) GO TO 130 IF(KL .NE. KLP) GO TO 130 IF(SIZ.NE.SIZP) GO TO 130 IF(IAND(ITYP,15).NE.0) GO TO 130 C .. check for same operator except power of N*(N+1) IF(IDVAL.EQ.IDTMP) THEN ISWAP=KK NJQT =-NJQT GO TO 140 ELSE IF(ISWAP.EQ.NPARX) THEN ISWAP=KK ENDIF 130 CONTINUE 140 CONTINUE ENDIF IF(ISWAP.LT.NPARX) THEN IF(ISWAP.GT.0) IDVAL=-IDVAL ISWAP=ISWAP+1 DO 145 K=NPARQ,ISWAP,-1 IP(K+1) =IP(K) IPTYP(K+1)=IPTYP(K) IPSYM(K+1)=IPSYM(K) ZFAC(K+1) =ZFAC(K) 145 NJQ(K+1) =NJQ(K) ENDIF IP(ISWAP) =IPAR IPTYP(ISWAP)=IDVAL IPSYM(ISWAP)=KSY ZFAC(ISWAP) =ZVAL NJQ(ISWAP) =NJQT IF(ITYI.GT.1) GO TO 111 110 IF(IV12Q.EQ.9999) THEN IV12D=IV12D+101 IF(IV12D.GT.NNVIB) THEN IV12D=0 IPAR=IPAR-1 ENDIF ELSE IV12D=0 IPAR=IPAR-1 ENDIF GO TO 112 ENDIF c write(*,'(5i10)') (ipar,ip(ipar),iptyp(ipar),ipsym(ipar), c & njq(ipar),ipar=1,nparx) RETURN END SUBROUTINE IDPARI(IDVAL,ITP,NSQ,ISYM,ZVAL) C formal parameters: INTEGER*4 IDVAL INTEGER ITP,NSQ,ISYM REAL*8 ZVAL C subroutine to parse parameter identifier for complex parameter type C on entry: C IDVAL= parameter ID / vibrational field C ITP= value of sub-parameter type for hybrid operator C on return: C IDVAL= modified IDVAL + ITP - NSQ C ITP= new value C NSQ= power of N*(N+1) +1 C ISYM= symmetry C ZVAL= factor for operator INTEGER*4 L10K,ND,NDP,MXSPIN,NDSPIN PARAMETER (L10K=10000,ND=9,NDP=ND+1,MXSPIN=5,NDSPIN=MXSPIN*9) INTEGER ITPNXT(ND),SI1,SI2,ITY,ITMP,IDHIGH,IDLOW,ISIN,INS,I,ITPIN, & ITPOP(ND),IPHAZ REAL*8 ZFAC(ND) INTEGER KAPPA,NSPIN,NTSPIN,IISPN(NDSPIN),ITPTR,ITPTRM,IDIAG,NVIB COMMON /SQN/ KAPPA,NSPIN,NTSPIN,IISPN,ITPTR,ITPTRM,IDIAG,NVIB SAVE /SQN/ C 1 2 3 4 5 6 7 8 9 DATA ITPNXT/0, 0, 2, 6, 6, 2, 0, 9, 0/ DATA ITPOP /0, 0, 100, 400, 400, 100, 10, 0, 10/ DATA ZFAC /1.D0,3.D0,1.5D0,-4.D0,4.D0,-3.D0,1.D0,2.D0,-2.D0/ IPHAZ=0 IF(ITP.GT.0) THEN C .. get current sub-parameter type from previous IF(ITP.LE.ND) THEN ITPIN=ITP ITP=ITPNXT(ITP) IF(ITP.EQ.0) RETURN ELSE IPHAZ=ITP/NDP ITPIN=ITP-IPHAZ*NDP IF(ITPIN.NE.0) ITP=ITPNXT(ITPIN) ENDIF ENDIF ZVAL=ZFAC(1) IDHIGH=IDVAL/L10K ISIN=-1 INS=0 IF(IDHIGH.EQ.0) THEN IDLOW=IDVAL SI2=0 SI1=0 ELSE IDLOW=IDVAL-L10K*IDHIGH IF(IDHIGH.GE.1000) THEN ISIN=IDHIGH/1000 IDHIGH=IDHIGH-1000*ISIN ENDIF SI2=IDHIGH/100 INS=IDHIGH-100*SI2 SI1=INS/10 INS=INS-SI1*10 IF(INS.GT.0.AND.NSPIN.EQ.0) RETURN IF(SI2.GT.SI1) THEN C make sure SI1 > SI2 I=SI1 SI1=SI2 SI2=I IDVAL=IDVAL+(SI2-SI1)*900000 ENDIF IF(SI1.GT.NSPIN) RETURN ENDIF ITY=IDLOW/100 NSQ=MOD(IDLOW,10) IDVAL=IDVAL-NSQ NSQ=NSQ+1 C exclude S * S and I * S from aa, bb, cc operators IF(ITP.EQ.2) THEN IF(SI2.NE.0) THEN ITP=ITPNXT(ITP) IF(ITP.EQ.0.AND.IPHAZ.EQ.0) RETURN ENDIF ENDIF IF(ITP.NE.0) THEN IF(ITP.EQ.6) THEN IF(SI1.EQ.0) THEN ITP=8 NSQ=NSQ+1 ENDIF ENDIF ELSE IF(IPHAZ.NE.0) THEN IPHAZ=IPHAZ-1 ELSE IF(INS.GE.5) THEN IPHAZ=4 ENDIF IF(ITY.EQ.0) THEN ITP=1 ELSE IF(ITY.LE.3) THEN ITP=ITY+2 IF(KAPPA.GT.0) ITP=8-ITP C ITP=3,4,5 for A,B,C IF(ITP.EQ.3) THEN IF(SI1.EQ.0) ITP=7 ENDIF ELSE ITP=1 ITMP=ITY IF(ITY.GE.92) ITMP=(ITY-92)/2 C negate odd powers of P+**2 IF(ITMP.LE.19) THEN IF(IAND(ITMP,1).EQ.0) ZVAL=-ZVAL ENDIF ENDIF ENDIF ISYM=ITY/20 IF(ISYM.GT.3) THEN ISYM=0 ELSE IF(ISYM.NE.0) THEN IF(KAPPA.GT.0) THEN IDVAL=IDVAL+4000*(2-ISYM) ELSE ISYM=4-ISYM ENDIF ENDIF IF(ISIN.GT.10) ISYM=3-ISYM IF(ITP.GT.1) THEN IDVAL=IDVAL-ITY*100+ITPOP(ITP) ZVAL=ZVAL/ZFAC(ITP) ENDIF IF(INS.GE.5) THEN IDVAL=IDVAL+IPHAZ ITP=ITP+IPHAZ*NDP ENDIF IF(IDVAL.EQ.0) IDVAL=1 RETURN END SUBROUTINE SETWT(IVIB,IWT,IAX,IWTPL,IWTMN,VSYM) C formal parameters: INTEGER IVIB,IAX,IWTPL,IWTMN INTEGER*2 IWT(*) INTEGER*4 VSYM C C set weights C on entry: C NVIB= number of vibrational states * number of wang states C IAX= axis of symmetry C IWTPL= weight for even rotational states C IWTMN= weight for odd rotational states C VSYM= symmetry code for vibrational states C NCOPY= number of spin states C on return: C IWT= weights for vibration rotation states INTEGER*4 ITMP,LVSYM INTEGER I,ISYM,II,JJ,IBGN,IEND,IAXX,IX LVSYM=0 IF(IVIB.LT.0) THEN IBGN=1 IEND=1-4*IVIB IF(VSYM.GT.0) LVSYM=VSYM ELSE IBGN=IVIB*4+1 IEND=IBGN ENDIF IX=(IAX-1)/3 IAXX=IAX-IX*3 DO 5 I=IBGN,IEND,4 ITMP=LVSYM LVSYM=LVSYM/10 ISYM=ITMP-10*LVSYM IF (IAND(ISYM,1).EQ.0) THEN II=IWTPL JJ=IWTMN ELSE II=IWTMN JJ=IWTPL ENDIF IWT(I)=II IF (IAXX.EQ.3) THEN C .. WEIGHTS FOR X IWT(I+1)=II IWT(I+2)=JJ IWT(I+3)=JJ ELSE IF (IAXX.EQ.2) THEN C .. WEIGHTS FOR B IWT(I+1)=JJ IWT(I+2)=II IWT(I+3)=JJ ELSE C .. WEIGHTS FOR Z IWT(I+1)=JJ IWT(I+2)=JJ IWT(I+3)=II ENDIF IF (IX.EQ.1) THEN IWT(I+1)=0 IWT(I+2)=0 ELSE IF (IX.EQ.2) THEN IWT(I)=0 IWT(I+3)=0 ENDIF 5 CONTINUE RETURN END SUBROUTINE SETSP(ISPNX,IVIB,NT,NSSV) C formal parameters: INTEGER*4 ISPNX INTEGER IVIB,NT,NSSV(0:*) C subroutine to set up spin structure C on entry: C ISPNX = spin code C IVIB = vibrational quantum number C on return: C NT = total number of spin combinations C NSSV = number of spin combinations in high F limit C /MOLDS/ initialized with spins*2 - F*2 INTEGER NX,IC1,IC2,IBGN,IEND,ITOT,MXSPIN,NDSPIN,I10,IORG, : NSPINL,MINFT,MINF0 PARAMETER (NX=170,MXSPIN=5,NDSPIN=MXSPIN*9) INTEGER*2 NVSM,JJB(0:NX,MXSPIN),MINF(0:NX),IIVPT(0:98),IISPT(0:98) : ,LVQN(0:98) INTEGER ISPIN,I,II,IX,NN,IVAL,IV,NV,NVZ,IBASE,JJT(MXSPIN),JJ, : ITMP,NSSTAT,ITOTPT,IFHALF INTEGER KAPPA,NSPIN,NTSPIN,IISPN(NDSPIN),ITPTR,ITPTRM,IDIAG,NVIB COMMON /SQN/ KAPPA,NSPIN,NTSPIN,IISPN,ITPTR,ITPTRM,IDIAG,NVIB COMMON /MOLDS/ NVSM,JJB,MINF,IIVPT,IISPT,LVQN SAVE /SQN/,/MOLDS/ IV=IVIB IF(IVIB.LT.0) THEN C initialize NN=-1-IVIB DO 1 I=0,NN IIVPT(I)=0 IISPT(I)=0 1 CONTINUE IV=0 NSPIN=0 NTSPIN=0 ELSE IF(NTSPIN.GE.NDSPIN) THEN WRITE(*,*) 'Too many kinds of spins' STOP ENDIF I10=10 IF(MOD(ISPNX,10).EQ.0) THEN ISPNX=ISPNX/I10 I10=100 ENDIF IF(ISPNX.EQ.1) ISPNX=0 IORG=NTSPIN NSSTAT=1 NSPINL=0 ISPIN=0 II=NT C decode spins DO 5 I=1,MXSPIN JJB(II,I)=0 JJT(I)=0 IX=I+IORG IF(ISPNX.EQ.0) THEN IISPN(IX)=0 ELSE NSPINL=I IISPN(IX)=MOD(ISPNX,I10)-1 ISPNX=ISPNX/I10 ISPIN=ISPIN+IISPN(IX) NSSTAT=NSSTAT*(IISPN(IX)+1) ENDIF 5 CONTINUE IF(NSPIN.LT.NSPINL) THEN NSPIN=NSPINL ELSE C check previous combinations for match DO 6 IBASE=1,IORG,MXSPIN DO 7 I=1,NSPIN IF(IISPN(I+IBASE-1).NE.IISPN(I+IORG)) GO TO 6 7 CONTINUE IISPT(IV)=IBASE-1 C match obtained, find a matching V and return DO 8 I=0,98 IF(IISPT(IV).EQ.IISPT(I)) THEN IIVPT(IV)=IIVPT(I) NSSV(IV)=NSSV(I) RETURN ENDIF 8 CONTINUE 6 CONTINUE ENDIF NTSPIN=NTSPIN+MXSPIN NSSV(IV)=NSSTAT IIVPT(IV)=NT IISPT(IV)=IORG NT=NT+NSSTAT IFHALF=IAND(ISPIN,1) MINF(II)=IFHALF IF(NSPINL.EQ.0) RETURN IF(ITPTR.LE.0) THEN ITOTPT=NSPINL+IORG IBGN=IISPN(ITOTPT) IEND=IBGN NVZ=NSSTAT/(IBGN+1) ELSE C .. case for ITOT coupling IF(IVIB.LT.0) ITPTR=NSPIN IF(NSPINL.LT.ITPTR) NSPINL=ITPTR ITOTPT=ITPTR+IORG IC1=IISPN(ITOTPT) IC2=IISPN(ITOTPT-1) IISPN(ITOTPT-1)=0 IBGN=ABS(IC1-IC2) IEND=IC1+IC2 NVZ=NSSTAT/((IC1+1)*(IC2+1)) ENDIF C select N quanta DO 10 NN=-ISPIN,ISPIN,2 MINF0=IFHALF IF(NN.LT.0) MINF0=-NN C loop through ITOT DO 15 ITOT=IBGN,IEND,2 IISPN(ITOTPT)=ITOT NV=NVZ*(ITOT+1)-1 C try all possible spins DO 20 IV=0,NV JJ=NN IVAL=IV MINFT=MINF0 DO 25 I=1,NSPINL IX=I+IORG IBASE=IISPN(IX)+1 JJT(I)=JJ ITMP=IVAL IVAL=IVAL/IBASE ITMP=ITMP-IVAL*IBASE JJ=JJ+ITMP-IISPN(IX) C make (J+N) GE S IF(MINFT+JJ.LT.0) MINFT=-JJ JJ=JJ+ITMP 25 CONTINUE IF(JJ.EQ.0) THEN C this possibility gives F=0 IF(ITPTR.GT.0) JJT(ITPTR)=ITOT IF(II.LE.NX) THEN IF(IAND(MINFT,1).NE.IFHALF) MINFT=MINFT+1 MINF(II)=MINFT DO 30 I=1,MXSPIN 30 JJB(II,I)=JJT(I) II=II+1 ENDIF ENDIF 20 CONTINUE 15 CONTINUE 10 CONTINUE IF(II.NE.NT) THEN WRITE(*,*) & ' warning spin states truncated in SETSP ',II, NT IV=MAX(IVIB,0) NSSV(IV)=NSSV(IV)-(NT-II) NT=II ENDIF IF(ITPTR.GT.0) THEN ITPTRM=ITPTR-1 IISPN(ITOTPT)=IC1 IISPN(ITOTPT-1)=IC2 ENDIF RETURN END FUNCTION GETQQ(MXBLK,IBLK,F,IWTB,SBKPTR,KMIN,VS) C formal parameters: INTEGER GETQQ,IBLK,F,IWTB,MXBLK INTEGER*2 SBKPTR(*),KMIN(*),VS(*) C subroutine to get sub-block structure and major quantum numbers C on entry: C MXBLK=biggest block size C IBLK= block number C on return: C GETQQ= number of sub-blocks C F= F quantum number *2 C IWTB= statistical weight C SBKPTR= pointer to beginning index for sub-block C KMIN= minimum K for each sub-block C VS= vibration and symmetry for each sub-block INTEGER MAXQ,IBGN,IEND,MXSPIN,NDSPIN,NX PARAMETER (MAXQ=500,MXSPIN=5,NDSPIN=MXSPIN*9,NX=170) INTEGER ISBLK,N,I,II,MVS,MSS,NN,ISYM,KK,IV INTEGER NBKPJ,KNMIN(0:3,0:98),KNMAX(0:98),IXZ INTEGER*2 IWT(MAXQ),BLKPTR(MAXQ+1),MOLDVS(MAXQ) INTEGER*2 NVSM,JJB(0:NX,MXSPIN),MINF(0:NX),IIVPT(0:98),IISPT(0:98) : ,LVQN(0:98) COMMON /MOLDS/ NVSM,JJB,MINF,IIVPT,IISPT,LVQN COMMON /MOLD/ NBKPJ,KNMIN,KNMAX,IXZ,IWT,BLKPTR,MOLDVS SAVE /MOLD/,/MOLDS/ F=(IBLK-1)/NBKPJ ISBLK=IBLK-NBKPJ*F C get information for sub-blocks N=1 I=1 IBGN=BLKPTR(ISBLK) MSS=MOLDVS(IBGN)/NVSM F=F+F-IAND(MINF(MSS),1) IEND=BLKPTR(ISBLK+1)-1 DO 10 II=IBGN,IEND MVS=MOLDVS(II) MSS=MVS/NVSM IF(F.LT.MINF(MSS)) GO TO 10 IV=(MVS-MSS*NVSM)/4 ISYM=IAND(MVS,3) NN=(JJB(MSS,1)+F)/2 IF(IAND(NN,1).NE.0) ISYM=3-ISYM KK=MIN(KNMAX(IV),NN)-KNMIN(ISYM,IV) IF(KK.GE.0) THEN IF(I.GT.MXBLK) THEN WRITE(*,*) ' GETQQ: too many sub-blocks',MXBLK STOP ENDIF SBKPTR(I)=N KMIN(I)=KNMIN(ISYM,IV) VS(I)=MVS N=N+1+KK/2 I=I+1 ENDIF 10 CONTINUE IWTB=IWT(IBGN) SBKPTR(I)=N GETQQ=I-1 RETURN END SUBROUTINE GETQN (IBLK,INDX,IQN,IDGN,NCOD) C formal parameters: INTEGER IBLK,INDX,IQN(*),IDGN,NCOD C subroutine to get quantum numbers C on entry: C IBLK= block number C INDX= index within block C on return: C IQN= list of quantum numbers C IDGN= degeneracy C NCOD= seraching aid to find state within a Wang block INTEGER MAXSV,MAXS,MAXS2,MXSPIN,NDSPIN,I0,GETQQ,KINC PARAMETER (MAXS=64,MAXS2=MAXS+2,MAXSV=MAXS+MAXS+1,KINC=2) PARAMETER (MXSPIN=5,NDSPIN=MXSPIN*9,I0=0) INTEGER NSBL,JWT,JSV,NQNN,NQN,IFVIB,ESYMWT(0:98),I,IBGN,IEND, & IV,K,IXCOM(6),ISCOM(0:MXSPIN+1),IBLKSV(2),IBGNSV(2),IQSP, & IQF,LBLK,OLDBLK,NSBLK,KLSYM,ISYM,LOWER INTEGER*2 KMIN(MAXSV),VS(MAXSV),SBKPT(MAXSV+1) INTEGER KAPPA,NSPIN,NTSPIN,IISPN(NDSPIN),ITPTR,ITPTRM,IDIAG,NVIB COMMON /SQN/ KAPPA,NSPIN,NTSPIN,IISPN,ITPTR,ITPTRM,IDIAG,NVIB COMMON /CGETQ/ NSBL(2),JWT(2),JSV(2),SBKPT,KMIN,VS COMMON /SQFMT/NQNN,NQN,IFVIB,KLSYM,ESYMWT SAVE IBLKSV,IBGNSV,LBLK,OLDBLK SAVE /SQN/,/CGETQ/,/SQFMT/ DATA IBLKSV,LBLK,OLDBLK,IBGNSV/0,0,1,2,1,MAXS2/ C C get sub-block info from store or GETQQ C C look at last used then oldest IF(IBLK.NE.IBLKSV(LBLK)) THEN IF(IBLK.EQ.IBLKSV(OLDBLK)) THEN LBLK=OLDBLK OLDBLK=3-OLDBLK ELSE C use GETQQ because IBLK not stored IBGN=IBGNSV(OLDBLK) NSBLK=GETQQ(MAXS,IBLK,I,JWT(OLDBLK),SBKPT(IBGN), & KMIN(IBGN),VS(IBGN)) IF(NSBLK.EQ.0) THEN IDGN=0 RETURN ELSE LBLK=OLDBLK OLDBLK=3-LBLK IBLKSV(LBLK)=IBLK JSV(LBLK)=I NSBL(LBLK)=IBGN+NSBLK-1 ENDIF ENDIF ENDIF C C check for request of size C IEND=NSBL(LBLK) IF(INDX.LE.0) THEN IDGN=SBKPT(IEND+1)-1 RETURN ENDIF C binary search for sub-block IBGN=IBGNSV(LBLK) I=IBGN C DO WHILE(IBGN.LT.IEND) 10 IF(IBGN.LT.IEND) THEN I=(IBGN+IEND+1)/2 IF(INDX.GE.SBKPT(I)) THEN IBGN=I ELSE I=I-1 IEND=I ENDIF GO TO 10 ENDIF C END DO C assemble quanta NCOD=SBKPT(I+1)-SBKPT(I) IQF=JSV(LBLK) IQSP=VS(I) CALL GETQS(IQSP,IQF,I0,I0,IXCOM,ISCOM,IV) IDGN=JWT(LBLK) IV=IXCOM(6) IQN(1)=IXCOM(3) ISYM=IXCOM(2) LOWER=IQN(1)+ISYM K=INDX-SBKPT(I) IQN(2)=KMIN(I)+K*KINC IF(ESYMWT(IV).GE.0) THEN IF(MOD(IQN(2)-IXCOM(5),3).NE.0) THEN C .. E symmetry IDGN=ESYMWT(IV) ISYM=4 ENDIF IF(IXCOM(5).NE.0) THEN C .. l-doubled state NCOD=1 IF(KLSYM.EQ.0) ISYM=4 IF(ISYM.EQ.4) LOWER=IXCOM(5) ENDIF ENDIF IF(NQNN.EQ.2) THEN IF(ISYM.EQ.1.OR.ISYM.EQ.3) IQN(2)=-IQN(2) NCOD=1 ELSE IF(KAPPA.GT.0) THEN IQN(3)=IQN(2) IQN(2)=IQN(1)-IQN(2)+IAND(LOWER,1) NCOD=-NCOD ELSE IQN(3)=IQN(1)-IQN(2)+IAND(LOWER,1) ENDIF IDGN=IDGN*(IQF+1) K=NQNN+1 IF(IFVIB.NE.0) THEN IQN(K)=IV K=K+1 ENDIF IF(NQN.GT.6) THEN IQN(K)=IQSP IQN(K+1)=(IQF+1)/2 ELSE DO 30 I=1,NSPIN IQN(K)=( ISCOM(I+1)+1 )/2 30 K=K+1 ENDIF IF(IDGN.LT.0) IDGN=0 RETURN END SUBROUTINE DIRCOS(XBRA,XKET,L,KD,NCASE,DIREL,IBRA,IKET,IFUP,LSYM) C formal parameters: INTEGER XBRA(6),XKET(6),L,KD,NCASE,IFUP,LSYM INTEGER*2 IBRA(*),IKET(*) REAL*8 DIREL(*) C C SUBROUTINE TO CALCULATE DIRECTION COSINE ELEMENTS C C XBRA,XKET INTEGER VECTORS CONTAINING: C 1: DIMENSION OF SUB-BLOCK C 2: SYMMETRY CODE FOR BLOCK (0=A,1=Bx,2=By,3=Bz) C 3: 'N' QUANTUM NUMBER C 4: BEGINNING K C 5: L VALUE C 6: VIBRATIONAL QUANTA C L IS TENSOR ORDER OF COSINE C KD IS TOTAL K QUANTIUM CHANGE C DIREL IS A VECTOR OF DIRECTION COSINE ELEMENTS C NCASE IS LENGTH OF DIREL (PRESET TO MAX LENGTH) C IBRA,IKET INDICES OF MATRIX ELEMENTS IN DIREL (START WITH ZERO) C IFUP =0 IF NO UPPER TRIANGLE TO BE CALC. C IFUP <0 IF ELECTRIC DIPOLE C ODD ORDER PARAMETERS ARE IMAGINARY AND EVEN ARE REAL IF NOT ELECTRIC DIPOLE C FOR ELECTRIC DIPOLE ODD ORDER IS REAL AND EVEN ARE IMAGINARY C LSYM NOT ZERO IF K-L SYMMETRY ENFORCED C C L=0 KD=0 DIREL=1. C L=1 KD=0 DIREL=PHI(Z) C L=1 KD=1 DIREL=PHI(X) OR PHI(Y) C L=2 KD=0 DIREL=0.5*(3* PHI(Z)**2 -1.) C L=2 KD=1 DIREL=PHI(Z)*PHI(X)+PHI(X)*PHI(Z) OR C =PHI(Z)*PHI(Y)+PHI(Y)*PHI(Z) C L=2 KD=2 DIREL=2*(PHI(X)**2 - PHI(Y)**2) OR C =PHI(X)*PHI(Y)+PHI(Y)*PHI(X) C C IF KD > L THEN (P-)**(KD-L) OPERATOR MULTIPLIED BY C L,L OPERATOR C WANG FUNCTIONS MULTIPLIED BY i**(SYMMETRY) TO MAKE ALL C ODD ANGULAR MOMENTA IMAG. AND ALL EVEN ONES REAL C BEST RESULTS ON SPEED IF N FOR BRA CHANGES SLOWEST C C ********************************************************************* INTEGER KINC,KINC2 PARAMETER (KINC=2,KINC2=KINC+KINC) INTEGER NCMAX,I,II,KDEL,KDEL2,KKDEL,IKDEL, & LL,LMAX,NSBRA,NSKET,ISBRA,ISKET,NNBRA,NNKET,KBRA0, & KKET0,ISGN,ITMP,ISREV,IL,IR,KBRA,KKET,N,NBGN,LBRA,LKET, & KKET2,KK,KBGN,K,ISUM,KS,IFSIN REAL*8 FAC(0:16),FF(359),ELE,C3JJ,ELEX,ELETMP LOGICAL MATSYM SAVE FAC,FF LBRA=-XBRA(5) LKET=-XKET(5) KBRA0=XBRA(4) KKET0=XKET(4) IF(NCASE.LE.0) THEN IF(NCASE.EQ.0) THEN C .. initialize FAC(0)=1. DO 10 LL=1,16 ELE=LL-0.5 10 FAC(LL)=FAC(LL-1)*SQRT(LL/ELE) RETURN ENDIF C special for adding additional odd K operators NCMAX=-NCASE IFSIN=1 ELSE NCMAX=NCASE IF (L.LE.0.AND.KD.EQ.0) THEN C ... SCALAR COSINE IF(LBRA.NE.LKET) THEN NCASE=0 RETURN ENDIF II=0 IR=0 KDEL=KBRA0-KKET0 IF(KDEL.NE.0) THEN KDEL=KDEL/KINC IF(KDEL.GT.0) THEN IR= KDEL ELSE II=-KDEL ENDIF ENDIF NCASE=MIN(XBRA(1)-II,XKET(1)-IR) IF(NCASE.GT.NCMAX) GO TO 900 DO 100 I=1,NCASE DIREL(I)=1 IBRA(I)=II IKET(I)=II+KDEL 100 II=II+1 RETURN ENDIF IFSIN=0 ENDIF C NCASE=0 NNBRA=XBRA(3) NNKET=XKET(3) IF(NNBRA+NNKET.LT.L) RETURN NSBRA=XBRA(1) KDEL=KD IKDEL=KDEL-L IF(IKDEL.LE.0) THEN LMAX=L KKDEL=KDEL ELSE IF (NNBRA+NNKET.LT.KDEL) RETURN K=KBRA0+IKDEL+KINC*(NSBRA-1) IF(K.GT.NNBRA) K=NNBRA CALL FFCAL(NNBRA,K,FF) LMAX=KDEL KKDEL=L ENDIF IKDEL=IKDEL-1 KDEL2=KKDEL+KKDEL LL=L+L ISBRA=XBRA(2) NSKET=XKET(1) ISKET=XKET(2) NNBRA=NNBRA+NNBRA NNKET=NNKET+NNKET IF(L.NE.0) THEN ELE=NNBRA+1 IF(NNKET.NE.NNBRA) ELE=SQRT(ELE*(NNKET+1)) IF(KKDEL.NE.0) THEN ELE=ELE*FAC(L) IF(ISBRA.NE.ISKET) ELE=0.5*ELE ENDIF ISGN=NNBRA+KKET0+KKET0+KDEL2 ELSE ELE=1 IF(ISBRA.NE.ISKET) ELE=0.5 ISGN=0 ENDIF ISREV=L+ISBRA/2+ISKET/2+IFSIN C correct for imaginary operators LMAX=LMAX+IFSIN IR=LMAX+ISKET+ISBRA C .. correct to make dipoles real IL=4 IF(IFUP.LT.0) IL=5 C .. find Lz operators IF(LBRA+LKET.NE.0)THEN IF(IAND(LMAX+IL,1).NE.0.AND.LSYM.EQ.0) THEN ISREV=ISREV+1 IR=IR+1 IF(LKET.EQ.-2) ISGN=ISGN+2 ENDIF ENDIF C .. correct for imaginary parameters ITMP=ISGN+IL+ISKET-ISBRA+IAND(IR,1) ISGN=ITMP/2 IF(ISGN+ISGN.NE.ITMP) THEN IF(XBRA(6).LT.XKET(6)) ISGN=ISGN+1 ENDIF IF(IAND(ISGN+ISREV,1).NE.0) ELE=-ELE MATSYM=.FALSE. DO 115 K=1,5 IF(XBRA(K).NE.XKET(K)) GO TO 120 115 CONTINUE MATSYM=.TRUE. C START CALCULATION OF C KBRA=KBRA0+ (IL+ i)*KINC IL >= 0 C KKET=KKET0+ (IR+ i)*KINC IR >= 0, i=0,... C KBRA=KKET+KDEL 120 IL=(KDEL+KKET0-KBRA0)/KINC IF(IL.GE.0) THEN IR=0 KBRA=KDEL+KKET0 KKET=KKET0 ELSE IR=-IL IL=0 KBRA=KBRA0 KKET=KBRA0-KDEL ENDIF N=MIN(NSKET-IR,NSBRA-IL) IF(LSYM.NE.0) THEN IF(MOD(KDEL+LBRA-LKET,3).NE.0) N=0 ELSE IF(LBRA.NE.LKET) THEN N=0 ENDIF IF(N.GT.0) THEN NBGN=NCASE+1 NCASE=NCASE+N IF(NCASE.GT.NCMAX) GO TO 900 KKET2=KKET+KKET DO 130 N=NBGN,NCASE IBRA(N)=IL IKET(N)=IR IF(L.EQ.0) THEN DIREL(N)=ELE ELSE KK=-KDEL2-KKET2 DIREL(N)=ELE*C3JJ(NNBRA,LL,NNKET,KK,KDEL2,KKET2) KKET2=KKET2+KINC2 ENDIF IL=IL+1 130 IR=IR+1 IF (KDEL.EQ.0) RETURN C .. CORRECT FOR K=0 IF (KBRA.EQ.LBRA.OR.KKET.EQ.LKET) THEN DIREL(NBGN)=DIREL(NBGN)*FAC(1) ENDIF IF (IKDEL.GE.0) THEN C ADD ON P- OR P+ OPERATOR TO LEFT IF(KDEL.GT.0) THEN KBGN=KBRA-IKDEL ELSE KBGN=KBRA+1 ENDIF DO 150 N=NBGN,NCASE ELETMP=DIREL(N) DO 155 K=KBGN,KBGN+IKDEL 155 ELETMP=ELETMP*FF(K) DIREL(N)=ELETMP 150 KBGN=KBGN+KINC ENDIF ENDIF C C DO QUASI-DIAGONAL PART C K=1,...,KDEL-1 C KBRA=KBRA0+ (IL - i)*KINC IL>=n C KKET=KKET0+ (IR + i)*KINC IR>=0 , i=0 .. n C KDEL=KKET+KBRA C IF(KDEL.LE.0) RETURN KDEL2=-KDEL2 IF (IAND(ISREV,1).NE.0) ELE=-ELE ISUM=KDEL-KBRA0-KKET0 IF(ISUM.LT.0) GO TO 200 IF(LSYM.NE.0) THEN IF(MOD(KDEL+LBRA+LKET,3).NE.0) GO TO 200 ELSE IF(LBRA.NE.0) THEN IF(LBRA.EQ.LKET) GO TO 200 ENDIF ISUM=ISUM/KINC C .. ISUM = IL+IR IL=ISUM IF(KKET0.EQ.LKET) IL=IL-1 IF(IL.GE.NSBRA) IL=NSBRA-1 IR=ISUM-IL IF(MATSYM) THEN N=(IL-IR+2)/2 ELSE N=ISUM IF(KBRA0.NE.LBRA) N=N+1 IF(N.GT.NSKET) N=NSKET N=N-IR ENDIF IF(N.LE.0) GO TO 200 KKET=KKET0+IR*KINC KBRA=KDEL-KKET NBGN =NCASE+1 NCASE=NCASE+N IF(NCASE.GT.NCMAX) GO TO 900 IF(IAND(ISBRA+NNBRA,3).LT.2) THEN ELEX=ELE ELSE ELEX=-ELE ENDIF KKET2=KKET+KKET DO 165 N=NBGN,NCASE IBRA(N)=IL IKET(N)=IR IF(L.EQ.0) THEN ELETMP=ELEX ELSE KS=-KDEL2-KKET2 ELETMP=ELEX*C3JJ(NNBRA,LL,NNKET,KS,KDEL2,KKET2) KKET2=KKET2+KINC2 ENDIF DO 220 K=KBRA-IKDEL,KBRA IF(K.LE.0) THEN ELETMP=ELETMP*FF(1-K) ELSE ELETMP=ELETMP*FF(K) ENDIF 220 CONTINUE DIREL(N)=ELETMP IR=IR+1 IL=IL-1 165 KBRA=KBRA-KINC 200 CONTINUE C C START CALCULATION OF C IF(MATSYM) THEN IF(IFUP.EQ.0) RETURN N=NCASE IF(IBRA(N).EQ.IKET(N)) N=N-1 IF(NCASE+N.GT.NCMAX) GO TO 900 DO 170 I=1,N NCASE=NCASE+1 DIREL(NCASE)=DIREL(I) IBRA(NCASE)=IKET(I) 170 IKET(NCASE)=IBRA(I) RETURN ENDIF C KDEL=-KDEL KKDEL=-KKDEL GO TO 120 900 WRITE(*,*) 'DIRCOS WORKING VECTOR TOO SHORT' STOP END SUBROUTINE FFCAL (NFF,KFF,FF) INTEGER NFF,KFF,NLAST,KLAST,KBGN,K REAL*8 FF(*),SQ SAVE NLAST,KLAST,SQ DATA NLAST /-1/ IF(NFF.NE.NLAST) THEN NLAST=NFF SQ=NFF+1 SQ=SQ*NFF KBGN=1 ELSE IF(KLAST.GE.KFF) RETURN KBGN=KLAST+1 ENDIF KLAST=KFF DO 20 K=KBGN,KLAST FF(K)=SQRT(SQ) SQ=SQ-(K+K) 20 CONTINUE RETURN END