C****************************************************************** C* KOBEPACK/1 VERSION 1.0 * C* MARCH 1992 * C* COPYRIGHT (C) M. KABURAGI, T. NISHINO AND T. TONEGAWA * C* WITH THE HELP OF TITPACK VERSION 2 * C* FEBRUARY 1991 * C* COPYRIGHT (C) HIDETOSHI NISHIMORI * C****************************************************************** C***** SMALL-SIZE ***** C****************************************************************** C C*************** CHECK OF THE EIGENVECTOR AND EIGENVALUE ********* C SUBROUTINE CHECKS(ID,IS,ELMNT,IDIM,NDCLRS,VEC,NDCLRV,MDCLRV, &HEXPC) C C ELMNT @ NONZERO ELEMENTS C IDIM @ MATRIX DIMENSION C NDCLRS @ DECLARED ARRAY SIZE IN THE MAIN PROGRAM C VEC(*,IS) @ EIGENVECTOR TO BE CHECKED C VEC(*,ID) # H*VEC(ISRC)/VEC(ISRC) C HEXPC # C VEC # EIGEN VECTOR AND WORK SPACE C NDCLRV @ SIZE OF VEC C & IMPLICIT REAL*8(A-H,O-Z) DIMENSION ELMNT(0:NDCLRS,0:IDIM-1) DIMENSION VEC(0:NDCLRV,0:MDCLRV) C CALL INPRO(IDIM,VEC(0,IS),VEC(0,IS),PROD) PROD=SQRT(PROD) IF(PROD.LT.1.D-30)THEN PRINT *,' #(W01)# NULL VECTOR GIVEN TO CHECKS.' RETURN ENDIF IF(ID.EQ.IS)THEN PRINT *,' #(W02)# ID EQUALS TO IS IN CHECKS.' RETURN ENDIF C CALL VEC01(IDIM,VEC(0,ID)) C DO 20 J=0,IDIM-1 DO 20 K=0,IDIM-1 VEC(J,ID)=VEC(J,ID)+ELMNT(J,K)*VEC(K,IS) 20 CONTINUE C CALL INPRO(IDIM,VEC(0,ID),VEC(0,IS),PROD) HEXPC=PROD PRINT * PRINT 200 200 FORMAT( &' --------------------------- INFORMATION FROM CHECKS') PRINT 210,PROD 210 FORMAT(' =',1PD16.8) DO 40 I=0,IDIM-1 IF(ABS(VEC(I,IS)).GT.1.D-40) THEN VEC(I,ID)=VEC(I,ID)/VEC(I,IS) ELSE VEC(I,ID)=1.D40 ENDIF 40 CONTINUE PRINT 220 220 FORMAT( &' H*X(J)/X(J) (J=MIN(IDIM/3,13),IDIM-1,MAX(1,IDIM/8))') PRINT 230,(VEC(I,ID),I=MIN(IDIM/3,13),IDIM-1,MAX(1,IDIM/8)) 230 FORMAT(1H ,1P4D18.9) PRINT 240 240 FORMAT( &' ---------------------------------------------------') RETURN END C C C SUBROUTINE SMALL(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,NVEC,NE,E,VEC,NDCLRV,MDCLRV,ELMNT,VECS,EIGN,IWK, &NDCLRS) C C****************************************************************** C* CORRESPONDENCE BETWEEN ISTAT1 AND SZ(I) * C* 0 === SZ(I)=-1 * C* 1 === SZ(I)= 0 * C* 2 === SZ(I)=+1 * C***************************************************************** C*** VARIABLES MARKED @ SHOULD BE GIVEN FROM THE MAIN PROGRAM C*** VARIABLES MARKED # ARE EVALUATED AND RETURNED C NS @ LATTICE SIZE C NSLMAX @ 2*NSLMAX DIMENSION OF THE MATRIX C ITSZ @ TOTAL BZ=SZ+1 C ISTAT1(I) # I-TH SPIN CONFIGURATION C ISTAT2 # INVERSE LIST OF ISTAT1 EXPRESSED BY THE C 2-DIM SEARCH ALGORITHM OF NISHINO C NVEC @ NUMBER OF EIGENVECTORS TO CALCULATE C E,EIGN # EIGENVALUES C ITR # NUMBER OF ITERATIONS REQUIRED FOR CONVERGENCE C & IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION E(NE) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION VEC(0:NDCLRV,0:MDCLRV) DIMENSION ELMNT(0:NDCLRS,0:NDCLRS) DIMENSION VECS(0:NDCLRS,0:NVEC),EIGN(0:NE-1) DIMENSION IWK(0:NDCLRS) CHARACTER*4 MSGBC(2) C I=NH(1) I=NSTAT(1,0) IDIM=IBASE(NSB) C CALL CRELMS(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC, &ELMNT,NDCLRS) C IDCLR=NDCLRS+1 C EPS=-1.D-13 C ISW=-1 C IF(IDIM.LT.2) THEN C WRITE (6,*) C & ' #(E01)# WRONG VALUE GIVEN TO IDIM IN SMALL.' C STOP C END IF C IF(NVEC.GT.0)THEN C IF((IDIM .GT.NDCLRS+1).OR.((MDCLRV+1)*(NDCLRV+1) C & .LT.8*IDIM)) THEN C WRITE(6,*) C & ' #(W03)# CANNOT DIAGONALIZE THE HAMILTONIAN IN SMALL.' C RETURN C END IF C CALL DCSMSS(ELMNT,IDCLR,IDIM,EPS,EIGN,NVEC,VECS,IDCLR,ISW, C & IWK,VEC,IERR) C IF(IERR.NE.0)THEN C WRITE(6,*) C & ' #(W04)# ERROR IN DCSMSS CALLED FROM SMALL. IERR=',IERR C RETURN C END IF C IF(NE.GT.NVEC) THEN C CALL CRELMS(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, C & ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC, C & ELMNT,NDCLRS) C CALL DCSMSN(ELMNT,IDCLR,IDIM,EPS,EIGN,NE,ISW,VEC,IERR) C IF(IERR.NE.0)THEN C WRITE(6,*) C & ' #(W05)# ERROR IN DCSMSN CALLED FROM SMALL. IERR=', C & IERR C RETURN C END IF C END IF C ELSE C IF((IDIM .GT.NDCLRS+1).OR.((MDCLRV+1)*(NDCLRV+1) C & .LT.5*IDIM)) THEN C WRITE(6,*) C & ' #(W06)# CANNOT DIAGONALIZE THE HAMILTONIAN IN SMALL. ' C RETURN C END IF C CALL DCSMSN(ELMNT,IDCLR,IDIM,EPS,EIGN,NE,ISW,VEC,IERR) C IF(IERR.NE.0)THEN C WRITE(6,*) C & ' #(W05)# ERROR IN DCSMSN CALLED FROM SMALL. IERR=',IERR C RETURN C ENDIF C ENDIF EPS=1.D-13 IF(IDIM.LT.3) THEN WRITE (6,*) & ' #(E01)# WRONG VALUE GIVEN TO IDIM IN SMALL.' STOP END IF IF(NVEC.GT.0)THEN IF((IDIM .GT.NDCLRS+1).OR.((MDCLRV+1)*(NDCLRV+1) & .LT.8*IDIM)) THEN WRITE(6,*) & ' #(W03)# CANNOT DIAGONALIZE THE HAMILTONIAN IN SMALL. ' RETURN END IF ELSE IF((IDIM .GT.NDCLRS+1).OR.((MDCLRV+1)*(NDCLRV+1) & .LT.6*IDIM)) THEN WRITE(6,*) & ' #(W06)# CANNOT DIAGONALIZE THE HAMILTONIAN IN SMALL. ' RETURN END IF ENDIF CALL DIAGNLZ(ELMNT,IDCLR,IDIM,EIGN,VECS,NE,NVEC,EPS,VEC,IWK) DO 140 I=1,NE E(I)=EIGN(I-1) 140 CONTINUE C IF(NVEC.GT.0)THEN IF(NDCLRV.LT.IDIM-1.OR.MDCLRV.LT.NVEC-1) THEN WRITE (6,*) & ' #(W07)# ARRAY SIZE FOR VEC IN SMALL IS TOO SMALL.' RETURN END IF DO 200 I=0,NVEC-1 DO 150 J=0,IDIM-1 VEC(J,I)=VECS(J,I) 150 CONTINUE 200 CONTINUE ENDIF C RETURN END C C CREATION OF SMALL-SIZE MATRIX ELEMENT C SUBROUTINE CRELMS(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,ELMNT,NDCLRS) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION ELMNT(0:NDCLRS,0:NDCLRS) CHARACTER*4 MSGBC(2) C CALL PAIRCHK(IPAIR,IBOND,NS) IDIM=IBASE(NSB) DO 10 J=0,IDIM-1 DO 10 I=0,IDIM-1 ELMNT(I,J)=0.D0 10 CONTINUE CALL ELMSDIA(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC, &ELMNT,NDCLRS) CALL ELMSOFD(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ELMNT,NDCLRS) C RETURN END C C CALCULATION OF DIAGONAL ELEMENT C SUBROUTINE ELMSDIA(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,ELMNT,NDCLRS) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) CHARACTER*4 MSGBC(2) DIMENSION ELMNT(0:NDCLRS,0:NDCLRS) C IDIM=IBASE(NSB) *VDIR LOOPCNT=5196627 DO 10 I=0,IDIM- 1 ELMNT(I,I)= 0.0D0 10 CONTINUE C IANFLD=0 DO 20 K=1,NS IF(ABS(ANISTRPY(K)).GT.1.D-60 & .OR.ABS(HFIELD(K)).GT.1.D-60) THEN IANFLD=IANFLD+1 ENDIF 20 CONTINUE C IF(IANFLD.NE.0)THEN CALL ELMSDA(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, & NDCLR1, & ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC, & ELMNT,NDCLRS) ENDIF C CALL ELMSDBX(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC, &ELMNT,NDCLRS) C RETURN END C C C SUBROUTINE ELMSDA(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,ELMNT,NDCLRS) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION ELMNT(0:NDCLRS,0:NDCLRS) CHARACTER*4 MSGBC(2) C MSGBC(1)=MSGBC(1) I=IPAIR(1) I=ISTAT2(1) A=AJX(1) A=AJZ(1) A=ABIQD(1) IHSHFT=3**NSLMAX C DO 90 I=0,NSB-1 LL=NSTAT(NL(I),0) IBX=IBASE(I) DO 80 KIS=1,NS IST=3**(KIS-1) AN=ANISTRPY(KIS) HF=HFIELD(KIS) IF(ABS(AN).LE.1.D-60.AND.ABS(HF).LE.1.D-60)THEN GOTO 80 ENDIF DO 30 J=0,NSTAT(NH(I),1)-1 IX=IBX+J*LL *VDIR LOOPCNT=73789 DO 20 K=0,LL -1 ITP=ISTAT1(J,NH(I))*IHSHFT+ISTAT1(K,NL(I)) ELMNT(IX+K,IX+K)=ELMNT(IX+K,IX+K) & -HF*DBLE(MOD(ITP/IST,3)-1) & +AN*DBLE((MOD(ITP/IST,3)-1)**2) 20 CONTINUE 30 CONTINUE 80 CONTINUE 90 CONTINUE C RETURN END C C C SUBROUTINE ELMSDBX(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,ELMNT,NDCLRS) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION ELMNT(0:NDCLRS,0:NDCLRS) CHARACTER*4 MSGBC(2) C MSGBC(1)=MSGBC(1) I=ISTAT2(1) A=AJX(1) A=ANISTRPY(1) A=HFIELD(1) IHSHFT=3**NSLMAX C DO 90 I=0,NSB-1 LL=NSTAT(NL(I),0) IBX=IBASE(I) DO 80 KIS=1,IBOND ISITE1=(IPAIR(KIS*2-1)-1) IS1=3**ISITE1 ISITE2=(IPAIR(KIS*2 )-1) IS2=3**ISITE2 XJZ=AJZ(KIS) XB=ABIQD(KIS) IF(ABS(XJZ).LE.1.D-60.AND.ABS(XB).LE.1.D-60)THEN GOTO 80 ENDIF IF(ABS(XB).GT.1.D-60)THEN XBXX=ABIQD(KIS) XBZZ=ABIQD(KIS) DO 30 J=0,NSTAT(NH(I),1)-1 IX=IBX+LL*J *VDIR LOOPCNT=73789 DO 20 K=0,LL -1 ITP=ISTAT1(J,NH(I))*IHSHFT+ISTAT1(K,NL(I)) IST1=(MOD(ITP/IS1,3)-1) IST2=(MOD(ITP/IS2,3)-1) ELMNT(IX+K,IX+K)=ELMNT(IX+K,IX+K) & +XJZ*DBLE(IST1*IST2) & +XB*DBLE((IST1*IST2)**2 & +1+(IST1*IST2+1)*(IST1*IST2+2)/2 & -(IST1+IST2)**2) CCC & +XBZZ*DBLE((IST1*IST2)**2) CCC & +XBXX*DBLE(1+(IST1*IST2+1)*(IST1*IST2+2)/2 CCC & -(IST1+IST2)**2) 20 CONTINUE 30 CONTINUE ELSE DO 70 J=0,NSTAT(NH(I),1)-1 IX=IBX+LL*J *VDIR LOOPCNT=73789 DO 60 K=0,LL -1 ITP=ISTAT1(J,NH(I))*IHSHFT+ISTAT1(K,NL(I)) ELMNT(IX+K,IX+K)=ELMNT(IX+K,IX+K) & +XJZ*DBLE((MOD(ITP/IS1,3)-1)*(MOD(ITP/IS2,3)- & 1)) 60 CONTINUE C 70 CONTINUE ENDIF 80 CONTINUE 90 CONTINUE C RETURN END C C MULTIPLICATION OF H(OFFDIAG) BY VEC C VEC(I,ID)=SUMJ(HOFFDIAG(I,J)*VEC(J,IS)) C SUBROUTINE ELMSOFD(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ELMNT,NDCLRS) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ELMNT(0:NDCLRS,0:NDCLRS) C CALL ELMSLLBX(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ELMNT,NDCLRS) C IF(NS.GE.NSLMAX+1 ) THEN CALL ELMSLHBX(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, & NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND, & ELMNT,NDCLRS) END IF C IF(NS.GE.NSLMAX+2) THEN CALL ELMSHHBX(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, & NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND, & ELMNT,NDCLRS) END IF C RETURN END C C MULTIPLICATION OF H(LOW,LOW) BY VEC WITH BIQUADRATIC TERM C VEC(ID)=H(LOW,LOW)*VEC(ISRC) C SUBROUTINE ELMSLLBX(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ELMNT,NDCLRS) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ELMNT(0:NDCLRS,0:NDCLRS) C I=NS A=AJZ(1) DO 90 I=0,NSB-1 ITL=NL(I) LL=NSTAT(ITL,0) IBX=IBASE(I) DO 80 KIS=1,IBOND ISITE1=(IPAIR(KIS*2-1)-1) ISITE2=(IPAIR(KIS*2 )-1) IF((ISITE1+1.LE.NSLMAX).AND.(ISITE2+1.LE.NSLMAX))THEN MSKH=3**ISITE1 MSKL=3**ISITE2 ELSE GOTO 80 ENDIF XJX=AJX(KIS) XB=ABIQD(KIS) IF(ABS(XJX).LE.1.D-60.AND.ABS(XB).LE.1.D-60) THEN GOTO 80 ENDIF IF(ABS(XB).GT.1.D-60) THEN XBXX=XB XBXZ=XB DO 30 J=0,NSTAT(NH(I),1)-1 IX=IBX+LL*J *VDIR NODEP *VDIR LOOPCNT=73789 DO 20 K=0,LL -1 ITP=ISTAT1(K,ITL) IF(MOD(ITP/MSKH,3).NE.0.AND.MOD(ITP/MSKL,3) & .NE.2) THEN XBJ=XJX-XBXZ*DBLE(1- & (MOD(ITP/MSKH,3)+MOD(ITP/MSKL,3)-2)**2) IDP=ISTAT2(ITP - MSKH+MSKL) ELMNT(IX+ K,IX+IDP)=ELMNT(IX+ K,IX+IDP)+ & XBJ ELMNT(IX+IDP,IX+K)=ELMNT(IX+IDP,IX+K)+XBJ END IF IF(MOD(ITP/MSKH,3).EQ.2.AND.MOD(ITP/MSKL,3) & .EQ.0) THEN IDP=ISTAT2(ITP -2*MSKH+2*MSKL) ELMNT(IX+ K,IX+IDP)=ELMNT(IX+ K,IX+IDP)+ & XBXX ELMNT(IX+IDP,IX+K)=ELMNT(IX+IDP,IX+K)+XBXX END IF 20 CONTINUE 30 CONTINUE C ELSE DO 70 J=0,NSTAT(NH(I),1)-1 IX=IBX+LL*J *VDIR NODEP *VDIR LOOPCNT=73789 DO 60 K=0,LL -1 ITP=ISTAT1(K,ITL) IF(MOD(ITP/MSKH,3).NE.0.AND.MOD(ITP/MSKL,3) & .NE.2) THEN IDP=ISTAT2(ITP - MSKH+MSKL) ELMNT(IX+ K,IX+IDP)=ELMNT(IX+ K,IX+IDP)+ & XJX ELMNT(IX+IDP,IX+K)=ELMNT(IX+IDP,IX+K)+XJX END IF 60 CONTINUE 70 CONTINUE ENDIF 80 CONTINUE 90 CONTINUE C RETURN END C C MULTIPLICATION OF H(LOW,HIGH) BY ELMNT WITH BIQUADRATIC TERM C ELMNT(ID)=H(LOW,HIGH)*ELMNT(ISRC) C SUBROUTINE ELMSLHBX(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ELMNT,NDCLRS) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ELMNT(0:NDCLRS,0:NDCLRS) C I=NS A=AJZ(1) DO 90 I=0,NSB-1 ITH=NH(I) ITL=NL(I) LH=NSTAT(ITH,1) LL=NSTAT(ITL,0) IBX=IBASE(I) IBX1=IBASE(I+1) LL1=NSTAT(ITL+1,0) DO 80 KIS=1,IBOND ISITE1=(IPAIR(KIS*2-1)-1) ISITE2=(IPAIR(KIS*2 )-1) IF((ISITE1+1.LE.NSLMAX).AND.(ISITE2+1.GT.NSLMAX))THEN MSKL=3**ISITE1 MSKH=3**(ISITE2-NSLMAX) ELSEIF((ISITE1+1.GT.NSLMAX).AND.(ISITE2+ & 1.LE.NSLMAX))THEN MSKH=3**(ISITE1-NSLMAX) MSKL=3**ISITE2 ELSE GOTO 80 ENDIF XJX=AJX(KIS) XB=ABIQD(KIS) IF(ABS(XJX).LE.1.D-60.AND.ABS(XB).LE.1.D-60) THEN GOTO 80 ENDIF IF(ABS(XB).GT.1.D-60) THEN XBXZ=XB XBXX=XB DO 40 J=0,NSTAT(ITH,1)-1 IX=IBX+J*LL ITPH=ISTAT1(J,ITH) IF(MOD(ITPH/MSKH,3) .NE. 0) THEN IY=IBX1+LL1*ISTAT2(ITPH-MSKH) *VDIR NODEP *VDIR LOOPCNT=73789 DO 20 K=0,LL - 1 ITPL=ISTAT1(K,ITL) IF(MOD(ITPL/MSKL,3).NE. 2)THEN XBJ=XJX-XBXZ*DBLE(1- & (MOD(ITPH/MSKH,3)+MOD(ITPL/MSKL,3)-2) & **2) IDP=ISTAT2(ITPL+MSKL) ELMNT(IY+IDP,IX+K)=ELMNT(IY+IDP,IX+K)+ & XBJ ELMNT(IX+ K,IY+IDP)=ELMNT(IX+ K,IY+IDP) & +XBJ END IF 20 CONTINUE END IF IF(MOD(ITPH/MSKH,3).EQ.2)THEN IBX2=IBASE(I+2) LL2=NSTAT(ITL+2,0) IY=IBX2+LL2*ISTAT2(ITPH-2*MSKH) *VDIR NODEP *VDIR LOOPCNT=73789 DO 30 K=0,LL - 1 ITPL=ISTAT1(K,ITL) IF(MOD(ITPL/MSKL,3).EQ.0)THEN IDP=ISTAT2(ITPL+2*MSKL) ELMNT(IY+IDP,IX+K)=ELMNT(IY+IDP,IX+K)+ & XBXX ELMNT(IX+ K,IY+IDP)=ELMNT(IX+ K,IY+IDP) & +XBXX END IF 30 CONTINUE END IF 40 CONTINUE C ELSE DO 70 J=0,NSTAT(ITH,1)-1 IX=IBX+J*LL ITPH=ISTAT1(J,ITH) IF(MOD(ITPH/MSKH,3) .NE. 0) THEN IY=IBX1+LL1*ISTAT2(ITPH-MSKH) *VDIR NODEP *VDIR LOOPCNT=73789 DO 60 K=0,LL - 1 ITPL=ISTAT1(K,ITL) IF(MOD(ITPL/MSKL,3).NE. 2)THEN IDP=ISTAT2(ITPL+MSKL) ELMNT(IY+IDP,IX+K)=ELMNT(IY+IDP,IX+K)+ & XJX ELMNT(IX+ K,IY+IDP)=ELMNT(IX+ K,IY+IDP) & +XJX END IF 60 CONTINUE END IF 70 CONTINUE ENDIF 80 CONTINUE 90 CONTINUE C RETURN END C C MULTIPLICATION OF H(HIGH,HIGH) BY ELMNT WITH BIQUADRATIC TERM C ELMNT(ID)=H(HIGH,HIGH)*ELMNT(IS) C SUBROUTINE ELMSHHBX(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ELMNT,NDCLRS) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ELMNT(0:NDCLRS,0:NDCLRS) C I=NS A=AJZ(1) DO 90 I=0,NSB-1 ITH=NH(I) ITL=NL(I) LH=NSTAT(ITH,1) LL=NSTAT(ITL,0) IBX=IBASE(I) DO 80 KIS=1,IBOND ISITE1=(IPAIR(KIS*2-1)-1) ISITE2=(IPAIR(KIS*2 )-1) IF((ISITE1+1.GT.NSLMAX).AND.(ISITE2+1.GT.NSLMAX))THEN MSKH=3**(ISITE1-NSLMAX) MSKL=3**(ISITE2-NSLMAX) ELSE GOTO 80 ENDIF XJX=AJX(KIS) XB=ABIQD(KIS) IF(ABS(XJX).LE.1.D-60.AND.ABS(XB).LE.1.D-60) THEN GOTO 80 ENDIF IF(ABS(XB).GT.1.D-60) THEN XBXZ=XB XBXX=XB DO 40 J=0,LH - 1 ITP=ISTAT1(J,ITH) IF(MOD(ITP/MSKH,3).NE.0 .AND. MOD(ITP/MSKL,3) & .NE.2) THEN XBJ=XJX-XBXZ*DBLE(1- & (MOD(ITP/MSKH,3)+MOD(ITP/MSKL,3)-2)**2) IX=IBX+LL*J IY=IBX+LL*ISTAT2(ITP - MSKH+MSKL) *VDIR NODEP *VDIR LOOPCNT=73789 DO 20 K=0,LL - 1 ELMNT(IY+K,IX+K)=ELMNT(IY+K,IX+K)+XBJ ELMNT(IX+K,IY+K)=ELMNT(IX+K,IY+K)+XBJ 20 CONTINUE END IF IF(MOD(ITP/MSKH,3).EQ.2.AND.MOD(ITP/MSKL,3) & .EQ.0)THEN IX=IBX+LL*J IY=IBX+LL*ISTAT2(ITP-2*MSKH+2*MSKL) *VDIR NODEP *VDIR LOOPCNT=73789 DO 30 K=0,LL - 1 ELMNT(IY+K,IX+K)=ELMNT(IY+K,IX+K)+XBXX ELMNT(IX+K,IY+K)=ELMNT(IX+K,IY+K)+XBXX 30 CONTINUE END IF 40 CONTINUE C ELSE DO 70 J=0,LH - 1 ITP=ISTAT1(J,ITH) IF(MOD(ITP/MSKH,3).NE.0 .AND. MOD(ITP/MSKL,3) & .NE.2) THEN IX=IBX+LL*J IY=IBX+LL*ISTAT2(ITP - MSKH+MSKL) *VDIR NODEP *VDIR LOOPCNT=73789 DO 60 K=0,LL - 1 ELMNT(IY+K,IX+K)=ELMNT(IY+K,IX+K)+XJX ELMNT(IX+K,IY+K)=ELMNT(IX+K,IY+K)+XJX 60 CONTINUE END IF 70 CONTINUE ENDIF 80 CONTINUE 90 CONTINUE C RETURN END C c****************************************************************** c* From TITPACK Version 2 * c* Copyright (C) Hidetoshi Nishimori * c* Thanks to Professor H. Nishimori * c****************************************************************** c c***************** eigenvalues of a small matrix ****************** c subroutine DIAGNLZ(elmnt,idclr,idim,E,v,ne,nvec,eps,wk,iwk) c c elmnt @ matrix elements c idclr @ declared array size in the main program c idim @ matrix dimension c E # eigenvalues c v # eigenvector c ne @ number of eigenvalues to calculate c nvec @ number of eigenvectors to calculate c eps @ limit of error c wk,iwk working areas c implicit real*8(a-h,o-z) dimension E(ne),v(idclr,nvec),elmnt(idclr,idim) dimension wk(idim,8),iwk(idim) c if(nvec.lt.0.or.nvec.gt.ne)then print *,' #(E02)# NVEC GIVEN TO DIAGNLZ, ', & 'CALLED FROM SMALL, OUT OF RANGE.' stop end if c call hshldr(elmnt,idclr,idim,wk(1,1),wk(1,2),wk(1,3),wk(1,4) &,wk(1,5),wk(1,6)) call bisec(wk(1,1),wk(1,2),idim,E,ne,eps) if(nvec.eq.0)return call vec3(E,elmnt,idclr,idim,ne,nvec,wk(1,4),wk(1,5),wk(1,6) &,wk(1,7),wk(1,8),iwk,wk(1,1),wk(1,2),wk(1,3),v) return end c c***************** Householder tridiagonalization ***************** c subroutine hshldr(elmnt,idclr,idim,alpha,beta,c,w,p,q) c c elmnt @ matrix c idclr @ declared array size in the main program c idim @ matrix dimension c alpha # diagonal element c beta # subdiagonal element c c # normalization factor for eigenvector calculation c w,p,q working areas c implicit real*8(a-h,o-z) dimension elmnt(idclr,idim),alpha(idim),beta(idim),c(idim) dimension w(idim),p(idim),q(idim) c do 10 k=1,idim-2 s=0.d0 do 20 i=k+1,idim 20 s=s+elmnt(i,k)**2 s=sqrt(s) if(elmnt(k+1,k).lt.0.d0)s=-s c alpha(k)=elmnt(k,k) beta(k)=-s if(s**2.lt.1.d-50)goto 10 c c(k)=1.d0/(s**2+elmnt(k+1,k)*s) w(k+1)=elmnt(k+1,k)+s do 30 i=k+2,idim 30 w(i)=elmnt(i,k) elmnt(k+1,k)=w(k+1) c do 40 i=k+1,idim t=0.d0 do 50 j=k+1,i 50 t=t+elmnt(i,j)*w(j) do 55 j=i+1,idim 55 t=t+elmnt(j,i)*w(j) p(i)=c(k)*t 40 continue c t=0.d0 do 60 i=k+1,idim 60 t=t+p(i)*w(i) s=0.5d0*c(k)*t do 70 i=k+1,idim 70 q(i)=p(i)-s*w(i) c do 80 j=k+1,idim do 80 i=j,idim 80 elmnt(i,j)=elmnt(i,j)-w(i)*q(j)-q(i)*w(j) c 10 continue c alpha(idim-1)=elmnt(idim-1,idim-1) alpha(idim)=elmnt(idim,idim) beta(idim-1)=elmnt(idim,idim-1) return end c c**** eigenvector of a tridiagonal matrix by inverse iteration *** c subroutine vec3(E,elmnt,idclr,idim,ne,nvec,di,bl,bu,bv,cm, &lex,alpha,beta,c,v) c c E @ eigenvalues c elmnt @ matrix elements c idclr @ declared array size in the main program c idim @ matrix dimension c ne @ number of eigenvalues c nvec @ number of eigenvectors c di - lex working areas c alpha @ diagonal element in the triDIAGNLZonal form c beta @ subdiagonal element in the triDIAGNLZonal form c c @ normalization factor given from hshldr c v # eigenvectors c implicit real*8 (a-h,o-z) dimension E(ne),elmnt(idclr,idim) dimension di(idim),bl(idim),bu(idim),bv(idim),cm(idim), &lex(idim) dimension alpha(idim),beta(idim),c(idim),v(idclr,nvec) c if(nvec.gt.ne)then print *,' #(E03)# NVEC GIVEN TO VEC3, CALLED ORIGINALLY ', & 'FROM SMALL, OUT OF RANGE.' stop end if do 10 k=1,nvec c do 100 j=1,idim di(j)=E(k)-alpha(j) bl(j)=-beta(j) bu(j)=-beta(j) 100 continue c c*** LU decomposition do 110 j=1,idim-1 if(abs(di(j)).gt.abs(bl(j)))then c--- non pivoting lex(j)=0 if(abs(di(j)).lt.1.d-40)di(j)=1.d-40 cm(j+1)=bl(j)/di(j) di(j+1)=di(j+1)-cm(j+1)*bu(j) bv(j)=0.d0 else c--- pivoting lex(j)=1 cm(j+1)=di(j)/bl(j) di(j)=bl(j) s=bu(j) bu(j)=di(j+1) bv(j)=bu(j+1) di(j+1)=s-cm(j+1)*bu(j) bu(j+1)= -cm(j+1)*bv(j) end if 110 continue if(abs(di(idim)).lt.1.d-40)di(idim)=1.d-40 c c--- initial vector do 120 j=1,idim 120 v(j,k)=1.d0/dble(j*5) c c*** degeneracy check up if(k.eq.1)then km=k else if(abs(E(k)-E(km)).gt.1.d-13)then km=k else do 130 i=km,k-1 prd=0.d0 do 140 j=1,idim 140 prd=prd+v(j,i)*v(j,k) do 150 j=1,idim 150 v(j,k)=v(j,k)-prd*v(j,i) 130 continue end if c c*** inverse iteration do 160 l=1,k-km+3 if((l.ne.1).or.(k.ne.km))then c--- forward substitution do 170 j=1,idim-1 if(lex(j).eq.0)then v(j+1,k)=v(j+1,k)-cm(j+1)*v(j,k) else s=v(j,k) v(j,k)=v(j+1,k) v(j+1,k)=s-cm(j+1)*v(j,k) end if 170 continue end if c--- backward substitution do 180 j=idim,1,-1 s=v(j,k) if(j.le.idim-1)s=s-bu(j)*v(j+1,k) if(j.le.idim-2)s=s-bv(j)*v(j+2,k) v(j,k)=s/di(j) 180 continue c c*** normalization dnorm=0.d0 do 190 j=1,idim 190 dnorm=dnorm+v(j,k)**2 if(dnorm.gt.1.d-50)dnorm=1./sqrt(dnorm) do 200 j=1,idim 200 v(j,k)=v(j,k)*dnorm 160 continue c 10 continue c c*** back transformation to the original representation do 210 k=1,nvec c do 220 i=idim-2,1,-1 prd=0.d0 do 230 j=i+1,idim 230 prd=prd+elmnt(j,i)*v(j,k) s=prd*c(i) do 240 j=i+1,idim 240 v(j,k)=v(j,k)-s*elmnt(j,i) 220 continue 210 continue c c*** orthogonalization for degenerate case km=1 do 250 k=2,nvec if(abs(E(k)-E(km)).ge.1.0d-13)then km=k else do 260 i=km,k-1 prd=0.d0 do 270 j=1,idim 270 prd=prd+v(j,i)*v(j,k) do 280 j=1,idim 280 v(j,k)=v(j,k)-prd*v(j,i) 260 continue c dnorm=0.0d0 do 290 j=1,idim 290 dnorm=dnorm+v(j,k)**2 s=1.d0/sqrt(dnorm) do 300 j=1,idim 300 v(j,k)=v(j,k)*s c end if 250 continue return end C C****************************************************************** C* KOBEPACK/1 VERSION 1.0 * C* MARCH 1992 * C* COPYRIGHT (C) M. KABURAGI, T. NISHINO AND T. TONEGAWA * C* WITH THE HELP OF TITPACK VERSION 2 * C* FEBRUARY 1991 * C* COPYRIGHT (C) HIDETOSHI NISHIMORI * C****************************************************************** C***** MID-SIZE(FOR ASL) ***** C****************************************************************** C* CORRESPONDENCE BETWEEN ISTAT1 AND SZ(I) * C* 0 === SZ(I)=-1 * C* 1 === SZ(I)= 0 * C* 2 === SZ(I)=+1 * C****************************************************************** C*** VARIABLES MARKED @ SHOULD BE GIVEN FROM THE MAIN PROGRAM C*** VARIABLES MARKED # ARE EVALUATED AND RETURNED C*** THE FOLLOWING VARIABLES ARE COMMON TO ALL ROUTINES C NS @ LATTICE SIZE C NSLMAX @ MAXIMUM SIZE OF LOWER PART OF LATTICE C ITSZ @ TOTAL SZ C NSB,NH,NL, @ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C MSGBC,NSTAT, @ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C ISTAT1,NDCLR1,@ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C ISTAT2,IBASE @ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C ISTAT1(I,IBZ) @ I-TH SPIN CONFIGURATION FOR IBZL(OR IBZH)=IBZ C BZ=SZ+1 C ISTAT2 @ INVERSE LIST OF ISTAT1 EXPRESSED BY THE C 2-DIM SEARCH ALGORITHM OF NISHINO C IPAIR @ PAIRS OF SITES CONNECTED BY BONDS C AJX @ EXCHANGE INTERACTION OF EACH BOND JX,Y C AJZ @ EXCHANGE INTERACTION OF EACH BOND JZ C ABIQD @ BIQUADRATIC INTERCTION C IBOND @ NUMBER OF BONDS C ANISTRPY @ SINGLE ION ANISOTROPY ENERGY AT EACH SITE C HFIELD @ MAGNETIC FIELD AT EACH SITE C DIAG @ DIAGONAL ELEMENT C LBOND(I,1) @# BOND NUMBER OF I-TH L-L BONDS C NBOND(1) @# TOTAL NUMBER OF L-L BONDS C LBOND(I,2) @# BOND NUMBER OF I-TH H-H BONDS C NBOND(2) @# TOTAL NUMBER OF H-H BONDS C LBOND(I,3) @# BOND NUMBER OF I-TH L-H BONDS C NBOND(3) @# TOTAL NUMBER OF L-H BONDS C IBSPIN(I,1) @# SITE NUMBER OF I-TH BOUNDARY L-SPINS C NBSPIN(1) @# TOTAL NUMBER OF BOUNDARY L-SPINS C IBSPIN(I,2) @# SITE NUMBER OF I-TH BOUNDARY H-SPINS C NBSPIN(2) @# TOTAL NUMBER OF BOUNDARY H-SPINS C IBSPAIR(2K-1,2K) C # BOUNDARY SPIN NUMBER OF K-TH BOUNDARY PAIRS C (2K-1:L, 2K:H) C IPAIR(2K-1,2K) IS REPLACED WITH IPAIR(2K-1):L OR H C LBASE(I,1) @# BASE ADDRESS FOR L-L DIAGONAL ELEMENT C LBASE(I,2) @# BASE ADDRESS FOR H-H DIAGONAL ELEMENT C LBASE(I,2+J) @# BASE ADDRESS FOR BOUNDARY L/H-SPIN DIAGONAL ELEMENT C MBASE(I,K,J) @# BASE ADDRESS FOR K-BOND OFF-DIAGONAL ELEMENT C IFLGDELM(I,1) @# FLAG FOR L-L DIAGONAL ELEMENT C IFLGDELM(I,2) @# FLAG FOR H-H DIAGONAL ELEMENT C IFLGDELM(I,2+J)@# FLAG FOR BOUNDARY L/H-SPIN DIAGONAL ELEMENT C IFLGOELM(I,K,J)@# FLAG FOR K-BOND OFF-DIAGONAL ELEMENT C ELMNT @ OFF-DIAGONAL ELEMENT BY NISHIMORI'S EXPRESSION C LOC @ LOCATION OF ADDRESS BY NISHIMORI'S EXPRESSION C C****************************************************************** C C CHECK OF THE EIGENVECTOR AND EIGENVALUE C SUBROUTINE CHECKM(ID,IS,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY, &HFIELD,MSGBC,DIAG,NDCLRD,VEC,NDCLRV,MDCLRV,LBOND,NBOND, &LBASE,MBASE,IFLGDELM,IFLGOELM,ELMNT,LOC,NDCLRM,HEXPC) C C & C*** VARIABLES MARKED @ SHOULD BE GIVEN FROM THE MAIN PROGRAM C*** VARIABLES MARKED # ARE EVALUATED AND RETURNED C VEC(*,ISRC) @ EIGENVECTOR TO BE CHECKED C VEC(*,IDIS) # H*VEC(ISRC)/VEC(ISRC) C HEXPEC # C VEC # EIGEN VECTOR AND WORK SPACE C NDCLRV @ SIZE OF VEC C & IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION DIAG(0:NDCLRD),VEC(0:NDCLRV,0:MDCLRV) CHARACTER*4 MSGBC(2) DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION MBASE(0:2*NSLMAX,IBOND+NS,6) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) C IDIM=IBASE(NSB) CALL INPRO(IDIM,VEC(0,IS),VEC(0,IS),PROD) PROD=SQRT(PROD) IF(PROD.LT.1.D-30)THEN PRINT *,' #(W08)# NULL VECTOR GIVEN TO CHECKM.' RETURN ENDIF IF(ID.EQ.IS)THEN PRINT *,' #(W09)# ID EQUALS TO IS IN CHECKM.' RETURN ENDIF CALL VEC01(IDIM,VEC(0,ID)) CALL MLTPLY(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC,DIAG, &NDCLRD,VEC(0,ID),VEC(0,IS),IDIM,LBOND,NBOND,LBASE,MBASE, &IFLGDELM,IFLGOELM,ELMNT,LOC,NDCLRM) C CALL INPRO(IDIM,VEC(0,ID),VEC(0,IS),PROD) HEXPC=PROD PRINT * PRINT 200 200 FORMAT( &' --------------------------- INFORMATION FROM CHECKM') PRINT 210,PROD 210 FORMAT(' =',1PD16.8) DO 40 I=0,IDIM-1 IF(ABS(VEC(I,IS)).GT.1.D-40) THEN VEC(I,ID)=VEC(I,ID)/VEC(I,IS) ELSE VEC(I,ID)=1.D40 ENDIF 40 CONTINUE PRINT 220 220 FORMAT(' H*X(J)/X(J) (J=13,IDIM-1,IDIM/8)') PRINT 230,(VEC(I,ID),I=13,IDIM-1,IDIM/8) 230 FORMAT(1H ,1P4D18.9) PRINT 240 240 FORMAT( &' ---------------------------------------------------') RETURN END C C DUMMY ROUTINE FOR LANCZOS C SUBROUTINE LNCZM(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,DIAG,NDCLRD,VEC,NDCLRV,MDCLRV,LBOND,NBOND,LBASE,MBASE, &IFLGDELM,IFLGOELM,ELMNT,LOC,NDCLRM,NVEC,IV,E,ITR) C C*** VARIABLES MARKED @ SHOULD BE GIVEN FROM THE MAIN PROGRAM C*** VARIABLES MARKED # ARE EVALUATED AND RETURNED C NVEC @ NUMBER OF EIGENVECTORS TO CALCULATE IN LANCZ1 C IV @ LOCATION OF THE NONZERO ELEMENT OF THE INITIAL C E # EIGENVALUES C ITR # NUMBER OF ITERATIONS REQUIRED FOR CONVERGENCE C VEC # EIGEN VECTOR AND WORK SPACE C NDCLRV @ SIZE OF VEC C & IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION DIAG(0:NDCLRD),VEC(0:NDCLRV,0:MDCLRV) CHARACTER*4 MSGBC(2),MSG DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION MBASE(0:2*NSLMAX,IBOND+NS,6) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) DIMENSION E(4) DIMENSION IV(20) C MSG=MSGBC(1) IDIM=IBASE(NSB) CALL PAIRCHK(IPAIR,IBOND,NS) IF(MDCLRV.LT.1)THEN PRINT *,' #(W10)# WRONG VALUE GIVEN TO MDCLRV IN LNCZM.' RETURN END IF IF(NVEC.LT.0)THEN PRINT *,' #(W11)# WRONG VALUE GIVEN TO NVEC IN LNCZM.' PRINT *,' ONLY THE EIGENVALUES ARE CALCULATED.' NVEC=0 END IF IF(NVEC.GT.4)THEN PRINT *,' #(W12)# WRONG VALUE GIVEN TO NVEC IN LNCZM.' PRINT *,' NVEC IS REPLACED BY 4.' NVEC=4 END IF C CALL LNCZ2E(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC,DIAG, &NDCLRD,VEC,NDCLRV,MDCLRV,LBOND,NBOND,LBASE,MBASE,IFLGDELM, &IFLGOELM,ELMNT,LOC,NDCLRM,NVEC,IV,E,ITR) C RETURN END C C LANCZOS METHOD C SUBROUTINE LNCZ2E(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,DIAG,NDCLRD,VEC,NDCLRV,MDCLRV,LBOND,NBOND,LBASE,MBASE, &IFLGDELM,IFLGOELM,ELMNT,LOC,NDCLRM,NVEC,IV,E,ITR) C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (ITRMAX=1000,NROTAT=5,ITRLNCZ=ITRMAX/NROTAT) PARAMETER (ERRLNCZ=1.D-13) PARAMETER (NDCLRWK=5*(ITRLNCZ+2)) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION DIAG(0:NDCLRD),VEC(0:NDCLRV,0:MDCLRV) CHARACTER*4 MSGBC(2),MSG DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION MBASE(0:2*NSLMAX,IBOND+NS,6) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) DIMENSION E(4) DIMENSION IV(20),IDS(0:3) DIMENSION IWK1(NDCLRWK),WK1(NDCLRWK,6),EE(NDCLRWK) COMMON/VECDAT/COEF(NDCLRWK,5),ALPHA(NDCLRWK),BETA(NDCLRWK) C MSG=MSGBC(1) NDCLRW=NDCLRWK IDS(0)=1 IDS(1)=0 IDS(2)=2 IDS(3)=3 I0=IDS(0) I1=IDS(1) I2=IDS(2) I3=IDS(3) IVEC=2 CALL PAIRCHK(IPAIR,IBOND,NS) IDIM=IBASE(NSB) CALL VEC02(IDIM,VEC(0,I0),VEC(0,I1)) CALL VECINI(IDIM,VEC(0,I0),IV) CALL MLTPLY(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC,DIAG, &NDCLRD,VEC(0,I1),VEC(0,I0),IDIM,LBOND,NBOND,LBASE,MBASE, &IFLGDELM,IFLGOELM,ELMNT,LOC,NDCLRM) C LK IS AN INDEX FOR ALPHA AND BETA EPS=1.D-10 LK=1 E2P=1.D50 DO 50 LLANCZ=1,ITRLNCZ C REPORT PER 5 ROTATIONS DO 40 J=0,NROTAT-1 CALL MACRO2(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, & ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND, & ANISTRPY,HFIELD,MSGBC,DIAG,NDCLRD, & VEC(0,I0),VEC(0,I1),IDIM,VEC(0,IVEC),NDCLRV,NVEC, & LBOND,NBOND,LBASE,MBASE,IFLGDELM,IFLGOELM, & ELMNT,LOC,NDCLRM,LK,0,PRDA) IF(LK.LT.0)THEN ITR=-LLANCZ RETURN ENDIF 40 CONTINUE LL=LLANCZ C C CALL MALNCNEC(LL,LK,ALPHA,BETA,EE,E2P,COEF,NDCLRW,NVEC, C & ERRLNCZ,IWK1,WK1,INDRET) C CALL MALNCTIT(LL,LK,ALPHA,BETA,EE,E2P,COEF,NDCLRW,NVEC, & ERRLNCZ,IWK1,WK1,INDRET) C E(1)=EE(1) E(2)=EE(2) E(3)=EE(3) E(4)=EE(4) IF(INDRET.NE.0)THEN ITR=INDRET*LLANCZ RETURN ENDIF 50 CONTINUE WRITE(6,*) ' #(W13)# ITERATION IN LNCZ2E, CALLED FROM LNCZM, ', & 'EXCEEDS THE MAXIMUM (1000).' INDRET=1 IF(NVEC.GT.0)THEN LL=ITRLNCZ ERRLN=1.D10 C CALL MALNCNEC(LL,LK,ALPHA,BETA,EE,E2P,COEF,NDCLRW,NVEC, C & ERRLN,IWK1,WK1,INDRET) CALL MALNCTIT(LL,LK,ALPHA,BETA,EE,E2P,COEF,NDCLRW,NVEC, & ERRLN,IWK1,WK1,INDRET) ENDIF ITR=INDRET*ITRLNCZ RETURN END C C IVEC-TH EIGEN VECTOR, VEC(*,0), SHOULD BE THE SAME AS THAT IN LNCZL C SUBROUTINE LNCZMINV(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,DIAG,NDCLRD,VEC,NDCLRV,MDCLRV,LBOND,NBOND,LBASE,MBASE, &IFLGDELM,IFLGOELM,ELMNT,LOC,NDCLRM,NVEC,IV,E,ITR) C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (ITRMAX=1000,NROTAT=5,ITRLNCZ=ITRMAX/NROTAT) PARAMETER (ERRLNCZ=1.D-13) PARAMETER (NDCLRWK=5*(ITRLNCZ+2)) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION DIAG(0:NDCLRD),VEC(0:NDCLRV,0:MDCLRV) CHARACTER*4 MSGBC(2),MSG DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION MBASE(0:2*NSLMAX,IBOND+NS,6) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) DIMENSION E(4) DIMENSION IV(20),IDS(0:3) COMMON/VECDAT/COEF(NDCLRWK,5),ALPHA(NDCLRWK),BETA(NDCLRWK) C MSG=MSGBC(1) IF(ITR.LT.0) THEN PRINT *,' #(W14)# NEGATIVE VALUE GIVEN TO ITR ', & 'IN LNCZMINV. ITR=',ITR RETURN END IF IF(NVEC+2.GT.MDCLRV+1)THEN NVEC=MDCLRV-1 PRINT *,' #(W15)# WRONG VALUE GIVEN TO NVEC IN LNCZMINV.' PRINT *,' NVEC IS REPLACED BY ',NVEC,'.' END IF IF(NVEC.LE.0)THEN PRINT *,' #(W16)# WRONG VALUE GIVEN TO NVEC IN LNCZMINV.' NVEC=0 RETURN END IF IF(NVEC.GT.4)THEN NVEC=4 PRINT *,' #(W17)# WRONG VALUE GIVEN TO NVEC IN LNCZMINV.' PRINT *,' NVEC IS REPLACED BY ',NVEC,'.' END IF C IDS(0)=NVEC IDS(1)=NVEC+1 IDS(2)=NVEC+2 IDS(3)=NVEC+3 IVEC=0 I0=IDS(0) I1=IDS(1) I2=IDS(2) I3=IDS(3) LLANCZ=ITR NDCLRW=NDCLRWK CALL PAIRCHK(IPAIR,IBOND,NS) IDIM=IBASE(NSB) CALL VEC02(IDIM,VEC(0,I0),VEC(0,I1)) CALL VECINI(IDIM,VEC(0,I0),IV) C (I) LK=1 DO 10 I=1,NVEC CALL ADDER(IDIM,VEC(0,I-1),VEC(0,I-1),VEC(0,I0),0.0D0, & COEF(LK,I)) 10 CONTINUE CALL MLTPLY(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC,DIAG, &NDCLRD,VEC(0,I1),VEC(0,I0),IDIM,LBOND,NBOND,LBASE,MBASE, &IFLGDELM,IFLGOELM,ELMNT,LOC,NDCLRM) DO 150 I=1,LLANCZ DO 140 J=0,4 CALL MACRO2(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, & ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND, & ANISTRPY,HFIELD,MSGBC,DIAG,NDCLRD, & VEC(0,I0),VEC(0,I1),IDIM,VEC(0,IVEC),NDCLRV,NVEC, & LBOND,NBOND,LBASE,MBASE,IFLGDELM,IFLGOELM, & ELMNT,LOC,NDCLRM, & LK,1,PRDA) 140 CONTINUE 150 CONTINUE IF(MDCLRV.LT.NVEC+2)THEN RETURN ENDIF IDS(1)=NVEC IDS(2)=NVEC+1 IDS(3)=NVEC+2 I1=IDS(1) I2=IDS(2) I3=IDS(3) DO 160 I=1,NVEC IDS(0)=I-1 I0=IDS(0) ENERGY=E(I) CALL INVERS2(NS,NSLMAX,NSB,NH,NL,IBASE, & NSTAT,ISTAT1,NDCLR1,ISTAT2, & IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC, & DIAG,NDCLRD,VEC,NDCLRV,MDCLRV, & LBOND,NBOND,LBASE,MBASE,IFLGDELM,IFLGOELM, & ELMNT,LOC,NDCLRM,NVEC,ENERGY,IDS) CALL VMOVE(IDIM,VEC(0,I-1),VEC(0,I1)) 160 CONTINUE RETURN END C C MACRO OF THE CALCULATION OF ALPHA AND BETA IN LANCZOS METHOD C SUBROUTINE MACRO2(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,DIAG,NDCLRD,VEC0,VEC1,IDIM,VEC,NDCLRV,NVEC,LBOND, &NBOND,LBASE,MBASE,IFLGDELM,IFLGOELM,ELMNT,LOC,NDCLRM,LK,ISW, &PRDA) C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (ITRMAX=1000,NROTAT=5,ITRLNCZ=ITRMAX/NROTAT) PARAMETER (ERRLNCZ=1.D-13) PARAMETER (NDCLRWK=5*(ITRLNCZ+2)) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION DIAG(0:NDCLRD) DIMENSION VEC0(0:IDIM-1),VEC1(0:IDIM-1),VEC(0:NDCLRV,0:NVEC) CHARACTER*4 MSGBC(2),MSG DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION MBASE(0:2*NSLMAX,IBOND+NS,6) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) COMMON/VECDAT/COEF(NDCLRWK,5),ALPHA(NDCLRWK),BETA(NDCLRWK) C MSG=MSGBC(1) CALL INPRO(IDIM,VEC0,VEC1,PRDA) ALPHA(LK)=PRDA CALL NORMER(IDIM,PRDB,VEC1,VEC0,1.D0,-PRDA) BETA(LK)=PRDB IF(BETA(LK).LT.0.5D-30)THEN PRINT *, & ' #(W18)# TRIDIAGONALIZATION UNSUCCESSFUL IN MACRO2 ', & 'CALLED ORIGINALLY' PRINT *,' FROM LNCZM. BETA(LK) IS TOO SMALL ', & 'AT LK=',LK,'.' LK=-LK RETURN END IF A00=-ALPHA(LK)/BETA(LK) A01=1.D0/BETA(LK) A10=-BETA(LK) A11=0.D0 CALL SWAPPER(IDIM,VEC0,VEC1,A00,A01,A10,A11) CALL MLTPLY(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC,DIAG, &NDCLRD,VEC1,VEC0,IDIM,LBOND,NBOND,LBASE,MBASE,IFLGDELM, &IFLGOELM,ELMNT,LOC,NDCLRM) LK=LK+1 IF(ISW.EQ.1) THEN DO 10 I=1,NVEC CALL ADDER(IDIM,VEC(0,I-1),VEC(0,I-1),VEC0,1.0D0, & COEF(LK,I)) 10 CONTINUE ENDIF RETURN END C C MAIN OF LANCZOS FOR NEC SX2 C SUBROUTINE MALNCNEC(LLANCZ,LK,ALPHA,BETA,EE,E2P,COEF,NDCLRW, &NVEC,ERRLNCZ,IWK1,WK1,INDRET) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IWK1(NDCLRW) DIMENSION ALPHA(NDCLRW),BETA(NDCLRW) DIMENSION EE(NDCLRW),COEF(NDCLRW,5),WK1(NDCLRW,6) C INDRET=0 EPS=-1.D-10 IP2=2 N=LK-1 M=4 LNV=NDCLRW ISW=-1 IF(LLANCZ.EQ.4)THEN CALL DCSTSN(ALPHA,N,BETA,EPS,EE,M,ISW,WK1,IERR) E2=EE(IP2) IF(IERR.NE.0) THEN WRITE(6,100)IERR INDRET=-1 RETURN END IF END IF IF(LLANCZ.GT.4)THEN CALL DCSTSN(ALPHA,N,BETA,EPS,EE,M,ISW,WK1,IERR) E2=EE(IP2) IF(IERR.NE.0) THEN WRITE(6,100)IERR INDRET=-1 RETURN ENDIF IF(ABS(E2P-E2).LT.MAX(ABS(E2*ERRLNCZ),ERRLNCZ))THEN IF(NVEC.GT.0)THEN CALL DCSTSS(ALPHA,N,BETA,EPS,EE,M,COEF,LNV, & ISW,IWK1,WK1,IERR) ENDIF INDRET=1 IF(IERR.NE.0) THEN WRITE(6,101)IERR INDRET=-1 ENDIF RETURN END IF ENDIF E2P=E2 100 FORMAT(1H , &' #(W19)# ERROR IN DCSTSN CALLED ORIGINALLY FROM LNCZM. ', & 'IERR=',I10) 101 FORMAT(1H , &' #(W20)# ERROR IN DCSTSS CALLED ORIGINALLY FROM LNCZM. ', & 'IERR=',I10) RETURN END C C MAIN OF LANCZOS FOR TITPACK C SUBROUTINE MALNCTIT(LLANCZ,LK,ALPHA,BETA,EE,E2P,COEF,NDCLRW, &NVEC,ERRLNCZ,IWK1,WK1,INDRET) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IWK1(NDCLRW) DIMENSION ALPHA(NDCLRW),BETA(NDCLRW) DIMENSION EE(NDCLRW),COEF(NDCLRW,5),WK1(NDCLRW,6) C A=COEF(1,1) INDRET=0 EPS=1.D-10 IP2=2 M=4 IF(LLANCZ.EQ.4)THEN CALL BISEC(ALPHA,BETA,LK-1,EE,4,EPS) E2=EE(IP2) END IF IF(LLANCZ.GT.4)THEN CALL BISEC(ALPHA,BETA,LK-1,EE,4,EPS) E2=EE(IP2) IF(ABS(E2P-E2).LT.MAX(ABS(E2*ERRLNCZ),ERRLNCZ))THEN IF(NVEC.GT.0)THEN CALL VEC12(EE,LK-1,M,WK1(1,1),WK1(1,2), & WK1(1,3),WK1(1,4),WK1(1,5),IWK1) ENDIF INDRET=1 RETURN END IF ENDIF E2P=E2 RETURN END C C INVERSE ITERATION C SUBROUTINE INVERS2(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,DIAG,NDCLRD,VEC,NDCLRV,MDCLRV,LBOND,NBOND,LBASE,MBASE, &IFLGDELM,IFLGOELM,ELMNT,LOC,NDCLRM,NVEC,ALAMBDA,IDS) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION DIAG(0:NDCLRD),VEC(0:NDCLRV,0:MDCLRV) CHARACTER*4 MSGBC(2),MSG DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION MBASE(0:2*NSLMAX,IBOND+NS,6) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) DIMENSION IDS(0:3) C I=NVEC MSG=MSGBC(1) I0=IDS(0) I1=IDS(1) I2=IDS(2) I3=IDS(3) IDIM=IBASE(NSB) LROT=0 GOTO 11 10 CONTINUE CALL CGMETH2(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC,DIAG, &NDCLRD,VEC,NDCLRV,MDCLRV,LBOND,NBOND,LBASE,MBASE,IFLGDELM, &IFLGOELM,ELMNT,LOC,NDCLRM,NVEC,ALAMBDA,IDS) 11 CALL VMOVE(IDIM,VEC(0,I1),VEC(0,I0)) CALL INPRO(IDIM,VEC(0,I1),VEC(0,I1),DTP) DTP=SQRT(DTP) CALL DIVER(IDIM,VEC(0,I1),VEC(0,I1),DTP) CALL VEC01(IDIM,VEC(0,I2)) CALL MLTPLY(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC,DIAG, &NDCLRD,VEC(0,I2),VEC(0,I1),IDIM,LBOND,NBOND,LBASE,MBASE, &IFLGDELM,IFLGOELM,ELMNT,LOC,NDCLRM) CALL INPRO(IDIM,VEC(0,I2),VEC(0,I2),DTP) DTP=SQRT(DTP) CALL DIVER(IDIM,VEC(0,I2),VEC(0,I2),DTP) C IF(ALAMBDA .LT. 0.D0) THEN CALL ADDER(IDIM,VEC(0,I3),VEC(0,I2),VEC(0,I1), & 1.0D0,1.0D0) ELSE CALL ADDER(IDIM,VEC(0,I3),VEC(0,I2),VEC(0,I1), & -1.0D0,1.0D0) END IF C CALL INPRO(IDIM,VEC(0,I3),VEC(0,I3),DTR) LROT=LROT+1 IF(LROT .GE. 15) THEN WRITE(6,*) & ' #(W21)# DID NOT CONVERGE FOR THE ',I0+1,'-TH VECTOR ', & 'IN INVERS2 CALLED FROM' WRITE(6,*) & ' LNCZMINV.' RETURN END IF IF(SQRT(DTR) .GE. 0.5D-11) GOTO 10 RETURN END C C CG METHOD C SUBROUTINE CGMETH2(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,DIAG,NDCLRD,VEC,NDCLRV,MDCLRV,LBOND,NBOND,LBASE,MBASE, &IFLGDELM,IFLGOELM,ELMNT,LOC,NDCLRM,NVEC,ALAMBDA,IDS) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION DIAG(0:NDCLRD),VEC(0:NDCLRV,0:MDCLRV) CHARACTER*4 MSGBC(2),MSG DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION MBASE(0:2*NSLMAX,IBOND+NS,6) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) DIMENSION IDS(0:3) C I=NVEC A=ALAMBDA MSG=MSGBC(1) I0=IDS(0) I1=IDS(1) I2=IDS(2) I3=IDS(3) IDIM=IBASE(NSB) C (I) CALL INPRO(IDIM,VEC(0,I1),VEC(0,I1),DB) CALL VMOVE(IDIM,VEC(0,I2),VEC(0,I1)) CALL VEC01(IDIM,VEC(0,I0)) C (II) EPSRON=1.0D-8 LROT=0 CALL INPRO(IDIM,VEC(0,I1),VEC(0,I1),DR1) 5 CONTINUE CALL ADDER(IDIM,VEC(0,I3),VEC(0,I3),VEC(0,I2),0.D0,-ALAMBDA) CALL MLTPLY(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC,DIAG, &NDCLRD,VEC(0,I3),VEC(0,I2),IDIM,LBOND,NBOND,LBASE,MBASE, &IFLGDELM,IFLGOELM,ELMNT,LOC,NDCLRM) CALL INPRO(IDIM,VEC(0,I3),VEC(0,I2),DTP) ALPHA=DR1/DTP CALL ADDER(IDIM,VEC(0,I0),VEC(0,I0),VEC(0,I2),1.0D0,ALPHA) CALL ADDER(IDIM,VEC(0,I1),VEC(0,I1),VEC(0,I3),1.0D0,-ALPHA) CALL INPRO(IDIM,VEC(0,I1),VEC(0,I1),DR2) IF(DR2/DB .LE. EPSRON) GOTO 200 BETA=DR2 / DR1 DR1=DR2 LROT= LROT+1 CALL ADDER(IDIM,VEC(0,I2),VEC(0,I1),VEC(0,I2),1.0D0,BETA) IF(LROT.GE. 40) THEN C WRITE(6,132) C 132 FORMAT(1H ,' #(Wxx)# DID NOT CONVERGE IN CGMETH2 ', C & 'CALLED ORIGINALLY FROM LNCZMINV.' GOTO 154 END IF GOTO 5 200 CONTINUE 154 CONTINUE C RETURN END C c****************************************************************** c* From TITPACK Version 2 * c* Copyright (C) Hidetoshi Nishimori * c* Thanks to Professor H. Nishimori * c****************************************************************** c c**** eigenvector of a tridiagonal matrix by inverse iteration **** c for the large/medium routines c subroutine vec12(E,ndim,nvec,di,bl,bu,bv,cm,lex) c c E(4) @ 4 lowest eigenvalues c ndim @ matrix dimension c nvec @ number of vectors to calculate c di - lex working areas c implicit real*8 (a-h,o-z) PARAMETER (ITRMAX=1000,NROTAT=5,ITRLNCZ=ITRMAX/NROTAT) PARAMETER (NDCLRWK=5*(ITRLNCZ+2)) dimension E(4) dimension di(ndim),bl(ndim),bu(ndim),bv(ndim),cm(ndim), &lex(ndim) COMMON/VECDAT/V(NDCLRWK,5),ALPHA(NDCLRWK),BETA(NDCLRWK) c do 10 k=1,nvec c do 100 j=1,ndim di(j)=E(k)-alpha(j) bl(j)=-beta(j) bu(j)=-beta(j) 100 continue c c*** LU decomposition do 110 j=1,ndim-1 if(abs(di(j)).gt.abs(bl(j)))then c--- non pivoting lex(j)=0 if(abs(di(j)).lt.1.d-40)di(j)=1.d-40 cm(j+1)=bl(j)/di(j) di(j+1)=di(j+1)-cm(j+1)*bu(j) bv(j)=0.d0 else c--- pivoting lex(j)=1 cm(j+1)=di(j)/bl(j) di(j)=bl(j) s=bu(j) bu(j)=di(j+1) bv(j)=bu(j+1) di(j+1)=s-cm(j+1)*bu(j) bu(j+1)= -cm(j+1)*bv(j) end if 110 continue if(abs(di(ndim)).lt.1.d-40)di(ndim)=1.d-40 c c--- initial vector do 120 j=1,ndim 120 v(j,k)=1.d0/dble(j*5) c c*** degeneracy check up if(k.eq.1)then km=k else if(abs(E(k)-E(km)).gt.1.d-13)then km=k else do 130 i=km,k-1 prd=0.d0 do 140 j=1,ndim 140 prd=prd+v(j,i)*v(j,k) do 150 j=1,ndim 150 v(j,k)=v(j,k)-prd*v(j,i) 130 continue end if c c*** inverse iteration do 160 l=1,k-km+3 if((l.ne.1).or.(k.ne.km))then c--- forward substitution do 170 j=1,ndim-1 if(lex(j).eq.0)then v(j+1,k)=v(j+1,k)-cm(j+1)*v(j,k) else s=v(j,k) v(j,k)=v(j+1,k) v(j+1,k)=s-cm(j+1)*v(j,k) end if 170 continue end if c--- backward substitution do 180 j=ndim,1,-1 s=v(j,k) if(j.le.ndim-1)s=s-bu(j)*v(j+1,k) if(j.le.ndim-2)s=s-bv(j)*v(j+2,k) v(j,k)=s/di(j) 180 continue c c*** normalization dnorm=0.d0 do 190 j=1,ndim 190 dnorm=dnorm+v(j,k)**2 if(dnorm.gt.1.d-50)dnorm=1./sqrt(dnorm) do 200 j=1,ndim 200 v(j,k)=v(j,k)*dnorm 160 continue c 10 continue return end C C****************************************************************** C MULTIPLICATION OF H BY VEC FOR MID-SIZE C VEC(I,ID)=VEC(I,ID)+SUMJ(H(I,J)*VEC(J,IS)) C SUBROUTINE MLTPLY(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,DIAG,NDCLRD,VD,VS,IDIM,LBOND,NBOND,LBASE,MBASE, &IFLGDELM,IFLGOELM,ELMNT,LOC,NDCLRM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION DIAG(0:NDCLRD),VD(0:IDIM-1),VS(0:IDIM-1) CHARACTER*4 MSGBC(2) DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION MBASE(0:2*NSLMAX,IBOND+NS,6) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) C C MULTIPLICATION OF DIAG BY VEC C VEC(ID)=VEC(ID)+DIAG*VEC(IS) C CALL MLTD(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC,DIAG, &NDCLRD,VD,VS,IDIM,LBOND,NBOND,LBASE,IFLGDELM) C CALL MLTO(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VD,VS,IDIM,LBOND,NBOND, &MBASE,IFLGOELM,ELMNT,LOC,NDCLRM) C RETURN END C C MULTIPLICATION OF DAIG BY VEC C VD=DIAG*VS C SUBROUTINE MLTD(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,DIAG,NDCLRD,VD,VS,IDIM,LBOND,NBOND,LBASE,IFLGDELM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION DIAG(0:NDCLRD),VD(0:IDIM-1),VS(0:IDIM-1) CHARACTER*4 MSGBC(2),MSG DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) C MSG=MSGBC(1) IDIM=IBASE(NSB) C CALL VEC01(IDIM,VD) C CALL MLTDLL(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC,DIAG, &NDCLRD,VD,VS,IDIM,LBOND,NBOND,LBASE,IFLGDELM) IF(NS.GE.NSLMAX+1)THEN CALL MLTDLH(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, & NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND, & DIAG,NDCLRD,VD,VS,IDIM,LBOND,NBOND,LBASE,IFLGDELM) C CALL MLTDHH(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, & NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, & MSGBC, & DIAG,NDCLRD,VD,VS,IDIM,LBOND,NBOND,LBASE,IFLGDELM) ENDIF C RETURN END C C MULTIPLICATION OF DIAGONAL ELEMENT FOR LL C SUBROUTINE MLTDLL(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,DIAG,NDCLRD,VD,VS,IDIM,LBOND,NBOND,LBASE,IFLGDELM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION DIAG(0:NDCLRD),VD(0:IDIM-1),VS(0:IDIM-1) CHARACTER*4 MSGBC(2),MSG DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) C MSG=MSGBC(1) DO 100 II=0,NSB-1 I=II ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) LBS1=LBASE(I,1,1) LBS2=LBASE(I,2,1) IF(NS.GE.NSLMAX+1.AND.IFLGDELM(I,2).GT.0 & .AND.IFLGDELM(I,1).GT.0)THEN CALL MLTDEL(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, & NDCLR1,ISTAT2,DIAG(LBS1),LL,DIAG(LBS2),LH,VD,VS,IDIM) ELSEIF(IFLGDELM(I,1).GT.0)THEN CALL MLTDLLEL(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, & ISTAT1,NDCLR1,ISTAT2,DIAG(LBS1),LL,VD,VS,IDIM) ELSE CALL MLTDLLBX(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, & ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND, & ANISTRPY,HFIELD,MSGBC,VD,VS,IDIM, & LBOND,NBOND) ENDIF 100 CONTINUE RETURN END C C MULTIPLICATION OF DIAGONAL ELEMENT FOR LL AND HH C SUBROUTINE MLTDEL(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,DIAGL,LML,DIAGH,LMH,VD,VS,IDIM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION DIAGL(0:LML-1),DIAGH(0:LMH-1) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) C J=NS J=NSB J=ISTAT2(1) J=ISTAT1(1,1) ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) C DO 30 J=0,LH-1 IX=IBX+J*LL *VDIR LOOPCNT=73789 DO 20 K=0,LL -1 VD(IX+K)=VD(IX+K) & +(DIAGH(J)+DIAGL(K))*VS(IX+K) 20 CONTINUE 30 CONTINUE C RETURN END C C MULTIPLICATION OF DIAGONAL ELEMENT FOR LL C SUBROUTINE MLTDLLEL(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,DIAG,LM,VD,VS,IDIM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION DIAG(0:LM-1) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) C J=NS J=NSB J=ISTAT2(1) J=ISTAT1(1,1) ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) C DO 30 J=0,LH-1 IX=IBX+LL*J *VDIR LOOPCNT=73789 DO 20 K=0,LL -1 VD(IX+K)=VD(IX+K) & +DIAG(K)*VS(IX+K) 20 CONTINUE 30 CONTINUE C RETURN END C C MULTIPLICATION OF DIAGONAL ELEMENT FOR LL FOR BIQD.NE.0 C SUBROUTINE MLTDLLBX(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY, &HFIELD,MSGBC,VD,VS,IDIM,LBOND,NBOND) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) CHARACTER*4 MSGBC(2) DIMENSION LBOND(IBOND+6,0:5),NBOND(5) C MSGBC(1)=MSGBC(1) J=NS J=NSB J=ISTAT2(1) A=AJX(1) ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) C DO 90 KIS=1,MIN(NS,NSLMAX) IST=3**(KIS-1) AN=ANISTRPY(KIS) HF=HFIELD(KIS) IF(ABS(ANISTRPY(KIS)).LE.1.D-60 & .AND.ABS(HFIELD(KIS)).LE.1.D-60) THEN GOTO 90 ENDIF DO 80 J=0,LH-1 IX=IBX+LL*J *VDIR LOOPCNT=73789 DO 20 K=0,LL -1 ITP=ISTAT1(K,NL(I)) VD(IX+K)=VD(IX+K) & -(HF*DBLE(MOD(ITP/IST,3)-1) & -AN*DBLE((MOD(ITP/IST,3)-1)**2))*VS(IX+K) 20 CONTINUE C 80 CONTINUE 90 CONTINUE C DO 140 KS=1,NBOND(1) KIS=LBOND(KS,1) ISITE1=(IPAIR(KIS*2-1)-1) IS1=3**ISITE1 ISITE2=(IPAIR(KIS*2 )-1) IS2=3**ISITE2 XJZ=AJZ(KIS) XB=ABIQD(KIS) XBXX=ABIQD(KIS) XBZZ=ABIQD(KIS) IF(ABS(AJZ(KIS)).LE.1.D-60 & .AND.ABS(ABIQD(KIS)).LE.1.D-60) THEN GOTO 140 ENDIF IF(ABS(ABIQD(KIS)).GT.1.D-60) THEN DO 110 J=0,LH -1 IX=IBX+LL*J *VDIR LOOPCNT=73789 DO 100 K=0,LL -1 ITP=ISTAT1(K,NL(I)) IST1=(MOD(ITP/IS1,3)-1) IST2=(MOD(ITP/IS2,3)-1) VD(IX+K)=VD(IX+K) & +(XJZ*DBLE(IST1*IST2) & +XB*DBLE((IST1*IST2)**2 & +1+(IST1*IST2+1)*(IST1*IST2+2)/2 & -(IST1+IST2)**2))*VS(IX+K) CCC & +XBZZ*DBLE((IST1*IST2)**2) CCC & +XBXX*DBLE(1+(IST1*IST2+1)*(IST1*IST2+2)/2 CCC & -(IST1+IST2)**2) 100 CONTINUE 110 CONTINUE ELSE DO 130 J=0,LH -1 IX=IBX+LL*J *VDIR LOOPCNT=73789 DO 120 K=0,LL -1 ITP=ISTAT1(K,NL(I)) IST1=(MOD(ITP/IS1,3)-1) IST2=(MOD(ITP/IS2,3)-1) VD(IX+K)=VD(IX+K) & +XJZ*DBLE(IST1*IST2)*VS(IX+K) 120 CONTINUE 130 CONTINUE ENDIF 140 CONTINUE C RETURN END C C MULTIPLICATION OF DIAGONAL ELEMENT FOR HH C SUBROUTINE MLTDHH(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,DIAG,NDCLRD,VD,VS,IDIM,LBOND,NBOND,LBASE,IFLGDELM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION DIAG(0:NDCLRD),VD(0:IDIM-1),VS(0:IDIM-1) CHARACTER*4 MSGBC(2) DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) C DO 100 II=0,NSB-1 I=II ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) LBS2=LBASE(I,2,1) IF(NS.GE.NSLMAX+1.AND.IFLGDELM(I,2).GT.0 & .AND.IFLGDELM(I,1).GT.0)THEN GOTO 100 ELSEIF(IFLGDELM(I,2).GT.0)THEN CALL MLTDHHEL(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, & ISTAT1,NDCLR1,ISTAT2,DIAG(LBS2),LH,VD,VS,IDIM) ELSE CALL MLTDHHBX(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, & ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND, & ANISTRPY,HFIELD,MSGBC,VD,VS,IDIM, & LBOND,NBOND) ENDIF 100 CONTINUE RETURN END C C MULTIPLICATION OF DIAGONAL ELEMENT FOR LL C SUBROUTINE MLTDHHEL(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,DIAG,LM,VD,VS,IDIM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION DIAG(0:LM-1),VD(0:IDIM-1),VS(0:IDIM-1) C J=NS J=NSB J=ISTAT2(1) J=ISTAT1(1,1) ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) C DO 30 J=0,LH-1 IX=IBX+J*LL *VDIR LOOPCNT=73789 DO 20 K=0,LL -1 VD(IX+K)=VD(IX+K) & +DIAG(J)*VS(IX+K) 20 CONTINUE 30 CONTINUE C RETURN END C C MULTIPLICATION OF DIAGONAL ELEMENT FOR LL FOR BIQD.NE.0 C SUBROUTINE MLTDHHBX(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY, &HFIELD,MSGBC,VD,VS,IDIM,LBOND,NBOND) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) CHARACTER*4 MSGBC(2) DIMENSION LBOND(IBOND+6,0:5),NBOND(5) C MSGBC(1)=MSGBC(1) J=NS J=NSB J=IPAIR(1) J=ISTAT2(1) J=LBOND(1,1) J=NBOND(1) J=ISTAT1(1,1) A=AJX(1) ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) C DO 40 KIS=MIN(NS,NSLMAX)+1,NS IST=3**(KIS-1-NSLMAX) AN=ANISTRPY(KIS) HF=HFIELD(KIS) IF(ABS(ANISTRPY(KIS)).LE.1.D-60 & .AND.ABS(HFIELD(KIS)).LE.1.D-60) THEN GOTO 40 ENDIF DO 30 J=0,LH -1 IX=IBX+J*LL ITP=ISTAT1(J,NH(I)) DGH=-HF*DBLE(MOD(ITP/IST,3)-1) & +AN*DBLE((MOD(ITP/IST,3)-1)**2) *VDIR LOOPCNT=73789 DO 20 K=0,LL-1 VD(IX+K)=VD(IX+K)+DGH*VS(IX+K) 20 CONTINUE 30 CONTINUE C 40 CONTINUE C DO 110 KS=1,NBOND(2) KIS=LBOND(KS,2) ISITE1=(IPAIR(KIS*2-1)-1) IS1=3**(ISITE1-NSLMAX) ISITE2=(IPAIR(KIS*2 )-1) IS2=3**(ISITE2-NSLMAX) XJZ=AJZ(KIS) XB=ABIQD(KIS) XBXX=ABIQD(KIS) XBZZ=ABIQD(KIS) IF(ABS(AJZ(KIS)).LE.1.D-60 & .AND.ABS(ABIQD(KIS)).LE.1.D-60) THEN GOTO 110 ENDIF IF(ABS(ABIQD(KIS)).GT.1.D-60) THEN DO 80 J=0,LH -1 IX=IBX+J*LL ITP=ISTAT1(J,NH(I)) IST1=(MOD(ITP/IS1,3)-1) IST2=(MOD(ITP/IS2,3)-1) DGH=XJZ*DBLE(IST1*IST2) & +XB*DBLE((IST1*IST2)**2 & +1+(IST1*IST2+1)*(IST1*IST2+2)/2 & -(IST1+IST2)**2) CCC & +XBZZ*DBLE((IST1*IST2)**2) CCC & +XBXX*DBLE(1+(IST1*IST2+1)*(IST1*IST2+2)/2 CCC & -(IST1+IST2)**2) *VDIR LOOPCNT=73789 DO 70 K=0,LL-1 VD(IX+K)=VD(IX+K)+DGH*VS(IX+K) 70 CONTINUE 80 CONTINUE ELSE DO 100 J=0,LH -1 IX=IBX+J*LL ITP=ISTAT1(J,NH(I)) IST1=(MOD(ITP/IS1,3)-1) IST2=(MOD(ITP/IS2,3)-1) DGH=XJZ*DBLE(IST1*IST2) *VDIR LOOPCNT=73789 DO 90 K=0,LL-1 VD(IX+K)=VD(IX+K)+DGH*VS(IX+K) 90 CONTINUE 100 CONTINUE ENDIF 110 CONTINUE C RETURN END C C MULTIPLICATION OF DIAGONAL H(LOW,HIGH) BY VEC C VEC(ID)=H(LOW,HIGH)*VEC(IS) C SUBROUTINE MLTDLH(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,DIAG,NDCLRD,VD,VS, &IDIM,LBOND,NBOND,LBASE,IFLGDELM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION DIAG(0:NDCLRD),VD(0:IDIM-1),VS(0:IDIM-1) DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) C C L-H BOND DO 100 KS=1,NBOND(3) KIS=LBOND(KS,3) DO 90 II=0,NSB-1 I=II ITH=NH(I) ITL=NL(I) LH=NSTAT(ITH,1) LL=NSTAT(ITL,0) IF(IFLGDELM(I,2+KS).GT.0)THEN LBS3=LBASE(I,2+KS,1) LBS4=LBASE(I,2+KS,2) CALL MLTDLHEL(I,KIS,NS,NSLMAX,NSB,NH,NL,IBASE, & NSTAT,ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD, & IBOND, & DIAG(LBS3),LL,DIAG(LBS4),LH, & VD,VS,IDIM) ELSE CALL MLTDLHBX(I,KIS,NS,NSLMAX,NSB,NH,NL,IBASE, & NSTAT,ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD, & IBOND, & VD,VS,IDIM) ENDIF 90 CONTINUE 100 CONTINUE RETURN END C C MULTIPLICATION OF LH ELEMENT C SUBROUTINE MLTDLHEL(I,KBOND,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,DIAGL,LML, &DIAGH,LMH,VD,VS,IDIM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) DIMENSION DIAGL(0:LML-1),DIAGH(0:LMH-1) C J=NS J=NSB J=ISTAT2(1) J=ISTAT1(1,1) J=IPAIR(2*KBOND) A=AJX(1) C KIS=KBOND ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) XJZ=AJZ(KIS) XB=ABIQD(KIS) IF(ABS(XJZ).LE.1.D-60.AND.ABS(XB).LE.1.D-60) THEN RETURN ENDIF IF(ABS(XB).GT.1.D-60) THEN XBXX=XB XBZZ=XB DO 50 J=0,LH -1 IX=IBASE(I)+J*LL *VDIR LOOPCNT=73789 DO 40 K=0,LL -1 ST1=DIAGL(K) VD(IX+K)=VD(IX+K) & +(XJZ*DIAGL(K)*DIAGH(J) & +XB*((DIAGL(K)*DIAGH(J))**2 & +1.D0+(DIAGL(K)*DIAGH(J)+1.D0)* & (DIAGL(K)*DIAGH(J)+2.D0)/2.D0 & -(DIAGL(K)+DIAGH(J))**2))*VS(IX+K) CCC & +XBZZ*((DIAGL(K)*DIAGH(J))**2) CCC & +XBXX*(1.D0+(DIAGL(K)*DIAGH(J)+1.D0) CCC & *(DIAGL(K)*DIAGH(J)+2.D0)/2.D0 CCC & -(DIAGL(K)+DIAGH(J))**2)*VS(IX+K) 40 CONTINUE 50 CONTINUE ELSE DO 70 J=0,LH -1 IX=IBASE(I)+J*LL *VDIR LOOPCNT=73789 DO 60 K=0,LL -1 VD(IX+K)=VD(IX+K)+ & XJZ*DIAGL(K)*DIAGH(J)*VS(IX+K) 60 CONTINUE 70 CONTINUE ENDIF C RETURN END C C MULTIPLICATION OF LH ELEMENT FOR BIQD.NE.0 C SUBROUTINE MLTDLHBX(I,KBOND,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VD,VS,IDIM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) C J=NS J=NSB J=ISTAT2(1) J=ISTAT1(1,1) A=AJX(1) C KIS=KBOND ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) ISITE1=(IPAIR(KIS*2-1)-1) ISITE2=(IPAIR(KIS*2 )-1) IF((ISITE1+1.LE.NSLMAX).AND.(ISITE2+1.GT.NSLMAX))THEN MSKL=3**ISITE1 MSKH=3**(ISITE2-NSLMAX) ELSEIF((ISITE1+1.GT.NSLMAX).AND.(ISITE2+1.LE.NSLMAX))THEN MSKH=3**(ISITE1-NSLMAX) MSKL=3**ISITE2 ELSE GOTO 100 ENDIF XJZ=AJZ(KIS) XB=ABIQD(KIS) IF(ABS(XJZ).LE.1.D-60.AND.ABS(XB).LE.1.D-60) THEN GOTO 100 ENDIF IF(ABS(XB).GT.1.D-60)THEN XBXZ=XB XBXX=XB DO 70 J=0,LH-1 IX=IBX+J*LL ITPH=ISTAT1(J,ITH) IST2=MOD(ITPH/MSKH,3)-1 *VDIR NODEP *VDIR LOOPCNT=73789 DO 60 K=0,LL - 1 ITPL=ISTAT1(K,ITL) IST1=MOD(ITPL/MSKL,3)-1 VD(IX+K)=VD(IX+K) & +(XJZ*DBLE(IST1*IST2) & +XB*DBLE((IST1*IST2)**2 & +1+(IST1*IST2+1)*(IST1*IST2+2)/2 & -(IST1+IST2)**2))*VS(IX+K) CCC & +XBZZ*DBLE((IST1*IST2)**2) CCC & +XBXX*DBLE(1+(IST1*IST2+1)*(IST1*IST2+2)/2 CCC & -(IST1+IST2)**2) 60 CONTINUE 70 CONTINUE ELSE DO 90 J=0,LH-1 IX=IBX+J*LL ITPH=ISTAT1(J,ITH) IST2=MOD(ITPH/MSKH,3)-1 *VDIR NODEP *VDIR LOOPCNT=73789 DO 80 K=0,LL - 1 ITPL=ISTAT1(K,ITL) IST1=MOD(ITPL/MSKL,3)-1 VD(IX+K)=VD(IX+K) & +XJZ*DBLE(IST1*IST2)*VS(IX+K) 80 CONTINUE 90 CONTINUE ENDIF C 100 CONTINUE C RETURN END C C MULTIPLICATION OF H(OFFDIAG) BY VEC C VEC(I,ID)=SUMJ(HOFFDIAG(I,J)*VEC(J,IS)) C SUBROUTINE MLTO(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VD,VS,IDIM,LBOND, &NBOND,MBASE,IFLGOELM,ELMNT,LOC,NDCLRM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) DIMENSION LBOND(IBOND+6,0:5),NBOND(5) DIMENSION MBASE(0:2*NSLMAX,IBOND+NS,6) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) C CALL MLTOLL(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VD,VS,IDIM,LBOND,NBOND, &MBASE,IFLGOELM,ELMNT,LOC,NDCLRM) C IF(NS.GE.NSLMAX+1 ) THEN CALL MLTOLH(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, & NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VD,VS,IDIM, & LBOND,NBOND,MBASE,IFLGOELM, & ELMNT,LOC,NDCLRM) END IF C IF(NS.GE.NSLMAX+2) THEN CALL MLTOHH(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, & NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VD,VS,IDIM, & LBOND,NBOND,MBASE,IFLGOELM, & ELMNT,LOC,NDCLRM) END IF C RETURN END C C MULTIPLICATION OF OFF-DIAGONAL H(LOW,LOW) BY VEC C VEC(ID)=H'(LOW,LOW)*VEC(IS) C SUBROUTINE MLTOLL(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VD,VS,IDIM,LBOND, &NBOND,MBASE,IFLGOELM,ELMNT,LOC,NDCLRM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) DIMENSION LBOND(IBOND+6,0:5),NBOND(5) DIMENSION MBASE(0:2*NSLMAX,IBOND+NS,6) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) C C L-L PAIR DO 100 KS=1,NBOND(1) KIS=LBOND(KS,1) DO 100 II=0,NSB-1 I=II ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) MBS1E=MBASE(I,KIS,1) MBS1L=MBASE(I,KIS,2) IC1=MBASE(I,KIS,3) IF(IC1.GT.0)THEN IF(IFLGOELM(I,KIS).GT.0)THEN CALL MLTOLLEL(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, & ISTAT1,NDCLR1,ISTAT2,VD,VS,IDIM, & ELMNT(MBS1E),LOC(MBS1L),LL,IC1) ELSE CALL MLTOLLBX(I,KIS,NS,NSLMAX,NSB,NH,NL,IBASE, & NSTAT, & ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND, & VD,VS,IDIM) ENDIF ENDIF 100 CONTINUE RETURN END C C MULTIPLICATION OF OFF-DIAGONAL LL ELEMENT C SUBROUTINE MLTOLLEL(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,VD,VS,IDIM,ELMNT,LOC,LM,ICC) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) DIMENSION ELMNT(0:LM-1,ICC),LOC(0:LM-1,ICC) C J=NS J=NSB J=ISTAT2(1) J=ISTAT1(1,1) ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) DO 80 KS=1,ICC DO 70 J=0,LH-1 IX=IBASE(I)+LL*J *VDIR NODEP *VDIR LOOPCNT=73789 DO 60 K=0,LL -1 IF(LOC(K,KS).NE.K)THEN IDP=LOC(K,KS) VD(IX+K)=VD(IX+K) & +ELMNT(K,KS)*VS(IX+IDP) VD(IX+IDP)=VD(IX+IDP) & +ELMNT(K,KS)*VS(IX+ K) ENDIF 60 CONTINUE 70 CONTINUE 80 CONTINUE C RETURN END C C C SUBROUTINE MLTOLLBX(I,KIS,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VD,VS,IDIM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) C J=NS J=NSB J=ISTAT2(1) A=AJZ(1) ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) C ISITE1=(IPAIR(KIS*2-1)-1) ISITE2=(IPAIR(KIS*2 )-1) IF((ISITE1+1.LE.NSLMAX).AND.(ISITE2+1.LE.NSLMAX))THEN MSKH=3**ISITE1 MSKL=3**ISITE2 ELSE GOTO 100 ENDIF XJX=AJX(KIS) XB=ABIQD(KIS) IF(ABS(XJX).LE.1.D-60.AND.ABS(XB).LE.1.D-60) THEN GOTO 100 ENDIF IF(ABS(XB).GT.1.D-60) THEN XBXX=XB XBXZ=XB DO 70 J=0,LH-1 IX=IBX+LL*J *VDIR NODEP *VDIR LOOPCNT=73789 DO 60 K=0,LL -1 ITP=ISTAT1(K,ITL) IF(MOD(ITP/MSKH,3).NE.0.AND.MOD(ITP/MSKL,3).NE.2) & THEN XBJ=XJX-XBXZ*DBLE(1- & (MOD(ITP/MSKH,3)+MOD(ITP/MSKL,3)-2)**2) IDP=ISTAT2(ITP - MSKH+MSKL) VD(IX+K)=VD(IX+K)+XBJ*VS(IX+IDP) VD(IX+IDP)=VD(IX+IDP)+XBJ*VS(IX+K) END IF IF(MOD(ITP/MSKH,3).EQ.2.AND.MOD(ITP/MSKL,3).EQ.0) & THEN IDP=ISTAT2(ITP -2*MSKH+2*MSKL) VD(IX+K)=VD(IX+K)+XBXX*VS(IX+IDP) VD(IX+IDP)=VD(IX+IDP)+XBXX*VS(IX+K) END IF 60 CONTINUE 70 CONTINUE C ELSE DO 90 J=0,LH-1 IX=IBX+LL*J *VDIR NODEP *VDIR LOOPCNT=73789 DO 80 K=0,LL -1 ITP=ISTAT1(K,ITL) IF(MOD(ITP/MSKH,3).NE.0.AND.MOD(ITP/MSKL,3).NE.2) & THEN IDP=ISTAT2(ITP - MSKH+MSKL) VD(IX+K)=VD(IX+K)+XJX*VS(IX+IDP) VD(IX+IDP)=VD(IX+IDP)+XJX*VS(IX+K) END IF 80 CONTINUE 90 CONTINUE ENDIF C 100 CONTINUE C RETURN END C C MULTIPLICATION OF OFF-DIAGONAL H(HIGH,HIGH) BY VEC C VEC(ID)=H'(HIGH,HIGH)*VEC(IS) C SUBROUTINE MLTOHH(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VD,VS,IDIM,LBOND, &NBOND,MBASE,IFLGOELM,ELMNT,LOC,NDCLRM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) DIMENSION LBOND(IBOND+6,0:5),NBOND(5),MBASE(0:2*NSLMAX, &IBOND+NS,6) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) C DO 100 KS=1,NBOND(2) KIS=LBOND(KS,2) DO 100 II=0,NSB-1 I=II ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) MBS2E=MBASE(I,KIS,1) MBS2L=MBASE(I,KIS,2) IC2=MBASE(I,KIS,3) IF(IC2.GT.0)THEN IF(IFLGOELM(I,KIS).GT.0)THEN CALL MLTOHHEL(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, & ISTAT1,NDCLR1,ISTAT2,VD,VS,IDIM, & ELMNT(MBS2E),LOC(MBS2L),LH,IC2) ELSE CALL MLTOHHBX(I,KIS,NS,NSLMAX,NSB,NH,NL,IBASE, & NSTAT, & ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND, & VD,VS,IDIM) ENDIF ENDIF 100 CONTINUE RETURN END C C MULTIPLICATION OF OFF-DIAGONAL HH ELEMENT C SUBROUTINE MLTOHHEL(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,VD,VS,IDIM,ELMNT,LOC,LM,ICC) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) DIMENSION ELMNT(0:LM-1,ICC),LOC(0:LM-1,ICC) C J=NS J=NSB J=ISTAT2(1) J=ISTAT1(1,1) ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) DO 80 KS=1,ICC DO 70 J=0,LH-1 IF(LOC(J,KS).NE.J)THEN IX=IBX+LL*J IY=IBX+LL*LOC(J,KS) *VDIR NODEP *VDIR LOOPCNT=73789 DO 60 K=0,LL - 1 VD(IY+K)=VD(IY+K) & +ELMNT(J,KS)*VS(IX+K) VD(IX+K)=VD(IX+K) & +ELMNT(J,KS)*VS(IY+K) 60 CONTINUE ENDIF 70 CONTINUE 80 CONTINUE C RETURN END C C C SUBROUTINE MLTOHHBX(I,KIS,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VD,VS,IDIM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) C J=NS J=NSB J=ISTAT2(1) J=ISTAT1(1,1) A=AJZ(1) ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) C ISITE1=(IPAIR(KIS*2-1)-1) ISITE2=(IPAIR(KIS*2 )-1) IF((ISITE1+1.GT.NSLMAX).AND.(ISITE2+1.GT.NSLMAX))THEN MSKH=3**(ISITE1-NSLMAX) MSKL=3**(ISITE2-NSLMAX) ELSE GOTO 100 ENDIF XJX=AJX(KIS) XB=ABIQD(KIS) IF(ABS(XJX).LE.1.D-60.AND.ABS(XB).LE.1.D-60) THEN GOTO 100 ENDIF IF(ABS(XB).GT.1.D-60) THEN XBXZ=XB XBXX=XB DO 70 J=0,LH - 1 ITP=ISTAT1(J,ITH) IF(MOD(ITP/MSKH,3).NE.0 .AND. MOD(ITP/MSKL,3).NE.2) & THEN XBJ=XJX-XBXZ*DBLE(1- & (MOD(ITP/MSKH,3)+MOD(ITP/MSKL,3)-2)**2) IX=IBX+LL*J IY=IBX+LL*ISTAT2(ITP - MSKH+MSKL) *VDIR NODEP *VDIR LOOPCNT=73789 DO 60 K=0,LL - 1 VD(IY+K)=VD(IY+K)+XBJ*VS(IX+K) VD(IX+K)=VD(IX+K)+XBJ*VS(IY+K) 60 CONTINUE END IF IF(MOD(ITP/MSKH,3).EQ.2.AND.MOD(ITP/MSKL,3).EQ.0)THEN IX=IBX+LL*J IY=IBX+LL*ISTAT2(ITP-2*MSKH+2*MSKL) *VDIR NODEP *VDIR LOOPCNT=73789 DO 65 K=0,LL - 1 VD(IY+K)=VD(IY+K)+XBXX*VS(IX+K) VD(IX+K)=VD(IX+K)+XBXX*VS(IY+K) 65 CONTINUE END IF 70 CONTINUE C ELSE DO 90 J=0,LH - 1 ITP=ISTAT1(J,ITH) IF(MOD(ITP/MSKH,3).NE.0 .AND. MOD(ITP/MSKL,3).NE.2) & THEN IX=IBX+LL*J IY=IBX+LL*ISTAT2(ITP-MSKH+MSKL) *VDIR NODEP *VDIR LOOPCNT=73789 DO 80 K=0,LL - 1 VD(IY+K)=VD(IY+K)+XJX*VS(IX+K) VD(IX+K)=VD(IX+K)+XJX*VS(IY+K) 80 CONTINUE END IF 90 CONTINUE ENDIF C 100 CONTINUE C RETURN END C C MULTIPLICATION OF OFF-DIAGONAL H(HIGH,HIGH) BY VEC C VEC(ID)=H'(HIGH,HIGH)*VEC(IS) C SUBROUTINE MLTOLH(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VD,VS,IDIM,LBOND, &NBOND,MBASE,IFLGOELM,ELMNT,LOC,NDCLRM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) DIMENSION LBOND(IBOND+6,0:5),NBOND(5) DIMENSION MBASE(0:2*NSLMAX,IBOND+NS,6) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) C C L-H BOND C DO 100 KS=1,NBOND(3) KIS=LBOND(KS,3) DO 90 II=0,NSB-1 I=II ITH=NH(I) ITL=NL(I) LH=NSTAT(ITH,1) LL=NSTAT(ITL,0) IF(IFLGOELM(I,KIS).GT.0)THEN MBS3E=MBASE(I,KIS,1) MBS3L=MBASE(I,KIS,2) IC3=MBASE(I,KIS,3) MBS4E=MBASE(I,KIS,4) MBS4L=MBASE(I,KIS,5) IC4=MBASE(I,KIS,6) IC=IC3 CALL MLTOLHEL(I,KIS,NS,NSLMAX,NSB,NH,NL,IBASE, & NSTAT,ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD, & IBOND, & VD,VS,IDIM,ELMNT(MBS3E),LOC(MBS3L),LL, & ELMNT(MBS4E),LOC(MBS4L),LH,IC) ELSE CALL MLTOLHBX(I,KIS,NS,NSLMAX,NSB,NH,NL,IBASE, & NSTAT,ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD, & IBOND, & VD,VS,IDIM) ENDIF 90 CONTINUE 100 CONTINUE RETURN END C C MULTIPLICATION OF LH ELEMENT FOR BIQD.NE.0 C SUBROUTINE MLTOLHEL(I,KBOND,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VD,VS,IDIM, &ELMNTL,LOCL,LML,ELMNTH,LOCH,LMH,ICC) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) DIMENSION ELMNTL(0:LML-1,ICC),LOCL(0:LML-1,ICC) DIMENSION ELMNTH(0:LMH-1,ICC),LOCH(0:LMH-1,ICC) C J=NS J=NSB J=IPAIR(1) J=ISTAT2(1) J=ISTAT1(1,1) A=AJZ(1) C KIS=KBOND ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) XJX=AJX(KIS) XB=ABIQD(KIS) IF(ABS(XJX).LE.1.D-60.AND.ABS(XB).LE.1.D-60) THEN RETURN ENDIF XBXZ=XB XBXX=XB IF(ABS(XB).GT.1.D-60.AND.ICC.GE.4)THEN IF(I.LE.NSB-2)THEN IBX1=IBASE(I+1) LL1=NSTAT(ITL+1,0) DO 70 J=0,LH-1 IX=IBX+J*LL IF(LOCH(J,1).GE.0)THEN IY=IBX1+LL1*LOCH(J,1) *VDIR NODEP *VDIR LOOPCNT=73789 DO 60 K=0,LL - 1 IF(LOCL(K,1).GE.0)THEN IDP=LOCL(K,1) VD(IY+IDP)=VD(IY+IDP)+ & (XJX*ELMNTL(K,1)*ELMNTH(J,1) & +XBXZ*(ELMNTL(K,2)*ELMNTH(J,2)+ & ELMNTL(K,3)*ELMNTH(J,3)))*VS(IX+K) VD(IX+K)=VD(IX+K)+ & (XJX*ELMNTL(K,1)*ELMNTH(J,1) & +XBXZ*(ELMNTL(K,2)*ELMNTH(J,2)+ & ELMNTL(K,3)*ELMNTH(J,3)))*VS(IY+IDP) ENDIF 60 CONTINUE ENDIF 70 CONTINUE ENDIF IF(I.LE.NSB-3)THEN IBX2=IBASE(I+2) LL2=NSTAT(ITL+2,0) DO 90 J=0,LH-1 IF(LOCH(J,4).GE.0)THEN IX=IBX+J*LL IY=IBX2+LL2*LOCH(J,4) *VDIR NODEP *VDIR LOOPCNT=73789 DO 80 K=0,LL - 1 IF(LOCL(K,4).GE.0)THEN IDP=LOCL(K,4) VD(IY+IDP)=VD(IY+IDP)+ & XBXX*ELMNTL(K,4)*ELMNTH(J,4)*VS(IX+K) VD(IX+K)=VD(IX+K)+ & XBXX*ELMNTL(K,4)*ELMNTH(J,4)*VS(IY+IDP) ENDIF 80 CONTINUE ENDIF 90 CONTINUE ENDIF C ELSEIF(ABS(XB).LE.1.D-60.AND.ABS(XJX).GT.1.D- &60.AND.ICC.GE.1)THEN IF(I.LE.NSB-2)THEN IBX1=IBASE(I+1) LL1=NSTAT(ITL+1,0) DO 110 J=0,LH-1 IF(LOCH(J,1).GE.0)THEN IX=IBX+J*LL IY=IBX1+LL1*LOCH(J,1) *VDIR NODEP *VDIR LOOPCNT=73789 DO 100 K=0,LL - 1 IF(LOCL(K,1).GE.0)THEN IDP=LOCL(K,1) VD(IY+IDP)=VD(IY+IDP)+ & XJX*ELMNTL(K,1)*ELMNTH(J,1)*VS(IX+K) VD(IX+K)=VD(IX+K)+ & XJX*ELMNTL(K,1)*ELMNTH(J,1)*VS(IY+IDP) ENDIF 100 CONTINUE ENDIF 110 CONTINUE ENDIF ENDIF C RETURN END C C MULTIPLICATION OF LH ELEMENT FOR BIQD.NE.0 C SUBROUTINE MLTOLHBX(I,KBOND,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VD,VS,IDIM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) C J=NS J=NSB J=ISTAT1(1,1) A=AJZ(1) C KIS=KBOND ISITE1=(IPAIR(KIS*2-1)-1) ISITE2=(IPAIR(KIS*2 )-1) IF((ISITE1+1.LE.NSLMAX).AND.(ISITE2+1.GT.NSLMAX))THEN MSKL=3**ISITE1 MSKH=3**(ISITE2-NSLMAX) ELSEIF((ISITE1+1.GT.NSLMAX).AND.(ISITE2+1.LE.NSLMAX))THEN MSKH=3**(ISITE1-NSLMAX) MSKL=3**ISITE2 ELSE GOTO 100 ENDIF XJX=AJX(KIS) XB=ABIQD(KIS) XBXX=ABIQD(KIS) XBXZ=ABIQD(KIS) IF(ABS(XJX).LE.1.D-60.AND.ABS(XB).LE.1.D-60) THEN GOTO 100 ENDIF ITH=NH(I) ITL=NL(I) LH=NSTAT(ITH,1) LL=NSTAT(ITL,0) IBX=IBASE(I) IF(ABS(XB).GT.1.D-60)THEN DO 70 J=0,LH-1 IX=IBX+J*LL ITPH=ISTAT1(J,ITH) IF(MOD(ITPH/MSKH,3) .NE. 0) THEN IBX1=IBASE(I+1) LL1=NSTAT(ITL+1,0) IY=IBX1+LL1*ISTAT2(ITPH-MSKH) *VDIR NODEP *VDIR LOOPCNT=73789 DO 60 K=0,LL - 1 ITPL=ISTAT1(K,ITL) IF(MOD(ITPL/MSKL,3).NE. 2)THEN XBJ=XJX-XBXZ*DBLE(1- & (MOD(ITPH/MSKH,3)+MOD(ITPL/MSKL,3)-2)**2) IDP=ISTAT2(ITPL+MSKL) VD(IY+IDP)=VD(IY+IDP)+XBJ*VS(IX+K) VD(IX+K)=VD(IX+K)+XBJ*VS(IY+IDP) ENDIF 60 CONTINUE ENDIF IF(MOD(ITPH/MSKH,3).EQ.2)THEN IBX2=IBASE(I+2) LL2=NSTAT(ITL+2,0) IY=IBX2+LL2*ISTAT2(ITPH-2*MSKH) *VDIR NODEP *VDIR LOOPCNT=73789 DO 65 K=0,LL - 1 ITPL=ISTAT1(K,ITL) IF(MOD(ITPL/MSKL,3).EQ.0)THEN IDP=ISTAT2(ITPL+2*MSKL) VD(IY+IDP)=VD(IY+IDP)+XBXX*VS(IX+K) VD(IX+K)=VD(IX+K)+XBXX*VS(IY+IDP) ENDIF 65 CONTINUE ENDIF 70 CONTINUE C ELSE DO 90 J=0,LH-1 IX=IBX+J*LL ITPH=ISTAT1(J,ITH) IF(MOD(ITPH/MSKH,3) .NE. 0) THEN IBX1=IBASE(I+1) LL1=NSTAT(ITL+1,0) IY=IBX1+LL1*ISTAT2(ITPH-MSKH) *VDIR NODEP *VDIR LOOPCNT=73789 DO 80 K=0,LL - 1 ITPL=ISTAT1(K,ITL) IF(MOD(ITPL/MSKL,3).NE. 2)THEN IDP=ISTAT2(ITPL+MSKL) VD(IY+IDP)=VD(IY+IDP)+XJX*VS(IX+K) VD(IX+K)=VD(IX+K)+XJX*VS(IY+IDP) ENDIF 80 CONTINUE ENDIF 90 CONTINUE ENDIF C 100 CONTINUE C RETURN END C C****************************************************************** C* KOBEPACK/1 VERSION 1.0 * C* MARCH 1992 * C* COPYRIGHT (C) M. KABURAGI, T. NISHINO AND T. TONEGAWA * C* WITH THE HELP OF TITPACK VERSION 2 * C* FEBRUARY 1991 * C* COPYRIGHT (C) HIDETOSHI NISHIMORI * C****************************************************************** C***** CALCULATION OF MATRIX ELEMENT FOR MID-SIZE ***** C****************************************************************** C* CORRESPONDENCE BETWEEN ISTAT1 AND SZ(I) * C* 0 === SZ(I)=-1 * C* 1 === SZ(I)= 0 * C* 2 === SZ(I)=+1 * C****************************************************************** C*** VARIABLES MARKED @ SHOULD BE GIVEN FROM THE MAIN PROGRAM C*** VARIABLES MARKED # ARE EVALUATED AND RETURNED C*** THE FOLLOWING VARIABLES ARE COMMON TO ALL ROUTINES C NS @ LATTICE SIZE C NSLMAX @ MAXIMUM SIZE OF LOWER PART OF LATTICE C ITSZ @ TOTAL SZ C NSB,NH,NL, @ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C MSGBC,NSTAT, @ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C ISTAT1,NDCLR1,@ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C ISTAT2,IBASE @ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C ISTAT1(I,IBZ) @ I-TH SPIN CONFIGURATION FOR IBZL(OR IBZH)=IBZ C BZ=SZ+1 C ISTAT2 @ INVERSE LIST OF ISTAT1 EXPRESSED BY THE C 2-DIM SEARCH ALGORITHM OF NISHINO C IPAIR @ PAIRS OF SITES CONNECTED BY BONDS C AJX @ EXCHANGE INTERACTION OF EACH BOND JX,Y C AJZ @ EXCHANGE INTERACTION OF EACH BOND JZ C ABIQD @ BIQUADRATIC INTERCTION C IBOND @ NUMBER OF BONDS C ANISTRPY @ SINGLE ION ANISOTROPY ENERGY AT EACH SITE C HFIELD @ MAGNETIC FIELD AT EACH SITE C DIAG @ DIAGONAL ELEMENT C*** CALCULATION OF ELEMENT (V1.2) C*** VARIABLES MARKED @ SHOULD BE GIVEN FROM THE MAIN PROGRAM C*** VARIABLES MARKED # ARE EVALUATED AND RETURNED C LBOND(I,1) # BOND NUMBER OF I-TH L-L BONDS C NBOND(1) # TOTAL NUMBER OF L-L BONDS C LBOND(I,2) # BOND NUMBER OF I-TH H-H BONDS C NBOND(2) # TOTAL NUMBER OF H-H BONDS C LBOND(I,3) # BOND NUMBER OF I-TH L-H BONDS C NBOND(3) # TOTAL NUMBER OF L-H BONDS C IBSPIN(I,1) # SITE NUMBER OF I-TH BOUNDARY L-SPINS C NBSPIN(1) # TOTAL NUMBER OF BOUNDARY L-SPINS C IBSPIN(I,2) # SITE NUMBER OF I-TH BOUNDARY H-SPINS C NBSPIN(2) # TOTAL NUMBER OF BOUNDARY H-SPINS C IBSPAIR(2K-1,2K) C # BOUNDARY SPIN NUMBER OF K-TH BOUNDARY PAIRS C (2K-1:L, 2K:H) C IPAIR(2K-1,2K) IS REPLACED WITH IPAIR(2K-1):L OR H C LBASE(I,J,1/2) @# BASE ADDRESS FOR J(LL/HH)-BOND DIAGONAL ELEMENT C LBASE(I,K,1/2) @# BASE ADDRESS FOR BOUNDARY L/H-SPIN DIAGONAL ELEMENT C MBASE(I,J,1/2/3) C @# BASE ADDRESS FOR J(LL/HH)-BOND OFF-DIAGONAL ELEMENT C MBASE(I,J,1/2/3) C @# BASE ADDRESS FOR BOUNDARY L/H-SPIN OFF-DIAGONAL C ELEMENT C IFLGDELM(I,J) @# FLAG FOR J-BOND DIAGONAL ELEMENT C IFLGOELM(I,J) @# FLAG FOR J-BOND OFF-DIAGONAL ELEMENT C ELMNT @ OFF-DIAGONAL ELEMENT BY NISHIMORI'S EXPRESSION C LOC @ LOCATION OF ADDRESS BY NISHIMORI'S EXPRESSION C C****************************************************************** C SUBROUTINE ELMM(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,DIAG,NDCLRD,LBOND,NBOND,LBASE,MBASE,IFLGDELM,IFLGOELM, &IBSPAIR,IBSPIN,NBSPIN,ELMNT,LOC,NDCLRM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION MBASE(0:2*NSLMAX,IBOND+NS,6) CHARACTER*4 MSGBC(2),MSG DIMENSION DIAG(0:NDCLRD) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) DIMENSION IBSPAIR(2*IBOND+12),IBSPIN(IBOND+6,2),NBSPIN(2) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) C MSG=MSGBC(1) C CALCULATION OF DIAGONAL ELEMENT C CALL CLSBND(NS,NSLMAX,IPAIR,IBOND,LBOND,NBOND,IBSPAIR, &IBSPIN,NBSPIN,MSGBC) C CALL ELMD(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD,MSGBC,DIAG, &NDCLRD,LBOND,NBOND,LBASE,IFLGDELM,IBSPAIR,IBSPIN,NBSPIN) C CALL ELMO(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,MSGBC,LBOND,NBOND,MBASE, &IFLGOELM,IBSPAIR,IBSPIN,NBSPIN,ELMNT,LOC,NDCLRM) RETURN END C C CLASSIFICATION OF THE BONDS INTO LL, LH AND HH **** C SUBROUTINE CLSBND(NS,NSLMAX,IPAIR,IBOND,LBOND,NBOND,IBSPAIR, &IBSPIN,NBSPIN,MSGBC) C DIMENSION IPAIR(2*IBOND+12),LBOND(IBOND+6,0:5),NBOND(5) DIMENSION IBSPAIR(2*IBOND+12),IBSPIN(IBOND+6,2),NBSPIN(2) CHARACTER*4 MSGBC(2) C MSGBC(1)=MSGBC(1) N=NS NBOND(1)=0 NBOND(2)=0 NBOND(3)=0 NBOND(4)=0 NBOND(5)=0 NBSPIN(1)=0 NBSPIN(2)=0 DO 100 KIS=1,IBOND ISITE1=(IPAIR(KIS*2-1)-1) ISITE2=(IPAIR(KIS*2 )-1) C L-L BOND IF((ISITE1+1.LE.NSLMAX).AND.(ISITE2+1.LE.NSLMAX))THEN NBOND(1)=NBOND(1)+1 LBOND(NBOND(1),1)=KIS LBOND(KIS,0)=1 C H-H BOND ELSEIF((ISITE1+1.GT.NSLMAX).AND.(ISITE2+1.GT.NSLMAX)) & THEN NBOND(2)=NBOND(2)+1 LBOND(NBOND(2),2)=KIS LBOND(KIS,0)=2 ENDIF C L-H BOND IF(((ISITE1+1.LE.NSLMAX).AND.(ISITE2+1.GT.NSLMAX)).OR. & ((ISITE1+1.GT.NSLMAX).AND.(ISITE2+1.LE.NSLMAX)))THEN NBOND(3)=NBOND(3)+1 LBOND(NBOND(3),3)=KIS LBOND(KIS,0)=3 IF((ISITE1+1.GT.NSLMAX).AND.(ISITE2+1.LE.NSLMAX))THEN ISITE2=(IPAIR(KIS*2-1)-1) ISITE1=(IPAIR(KIS*2 )-1) IPAIR(KIS*2-1)=ISITE1+1 IPAIR(KIS*2 )=ISITE2+1 ENDIF DO 10 I=1,NBSPIN(1) IF(IBSPIN(I,1).EQ.ISITE1+1)THEN GOTO 20 ENDIF 10 CONTINUE NBSPIN(1)=NBSPIN(1)+1 IBSPIN(NBSPIN(1),1)=ISITE1+1 I=NBSPIN(1) 20 CONTINUE IBSPAIR(2*NBOND(3)-1)=I C DO 30 I=1,NBSPIN(2) IF(IBSPIN(I,2).EQ.ISITE2+1)THEN GOTO 40 ENDIF 30 CONTINUE NBSPIN(2)=NBSPIN(2)+1 IBSPIN(NBSPIN(2),2)=ISITE2+1 I=NBSPIN(2) 40 CONTINUE IBSPAIR(2*NBOND(3))=I ENDIF 100 CONTINUE C RETURN END C C CALCULATION OF DIAGONAL ELEMENT C SUBROUTINE ELMD(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ANISTRPY,HFIELD, &MSGBC,DIAG,NDCLRD,LBOND,NBOND,LBASE,IFLGDELM,IBSPAIR,IBSPIN, &NBSPIN) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION LBOND(IBOND+6,0:5),NBOND(5),LBASE(0:2*NSLMAX, &IBOND+NS,2) DIMENSION IBSPAIR(2*IBOND+12),IBSPIN(IBOND+6,2),NBSPIN(2) CHARACTER*4 MSGBC(2) DIMENSION DIAG(0:NDCLRD) DIMENSION IFLGDELM(0:2*NSLMAX,IBOND+NS) DIMENSION IODER(0:24,3) C I=IBSPAIR(1) C DO 10 I=0,NSB-1 IODER(I,1)=I IODER(I,2)=I IODER(I,3)=I 10 CONTINUE DO 20 I=0,NSB-1 IL=IODER(I,1) IH=IODER(I,2) ILH=IODER(I,3) NSLMX=NSTAT(NL(IL),0) NSHMX=NSTAT(NH(IH),1) NSLHMX=NSTAT(NL(ILH),0)*NSTAT(NH(ILH),1) DO 20 J=I,NSB-1 JLP=IODER(J,1) NSLJ=NSTAT(NL(JLP),0) IF(NSLJ.GT.NSLMX)THEN IODER(J,1)=IODER(I,1) IODER(I,1)=JLP NSLMX=NSLJ ENDIF JHP=IODER(J,2) NSHJ=NSTAT(NH(JHP),1) IF(NSHJ.GT.NSHMX)THEN IODER(J,2)=IODER(I,2) IODER(I,2)=JHP NSHMX=NSHJ ENDIF JLHP=IODER(J,3) NSLHJ=NSTAT(NL(JLHP),0)*NSTAT(NH(JLHP),1) IF(NSLHJ.GT.NSLHMX)THEN IODER(J,3)=IODER(I,3) IODER(I,3)=JLHP NSLHMX=NSLHJ ENDIF 20 CONTINUE C DO 30 I=0,NSB-1 IFLGDELM(I,1)=0 IFLGDELM(I,2)=0 DO 30 J=1,NBSPIN(1)+NBSPIN(2) IFLGDELM(I,2+J)=0 30 CONTINUE C IANFLD=0 DO 50 K=1,NS IF(ABS(ANISTRPY(K)).GT.1.D-60 & .OR.ABS(HFIELD(K)).GT.1.D-60) THEN IANFLD=IANFLD+1 ENDIF 50 CONTINUE C DO 60 I=0,NDCLRD DIAG(I)= 0.0D0 60 CONTINUE C C L-L BOND C LBS=0 DO 70 IST=0,NSB-1 I=IODER(IST,1) LL=NSTAT(NL(I),0) LH=NSTAT(NH(I),1) LBASE(I,1,1)=LBS LBS1=LBASE(I,1,1) LBS=LBS+LL IF(LBS.LT.NDCLRD)THEN IF(IANFLD.NE.0)THEN CALL ELMDLLA(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, & ISTAT1,NDCLR1,ISTAT2,ANISTRPY,HFIELD,MSGBC, & DIAG(LBS1),LL) ENDIF CALL ELMDLLBX(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, & ISTAT1, & NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,MSGBC, & DIAG(LBS1),LL,LBOND,NBOND) IFLGDELM(I,1)=1 ELSE LBASE(I,1,1)=-1 IFLGDELM(I,1)=0 LBS=LBS-LL ENDIF 70 CONTINUE C IF(NS.GE.NSLMAX+1) THEN IF(NS.EQ.2*NSLMAX.AND.MSGBC(2).EQ.'SYMM')THEN DO 75 J=0,NSB-1 I=NSB-1-J LBASE(I,2,1)=LBASE(J,1,1) IFLGDELM(I,2)=IFLGDELM(J,1) 75 CONTINUE ELSE C H-H PAIR DO 80 IST=0,NSB-1 I=IODER(IST,2) LL=NSTAT(NL(I),0) LH=NSTAT(NH(I),1) LBASE(I,2,1)=LBS LBS2=LBASE(I,2,1) LBS=LBS+LH IF(LBS.LT.NDCLRD)THEN IF(IANFLD.NE.0)THEN CALL ELMDHHA(I,NS,NSLMAX,NSB,NH,NL,IBASE, & NSTAT, & ISTAT1,NDCLR1,ISTAT2,ANISTRPY,HFIELD,MSGBC, & DIAG(LBS2),LH) ENDIF CALL ELMDHHBX(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, & ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND, & MSGBC, & DIAG(LBS2),LH,LBOND,NBOND) IFLGDELM(I,2)=1 ELSE LBASE(I,2,1)=-1 IFLGDELM(I,2)=0 LBS=LBS-LH ENDIF 80 CONTINUE ENDIF C C C C L-H PAIR DO 100 ISS=0,NSB-1 I=IODER(ISS,3) LL=NSTAT(NL(I),0) LH=NSTAT(NH(I),1) KIS=1 90 CONTINUE C L-SPIN ELEMENT IF(KIS.LE.NBSPIN(1))THEN IST=IBSPIN(KIS,1) IBSP=2+NBOND(3)+KIS LBASE(I,IBSP,1)=LBS LBS3=LBASE(I,IBSP,1) LBS=LBS+LL IF(LBS.LT.NDCLRD)THEN CALL ELMDLS(I,IST,NS,NSLMAX,NSB,NH,NL,IBASE, & NSTAT, & ISTAT1,NDCLR1,ISTAT2,DIAG(LBS3),LL) IFLGDELM(I,IBSP)=1 ELSE LBASE(I,IBSP,1)=-1 IFLGDELM(I,IBSP)=0 LBS=LBS-LL ENDIF ENDIF C C H-SPIN IF(KIS.LE.NBSPIN(2))THEN IST=IBSPIN(KIS,2) IBSP=2+NBOND(3)+NBSPIN(1)+KIS LBASE(I,IBSP,1)=LBS LBS4=LBASE(I,IBSP,1) LBS=LBS+LH IF(LBS.LT.NDCLRD)THEN CALL ELMDHS(I,IST,NS,NSLMAX,NSB,NH,NL,IBASE, & NSTAT, & ISTAT1,NDCLR1,ISTAT2,DIAG(LBS4),LH) IFLGDELM(I,IBSP)=1 ELSE LBASE(I,IBSP,1)=-1 IFLGDELM(I,IBSP)=0 LBS=LBS-LH ENDIF ENDIF KIS=KIS+1 IF(KIS.LE.NBSPIN(1).OR.KIS.LE.NBSPIN(2))THEN GOTO 90 ENDIF 100 CONTINUE C DO 120 ISS=0,NSB-1 I=IODER(ISS,3) DO 110 KS=1,NBOND(3) ISPINL=IBSPAIR(KS*2-1) ISPINH=IBSPAIR(KS*2 ) IBSPL=2+NBOND(3)+ISPINL IBSPH=2+NBOND(3)+NBSPIN(1)+ISPINH IFLGSBL=IFLGDELM(I,IBSPL) IFLGSBH=IFLGDELM(I,IBSPH) IF(IFLGSBL.GT.0.AND.IFLGSBH.GT.0)THEN IFLGDELM(I,2+KS)=1 LBASE(I,2+KS,1)=LBASE(I,IBSPL,1) LBASE(I,2+KS,2)=LBASE(I,IBSPH,1) ELSE IFLGDELM(I,2+KS)=0 LBASE(I,2+KS,1)=-1 LBASE(I,2+KS,2)=-1 ENDIF 110 CONTINUE 120 CONTINUE C ENDIF C RETURN END C C CALCULATION OF L-L DIAGONAL ELEMENT FOR H AND D C SUBROUTINE ELMDLLA(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,ANISTRPY,HFIELD,MSGBC,DIAG,LM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION DIAG(0:LM-1) CHARACTER*4 MSGBC(2) C MSGBC(1)=MSGBC(1) J=NSB J=ISTAT2(1) J=IBASE(1) LL=NSTAT(NL(I),0) LH=NSTAT(NH(I),1) C DO 90 KIS=1,MIN(NS,NSLMAX) IST=3**(KIS-1) AN=ANISTRPY(KIS) HF=HFIELD(KIS) *VDIR LOOPCNT=73789 DO 20 K=0,LL -1 ITP=ISTAT1(K,NL(I)) DIAG(K)=DIAG(K) & -HF*DBLE(MOD(ITP/IST,3)-1) & +AN*DBLE((MOD(ITP/IST,3)-1)**2) 20 CONTINUE C 90 CONTINUE C 100 CONTINUE C RETURN END C C CALCULATION OF DIAGONAL L-L ELEMENT FOR BIQD.NE.0 C SUBROUTINE ELMDLLBX(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,MSGBC,DIAG, &LM,LBOND,NBOND) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION DIAG(0:LM-1) CHARACTER*4 MSGBC(2) DIMENSION LBOND(IBOND+6,0:5),NBOND(5) C MSGBC(1)=MSGBC(1) J=NS J=NSB J=ISTAT2(1) J=IBASE(1) A=AJX(1) LL=NSTAT(NL(I),0) LH=NSTAT(NH(I),1) C DO 90 KS=1,NBOND(1) KIS=LBOND(KS,1) ISITE1=(IPAIR(KIS*2-1)-1) IS1=3**ISITE1 ISITE2=(IPAIR(KIS*2 )-1) IS2=3**ISITE2 XJZ=AJZ(KIS) XB=ABIQD(KIS) XBXX=ABIQD(KIS) XBZZ=ABIQD(KIS) IF(ABS(XJZ).LE.1.D-60 & .AND.ABS(XB).LE.1.D-60) THEN GOTO 90 ENDIF IF(ABS(XB).GT.1.D-60) THEN *VDIR LOOPCNT=73789 DO 20 K=0,LL -1 ITP=ISTAT1(K,NL(I)) IST1=(MOD(ITP/IS1,3)-1) IST2=(MOD(ITP/IS2,3)-1) DIAG(K)=DIAG(K) & +XJZ*DBLE(IST1*IST2) & +XB*DBLE((IST1*IST2)**2 & +1+(IST1*IST2+1)*(IST1*IST2+2)/2 & -(IST1+IST2)**2) CCC & +XBZZ*DBLE((IST1*IST2)**2) CCC & +XBXX*DBLE(1+(IST1*IST2+1)*(IST1*IST2+2)/2 CCC & -(IST1+IST2)**2) 20 CONTINUE ELSE *VDIR LOOPCNT=73789 DO 30 K=0,LL -1 ITP=ISTAT1(K,NL(I)) IST1=(MOD(ITP/IS1,3)-1) IST2=(MOD(ITP/IS2,3)-1) DIAG(K)=DIAG(K) & +XJZ*DBLE(IST1*IST2) 30 CONTINUE ENDIF 90 CONTINUE C RETURN END C C CALCULATION OF H-H DIAGONAL ELEMENT FOR H AND D C SUBROUTINE ELMDHHA(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,ANISTRPY,HFIELD,MSGBC,DIAG,LM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION ANISTRPY(0:NS+1),HFIELD(0:NS+1) DIMENSION DIAG(0:LM-1) CHARACTER*4 MSGBC(2) C MSGBC(1)=MSGBC(1) J=NSB J=ISTAT2(1) J=IBASE(1) C LL=NSTAT(NL(I),0) LH=NSTAT(NH(I),1) DO 90 KIS=MIN(NS,NSLMAX)+1,NS IST=3**(KIS-NSLMAX-1) AN=ANISTRPY(KIS) HF=HFIELD(KIS) *VDIR LOOPCNT=73789 DO 20 J=0,LH -1 ITP=ISTAT1(J,NH(I)) DIAG(J)=DIAG(J) & -HF*DBLE(MOD(ITP/IST,3)-1) & +AN*DBLE((MOD(ITP/IST,3)-1)**2) 20 CONTINUE C 90 CONTINUE C RETURN END C C CALCULATION OF DIAGONAL H-H ELEMENT C SUBROUTINE ELMDHHBX(I,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,MSGBC,DIAG, &LM,LBOND,NBOND) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION DIAG(0:LM-1) CHARACTER*4 MSGBC(2) DIMENSION LBOND(IBOND+6,0:5),NBOND(5) C MSGBC(1)=MSGBC(1) J=NS J=NSB J=ISTAT2(1) J=IBASE(1) A=AJX(1) C LL=NSTAT(NL(I),0) LH=NSTAT(NH(I),1) DO 90 KS=1,NBOND(2) KIS=LBOND(KS,2) ISITE1=(IPAIR(KIS*2-1)-1) IS1=3**(ISITE1-NSLMAX) ISITE2=(IPAIR(KIS*2 )-1) IS2=3**(ISITE2-NSLMAX) XJZ=AJZ(KIS) XB=ABIQD(KIS) XBXX=ABIQD(KIS) XBZZ=ABIQD(KIS) IF(ABS(XJZ).LE.1.D-60 & .AND.ABS(XB).LE.1.D-60) THEN GOTO 90 ENDIF IF(ABS(XB).GT.1.D-60) THEN *VDIR LOOPCNT=73789 DO 20 J=0,LH -1 ITP=ISTAT1(J,NH(I)) IST1=(MOD(ITP/IS1,3)-1) IST2=(MOD(ITP/IS2,3)-1) DIAG(J)=DIAG(J) & +XJZ*DBLE(IST1*IST2) & +XB*DBLE((IST1*IST2)**2 & +1+(IST1*IST2+1)*(IST1*IST2+2)/2 & -(IST1+IST2)**2) CCC & +XBZZ*DBLE((IST1*IST2)**2) CCC & +XBXX*DBLE(1+(IST1*IST2+1)*(IST1*IST2+2)/2 CCC & -(IST1+IST2)**2) 20 CONTINUE ELSE DO 30 J=0,LH -1 ITP=ISTAT1(J,NH(I)) IST1=(MOD(ITP/IS1,3)-1) IST2=(MOD(ITP/IS2,3)-1) DIAG(J)=DIAG(J) & +XJZ*DBLE(IST1*IST2) 30 CONTINUE ENDIF 90 CONTINUE C RETURN END C C CALCULATION OF DIAGONAL L-SPIN ELEMENT C SUBROUTINE ELMDLS(I,IST,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,DIAG,LM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION DIAG(0:LM-1) C J=NS J=NSB J=ISTAT2(1) J=IBASE(1) LL=NSTAT(NL(I),0) LH=NSTAT(NH(I),1) C C FOR BOUNDARY L-SPIN ISITE1=IST-1 IS1=3**ISITE1 *VDIR LOOPCNT=73789 DO 20 K=0,LL -1 ITP=ISTAT1(K,NL(I)) IST1=(MOD(ITP/IS1,3)-1) DIAG(K)=DBLE(IST1) 20 CONTINUE C RETURN END C C CALCULATION OF DIAGONAL H SPIN C SUBROUTINE ELMDHS(I,IST,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,DIAG,LM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION DIAG(0:LM-1) C J=NS J=NSB J=ISTAT2(1) J=IBASE(1) LL=NSTAT(NL(I),0) LH=NSTAT(NH(I),1) C C FOR BOUNDARY H-SPIN ISITE2=IST-1 IS2=3**(ISITE2-NSLMAX) LH=NSTAT(NH(I),1) *VDIR LOOPCNT=73789 DO 180 J=0,LH-1 ITP=ISTAT1(J,NH(I)) IST2=(MOD(ITP/IS2,3)-1) DIAG(J)=DBLE(IST2) C 180 CONTINUE C RETURN END C C CALCULATION OF OFF-DIAGONAL ELEMENT C SUBROUTINE ELMO(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,MSGBC,LBOND,NBOND, &MBASE,IFLGOELM,IBSPAIR,IBSPIN,NBSPIN,ELMNT,LOC,NDCLRM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) CHARACTER*4 MSGBC(2) DIMENSION LBOND(IBOND+6,0:5),NBOND(5),MBASE(0:2*NSLMAX, &IBOND+NS,6) DIMENSION IBSPAIR(2*IBOND+12),IBSPIN(IBOND+6,2),NBSPIN(2) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) DIMENSION IFLGOELM(0:2*NSLMAX,IBOND+NS) DIMENSION IODER(0:24,3) C J=NSTAT(1,1) J=ISTAT1(1,1) J=ISTAT2(1) J=LBOND(1,1) J=IBSPAIR(1) DO 10 I=0,NSB-1 IODER(I,1)=I IODER(I,2)=I IODER(I,3)=I 10 CONTINUE DO 20 I=0,NSB-1 IL=IODER(I,1) IH=IODER(I,2) ILH=IODER(I,3) NSLMX=NSTAT(NL(IL),0) NSHMX=NSTAT(NH(IH),1) NSLHMX=NSTAT(NL(ILH),0)*NSTAT(NH(ILH),1) DO 20 J=I,NSB-1 JLP=IODER(J,1) NSLJ=NSTAT(NL(JLP),0) IF(NSLJ.GT.NSLMX)THEN IODER(J,1)=IODER(I,1) IODER(I,1)=JLP NSLMX=NSLJ ENDIF JHP=IODER(J,2) NSHJ=NSTAT(NH(JHP),1) IF(NSHJ.GT.NSHMX)THEN IODER(J,2)=IODER(I,2) IODER(I,2)=JHP NSHMX=NSHJ ENDIF JLHP=IODER(J,3) NSLHJ=NSTAT(NL(JLHP),0)*NSTAT(NH(JLHP),1) IF(NSLHJ.GT.NSLHMX)THEN IODER(J,3)=IODER(I,3) IODER(I,3)=JLHP NSLHMX=NSLHJ ENDIF 20 CONTINUE C DO 50 I=0,NSB-1 DO 50 J=1,IBOND+NBSPIN(1)+NBSPIN(2) IFLGOELM(I,J)=0 50 CONTINUE C C L-L BOND C MBS=0 DO 70 IST=0,NSB-1 I=IODER(IST,1) LL=NSTAT(NL(I),0) LH=NSTAT(NH(I),1) DO 70 KS=1,NBOND(1) KIS=LBOND(KS,1) XJX=AJX(KIS) XB=ABIQD(KIS) IC=2 IF(ABS(XB).LE.1.D-60.AND.ABS(XJX).GT.1.D-60)THEN IC=1 ELSEIF(ABS(XB).LE.1.D-60.AND.ABS(XJX).LE.1.D-60) & THEN IC=0 ENDIF MBSE=MBS MBSL=MBS MBASE(I,KIS,1)=MBSE MBASE(I,KIS,2)=MBSL MBASE(I,KIS,3)=IC MBS1E=MBASE(I,KIS,1) MBS1L=MBASE(I,KIS,2) IC1=MBASE(I,KIS,3) MBS=MBS+LL*IC MBSE=MBS MBSL=MBS IF(MBS.LT.NDCLRM.AND.IC.GT.0)THEN CALL ELMOLLBX(I,KIS,NS,NSLMAX,NSB,NH,NL,IBASE, & NSTAT, & ISTAT1, & NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND, & ELMNT(MBS1E),LOC(MBS1L),LL,IC1) IFLGOELM(I,KIS)=1 ELSE MBASE(I,KIS,1)=-1 MBASE(I,KIS,2)=-1 IFLGOELM(I,KIS)=0 MBS=MBS-LL*IC MBSE=MBS MBSL=MBS ENDIF 70 CONTINUE C IF(NS.GE.NSLMAX+2)THEN C IF(NS.EQ.2*NSLMAX.AND.MSGBC(2).EQ.'SYMM')THEN DO 75 J=0,NSB-1 I=NSB-1-J DO 75 KS=1,NBOND(2) KISL=LBOND(KS,1) KISH=LBOND(KS,2) MBASE(I,KISH,1)=MBASE(J,KISL,1) MBASE(I,KISH,2)=MBASE(J,KISL,2) MBASE(I,KISH,3)=MBASE(J,KISL,3) IFLGOELM(I,KISH)=IFLGOELM(J,KISL) 75 CONTINUE ELSE C H-H BOND DO 80 IST=0,NSB-1 I=IODER(IST,2) LL=NSTAT(NL(I),0) LH=NSTAT(NH(I),1) DO 80 KS=1,NBOND(2) KIS=LBOND(KS,2) XJX=AJX(KIS) XB=ABIQD(KIS) IC=2 IF(ABS(XB).LE.1.D-60.AND.ABS(XJX).GT.1.D-60)THEN IC=1 ELSEIF(ABS(XB).LE.1.D-60.AND.ABS(XJX) & .LE.1.D-60) THEN IC=0 ENDIF MBSE=MBS MBSL=MBS MBASE(I,KIS,1)=MBSE MBASE(I,KIS,2)=MBSL MBASE(I,KIS,3)=IC MBS2E=MBASE(I,KIS,1) MBS2L=MBASE(I,KIS,2) IC2=MBASE(I,KIS,3) MBS=MBS+LH*IC MBSE=MBS MBSL=MBS IF(MBS.LT.NDCLRM.AND.IC2.GT.0)THEN CALL ELMOHHBX(I,KIS,NS,NSLMAX,NSB,NH,NL, & IBASE, & NSTAT, & ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD, & IBOND, & ELMNT(MBS2E),LOC(MBS2L),LH,IC) IFLGOELM(I,KIS)=1 ELSE MBASE(I,KIS,1)=-1 MBASE(I,KIS,2)=-1 IFLGOELM(I,KIS)=0 MBS=MBS-LH*IC MBSE=MBS MBSL=MBS ENDIF 80 CONTINUE C ENDIF ENDIF C IF(NS.GE.NSLMAX+1)THEN C L-H SPIN C L-SPIN IBIQD=0 DO 85 KS=1,NBOND(3) KIS=LBOND(KS,3) IF(ABS(ABIQD(KIS)).GT.1.D-60) THEN IBIQD=IBIQD+1 ENDIF 85 CONTINUE C MBND=1 IF(IBIQD.GT.0)THEN MBND=4 ENDIF IC=MBND C DO 100 ISS=0,NSB-1 I=IODER(ISS,3) LL=NSTAT(NL(I),0) LH=NSTAT(NH(I),1) KIS=1 90 CONTINUE IF(KIS.LE.NBSPIN(1))THEN IST=IBSPIN(KIS,1) KBSP=IBOND+KIS MBSE=MBS MBSL=MBS MBASE(I,KBSP,1)=MBSE MBASE(I,KBSP,2)=MBSL MBASE(I,KBSP,3)=IC MBS3E=MBASE(I,KBSP,1) MBS3L=MBASE(I,KBSP,2) MBS=MBS+LL*IC MBSE=MBS MBSL=MBS IF(MBS.LT.NDCLRM)THEN CALL ELMOLSBX(I,IST,NS,NSLMAX,NSB,NH,NL,IBASE, & NSTAT, & ISTAT1,NDCLR1,ISTAT2,ELMNT(MBS3E),LOC(MBS3L),LL, & IC) IFLGOELM(I,KBSP)=1 ELSE MBASE(I,KBSP,1)=-1 MBASE(I,KBSP,2)=-1 MBASE(I,KBSP,3)=-1 IFLGOELM(I,KBSP)=0 MBS=MBS-LL*IC MBSE=MBS MBSL=MBS ENDIF ENDIF C H-SPIN IF(KIS.LE.NBSPIN(2))THEN IST=IBSPIN(KIS,2) KBSP=IBOND+NBSPIN(1)+KIS MBSE=MBS MBSL=MBS MBASE(I,KBSP,1)=MBSE MBASE(I,KBSP,2)=MBSL MBASE(I,KBSP,3)=IC MBS4E=MBASE(I,KBSP,1) MBS4L=MBASE(I,KBSP,2) MBS=MBS+LH*IC MBSE=MBS MBSL=MBS IF(MBS.LT.NDCLRM)THEN CALL ELMOHSBX(I,IST,NS,NSLMAX,NSB,NH,NL,IBASE, & NSTAT, & ISTAT1,NDCLR1,ISTAT2,ELMNT(MBS4E),LOC(MBS4L),LH, & IC) IFLGOELM(I,KBSP)=1 ELSE MBASE(I,KBSP,1)=-1 MBASE(I,KBSP,2)=-1 MBASE(I,KBSP,3)=-1 IFLGOELM(I,KBSP)=0 MBS=MBS-LH*IC MBSE=MBS MBSL=MBS ENDIF ENDIF KIS=KIS+1 IF(KIS.LE.NBSPIN(1).OR.KIS.LE.NBSPIN(2))THEN GOTO 90 ENDIF 100 CONTINUE C DO 120 ISS=0,NSB-1 I=IODER(ISS,3) DO 110 KS=1,NBOND(3) KIS=LBOND(KS,3) ISPINL=IBSPAIR(KS*2-1) ISPINH=IBSPAIR(KS*2 ) IFLGSBL=IFLGOELM(I,IBOND+ISPINL) IFLGSBH=IFLGOELM(I,IBOND+NBSPIN(1)+ISPINH) IF(IFLGSBL.GT.0.AND.IFLGSBH.GT.0)THEN IFLGOELM(I,KIS)=1 IBSP=IBOND+ISPINL MBASE(I,KIS,1)=MBASE(I,IBSP,1) MBASE(I,KIS,2)=MBASE(I,IBSP,2) MBASE(I,KIS,3)=MBASE(I,IBSP,3) IBSP=IBOND+NBSPIN(1)+ISPINH MBASE(I,KIS,4)=MBASE(I,IBSP,1) MBASE(I,KIS,5)=MBASE(I,IBSP,2) MBASE(I,KIS,6)=MBASE(I,IBSP,3) ELSE IFLGOELM(I,KIS)=0 MBASE(I,KIS,1)=-1 MBASE(I,KIS,2)=-1 MBASE(I,KIS,3)=-1 MBASE(I,KIS,4)=-1 MBASE(I,KIS,5)=-1 MBASE(I,KIS,6)=-1 ENDIF 110 CONTINUE 120 CONTINUE C ENDIF C RETURN END C C CALCULATION OF OFF-DIAGONAL LL ELEMENT FOR BIQD.NE.0 C SUBROUTINE ELMOLLBX(I,KBOND,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ELMNT,LOC,LM, &ICC) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ELMNT(0:LM-1,ICC),LOC(0:LM-1,ICC) C J=NS J=NSB J=ISTAT2(1) A=AJZ(1) C ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) C KIS=KBOND ISITE1=(IPAIR(KIS*2-1)-1) ISITE2=(IPAIR(KIS*2 )-1) MSKH=3**ISITE1 MSKL=3**ISITE2 XJX=AJX(KIS) XB=ABIQD(KIS) C IF(ICC.LE.0)THEN GOTO 80 ENDIF IF(ICC.EQ.2)THEN XBXX=XB XBXZ=XB *VDIR NODEP *VDIR LOOPCNT=73789 DO 30 K=0,LL -1 ITP=ISTAT1(K,ITL) IF(MOD(ITP/MSKH,3).NE.0.AND.MOD(ITP/MSKL,3).NE.2)THEN ELMNT(K,1)=XJX-XBXZ*DBLE(1- & (MOD(ITP/MSKH,3)+MOD(ITP/MSKL,3)-2)**2) LOC(K,1)=ISTAT2(ITP - MSKH+MSKL) ELSE ELMNT(K,1)=0.D0 LOC(K,1)=K ENDIF C IF(MOD(ITP/MSKH,3).EQ.2.AND.MOD(ITP/MSKL,3).EQ.0)THEN LOC(K,2)=ISTAT2(ITP -2*MSKH+2*MSKL) ELMNT(K,2)=XBXX ELSE ELMNT(K,2)=0.D0 LOC(K,2)=K ENDIF 30 CONTINUE C ELSEIF(ICC.EQ.1) THEN *VDIR NODEP *VDIR LOOPCNT=73789 DO 40 K=0,LL -1 ITP=ISTAT1(K,ITL) IF(MOD(ITP/MSKH,3).NE.0.AND.MOD(ITP/MSKL,3).NE.2)THEN ELMNT(K,1)=XJX LOC(K,1)=ISTAT2(ITP - MSKH+MSKL) ELSE ELMNT(K,1)=0.D0 LOC(K,1)=K ENDIF 40 CONTINUE ENDIF 80 CONTINUE C RETURN END C C CALCULATION OF OFF-DIAGONAL HH ELEMENT FOR BIQD.NE.0 C SUBROUTINE ELMOHHBX(I,KBOND,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,ELMNT,LOC,LM, &ICC) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBOND+12) DIMENSION AJX(IBOND+6),AJZ(IBOND+6),ABIQD(IBOND+6) DIMENSION ELMNT(0:LM-1,ICC),LOC(0:LM-1,ICC) C J=NS J=NSB J=IPAIR(1) J=ISTAT2(1) A=AJZ(1) IHSHFT=3**NSLMAX ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) IBX=IBASE(I) C KIS=KBOND ISITE1=(IPAIR(KIS*2-1)-1) ISITE2=(IPAIR(KIS*2 )-1) MSKH=3**(ISITE1-NSLMAX) MSKL=3**(ISITE2-NSLMAX) XJX=AJX(KIS) XB=ABIQD(KIS) XBXX=XB XBXZ=XB C IF(ICC.LE.0)THEN GOTO 80 ENDIF IF(ICC.EQ.2) THEN *VDIR NODEP *VDIR LOOPCNT=73789 DO 30 J=0,LH -1 ITP=ISTAT1(J,ITH) IF(MOD(ITP/MSKH,3).NE.0.AND.MOD(ITP/MSKL,3).NE.2)THEN ELMNT(J,1)=XJX-XBXZ*DBLE(1- & (MOD(ITP/MSKH,3)+MOD(ITP/MSKL,3)-2)**2) LOC(J,1)=ISTAT2(ITP - MSKH+MSKL) ELSE ELMNT(J,1)=0.D0 LOC(J,1)=J ENDIF IF(MOD(ITP/MSKH,3).EQ.2.AND.MOD(ITP/MSKL,3).EQ.0)THEN LOC(J,2)=ISTAT2(ITP -2*MSKH+2*MSKL) ELMNT(J,2)=XBXX ELSE LOC(J,2)=J ELMNT(J,2)=0.D0 ENDIF 30 CONTINUE C ELSEIF(ICC.EQ.1) THEN *VDIR NODEP *VDIR LOOPCNT=73789 DO 40 J=0,LH -1 ITP=ISTAT1(J,ITH) IF(MOD(ITP/MSKH,3).NE.0.AND.MOD(ITP/MSKL,3).NE.2)THEN ELMNT(J,1)=XJX LOC(J,1)=ISTAT2(ITP - MSKH+MSKL) ELSE ELMNT(J,1)=0.D0 LOC(J,1)=J ENDIF 40 CONTINUE ENDIF 80 CONTINUE C RETURN END C C CALCULATION OF OFF-DIAGONAL L-SPIN ELEMENT FOR BIQD.NE.0 C SUBROUTINE ELMOLSBX(I,IST,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,ELMNT,LOC,LM,ICC) C C ELMNT(K,1)=(Sx+iSy)* K)/SQRT(2) C ELMNT(K,2)=(Sx+iSy)Sz* K)/SQRT(2) C ELMNT(K,3)=Sz*(Sx+iSy) K)/SQRT(2) C ELMNT(K,4)=(Sx+iSy)**2 K)/2 C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION ELMNT(0:LM-1,ICC),LOC(0:LM-1,ICC) C J=NS J=NSB J=ISTAT2(1) ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) MSKL=3**(IST-1) IBX=IBASE(I) IF(ICC.EQ.4)THEN *VDIR NODEP *VDIR LOOPCNT=73789 DO 60 K=0,LL - 1 ITPL=ISTAT1(K,ITL) IF(MOD(ITPL/MSKL,3).NE. 2)THEN ELMNT(K,1)=1.D0 LOC(K,1)=ISTAT2(ITPL+MSKL) ELSE ELMNT(K,1)=0.D0 LOC(K,1)=-1 ENDIF IF(MOD(ITPL/MSKL,3).EQ.0)THEN ELMNT(K,2)=-1.D0 LOC(K,2)=ISTAT2(ITPL+MSKL) ELMNT(K,4)=1.D0 LOC(K,4)=ISTAT2(ITPL+2*MSKL) ELSE ELMNT(K,2)=0.D0 LOC(K,2)=-1 ELMNT(K,4)=0.D0 LOC(K,4)=-1 ENDIF IF(MOD(ITPL/MSKL,3).EQ.1)THEN ELMNT(K,3)=1.D0 LOC(K,3)=ISTAT2(ITPL+MSKL) ELSE ELMNT(K,3)=0.D0 LOC(K,3)=-1 ENDIF 60 CONTINUE C ELSEIF(ICC.EQ.1)THEN *VDIR NODEP *VDIR LOOPCNT=73789 DO 70 K=0,LL - 1 ITPL=ISTAT1(K,ITL) IF(MOD(ITPL/MSKL,3).NE. 2)THEN ELMNT(K,1)=1.D0 LOC(K,1)=ISTAT2(ITPL+MSKL) ELSE ELMNT(K,1)=0.D0 LOC(K,1)=-1 ENDIF 70 CONTINUE C ENDIF RETURN END C C CALCULATION OF OFF-DIAGONAL H-SPIN ELEMENT FOR BIQD.NE.0 C SUBROUTINE ELMOHSBX(I,IST,NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2,ELMNT,LOC,LM,ICC) C C ELMNT(J,1)=(Sx-iSy)* J)/SQRT(2) C ELMNT(J,2)=(Sx-iSy)Sz* J)/SQRT(2) C ELMNT(J,3)=Sz(Sx-iSy)* J)/SQRT(2) C ELMNT(J,4)=(Sx-iSy)**2 J)/2 C IMPLICIT REAL*8 (A-H,O-Z) INTEGER IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) INTEGER NSTAT(0:2*NSLMAX,0:1) INTEGER ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION ELMNT(0:LM-1,ICC),LOC(0:LM-1,ICC) C J=NS J=NSB J=ISTAT2(1) ITL=NL(I) ITH=NH(I) LL=NSTAT(ITL,0) LH=NSTAT(ITH,1) MSKH=3**(IST-1-NSLMAX) IBX=IBASE(I) IF(ICC.EQ.4)THEN *VDIR NODEP *VDIR LOOPCNT=73789 DO 60 J=0,LH - 1 ITPH=ISTAT1(J,ITH) IF(MOD(ITPH/MSKH,3).NE.0)THEN ELMNT(J,1)=1.D0 LOC(J,1)=ISTAT2(ITPH-MSKH) ELSE ELMNT(J,1)=0.D0 LOC(J,1)=-1 ENDIF IF(MOD(ITPH/MSKH,3).EQ.2)THEN ELMNT(J,2)=1.D0 LOC(J,2)=ISTAT2(ITPH-MSKH) ELMNT(J,4)=1.D0 LOC(J,4)=ISTAT2(ITPH-2*MSKH) ELSE ELMNT(J,2)=0.D0 LOC(J,2)=-1 ELMNT(J,4)=0.D0 LOC(J,4)=-1 ENDIF IF(MOD(ITPH/MSKH,3).EQ.1)THEN ELMNT(J,3)=-1.D0 LOC(J,3)=ISTAT2(ITPH-MSKH) ELSE ELMNT(J,3)=0.D0 LOC(J,3)=-1 ENDIF 60 CONTINUE C ELSEIF(ICC.EQ.1)THEN *VDIR NODEP *VDIR LOOPCNT=73789 DO 70 J=0,LH - 1 ITPH=ISTAT1(J,ITH) IF(MOD(ITPH/MSKH,3).NE.0)THEN ELMNT(J,1)=1.D0 LOC(J,1)=ISTAT2(ITPH-MSKH) ELSE ELMNT(J,1)=0.D0 LOC(J,1)=-1 ENDIF 70 CONTINUE ENDIF C RETURN END C C****************************************************************** C* KOBEPACK/1 VERSION 1.0 * C* MARCH 1992 * C* COPYRIGHT (C) M. KABURAGI, T. NISHINO AND T. TONEGAWA * C* WITH THE HELP OF TITPACK VERSION 2 * C* FEBRUARY 1991 * C* COPYRIGHT (C) HIDETOSHI NISHIMORI * C****************************************************************** C***** COMMON ROUTINE ***** C****************************************************************** C* CORRESPONDENCE BETWEEN ISTAT1 AND SZ(I) * C* 0 === SZ(I)=-1 * C* 1 === SZ(I)= 0 * C* 2 === SZ(I)=+1 * C****************************************************************** C*** VARIABLES MARKED @ SHOULD BE GIVEN FROM THE MAIN PROGRAM C*** VARIABLES MARKED # ARE EVALUATED AND RETURNED C*** THE FOLLOWING VARIABLES ARE COMMON TO ALL ROUTINES C NS @ LATTICE SIZE C NSLMAX @ MAXIMUM SIZE OF LOWER PART OF LATTICE C ITSZ @ TOTAL SZ C MSGBC @ MESSAGE OF BOUNDARY CONDITION C NSB,NH,NL, @ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C NSTAT, @ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C ISTAT1,NDCLR1,@ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C ISTAT2,IBASE @ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C ISTAT1(I,IBZ) @ I-TH SPIN CONFIGURATION FOR IBZL(OR IBZH)=IB C BZ=SZ+1 C ISTAT2 @ INVERSE LIST OF ISTAT1 EXPRESSED BY THE C 2-DIM SEARCH ALGORITHM OF NISHINO C IPAIR @ PAIRS OF SITES CONNECTED BY BONDS C AJX @ EXCHANGE INTERACTION OF EACH BOND JX,Y C AJZ @ EXCHANGE INTERACTION OF EACH BOND JZ C ABIQD @ BIQUADRATIC INTERCTION C IBOND @ NUMBER OF BONDS C ANISTRPY @ SINGLE ION ANISOTROPY ENERGY AT EACH SITE C HFIELD @ MAGNETIC FIELD AT EACH SITE C DIAG @ DIAGONAL ELEMENT C C NUMBER OF CONFIGURATIONS FOR GIVEN NS AND ISZ C FUNCTION NCONF(NS,ITSZ) N=NS ISZ=ITSZ NCONF=0 IASZ=IABS(ISZ) KMAX=(N-IASZ)/2 IF(KMAX.LT.0)RETURN I1=1 I2=1 IF(IASZ.GT.0) THEN DO 10 I=1,IASZ I2=I2*(N-I+1)/I 10 CONTINUE ENDIF NCONF=I1*I2 IF(KMAX.EQ.0)RETURN DO 20 K=1,KMAX I1=I1*(N-K+1)/K I2=I2*(N-IASZ-2*K+2)*(N-IASZ-2*K+1)/(N-K+1)/(IASZ+K) NCONF=NCONF+I1*I2 20 CONTINUE RETURN END C C DATA CHECK OF PAIRS OF SITES C SUBROUTINE PAIRCHK(IPAIR,IBOND,NS) DIMENSION IPAIR(IBOND*2) C N=NS DO 10 K=1,IBOND ISITE1=IPAIR(K*2-1) ISITE2=IPAIR(K*2 ) IF(ISITE1.LE.0.OR.ISITE2.LE.0.OR. & ISITE1.GT.N.OR.ISITE2.GT.N)THEN PRINT *,' #(E04)# INCORRECT DATA IN IPAIR.' PRINT *,' LOCATION: ',K*2-1,', ',K*2 STOP END IF 10 CONTINUE RETURN END C C CONFIGURATIONS WITH THE SPECIFIED SZ C SUBROUTINE CRESZ(NS,NSLMAX,ITSZ,NSB,NH,NL,IBASE,NSTAT, &ISTAT1,NDCLR1,ISTAT2) C C*** VARIABLES MARKED @ SHOULD BE GIVEN FROM THE MAIN PROGRAM C*** VARIABLES MARKED # ARE EVALUATED AND RETURNED C ITSZ @ TOTAL SZ C ISTAT1(I,IBZ) # I-TH SPIN CONFIGURATION FOR IBZL(OR IBZH)=IB C BZ=SZ+1 C ISTAT2 # INVERSE LIST OF ISTAT1 EXPRESSED BY THE C 2-DIM SEARCH ALGORITHM OF NISHINO C & DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) C CALL PRMCHK(NS,NSLMAX,ITSZ) CALL CRENUM(NS,NSLMAX,NSB,ITSZ,NH,NL) CALL CRESTA(NS,NSLMAX,NSTAT,ISTAT1,NDCLR1,ISTAT2) CALL CRBASE(NSB,NSLMAX,NH,NL,IBASE,NSTAT) RETURN END C SUBROUTINE PRMCHK(NS,NSLMAX,ITSZ) IMPLICIT REAL*8 (A-H,O-Z) C NSL12=12 IF(NS.LT.2.OR.NS.GT.2*NSL12) THEN WRITE (6,100) NS STOP 100 FORMAT(1H , & ' #(E05)# NS IS SMALLER THAN 2 OR LARGER THAN 24. ', & 'NS=',I5) END IF IF(ABS(ITSZ).GT.NS) THEN WRITE(6,101) NS,ITSZ STOP 101 FORMAT(1H , & ' #(E06)# ITSZ IS LARGER THAN NS OR SMALLER THAN -NS. ', & 'NS=',I5,', ITSZ=',I5) END IF IF(NSLMAX.LT.(NS+1)/2.OR.NSLMAX.GT.NS) THEN WRITE(6,102) NS,NSLMAX STOP 102 FORMAT(1H , & ' #(E07)# NSLMAX IS SMALLER THAN (NS+1)/2 OR LARGER ', & 'THAN NS.'/1H ,9X, & 'NS=',I5,', NSLMAX=',I5) END IF C RETURN END C C CREATION OF COMBINATION OF HIGHER-PART SZ AND LOWER-PART SZ C SUBROUTINE CRENUM(NS,NSLMAX,NSB,ITSZ,NH,NL) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION NH(0:2*NSLMAX),NL(0:2*NSLMAX) C NSB=1 IF(NS .LE. NSLMAX) THEN NH(0)= 0 NL(0)= NS+ITSZ GOTO 30 END IF C NH(0)= 2*(NS-NSLMAX) NL(0)= 0 10 CONTINUE C IF(NH(0)+ NL(0) .GT. NS+ITSZ) THEN NH(0)= NH(0)- 1 GOTO 10 END IF C IF(NH(0)+ NL(0) .LT. NS+ITSZ) THEN NL(0)= NL(0)+ 1 GOTO 10 END IF C 20 CONTINUE IF(NH(NSB-1).EQ.0 .OR. NL(NSB-1).EQ.2*NSLMAX) GOTO 30 NH(NSB)= NH(NSB-1) - 1 NL(NSB)= NL(NSB-1)+1 NSB=NSB+1 GOTO 20 C 30 CONTINUE RETURN END C C CREATION OF STATES C SUBROUTINE CRESTA(NS,NSLMAX,NSTAT,ISTAT1,NDCLR1,ISTAT2) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) C NSL12=12 IF(NSLMAX.GT.NSL12)THEN WRITE(6,*)' #(E08)# WRONG VALUE GIVEN TO NSLMAX IN ', & 'CRESTA CALLED FROM CRESZ.' WRITE(6,*)' NSLMAX=',NSLMAX STOP ENDIF C DO 20 I=0,2*NSLMAX NSTAT(I,0)= 0 NSTAT(I,1)= 0 20 CONTINUE C *VDIR LOOPCNT=531441 DO 30 I=0,3**MIN(NS,NSLMAX)-1 IBIT2=MOD(I/3** 0,3)+MOD(I/3** 1,3)+MOD(I/3** 2,3) & +MOD(I/3** 3,3)+MOD(I/3** 4,3)+MOD(I/3** 5,3) & +MOD(I/3** 6,3)+MOD(I/3** 7,3)+MOD(I/3** 8,3) & +MOD(I/3** 9,3)+MOD(I/3**10,3)+MOD(I/3**11,3) ISTAT1(NSTAT(IBIT2,0),IBIT2)= I ISTAT2(I)=NSTAT(IBIT2,0) NSTAT(IBIT2,0)=NSTAT(IBIT2,0)+1 30 CONTINUE C *VDIR LOOPCNT=531441 DO 40 I= 0,3**MAX(0,NS-NSLMAX)-1 IBIT2=MOD(I/3** 0,3)+MOD(I/3** 1,3)+MOD(I/3** 2,3) & +MOD(I/3** 3,3)+MOD(I/3** 4,3)+MOD(I/3** 5,3) & +MOD(I/3** 6,3)+MOD(I/3** 7,3)+MOD(I/3** 8,3) & +MOD(I/3** 9,3)+MOD(I/3**10,3)+MOD(I/3**11,3) ISTAT1(NSTAT(IBIT2,1),IBIT2)= I NSTAT(IBIT2,1)=NSTAT(IBIT2,1)+1 40 CONTINUE C RETURN END C C C SUBROUTINE CRBASE(NSB,NSLMAX,NH,NL,IBASE,NSTAT) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) C ****** IBASE(NSB) IS DIMENSION OF TOTAL SPACE ** IBASE(0)= 0 DO 40 I=1,NSB IBASE(I)= IBASE(I-1)+NSTAT(NH(I-1),1)* NSTAT(NL(I-1),0) 40 CONTINUE C RETURN END C C C FUNCTION ITOTBZ(I) IF(I.LT.0.OR.I.GT.3**12-1) THEN WRITE (6,*) & ' #(E09)# WRONG VALUE GIVEN TO I IN ITOTBZ. ', & 'I=',I STOP END IF ITOTBZ=MOD(I/3** 0,3)+MOD(I/3** 1,3)+MOD(I/3** 2,3)+ &MOD(I/3** 3,3)+MOD(I/3** 4,3)+MOD(I/3** 5,3)+MOD(I/3** 6,3)+ &MOD(I/3** 7,3)+MOD(I/3** 8,3)+MOD(I/3** 9,3)+MOD(I/3**10,3)+ &MOD(I/3**11,3) RETURN END C C INNER PRODUCT OF VEC C SUBROUTINE INPRO(IDIM,V1,V0,PROD) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION V1(0:IDIM-1),V0(0:IDIM-1) PROD=0.0D0 *VDIR LOOPCNT=5196627 DO 10 I=0,IDIM-1 PROD=PROD+V1(I)*V0(I) 10 CONTINUE RETURN END C C LINEAR COMBINATION: VD=A*VECS1+B*VECS2 C SUBROUTINE ADDER(IDIM,VD,VS1,VS2,A1,A2) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION VD(0:IDIM-1),VS1(0:IDIM-1),VS2(0:IDIM-1) *VDIR LOOPCNT=5196627 DO 10 I=0,IDIM-1 VD(I)=A1*VS1(I)+A2*VS2(I) 10 CONTINUE RETURN END C C DIVISION OF VEC BY D C SUBROUTINE DIVER(IDIM,VD,VS,D) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) *VDIR LOOPCNT=5196627 DO 10 I=0,IDIM-1 VD(I)=VS(I)/D 10 CONTINUE RETURN END C C CLEAR VEC C SUBROUTINE VEC01(IDIM,V) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION V(0:IDIM-1) *VDIR LOOPCNT=5196627 DO 5 J=0,IDIM-1 V(J)=0.0D0 5 CONTINUE RETURN END C SUBROUTINE VEC02(IDIM,V1,V2) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION V1(0:IDIM-1),V2(0:IDIM-1) *VDIR LOOPCNT=5196627 DO 5 J=0,IDIM-1 V1(J)=0.0D0 V2(J)=0.0D0 5 CONTINUE RETURN END C C CALCULATION OF NORM OF LINEAR COMBINATION: PROD=(A1*V1+A2*V2)**2 C SUBROUTINE NORMER(IDIM,PROD,V1,V2,A1,A2) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION V1(0:IDIM-1),V2(0:IDIM-1) PROD=0.D0 *VDIR LOOPCNT=5196627 DO 10 I=0,IDIM-1 PROD=PROD+(A1*V1(I)+A2*V2(I))**2 10 CONTINUE PROD=SQRT(PROD) RETURN END C C SWAP OF (VEC1,VEC2)=(A11*VEC1+A12*VEC2, A21*VEC1+A22*VEC2) C SUBROUTINE SWAPPER(IDIM,V1,V2,A11,A12,A21,A22) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION V1(0:IDIM-1),V2(0:IDIM-1) *VDIR LOOPCNT=5196627 DO 10 I=0,IDIM-1 TEMP1=A11*V1(I)+A12*V2(I) TEMP2=A21*V1(I)+A22*V2(I) V1(I)=TEMP1 V2(I)=TEMP2 10 CONTINUE RETURN END C C BLOCK MOVE C SUBROUTINE VMOVE(IDIM,VD,VS) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION VD(0:IDIM-1),VS(0:IDIM-1) *VDIR LOOPCNT=5196627 DO 10 I=0,IDIM-1 VD(I)=VS(I) 10 CONTINUE RETURN END C C VEC INITIALIZATION C SUBROUTINE VECINI(IDIM,VS,IV) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IV(20) DIMENSION VS(0:IDIM-1) CALL VEC01(IDIM,VS) IVM=0 DO 20 I=1,20 IF(IV(I).GE.0.AND.IV(I).LE.IDIM-1)THEN IVM=IVM+1 ENDIF 20 CONTINUE IF(IVM.EQ.0)THEN DO 30 I=1,20 C IV(I)=INT(RAND(FLOAT(IDIM-1))) IV(I)=MOD(IRAND(0),IDIM) 30 CONTINUE ENDIF DO 40 I=1,20 IF(IV(I).GE.0.AND.IV(I).LE.IDIM-1)THEN VS(IV(I))=1.0D0 ENDIF 40 CONTINUE CALL INPRO(IDIM,VS,VS,PROD) PROD=SQRT(PROD) CALL DIVER(IDIM,VS,VS,PROD) RETURN END C C C SUBROUTINE DISPNNF(NS,ITP,SP) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*1 SP(0:NS-1) C DO 10 I=0,NS-1 IF(MOD(ITP,3) .EQ. 0) SP(I)='-' IF(MOD(ITP,3) .EQ. 1) SP(I)='O' IF(MOD(ITP,3) .EQ. 2) SP(I)='+' ITP=ITP / 3 10 CONTINUE C RETURN END C C C SUBROUTINE CHKWV(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,VEC,NDCLRV,MDCLRV,ISS,NV,ICYCLST,ICYCLNO, &VCOEF,IDCLR) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION ICYCLST(IDCLR,NV),ICYCLNO(IDCLR,NV),VCOEF(IDCLR, &NV) DIMENSION VEC(0:NDCLRV,0:MDCLRV) C N=NS DO 7 IV=1,NV DO 7 I=1,NS ICYCLST(I,IV)=-1 ICYCLNO(I,IV)=-1 VCOEF(I,IV)=0.D0 7 CONTINUE C DO 100 IV=1,NV DO 30 I=0,NSB -1 IBX=IBASE(I) ITH=NH(I) ITL=NL(I) LH=NSTAT(ITH,1) LL=NSTAT(ITL,0) C *VDIR LOOPCNT=5196627 DO 20 J=0,LH*LL -1 IF(ABS(VEC(IBX+J,ISS)).LE.1.D-9) THEN GOTO 20 ENDIF DO 10 IIVV=1,IV-1 DO 10 IIV=1,N IF(IBX+J.EQ.ICYCLNO(IIV,IIVV))GOTO 20 10 CONTINUE ITP=ISTAT1(J/LL,ITH)*3**NSLMAX+ISTAT1(MOD(J,LL), & ITL) ICYCL=ITP*3 ICYCL=MOD(ICYCL,3**NS)+ICYCL/3**NS IF(ICYCL.EQ.ITP)THEN GOTO 20 ELSE JBX=IBX+J GOTO 40 ENDIF 20 CONTINUE 30 CONTINUE WRITE(6,*) ' #(W22)# I EXCEEDS IDIM IN CHKWV.' RETURN 40 CONTINUE ICYCLNO(1,IV)=JBX ICYCLST(1,IV)=ITP VCOEF(1,IV)=VEC(JBX,ISS) DO 70 K=1,NS-1 ICYCL=ITP*3 ICYCL=MOD(ICYCL,3**NS)+ICYCL/3**NS ISTL=MOD(ICYCL,3**NSLMAX) ISTH=ICYCL/3**NSLMAX ITCH=ITOTBZ(ISTH) ITCL=ITOTBZ(ISTL) DO 50 I=0,NSB-1 IF(NL(I).EQ.ITCL.AND.NH(I).EQ.ITCH)THEN GOTO 60 ENDIF 50 CONTINUE WRITE(6,*) ' #(W23)# ICYCL IS NOT FOUND IN CHKWV.' GOTO 100 60 CONTINUE IBX=IBASE(I) ITH=NH(I) ITL=NL(I) LH=NSTAT(ITH,1) LL=NSTAT(ITL,0) NEWCFG=IBX+LL*ISTAT2(ISTH)+ISTAT2(ISTL) ICYCLNO(K+1,IV)=NEWCFG ICYCLST(K+1,IV)=ICYCL VCOEF(K+1,IV)=VEC(NEWCFG,ISS) ITP=ICYCL 70 CONTINUE 100 CONTINUE C RETURN END C C SUBROUTIN FOR FINDING THE STATE NUMBER OF Siz=0 FOR ALL i C SUBROUTINE FIND0(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,IST0,ISTN0) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) C IF(NH(0)+NL(0).NE.NS)THEN WRITE(6,*) ' #(W24)# TOTAL SZ IS NOT ZERO IN FIND0.' RETURN ENDIF IHSHFT=3**NSLMAX C IST0=0 DO 20 I=0,NS-1 IST0=1+3*IST0 20 CONTINUE ISTL=MOD(IST0,IHSHFT) ISTH=IST0/IHSHFT CALL FINDHL(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,ISTH,ISTL,ISTN0) C RETURN END C C SUBROUTIN FOR FINDING STATE NUMBER OF (H=ISTH, L=ISTL) C SUBROUTINE FINDHL(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,ISTH,ISTL,ISTN) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) C I=NS I=ISTAT2(1) IF(ISTH.LT.0.OR.ISTH.GT.3**(NS-NSLMAX)-1.OR. & ISTL.LT.0.OR.ISTL.GT.3**NSLMAX-1) THEN WRITE (6,*) & ' #(W25)# WRONG VALUE(S) GIVEN TO ISTH AND/OR ISTL ', & 'IN FINDHL.' WRITE (6,*) & ' ISTH=',ISTH,' ISTL=',ISTL RETURN END IF IF(NH(0)+NL(0).NE.ITOTBZ(ISTH)+ITOTBZ(ISTL))THEN WRITE(6,*) ' #(W26)# TOTAL SZ IS INCONSISTENT IN FINDHL.' RETURN ENDIF C DO 70 I=0,NSB-1 IBH=NH(I) IBL=NL(I) LH=NSTAT(IBH,1) LL=NSTAT(IBL,0) IBX=IBASE(I) C DO 60 J=0,LH - 1 IX=IBX+J*LL ITH=ISTAT1(J,IBH) IF(ISTH.NE.ITH) THEN GOTO 60 ENDIF C *VDIR LOOPCNT=73789 DO 50 K=0,LL - 1 ITL=ISTAT1(K,IBL) IF(ITL.NE.ISTL)THEN GOTO 50 ELSE ISTN=IX+K RETURN ENDIF 50 CONTINUE 60 CONTINUE 70 CONTINUE WRITE(6,*) ' #(W27)# ISTN IS NOT FOUND IN FINDHL.' C RETURN END C c****************************************************************** c* From TITPACK Version 2 * c* Copyright (C) Hidetoshi Nishimori * c* Thanks to Professor H. Nishimori * c****************************************************************** c c************** eigenvalues by the bisection method *************** c subroutine bisec(alpha,beta,ndim,E,ne,eps) c c alpha @ diagonal element c beta @ subdiagonal element c ndim @ matrix dimension c E # eigenvalues c ne @ number of eigenvalues to calculate c eps @ limit of error implicit real*8 (a-h,o-z) dimension alpha(ndim),beta(ndim),E(ne),b2(2000) c if(ndim.gt.2000)then print *,' #(E10)# NDIM GIVEN TO BISEC, CALLED ORIGINALLY ', & 'FROM SMALL, EXCEEDS 2000.' stop end if if(ne.gt.ndim.or.ne.le.0)then print *,' #(E11)# NE GIVEN TO BISEC, CALLED ORIGINALLY ', & 'FROM SMALL, OUT OF RANGE.' stop end if c c*** initial bound range=abs(alpha(1))+abs(beta(1)) do 10 k=2,ndim-1 10 range=max(range,abs(beta(k-1))+abs(alpha(k))+abs(beta(k)) & ) range=max(range,abs(beta(ndim-1))+abs(alpha(ndim))) range=-range c b2(1)=0.d0 do 20 i=2,ndim 20 b2(i)=beta(i-1)**2 c epsabs=abs(range)*eps do 30 i=1,ne 30 E(i)=-range b=range c c*** bisection method do 100 k=1,ne a=E(k) do 110 j=1,100 c=(a+b)/2.d0 if(abs(a-b).lt.epsabs)goto 100 numneg=0 g=1.d0 ipass=0 do 120 i=1,ndim if(ipass.eq.0)then g=c-alpha(i)-b2(i)/g else if(ipass.eq.1)then ipass=2 else g=c-alpha(i) ipass=0 end if c if(ipass.eq.0)then if(g.le.0.d0)numneg=numneg+1 if(abs(g).le.abs(b2(i)*epsabs*eps))ipass=1 end if 120 continue numneg=ndim-numneg if(numneg.lt.k)then b=c else a=c do 130 i=k,min(numneg,ne) 130 E(i)=c end if 110 continue 100 continue return end C C****************************************************************** C* KOBEPACK/1 VERSION 1.0 * C* MARCH 1992 * C* COPYRIGHT (C) M. KABURAGI, T. NISHINO AND T. TONEGAWA * C* WITH THE HELP OF TITPACK VERSION 2 * C* FEBRUARY 1991 * C* COPYRIGHT (C) HIDETOSHI NISHIMORI * C****************************************************************** C***** CORRELATION FUNCTION ***** C****************************************************************** C* CORRESPONDENCE BETWEEN ISTAT1 AND SZ(I) * C* 0 === SZ(I)=-1 * C* 1 === SZ(I)= 0 * C* 2 === SZ(I)=+1 * C****************************************************************** C*** VARIABLES MARKED @ SHOULD BE GIVEN FROM THE MAIN PROGRAM C*** VARIABLES MARKED # ARE EVALUATED AND RETURNED C*** THE FOLLOWING VARIABLES ARE COMMON TO ALL ROUTINES C NS @ LATTICE SIZE C NSLMAX @ MAXIMUM SIZE OF LOWER PART OF LATTICE C ITSZ @ TOTAL SZ C MSGBC @ MESSAGE OF BOUNDARY CONDITION C NSB,NH,NL, @ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C NSTAT, @ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C ISTAT1,NDCLR1,@ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C ISTAT2,IBASE @ LISTS OF HIGH-LOW EXPRESSION BY NISHINO C ISTAT1(I,IBZ) @ I-TH SPIN CONFIGURATION FOR IBZL(OR IBZH)=IB C BZ=SZ+1 C ISTAT2 @ INVERSE LIST OF ISTAT1 EXPRESSED BY THE C 2-DIM SEARCH ALGORITHM OF NISHINO C IPAIR @ PAIRS OF SITES CONNECTED BY BONDS C AJX @ EXCHANGE INTERACTION OF EACH BOND JX,Y C AJZ @ EXCHANGE INTERACTION OF EACH BOND JZ C ABIQD @ BIQUADRATIC INTERCTION C IBOND @ NUMBER OF BONDS C ANISTRPY @ SINGLE ION ANISOTROPY ENERGY AT EACH SITE C HFIELD @ MAGNETIC FIELD AT EACH SITE C DIAG @ DIAGONAL ELEMENT C C************** SPIN AND T-T CORRELATION FUNCTION 1 *************** C*** CALCULATION OF (SZIAV,SZJAV,TPIAV,TPJAV,CFSXIJ,CFSZIJ,CFTPPIJ C*** FOR GIVEN (ISITE,JSITE) C SUBROUTINE CRFNST1(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,SZIAV,SZJAV, &TPIAV,TPJAV,CFSXIJ,CFSZIJ,CFTPPIJ) C C NS @ LATTICE SIZE C NSLMAX @ MAXIMUM SIZE OF LOWER PART OF LATTICE C ISITE @ 1ST SITE TO BE CALCULTED C JSITE @ 2ND SITE TO BE CALCULTED C VEC @ EIGENVETOR C SZIAV # C SZJAV # C TPIAV # = C TPJAV # = C CFSXIJ # -* CORRELATION FUNCTION C CFSZIJ # -* CORRELATION FUNCTION C CFTPPIJ # -* CORRELATION FUNCTION C ISTAT1,ISTAT2 @ SPIN CONFIGURATIONS GENERATED IN 'CRESZ' C & IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION VEC(0:NDCLRV,0:MDCLRV) C I1=JSITE-1 I3=ISITE-1 IF(I1.LT.0.OR.I1.GE.NS.OR.I3.LT.0.OR.I3.GE.NS)THEN PRINT *,' #(W28)# WRONG SITE NUMBER GIVEN TO CRFNST1.' RETURN END IF C C CALL AVRGSI(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &VEC,NDCLRV,MDCLRV,ISS,ISITE,SZIAV) C CALL AVRGSI(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &VEC,NDCLRV,MDCLRV,ISS,JSITE,SZJAV) C CALL AVRGSISJ(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,SXIJAV,SZIJAV) CFSXIJ=SXIJAV CFSZIJ=SZIJAV-SZIAV*SZJAV C CALL AVRGTI(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,TPIAV) C CALL AVRGTI(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,VEC,NDCLRV,MDCLRV,ISS,JSITE,TPJAV) C CALL AVRGTITJ(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,TPIJAV) CFTPPIJ=TPIJAV-TPIAV*TPJAV C RETURN END C C************** SPIN AND T-T CORRELATION FUNCTION 2 *************** C*** CALCULATION OF (CFSXIJ,CFSZIJ,CFTPPIJ) C*** FOR GIVEN (SZIAV,SZJAV,TPIAV,TPJAV,ISITE,JSITE) C SUBROUTINE CRFNST2(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,SZIAV,SZJAV, &TPIAV,TPJAV,CFSXIJ,CFSZIJ,CFTPPIJ) C C NS @ LATTICE SIZE C NSLMAX @ MAXIMUM SIZE OF LOWER PART OF LATTICE C ISITE @ 1ST SITE TO BE CALCULTED C JSITE @ 2ND SITE TO BE CALCULTED C VEC @ EIGENVETOR C SZIAV @ C SZJAV @ C TPIAV @ = C TPJAV @ = C CFSXIJ # -* CORRELATION FUNCTION C CFSZIJ # -* CORRELATION FUNCTION C CFTPPIJ # -* CORRELATION FUNCTION C ISTAT1,ISTAT2 @ SPIN CONFIGURATIONS GENERATED IN 'CRESZ' C & IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION VEC(0:NDCLRV,0:MDCLRV) C I1=JSITE-1 I3=ISITE-1 IF(I1.LT.0.OR.I1.GE.NS.OR.I3.LT.0.OR.I3.GE.NS)THEN PRINT *,' #(W28)# WRONG SITE NUMBER GIVEN TO CRFNST2.' RETURN END IF C CALL AVRGSISJ(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,SXIJAV,SZIJAV) CFSXIJ=SXIJAV CFSZIJ=SZIJAV-SZIAV*SZJAV C CALL AVRGTITJ(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,TPIJAV) CFTPPIJ=TPIJAV-TPIAV*TPJAV C RETURN END C C****************** SPIN CORRELATION FUNCTION 1 ******************* C*** CALCULATION OF (SZIAV,SZJAV,CFSXIJ, CFSZIJ) C*** FOR GIVEN (ISITE,JSITE) C SUBROUTINE CRFNS1(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,SZIAV,SZJAV, &CFSXIJ,CFSZIJ) C C NS @ LATTICE SIZE C NSLMAX @ MAXIMUM SIZE OF LOWER PART OF LATTICE C ISITE @ 1ST SITE TO BE CALCULTED C JSITE @ 2ND SITE TO BE CALCULTED C VEC @ EIGENVETOR C SZIAV # C SZJAV # C CFSXIJ # -* CORRELATION FUNCTION C CFSZIJ # -* CORRELATION FUNCTION C ISTAT1,ISTAT2 @ SPIN CONFIGURATIONS GENERATED IN 'CRESZ' C & IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION VEC(0:NDCLRV,0:MDCLRV) C I1=JSITE-1 I3=ISITE-1 IF(I1.LT.0.OR.I1.GE.NS.OR.I3.LT.0.OR.I3.GE.NS)THEN PRINT *,' #(W28)# WRONG SITE NUMBER GIVEN TO CRFNS1.' RETURN END IF C CALL AVRGSI(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &VEC,NDCLRV,MDCLRV,ISS,ISITE,SZIAV) C CALL AVRGSI(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &VEC,NDCLRV,MDCLRV,ISS,JSITE,SZJAV) C CALL AVRGSISJ(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,SXIJAV,SZIJAV) CFSXIJ=SXIJAV CFSZIJ=SZIJAV-SZIAV*SZJAV C RETURN END C C****************** SPIN CORRELATION FUNCTION 2 ******************* C*** CALCULATION OF (CFSXIJ,CFSZIJ) C*** FOR GIVEN (ISITE,JSITE,SZIAV,SZJAV) C SUBROUTINE CRFNS2(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,SZIAV,SZJAV, &CFSXIJ,CFSZIJ) C C NS @ LATTICE SIZE C NSLMAX @ MAXIMUM SIZE OF LOWER PART OF LATTICE C ISITE @ 1ST SITE TO BE CALCULTED C JSITE @ 2ND SITE TO BE CALCULTED C VEC @ EIGENVETOR C SZIAV @ C SZJAV @ C CFSXIJ # -* CORRELATION FUNCTION C CFSZIJ # -* CORRELATION FUNCTION C ISTAT1,ISTAT2 @ SPIN CONFIGURATIONS GENERATED IN 'SZ' C & IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION VEC(0:NDCLRV,0:MDCLRV) C I1=JSITE-1 I3=ISITE-1 IF(I1.LT.0.OR.I1.GE.NS.OR.I3.LT.0.OR.I3.GE.NS)THEN PRINT *,' #(W28)# WRONG SITE NUMBER GIVEN TO CRFNS2.' RETURN END IF C CALL AVRGSISJ(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,SXIJAV,SZIJAV) CFSXIJ=SXIJAV CFSZIJ=SZIJAV-SZIAV*SZJAV C RETURN END C C* CALCULATION AVRAGE(SZI) C SUBROUTINE AVRGSI(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,VEC,NDCLRV,MDCLRV,ISS,ISITE,SZIAV) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX) REAL*8 VEC(0:NDCLRV,0:MDCLRV) C I3=ISITE-1 IF(I3.LT.0.OR.I3.GE.NS)THEN PRINT *,' #(W28)# WRONG SITE NUMBER GIVEN TO AVRGSI.' RETURN END IF C SZI= 0.0D0 C DO 90 I=0,NSB -1 IBX=IBASE(I) ITH=NH(I) ITL=NL(I) LH=NSTAT(ITH,1) LL=NSTAT(ITL,0) MSK=3**I3 DO 70 J=0,LH -1 IX=IBX+J*LL *VDIR LOOPCNT=73789 DO 60 K=0,LL-1 ITP=ISTAT1(J,ITH)*3**NSLMAX+ISTAT1(K,ITL) SZI=SZI+(MOD(ITP/MSK,3)-1)*VEC(IX+K,ISS)*VEC(IX+K, & ISS) 60 CONTINUE 70 CONTINUE 90 CONTINUE SZIAV=SZI C RETURN END C C* CALCULATION OF AVERAGE(SI*SJ) C SUBROUTINE AVRGSISJ(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,SXIJAV, &SZIJAV) C PARAMETER(IBND=1,NDCLRM=0) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBND+12) DIMENSION AJX(IBND+6),AJZ(IBND+6),ABIQD(IBND+6) DIMENSION VEC(0:NDCLRV,0:MDCLRV) DIMENSION LBOND(IBND+6,0:5),NBOND(5) DIMENSION MBASE(0:24,IBND+16,6) DIMENSION IBSPAIR(2*IBND+12),IBSPIN(IBND+6,2),NBSPIN(2) DIMENSION IFLGOELM(0:24,IBND+16) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) CHARACTER*4 MSGBC(2) C I1=JSITE-1 I3=ISITE-1 IF(I1.LT.0.OR.I1.GE.NS.OR.I3.LT.0.OR.I3.GE.NS)THEN PRINT *,' #(W28)# WRONG SITE NUMBER GIVEN TO AVRGSISJ.' RETURN END IF ZIJ= 0.0D0 XIJ= 0.0D0 C IF(I1.NE.I3)THEN C IS=I1 ID=I3 C DO 70 I=0,NSB-1 ITH=NH(I) ITL=NL(I) LH=NSTAT(ITH,1) LL=NSTAT(ITL,0) IBX=IBASE(I) MSKH=3**ID MSKL=3**IS C DO 60 J=0,LH - 1 IX=IBX+J*LL *VDIR LOOPCNT=73789 DO 50 K=0,LL - 1 ITP=ISTAT1(J,ITH)*3**NSLMAX+ISTAT1(K,ITL) ZIJ= ZIJ+ & DBLE((MOD(ITP/MSKH,3)-1)*(MOD(ITP/MSKL,3)-1)) & *VEC(IX+K,ISS)*VEC(IX+K,ISS) 50 CONTINUE 60 CONTINUE 70 CONTINUE C MSGBC(1)='PERI' IDIM=IBASE(NSB) IBOND=IBND IPAIR(1)=I1+1 IPAIR(2)=I3+1 AJX(1)=1.D0 AJZ(1)=0.D0 ABIQD(1)=0.D0 ID=MDCLRV IS=0 C CALL CLSBND(NS,NSLMAX,IPAIR,IBOND,LBOND,NBOND,IBSPAIR, & IBSPIN, & NBSPIN,MSGBC) C CALL ELMO(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, & ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,MSGBC,LBOND,NBOND, & MBASE,IFLGOELM,IBSPAIR,IBSPIN,NBSPIN,ELMNT,LOC,NDCLRM) C CALL VEC01(IDIM,VEC(0,ID)) CALL MLTO(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, & ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND, & VEC(0,ID),VEC(0,IS),IDIM,LBOND,NBOND,MBASE,IFLGOELM, & ELMNT,LOC,NDCLRM) C CALL INPRO(IDIM,VEC(0,ID),VEC(0,IS),PROD) XIJ= PROD C ELSE IS=I1 C DO 130 I=0,NSB-1 ITH=NH(I) ITL=NL(I) LH=NSTAT(ITH,1) LL=NSTAT(ITL,0) IBX=IBASE(I) MSKL=3**IS C DO 120 J=0,LH - 1 IX=IBX+J*LL *VDIR LOOPCNT=73789 DO 110 K=0,LL - 1 ITP=ISTAT1(J,ITH)*3**NSLMAX+ISTAT1(K,ITL) ZIJ= ZIJ+DBLE((MOD(ITP/MSKL,3)-1)**2) & *VEC(IX+K,ISS)*VEC(IX+K,ISS) XIJ= XIJ+DBLE(2-(MOD(ITP/MSKL,3)-1)**2) & *VEC(IX+K,ISS)*VEC(IX+K,ISS) 110 CONTINUE 120 CONTINUE 130 CONTINUE C ENDIF C SZIJAV=ZIJ SXIJAV=XIJ/2.D0 C RETURN END C C****************** TP-TP CORRELATION FUNCTION 1 ****************** C*** CALCULATION OF (TPIAV,TPJAV,CFTPPIJ) C*** FOR GIVEN (ISITE,JSITE) C SUBROUTINE CRFNT1(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,TPIAV,TPJAV, &CFTPPIJ) C C NS @ LATTICE SIZE C NSLMAX @ MAXIMUM SIZE OF LOWER PART OF LATTICE C ISITE @ 1ST SITE TO BE CALCULTED C JSITE @ 2ND SITE TO BE CALCULTED C VEC @ EIGENVETOR C TPIAV # = C TPJAV # = C CFTPPIJ # -**2 CORRELATION FUNCTION C ISTAT1,ISTAT2 @ SPIN CONFIGURATIONS GENERATED IN 'CRESZ' C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION VEC(0:NDCLRV,0:MDCLRV) C I1=JSITE-1 I3=ISITE-1 IF(I1.LT.0.OR.I1.GE.NS.OR.I3.LT.0.OR.I3.GE.NS)THEN PRINT *,' #(W28)# WRONG SITE NUMBER GIVEN TO CRFNT1.' RETURN END IF C CALL AVRGTI(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &VEC,NDCLRV,MDCLRV,ISS,ISITE,TPIAV) C CALL AVRGTI(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &VEC,NDCLRV,MDCLRV,ISS,JSITE,TPJAV) C CALL AVRGTITJ(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,TPIJAV) CFTPPIJ=TPIJAV-TPIAV*TPJAV C RETURN END C C****************** TP-TP CORRELATION FUNCTION 2 ****************** C*** CALCULATION OF (CFTPPIJ) C*** FOR GIVEN (ISITE,JSITE,TPIAV,TPJAV) C SUBROUTINE CRFNT2(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,TPIAV,TPJAV, &CFTPPIJ) C C NS @ LATTICE SIZE C NSLMAX @ MAXIMUM SIZE OF LOWER PART OF LATTICE C ISITE @ 1ST SITE TO BE CALCULTED C JSITE @ 2ND SITE TO BE CALCULTED C VEC @ EIGENVETOR C TPIAV @ = C TPJAV @ = C CFTPPIJ # -**2 CORRELATION FUNCTION C ISTAT1,ISTAT2 @ SPIN CONFIGURATIONS GENERATED IN 'CRESZ' C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION VEC(0:NDCLRV,0:MDCLRV) C I1=JSITE-1 I3=ISITE-1 IF(I1.LT.0.OR.I1.GE.NS.OR.I3.LT.0.OR.I3.GE.NS)THEN PRINT *,' #(W28)# WRONG SITE NUMBER GIVEN TO CRFNT2.' RETURN END IF C CALL AVRGTITJ(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,TPIJAV) CFTPPIJ=TPIJAV-TPIAV*TPJAV C RETURN END C C CALCULATION AVRAGE(TPI) C SUBROUTINE AVRGTI(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,TPIAV) C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(IBND=1,NDCLRM=0) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBND+12) DIMENSION AJX(IBND+6),AJZ(IBND+6),ABIQD(IBND+6) DIMENSION VEC(0:NDCLRV,0:MDCLRV) DIMENSION LBOND(IBND+6,0:5),NBOND(5) DIMENSION MBASE(0:24,IBND+16,6) DIMENSION IBSPAIR(2*IBND+12),IBSPIN(IBND+6,2),NBSPIN(2) DIMENSION IFLGOELM(0:24,IBND+16) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) CHARACTER*4 MSGBC(2) C I3=ISITE-1 IF(I3.LT.0.OR.I3.GE.NS)THEN PRINT *,' #(W28)# WRONG SITE NUMBER GIVEN TO AVRGTI.' RETURN END IF JSITE=MOD(ISITE,NS)+1 C XIJ= 0.0D0 C MSGBC(1)='PERI' IDIM=IBASE(NSB) IBOND=IBND IPAIR(1)=ISITE IPAIR(2)=JSITE AJX(1)=1.D0 AJZ(1)=0.D0 ABIQD(1)=0.D0 IS=ISS ID=MDCLRV C CALL CLSBND(NS,NSLMAX,IPAIR,IBOND,LBOND,NBOND,IBSPAIR, &IBSPIN,NBSPIN,MSGBC) C CALL ELMO(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,MSGBC,LBOND,NBOND,MBASE, &IFLGOELM,IBSPAIR,IBSPIN,NBSPIN,ELMNT,LOC,NDCLRM) C CALL VEC01(IDIM,VEC(0,ID)) CALL MLTO(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VEC(0,ID),VEC(0,IS),IDIM, &LBOND,NBOND,MBASE,IFLGOELM,ELMNT,LOC,NDCLRM) C CALL INPRO(IDIM,VEC(0,ID),VEC(0,IS),PROD) XIJ=PROD C TPIAV=2.D0*XIJ C RETURN END C C************************ AVERAGE OF TITJ ************************* C*** CALCULATION OF AVRG(TPI*TPJ) FOR GIVEN (ISITE,JSITE) C SUBROUTINE AVRGTITJ(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,TPIJAV) C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(IBND=1,NDCLRM=0) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION IPAIR(2*IBND+12) DIMENSION AJX(IBND+6),AJZ(IBND+6),ABIQD(IBND+6) DIMENSION VEC(0:NDCLRV,0:MDCLRV) DIMENSION LBOND(IBND+6,0:5),NBOND(5) DIMENSION MBASE(0:24,IBND+16,6) DIMENSION IBSPAIR(2*IBND+12),IBSPIN(IBND+6,2),NBSPIN(2) DIMENSION IFLGOELM(0:24,IBND+16) DIMENSION ELMNT(0:NDCLRM),LOC(0:NDCLRM) CHARACTER*4 MSGBC(2) C I1=JSITE-1 I3=ISITE-1 IF(I1.LT.0.OR.I1.GE.NS.OR.I3.LT.0.OR.I3.GE.NS)THEN PRINT *,' #(W28)# WRONG SITE NUMBER GIVEN TO AVRGTITJ.' RETURN END IF PIJ= 0.0D0 C MSGBC(1)='PERI' IDIM=IBASE(NSB) IBOND=IBND IPAIR(1)=ISITE IPAIR(2)=MOD(ISITE,NS)+1 AJX(1)=1.D0 AJZ(1)=0.D0 ABIQD(1)=0.D0 IS=ISS ID=MDCLRV C C & CALL CLSBND(NS,NSLMAX,IPAIR,IBOND,LBOND,NBOND,IBSPAIR, &IBSPIN,NBSPIN,MSGBC) C CALL ELMO(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,MSGBC,LBOND,NBOND,MBASE, &IFLGOELM,IBSPAIR,IBSPIN,NBSPIN,ELMNT,LOC,NDCLRM) C CALL VEC01(IDIM,VEC(0,ID)) CALL MLTO(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VEC(0,ID),VEC(0,IS),IDIM, &LBOND,NBOND,MBASE,IFLGOELM,ELMNT,LOC,NDCLRM) C IPAIR(1)=JSITE IPAIR(2)=MOD(JSITE,NS)+1 IS=MDCLRV ID=MDCLRV-1 C CALL CLSBND(NS,NSLMAX,IPAIR,IBOND,LBOND,NBOND,IBSPAIR, &IBSPIN,NBSPIN,MSGBC) C CALL ELMO(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,MSGBC,LBOND,NBOND,MBASE, &IFLGOELM,IBSPAIR,IBSPIN,NBSPIN,ELMNT,LOC,NDCLRM) C CALL VEC01(IDIM,VEC(0,ID)) CALL MLTO(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1,NDCLR1, &ISTAT2,IPAIR,AJX,AJZ,ABIQD,IBOND,VEC(0,ID),VEC(0,IS),IDIM, &LBOND,NBOND,MBASE,IFLGOELM,ELMNT,LOC,NDCLRM) C IS=ISS ID=MDCLRV-1 CALL INPRO(IDIM,VEC(0,ID),VEC(0,IS),PROD) TPIJAV=PROD*4.D0 C RETURN END C C C SUBROUTINE STRODRZ(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,OSTZIJ) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION VEC(0:NDCLRV,0:MDCLRV) C NSL12=12 I=ISTAT2(1) I1=JSITE-1 I3=ISITE-1 IF(JSITE.GT.ISITE)THEN I1=ISITE-1 I3=JSITE-1 ENDIF IF(I1.LT.0.OR.I1.GE.NS.OR.I3.LT.0.OR.I3.GE.NS)THEN PRINT *,' #(W28)# WRONG SITE NUMBER GIVEN TO STRODRZ.' RETURN END IF C IHSHFT=3**NSLMAX IHSFT12=3**NSL12 STZIJ= 0.0D0 C IS=I1 ID=I3 IF(ID.GT.IS)THEN C DO 70 I=0,NSB-1 ITH=NH(I) ITL=NL(I) LH=NSTAT(ITH,1) LL=NSTAT(ITL,0) IBX=IBASE(I) MSKH=3**ID MSKL=3**IS C DO 60 J=0,LH - 1 IX=IBX+J*LL *VDIR LOOPCNT=73789 DO 50 K=0,LL - 1 ITP=ISTAT1(J,ITH)*IHSHFT+ISTAT1(K,ITL) IT=MOD((ITP/(3*MSKL))*(3*MSKL),MSKH) IB=MOD(IT/3** 0,3)+MOD(IT/3** 1,3) & +MOD(IT/3** 2,3)+MOD(IT/3** 3,3) & +MOD(IT/3** 4,3)+MOD(IT/3** 5,3) & +MOD(IT/3** 6,3)+MOD(IT/3** 7,3) & +MOD(IT/3** 8,3)+MOD(IT/3** 9,3) & +MOD(IT/3**10,3)+MOD(IT/3**11,3) IF(NS.LE.NSL12) THEN GOTO 40 ENDIF IT=IT/IHSFT12 IB=IB+MOD(IT/3** 0,3)+MOD(IT/3** 1,3) & +MOD(IT/3** 2,3)+MOD(IT/3** 3,3) & +MOD(IT/3** 4,3)+MOD(IT/3** 5,3) & +MOD(IT/3** 6,3)+MOD(IT/3** 7,3) & +MOD(IT/3** 8,3)+MOD(IT/3** 9,3) & +MOD(IT/3**10,3)+MOD(IT/3**11,3) 40 CONTINUE IB=IB+(ID-IS-1) STZIJ=STZIJ+DBLE(((-1)**IB) & *(MOD(ITP/MSKH,3)-1)*(MOD(ITP/MSKL,3)-1)) & *VEC(IX+K,ISS)*VEC(IX+K,ISS) 50 CONTINUE 60 CONTINUE 70 CONTINUE C ELSE C DO 100 I=0,NSB-1 ITH=NH(I) ITL=NL(I) LH=NSTAT(ITH,1) LL=NSTAT(ITL,0) IBX=IBASE(I) MSKH=3**ID MSKL=3**IS C DO 90 J=0,LH - 1 IX=IBX+J*LL *VDIR LOOPCNT=73789 DO 80 K=0,LL - 1 ITP=ISTAT1(J,ITH)*IHSHFT+ISTAT1(K,ITL) STZIJ=STZIJ+ & DBLE((MOD(ITP/MSKH,3)-1)*(MOD(ITP/MSKL,3)-1)) & *VEC(IX+K,ISS)*VEC(IX+K,ISS) 80 CONTINUE 90 CONTINUE 100 CONTINUE C ENDIF C OSTZIJ=STZIJ C RETURN END C C C SUBROUTINE STRODRX(NS,NSLMAX,NSB,NH,NL,IBASE,NSTAT,ISTAT1, &NDCLR1,ISTAT2,VEC,NDCLRV,MDCLRV,ISS,ISITE,JSITE,OSTXIJ) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IBASE(0:2*NSLMAX+1),NH(0:2*NSLMAX),NL(0:2*NSLMAX) DIMENSION NSTAT(0:2*NSLMAX,0:1) DIMENSION IA(0:24),IC(0:24),IQ(0:24),IR(0:24) DIMENSION ISTAT1(0:NDCLR1,0:2*NSLMAX),ISTAT2(0:3**NSLMAX) DIMENSION VEC(0:NDCLRV,0:MDCLRV) C NSL12=12 I=ISTAT2(1) I1=JSITE-1 I3=ISITE-1 IF(JSITE.GT.ISITE)THEN I1=ISITE-1 I3=JSITE-1 ENDIF IF(I1.LT.0.OR.I1.GE.NS.OR.I3.LT.0.OR.I3.GE.NS)THEN PRINT *,' #(W28)# WRONG SITE NUMBER GIVEN TO STRODRX.' RETURN END IF C IHSHFT=3**NSLMAX IHSFT12=3**NSL12 STXIJ= 0.0D0 C IS=I1 ID=I3 IF(ID.GT.IS)THEN C DO 20 I=0,22 IQ(I)=1 IR(I)=0 IC(I)=1 IF((I-ID)*(I-IS).LT.0)THEN IQ(I)=2 IR(I)=2 IC(I)=-1 ENDIF 20 CONTINUE ICC= IC(22)*IC(21)*IC(20)*IC(19)*IC(18)*IC(17) & *IC(16)*IC(15)*IC(14)*IC(13)*IC(12) ICC=ICC*IC(11)*IC(10)*IC( 9)*IC( 8)*IC( 7)*IC( 6)*IC( 5) & *IC( 4)*IC( 3)*IC( 2)*IC( 1)*IC( 0) C DO 70 I=0,NSB-1 IBH=NH(I) IBL=NL(I) LH=NSTAT(IBH,1) LL=NSTAT(IBL,0) IBX=IBASE(I) MSKH=3**ID MSKL=3**IS C DO 60 J=0,LH - 1 IX=IBX+J*LL ITH=ISTAT1(J,IBH) *VDIR LOOPCNT=73789 DO 50 K=0,LL - 1 ITL=ISTAT1(K,IBL) IP=ITH*IHSHFT+ITL IA( 0)=MOD(IQ( 0)*MOD(IP/3** 0,3)+IR( 0),3) IA( 1)=MOD(IQ( 1)*MOD(IP/3** 1,3)+IR( 1),3) IA( 2)=MOD(IQ( 2)*MOD(IP/3** 2,3)+IR( 2),3) IA( 3)=MOD(IQ( 3)*MOD(IP/3** 3,3)+IR( 3),3) IA( 4)=MOD(IQ( 4)*MOD(IP/3** 4,3)+IR( 4),3) IA( 5)=MOD(IQ( 5)*MOD(IP/3** 5,3)+IR( 5),3) IA( 6)=MOD(IQ( 6)*MOD(IP/3** 6,3)+IR( 6),3) IA( 7)=MOD(IQ( 7)*MOD(IP/3** 7,3)+IR( 7),3) IA( 8)=MOD(IQ( 8)*MOD(IP/3** 8,3)+IR( 8),3) IA( 9)=MOD(IQ( 9)*MOD(IP/3** 9,3)+IR( 9),3) IA(10)=MOD(IQ(10)*MOD(IP/3**10,3)+IR(10),3) IA(11)=MOD(IQ(11)*MOD(IP/3**11,3)+IR(11),3) IA(12)=MOD(IQ(12)*MOD(IP/3**12,3)+IR(12),3) IA(13)=MOD(IQ(13)*MOD(IP/3**13,3)+IR(13),3) IA(14)=MOD(IQ(14)*MOD(IP/3**14,3)+IR(14),3) IA(15)=MOD(IQ(15)*MOD(IP/3**15,3)+IR(15),3) IA(16)=MOD(IQ(16)*MOD(IP/3**16,3)+IR(16),3) IA(17)=MOD(IQ(17)*MOD(IP/3**17,3)+IR(17),3) IA(18)=MOD(IQ(18)*MOD(IP/3**18,3)+IR(18),3) IA(19)=MOD(IQ(19)*MOD(IP/3**19,3)+IR(19),3) IA(20)=MOD(IQ(20)*MOD(IP/3**20,3)+IR(20),3) IA(21)=MOD(IQ(21)*MOD(IP/3**21,3)+IR(21),3) IA(22)=MOD(IQ(22)*MOD(IP/3**22,3)+IR(22),3) INB=IA( 0)+IA( 1)+IA( 2)+IA( 3) & +IA( 4)+IA( 5)+IA( 6)+IA( 7) & +IA( 8)+IA( 9)+IA(10)+IA(11) & +IA(12)+IA(13)+IA(14)+IA(15) & +IA(16)+IA(17)+IA(18)+IA(19) & +IA(20)+IA(21)+IA(22) INP=IA(0)+3*(IA(1)+3*(IA(2)+3*(IA(3) & +3*(IA(4)+3*(IA(5)+3*(IA(6)+3*(IA(7) & +3*(IA(8)+3*(IA(9)+3*(IA(10)+3*(IA(11) & +3*(IA(12)+3*(IA(13)+3*(IA(14)+3*(IA(15) & +3*(IA(16)+3*(IA(17)+3*(IA(18)+3*(IA(19) & +3*(IA(20)+3*(IA(21)+3*IA(22)) & )))))))))))))))))))) IF((INB.EQ.IBH+IBL).AND. & (MOD(INP/MSKH,3).NE.0.AND.MOD(INP/MSKL,3).NE.2)) & THEN INP=INP-MSKH+MSKL ELSEIF((INB.EQ.IBH+IBL-2).AND. & (MOD(INP/MSKH,3).NE.2.AND.MOD(INP/MSKL,3) & .NE.2)) THEN INP=INP+MSKH+MSKL ELSE GOTO 50 ENDIF INPL=MOD(INP,IHSHFT) INPH=INP/IHSHFT INBL=MOD(INPL/3** 0,3)+MOD(INPL/3** 1,3) & +MOD(INPL/3** 2,3)+MOD(INPL/3** 3,3) & +MOD(INPL/3** 4,3)+MOD(INPL/3** 5,3) & +MOD(INPL/3** 6,3)+MOD(INPL/3** 7,3) & +MOD(INPL/3** 8,3)+MOD(INPL/3** 9,3) & +MOD(INPL/3**10,3)+MOD(INPL/3**11,3) INBH=MOD(INPH/3** 0,3)+MOD(INPH/3** 1,3) & +MOD(INPH/3** 2,3)+MOD(INPH/3** 3,3) & +MOD(INPH/3** 4,3)+MOD(INPH/3** 5,3) & +MOD(INPH/3** 6,3)+MOD(INPH/3** 7,3) & +MOD(INPH/3** 8,3)+MOD(INPH/3** 9,3) & +MOD(INPH/3**10,3)+MOD(INPH/3**11,3) IBS=INBL-NL(0) IF((IBS+1)*(IBS-NSB).GE.0)THEN WRITE(6,*) & ' #(W29)# IBASE ERROR IN STRODRX. ', & 'IBS=',IBS,', INBH=',INBH, & ', INBL=',INBL RETURN ENDIF IY=IBASE(IBS)+ISTAT2(INPH)*NSTAT(INBL,0)+ & ISTAT2(INPL) STXIJ=STXIJ+DBLE(2*ICC)* & VEC(IX+K,ISS)*VEC(IY,ISS) 50 CONTINUE 60 CONTINUE 70 CONTINUE C ELSE C DO 100 I=0,NSB-1 IBH=NH(I) IBL=NL(I) LH=NSTAT(IBH,1) LL=NSTAT(IBL,0) IBX=IBASE(I) MSKH=3**ID MSKL=3**IS C DO 90 J=0,LH - 1 IX=IBX+J*LL *VDIR LOOPCNT=73789 DO 80 K=0,LL - 1 ITP=ISTAT1(J,IBH)*IHSHFT+ISTAT1(K,IBL) STXIJ=STXIJ+ & DBLE(2-(MOD(ITP/MSKH,3)-1)*(MOD(ITP/MSKL,3)-1)) & *VEC(IX+K,ISS)*VEC(IX+K,ISS) 80 CONTINUE 90 CONTINUE 100 CONTINUE ENDIF C OSTXIJ=STXIJ/2.D0 C RETURN END