This is the source code of DMRG for S=1/2
Heisenberg Chain with Next Nearest
Neighbor Interaction.
This program runs with NEC's ASL library
(Please contact me if you wish to know about ASL.)
This program is really primitive.Use it just for tutorial purpouse.
*#RUN:L=LIB/ASL7 OPT=3 C*#RUN:L=LIB/ASL7 OPT=0 NIAP ARGCHK SUBCHK C123456789012345678901234567890123456789012345678901234567890123456789012 C C Sannomiya Kyuukou Version: as fast as the Hankyu Train. C Version 1.1 1995/1/15 by Nishino Tomotoshi: Free Soft. C C ***** Subroutines whose name begins with "DCS" belong to ***** C ***** Center Library of SX/ACOS (ASL) ***** C C +---- Density Matrix (Numerical) Renormalization Group ----- C ----- for the J-J' Heisenberg Chain. (Ladder mo dekiru!)----- C ----- ________ ________ ----- C ----- / Joint\ / \ ----- C ----- [L-Block===E====J]====O====O====[J====E===R-Block] ----- C ----- End \________/ \________/ ----- C ----- ----- C ----- === J(0), __/ J(1), Ising D(0), XY D(1) ----- C ----- {Internal Expression: J(2), J(3)} ----- C ----- ----- C +---- Special #: 0-1 -- Left/Right Sub System ----+ C | 100 -- Max.# of Site (0 to 99) | C | -1~101 -- Max.# of Sub Space = 100+1 | C | 10000 -- Matrix Stack Space (10Mbyte/Operator)| C | 40000 -- Extended Matrix Stack Space | C +------------ 100000 -- Max. Lanczos Vector Length ----------+ C INTEGER LV, NS, NSB, NL, MX, DX, LP, NP, RT INTEGER TZ, NR, MM(0:16), RP, IP, MRT INTEGER SZ(0:100,0:1), MP(0:100,0:5), X(-1:101,0:3,0:100) REAL*8 E0(0:3,0:16), BH(0:40000,0:1), SM(0:10000,0:5,0:100) REAL*8 T1,TIME(0:16), JH(0:4), E(0:3), HM(0:10000,0:100) REAL*8 FMX, FM(0:16), ER, ERE, V(0:100000,0:3) C C +------------------------------------------------------------+ C | Parameter Input and Initialization of Variables | C +------------------------------------------------------------+ C .......... Debug Level: 0-None, 1-Normal, 2-Detail, 3-Dump, 4-Debug C 5-FullDump-Debug, 6-Stop LV = 0 C .......... Desiable Accuracy of the Density Matrix (Desiable...) ER = 1.0E-15 C .......... Desiable Accuracy of the Ground State Energy ERE = 1.0E-7 C .......... Minimum Step Number in Lanczos Diagonalization C when we consider Excitations. LLS = 120 C .......... Exchange Parameters (N.N. J(0); N.N.N. J(1)) JH(0)= 1.00D0 JH(1)= 1.00D0 C .......... Ising/XY Strength (Ising J(2); XY J(3)) JH(2)= 1.00D0 JH(3)= 1.00D0 C .......... Number of Sites; it should be even number. (Ns = 100) DO 90 NS = 16, 16, 2 C .......... Total Sz (= an Integer, since Ns should be even.) DO 90 TZ = 0, 1 C .......... Maximum Number of Iteration (Max RT) MRT= 10 C .......... Initial "MAXIMUM Number" of Basis Kept MX = 20 C .......... Increase Rate of MX during the Finite Size Iteration DX = 10 C .......... Parameter Display/ Parameter Check CALL PRMDSP( LV, NS, TZ, MX, ER, JH) C .......... Make Initial Blocks: BL [O===O]---, and ---[O===O] BR CALL INIBLK( LV, JH, X(-1,0, 1), SM(0,0, 1), HM(0, 1)) CALL INIBLK( LV, JH, X(-1,0,NS-2), SM(0,0,NS-2), HM(0,NS-2)) C C +------------------------------------------------------------+ C | Iteration of Diagonalization and Renormalizatioin | C +------------------------------------------------------------+ C ...........Infinite Chain Iteration to NS-2: (Rotation=0) RT = 0 CALL GETSEC( T1 ) C ........... Dummy ..... IP = 0 FMX= 0.0D0 C ...........Initial Block Size = 2 for both L/R Blocks. NL = 2 NR = 2 10 CONTINUE LP = NL - 1 RP = NS - NR CALL MATSTR(LV, NL,NR,TZ,NSB, X(-1,0,LP), X(-1,0,RP),SZ, MP) CALL BLOCKH(LV,NSB, SZ, MP, X(-1,0,LP), X(-1,0,RP),JH, & HM(0,LP),HM(0,RP),SM( 0,0,LP),SM( 0,0,RP),BH ) IF(NL .EQ. NS/2-1 ) THEN LS = LLS ELSE LS = 10 END IF CALL LANZOS(LV,NSB,LS,SZ, MP, X(-1,0,LP), X(-1,0,RP),JH, & SM( 0,0,LP),SM( 0,0,RP),BH,E,V) IF( LV .GE. 1 ) WRITE(6,100) NL, NR,E(0), IP, FMX NL = NL + 1 NP = LP + 1 CALL RENRML(LV,NL, NSB,SZ(0,0),MP,X(-1,0,LP),X(-1,0,NP), & SM(0,0,LP),SM(0,0,NP), BH(0,0), HM(0,NP), & MX, IP, ER, FMX, V) NR = NR + 1 NP = RP - 1 CALL RENRMR(LV,NR, NSB,SZ(0,1),MP,X(-1,0,RP),X(-1,0,NP), & SM(0,0,RP),SM(0,0,NP), BH(0,1), HM(0,NP), & MX, IP, ER, FMX, V) IF( LV.GE. 6 ) STOP IF( NL+NR .LT. NS ) GOTO 10 IF( LV.GE. 5 ) STOP C ............ Second Step: Finite Chain Iteration for LEFT Block .... NR = NS - NL - 2 C E0(0,RT)= E(0) E0(1,RT)= E(1) E0(2,RT)= E(2) E0(3,RT)= E(3) MM(RT)= IP FM(RT)= FMX CALL EXTIME( T1, TIME(RT) ) WRITE(6,*) WRITE(6,110) RT, MM(RT), E0(0,RT), E0(1,RT), E0(2,RT), & FM(RT), TIME(RT) RT = RT + 1 MX = MX + DX 20 CONTINUE LP = NL - 1 RP = NS - NR IF( NR .EQ. 2 ) GOTO 30 CALL MATSTR(LV, NL,NR,TZ,NSB, X(-1,0,LP), X(-1,0,RP),SZ, MP) CALL BLOCKH(LV,NSB, SZ, MP, X(-1,0,LP), X(-1,0,RP),JH, & HM(0,LP),HM(0,RP),SM( 0,0,LP),SM( 0,0,RP),BH ) IF(NL .EQ. NS/2-1 ) THEN LS = LLS ELSE LS = 10 END IF CALL LANZOS(LV,NSB,LS,SZ, MP, X(-1,0,LP), X(-1,0,RP),JH, & SM( 0,0,LP),SM( 0,0,RP),BH,E,V) IF( LV .GE. 1 ) WRITE(6,100) NL, NR,E(0), IP, FMX IF( NL .EQ. NS/2-1 ) THEN E0(0,RT)= E(0) E0(1,RT)= E(1) E0(2,RT)= E(2) E0(3,RT)= E(3) MM(RT)= IP FM(RT)= FMX CALL EXTIME( T1, TIME(RT) ) WRITE(6,111) RT, MM(RT), E0(0,RT), E0(1,RT), E0(2,RT), & FM(RT), TIME(RT) C ( Please Do Not Increase MX here! ) IF((ABS(E0(0,RT)-E0(0,RT-1)).LT.ERE .AND. MM(RT).GT.MM(RT-1)) & .OR. RT.GT.MRT) GOTO 40 RT = RT + 1 END IF NL = NL + 1 NP = LP + 1 NR = NR - 1 CALL RENRML(LV,NL, NSB,SZ(0,0),MP,X(-1,0,LP),X(-1,0,NP), & SM(0,0,LP),SM(0,0,NP), BH(0,0), HM(0,NP), & MX, IP, ER, FMX, V) GOTO 20 C ............ Third Step: Finite Chain Iteration for LEFT Block .... 30 CONTINUE LP = NL - 1 RP = NS - NR IF( NL .EQ. 2 ) GOTO 20 CALL MATSTR(LV, NL,NR,TZ,NSB, X(-1,0,LP), X(-1,0,RP),SZ, MP) CALL BLOCKH(LV,NSB, SZ, MP, X(-1,0,LP), X(-1,0,RP),JH, & HM(0,LP),HM(0,RP),SM( 0,0,LP),SM( 0,0,RP),BH ) IF(NL .EQ. NS/2-1 ) THEN LS = LLS ELSE LS = 10 END IF CALL LANZOS(LV,NSB,LS,SZ, MP, X(-1,0,LP), X(-1,0,RP),JH, & SM( 0,0,LP),SM( 0,0,RP),BH,E,V) IF( LV .GE. 1 ) WRITE(6,100) NL, NR,E(0), IP, FMX IF( NL .EQ. NS/2 - 1 ) THEN E0(0,RT)= E(0) E0(1,RT)= E(1) E0(2,RT)= E(2) E0(3,RT)= E(3) MM(RT)= IP FM(RT)= FMX CALL EXTIME( T1, TIME(RT) ) WRITE(6,111) RT, MM(RT), E0(0,RT), E0(1,RT), E0(2,RT), & FM(RT), TIME(RT) IF((ABS(E0(0,RT)-E0(0,RT-1)).LT.ERE .AND. MM(RT).GT.MM(RT-1)) & .OR. RT.GT.MRT) GOTO 40 RT = RT + 1 MX = MX + DX END IF NR = NR + 1 NP = RP - 1 NL = NL - 1 CALL RENRMR(LV,NR, NSB,SZ(0,1),MP,X(-1,0,RP),X(-1,0,NP), & SM(0,0,RP),SM(0,0,NP), BH(0,1), HM(0,NP), & MX, IP, ER, FMX, V) GOTO 30 C ............ Observation Routines .................................. C Sweep Once More and Create SL/SR for Observation. 40 CONTINUE C WRITE(6,112) ERE C ............ Parameter Do-End .................................. 90 CONTINUE 100 FORMAT(1H ,I3,'-o-o-',I3,' E0:',F21.14,' m:',I3,' Min[DM]:',E10.2) 110 FORMAT(1H ,I2,' m:',I3,' E',3F16.11,' co.',E7.1,' ',F6.2,'s') 111 FORMAT(1H ,I2,' ',I3,' ',3F16.11,' co.',E7.1,' ',F6.2,'s') 112 FORMAT(1H , 'Conv.acc->',3F16.11) C Write(6,*) 'Observation Routines has not been coded yet.' CALL MEMDSP( LV, NL, NSB, MP, X(-1,0,LP) ) Write(6,*) 'Bye-Bye!' C STOP END CCC CC C SUBROUTINE GETSEC( TIME ) REAL*8 TIME C -NEC- CALL CLOCK( TIME ) RETURN END CCC SUBROUTINE EXTIME( T1, T2 ) REAL*8 T1, T2, T3 CALL GETSEC( T3 ) T2 = T3 - T1 T1 = T3 RETURN END CCC SUBROUTINE MEMDSP( LV, N, NSB, MP, X ) INTEGER LV, N, NSB, MP(0:100,0:5), X(-1:101,0:3) IF( LV .GE. 0 ) THEN WRITE(6,*) Write(6,*) 'Typical Memory Requirement for each Variable.' Write(6,*) 'V ------ Vector Length: ', MP(NSB,3) Write(6,*) 'HM,SM --- Matrix Stack: ', X(N+1,2) Write(6,*) 'BH(L/R) - Matrix Stack: ', X(N+1,2)*4 END IF RETURN END CCC SUBROUTINE PRMDSP( LV, NS, TZ, MX, ER, JH) INTEGER LV, NS, TZ, MX REAL*8 ER, JH(0:4) C WRITE(6,*) IF( LV .GE. 0 ) WRITE(6,100) NS, TZ, JH(0), JH(1) IF( LV .GE. 1 ) WRITE(6,*)'initial-m:', MX,' DM-cut-off:',ER IF( LV .GE. 1 ) WRITE(6,*) 'Ising:',JH(2), ' XY:',JH(3) 100 FORMAT(1H ,'NS=',I3,' SZ=',I2,' J1=',F5.2,' J2=',F5.2 ) IF( NS .LT. 8 ) THEN Write(6,*) 'Lattice Size is TOO SMALL. Set NS >= 8.' CALL ERRSTP( 'PRMDSP' ) END IF IF( TZ .GE. 6 ) THEN Write(6,*) 'Total Sz is Larger than (Initial System Size)/2.' Write(6,*) 'You have to modify MAIN-ROUTINE! (Sorry....) ' CALL ERRSTP( 'PRMDSP' ) END IF IF( NS .GT. 100 ) THEN Write(6,*) 'Lattice Size is TOO LARGE. Set NS = 100,', & ' or increase Memory Space.' CALL ERRSTP( 'PRMDSP' ) END IF IF( MOD(NS,2) .EQ. 1 ) THEN Write(6,*) 'Ns should be an even number.' CALL ERRSTP( 'PRMDSP' ) END IF IF( ABS(TZ) .GE. NS/2 ) THEN Write(6,*) 'Total Sz (Var.TZ) should be smaller than NU.' CALL ERRSTP( 'PRMDSP' ) END IF RETURN END CCC SUBROUTINE MATDSP( N, X, ME, CH1, CH2 ) INTEGER N, X(-1:101,0:3) REAL*8 ME(0:*) CHARACTER*3 CH1, CH2 CHARACTER*1 A(0:40,0:40),IMG C IF( X(N+1,1) .GE. 40 ) THEN Write(6,*)'Size (',X(N+1,1),') is too large to display.' RETURN END IF C DO 10 J = 0, X(N+1,1)-1 DO 10 I = 0, X(N+1,1)-1 A(I,J)= '.' 10 CONTINUE C IF( CH1 .EQ. 'BLD' ) THEN DO 20 K = 0, N DO 20 J = 0, X(K,0)-1 DO 20 I = 0, X(K,0)-1 A( X(K,1)+I, X(K,1)+J) = IMG( ME(X(K,2)+J*X(K,0)+I) ) 20 CONTINUE END IF IF( CH1 .EQ. 'BL+' ) THEN DO 30 K = 0, N-1 DO 30 J = 0, X(K+1,0)-1 DO 30 I = 0, X(K ,0)-1 A( X(K,1)+I, X(K+1,1)+J) = IMG( ME(X(K,3)+J*X(K,0)+I) ) 30 CONTINUE END IF IF( CH1 .EQ. 'BL-' ) THEN DO 40 K = 0, N-1 DO 40 J = 0, X(K ,0)-1 DO 40 I = 0, X(K+1,0)-1 A( X(K+1,1)+I,X(K,1)+J) = IMG( ME(X(K,3)+J*X(K+1,0)+I) ) 40 CONTINUE END IF C Write(6,*) '---',CH2,'---' DO 50 I = 0, X(N+1,1)-1 50 Write(6,100) (A(I,J), J=0, X(N+1,1)-1) 100 FORMAT(1H , 39A2 ) RETURN END CCC CHARACTER*1 FUNCTION IMG( F ) REAL*8 F IMG = ' ' IF( F .LT. -1.0E-8 ) IMG = '-' IF( F .LT. -1.0E-4 ) IMG = '=' IF( F .GT. 1.0E-8 ) IMG = '+' IF( F .GT. 1.0E-4 ) IMG = '#' RETURN END CCC SUBROUTINE MDUMP( N, X, H, S ) INTEGER N, X(-1:101,0:3) REAL*8 H(0:*), S(0:10000,0:5) CALL MATDSP( N, X, H, 'BLD', '-H-' ) CALL MATDSP( N, X, S(0,0), 'BLD', 'SzE' ) CALL MATDSP( N, X, S(0,1), 'BL+', 'S+E' ) CALL MATDSP( N, X, S(0,2), 'BL-', 'S-E' ) CALL MATDSP( N, X, S(0,3), 'BLD', 'SzJ' ) CALL MATDSP( N, X, S(0,4), 'BL+', 'S+J' ) CALL MATDSP( N, X, S(0,5), 'BL-', 'S-J' ) RETURN END CCC CC C SUBROUTINE INIBLK( LV, JH, X, S, H ) INTEGER LV, X(-1:101,0:3) REAL*8 H(0:*), JH(0:3), S(0:10000,0:5) C ....[Up,Up] -- [Up,Dw][Dw,Up] -- [Dw,Dw]: Set Sub-Block Size X(-1,0)= 0 X( 0,0)= 1 X( 1,0)= 2 X( 2,0)= 1 X( 3,0)= 0 CALL POINTR( LV, 2, X ) C ..................................Set Initial Hamiltonian H(0)= JH(0)*JH(2)/4.0D0 H(1)=-JH(0)*JH(2)/4.0D0 H(2)= JH(0)*JH(3)/2.0D0 H(3)= JH(0)*JH(3)/2.0D0 H(4)=-JH(0)*JH(2)/4.0D0 H(5)= JH(0)*JH(2)/4.0D0 C .............(Send^Z)..Set Spin Operator at the End/Joint point. C End Point (0,1,2), Joint Point (3,4,5): [E===J]===O===O S(0,0)= 1.0D0/2.0D0 S(1,0)= 1.0D0/2.0D0 S(2,0)= 0.0D0 S(3,0)= 0.0D0 S(4,0)=-1.0D0/2.0D0 S(5,0)=-1.0D0/2.0D0 C .........................(Send^Plus/Minus) S(0,1)= 0.0D0 S(1,1)= 1.0D0 S(2,1)= 1.0D0 S(3,1)= 0.0D0 CALL HCONJU( 2, X, S(0,1), S(0,2) ) C .........................(Sjoint^Z) S(0,3)= 1.0D0/2.0D0 S(1,3)=-1.0D0/2.0D0 S(2,3)= 0.0D0 S(3,3)= 0.0D0 S(4,3)= 1.0D0/2.0D0 S(5,3)=-1.0D0/2.0D0 C .........................(Sjoint^Plus/Minus) S(0,4)= 1.0D0 S(1,4)= 0.0D0 S(2,4)= 0.0D0 S(3,4)= 1.0D0 CALL HCONJU( 2, X, S(0,4), S(0,5) ) C IF( LV .GE. 4 ) CALL MDUMP( 2, X, H, S ) RETURN END CCC SUBROUTINE HCONJU( N, X, MP, MM ) INTEGER N, X(-1:101,0:3) REAL*8 MP(0:*), MM(0:*) DO 10 I = 0, N-1 DO 10 J = 0, X(I+1,0)-1 DO 10 K = 0, X(I ,0)-1 MM(X(I,3)+X(I+1,0)*K+J)= MP(X(I,3)+X(I,0)*J+K) 10 CONTINUE RETURN END CCC CC C SUBROUTINE POINTR( LV, N, X ) INTEGER LV, N, X(-1:101,0:3) C IF( N .GE. 100 ) THEN Write(6,*) 'System Size Over! ', & 'Increase Buffer Size [100].' CALL ERRSTP( 'POINTR' ) END IF C .....................Set Pointer to the Linear Matrix Address X(-1,1)= 0 DO 10 I = 0, N+1 X( I,1)= X(I-1,1) + X(I-1,0) 10 CONTINUE C .....................Set Pointer to Block-Diagonal Matrix-Element X(-1,2)= 0 DO 20 I = 0, N+1 X( I,2)= X(I-1,2) + X(I-1,0)**2 20 CONTINUE C .....................Set Pointer to Block-Offdiagonal Matrix-Element X(-1,3)= 0 DO 30 I = 0, N X( I,3)= X(I-1,3) + X(I-1,0)*X(I,0) 30 CONTINUE C .......................... Pointer Display/Check ... IF( LV .GE. 3 ) THEN Write(6,*) ' Pointers to Matrix Structure: L/R' DO 40 I = -1, N+1 IF( X(I,0) .NE. 0 ) Write(6,100) I,X(I,0),X(I,1),X(I,2),X(I,3) 100 FORMAT(1H ,5I8) 40 CONTINUE END IF IF( LV .GE. 2 ) THEN Write(6,*)'Pointer Maximums:' Write(6,100) N+1, X(N+1,0), X(N+1,1), X(N+1,2), X(N,3) END IF IF( X(N+1,2) .GT. 10000) THEN Write(6,*) 'Memory Over in Operator Matrices.' CALL ERRSTP( 'POINTR' ) END IF RETURN END CCC CC C SUBROUTINE MATSTR(LV, NL,NR, TZ,NSB, XL,XR, SZ, MP) INTEGER EZ, LV, NL,NR, TZ,NSB, XL(-1:101) INTEGER MP(0:100,0:5), XR(-1:101), SZ(0:100,0:1) C .......... Parameter Tansform (Fermion Reprisentation) EZ =(NL + NR + 2)/2 - TZ C............ Ns = NL + NR + 2 IF( EZ.LT.0 .OR. EZ.GT.NL+NR+2 .OR. NL+NR+2.GT.100 ) THEN Write(6,*)'Ez,NL,NR = ',EZ,NL,NR,' Something is wrong.' CALL ERRSTP( 'MATSTR' ) END IF C............................. Index: 0 --- Left, 1 --- Right. IF( EZ .LE. NR+1 ) THEN SZ(0,0)= 0 SZ(0,1)= EZ ELSE SZ(0,0)= EZ-(NR+1) SZ(0,1)= NR+1 END IF C............................. Sub-Space Scan NSB = 1 10 CONTINUE IF( SZ(NSB-1,0) .LT. NL+1 .AND. SZ(NSB-1,1) .GT. 0 ) THEN SZ(NSB,0)= SZ(NSB-1,0) +1 SZ(NSB,1)= SZ(NSB-1,1) -1 NSB = NSB+1 GOTO 10 END IF C ........................ Matrix Structure Pointers ... DO 20 I = 0, NSB-1 MP(I,0)= XL(SZ(I,0)-1) + XL(SZ(I,0)) MP(I,1)= XR(SZ(I,1)-1) + XR(SZ(I,1)) MP(I,2)= MP(I,0) * MP(I,1) 20 CONTINUE MP(0,3)= 0 MP(0,4)= 0 MP(0,5)= 0 DO 30 I = 1, NSB MP(I,3)= MP(I-1,3) + MP(I-1,2) MP(I,4)= MP(I-1,4) + MP(I-1,0)**2 MP(I,5)= MP(I-1,5) + MP(I-1,1)**2 30 CONTINUE C IF( LV .GE. 3 ) THEN Write(6,*) ' Szl, Szr, Dim[L], Dim[R], Dim[LR], Pointers:' DO 40 I = 0, NSB IF( MP(I,2) .NE. 0 ) & WRITE(6,100) SZ(I,0), SZ(I,1), MP(I,0), MP(I,1), & MP(I,2), MP(I,3), MP(I,4), MP(I,5) 40 CONTINUE 100 FORMAT(1H ,8I8) END IF C IF( MP(NSB,3) .GT. 100000 ) THEN Write(6,*) 'Memory over in Lanczos Vectors' CALL ERRSTP( 'MATSTR' ) END IF IF( MAX(MP(NSB,4),MP(NSB,5)) .GT. 40000 ) THEN Write(6,*) 'Memory over in Block Hamiltonian.' CALL ERRSTP( 'MATSTR' ) END IF C RETURN END CCC CC C SUBROUTINE RENRML(LV, N,NSB,SZ,MP,XO,XN,SO,SN, BH, H,MX,IP, & ER, FMAX, V) INTEGER BA,EP,BP,LV, N,NSB,MX, MP(0:100,0:5), XO(-1:101,0:3) INTEGER IP, J, FLG,SZ(0:100), Y(-1:101,0:3), XN(-1:101,0:3) REAL*8 ER, BH(0:40000), SO(0:10000,0:5), V(0:100000,0:3) REAL*8 FMAX, H(0:10000), SN(0:10000,0:5), TP(0:40000) C CALL EXINDX( LV, N, NSB, Y, SZ, MP(0,0) ) CALL VEC0( Y(N+1,2), V(0,2) ) CALL VEC0( Y(N+1,1), V(0,3) ) DO 10 I = 0, NSB-1 CALL GETDML( MP(I,0),MP(I,1), V(Y(SZ(I),2),2), V( MP(I,3),1)) CALL DIAGON( LV, MP(I,0), V(Y(SZ(I),2),2), V(Y(SZ(I),1),3)) 10 CONTINUE CALL SELECT( LV,FLG, MX, Y(N+1,1), ER, V, IP, FMAX) CALL NEWIDX( LV, N, Y, XN(-1,0), V(0,3) ) C DO 20 I = 0, N J = I-1 EP = XN(I,2) BP = Y(I,2) + Y(I,0)*(Y(I,0)-XN(I,0)) BA = Y(I,2) CALL ROTATD(FLG,Y(I,0),XN(I,0),BH(BA), H(EP), V(BP,2),V(0,3)) CALL CRESZE( Y(I,0),XO(J,0),TP, SO(XO(J,2),3)) CALL ROTATD(FLG,Y(I,0),XN(I,0),TP, SN(EP,0), V(BP,2),V(0,3)) CALL CRESZJ( Y(I,0),XO(J,0),TP ) CALL ROTATD(FLG,Y(I,0),XN(I,0),TP, SN(EP,3), V(BP,2),V(0,3)) IF( I .EQ. N ) GOTO 20 EP = XN(I,3) BA = Y(I+1,2) + Y(I+1,0)*(Y(I+1,0)-XN(I+1,0)) CALL CRESPE( Y(I,0),XO(J,0),TP, SO(XO(J,3),4)) CALL ROTATP(FLG,Y(I,0),XN(I,0),TP,SN(EP,1),V(BP,2),V(BA,2),V(0,3)) CALL CRESPJ( Y(I,0),XO(J,0),TP ) CALL ROTATP(FLG,Y(I,0),XN(I,0),TP,SN(EP,4),V(BP,2),V(BA,2),V(0,3)) 20 CONTINUE C CALL HCONJU( N, XN, SN(0,1), SN(0,2) ) CALL HCONJU( N, XN, SN(0,4), SN(0,5) ) IF( LV .GE. 4 ) CALL MDUMP( N, XN(-1,0), H, SN ) RETURN END CCC SUBROUTINE RENRMR(LV, N,NSB,SZ,MP,XO,XN,SO,SN, BH, H, MX,IP, & ER, FMAX, V) INTEGER BA,EP,BP,LV, N,NSB,MX, MP(0:100,0:5), XO(-1:101,0:3) INTEGER IP, J, FLG,SZ(0:100), Y(-1:101,0:3), XN(-1:101,0:3) REAL*8 ER, BH(0:40000), SO(0:10000,0:5), V(0:100000,0:3) REAL*8 FMAX, H(0:10000), SN(0:10000,0:5), TP(0:40000) C CALL EXINDX( LV, N, NSB, Y, SZ, MP(0,1) ) CALL VEC0( Y(N+1,2), V(0,2) ) CALL VEC0( Y(N+1,1), V(0,3) ) DO 10 I = 0, NSB-1 CALL GETDMR( MP(I,0),MP(I,1), V(Y(SZ(I),2),2), V( MP(I,3),1)) CALL DIAGON( LV, MP(I,1), V(Y(SZ(I),2),2), V(Y(SZ(I),1),3)) 10 CONTINUE CALL SELECT( LV,FLG, MX, Y(N+1,1), ER, V, IP, FMAX ) CALL NEWIDX( LV, N, Y, XN(-1,0), V(0,3) ) C DO 20 I = 0, N J = I-1 EP = XN(I,2) BP = Y(I,2) + Y(I,0)*(Y(I,0)-XN(I,0)) BA = Y(N+1,2) - Y(I+1,2) CALL ROTATD(FLG,Y(I,0),XN(I,0),BH(BA), H(EP), V(BP,2),V(0,3)) CALL CRESZE( Y(I,0),XO(J,0),TP, SO(XO(J,2),3)) CALL ROTATD(FLG,Y(I,0),XN(I,0),TP, SN(EP,0), V(BP,2),V(0,3)) CALL CRESZJ( Y(I,0),XO(J,0),TP ) CALL ROTATD(FLG,Y(I,0),XN(I,0),TP, SN(EP,3), V(BP,2),V(0,3)) IF( I .EQ. N ) GOTO 20 EP = XN(I,3) BA = Y(I+1,2) + Y(I+1,0)*(Y(I+1,0)-XN(I+1,0)) CALL CRESPE( Y(I,0),XO(J,0),TP, SO(XO(J,3),4)) CALL ROTATP(FLG,Y(I,0),XN(I,0),TP,SN(EP,1),V(BP,2),V(BA,2),V(0,3)) CALL CRESPJ( Y(I,0),XO(J,0),TP ) CALL ROTATP(FLG,Y(I,0),XN(I,0),TP,SN(EP,4),V(BP,2),V(BA,2),V(0,3)) 20 CONTINUE C CALL HCONJU( N, XN, SN(0,1), SN(0,2) ) CALL HCONJU( N, XN, SN(0,4), SN(0,5) ) IF( LV .GE. 4 ) CALL MDUMP( N, XN(-1,0), H, SN ) RETURN END CCC SUBROUTINE EXINDX( LV, N, NSB, Y, SZ, MP ) INTEGER MP(0:100), LV, N, NSB, Y(-1:101,0:3), SZ(0:100) DO 10 I =-1, N+1 Y(I,0)= 0 10 CONTINUE DO 20 I = 0, NSB-1 Y(SZ(I),0)= MP(I) 20 CONTINUE CALL POINTR( LV, N, Y ) IF( Y(N+1,2) .GT. 100000 ) THEN Write(6,*) 'Vector Overflow!' CALL ERRSTP( 'RENRML' ) END IF RETURN END CCC SUBROUTINE GETDML( WL, WR, V2, V1 ) INTEGER WL, WR REAL*8 V2(0:*), V1(0:*) IF( WL .GE. WR ) THEN DO 10 J = 0, WR-1 DO 10 K = 0, WL-1 DO 10 L = 0, WL-1 V2(WL*K+L)= V2(WL*K+L) + V1(WR*K+J) * V1(WR*L+J) 10 CONTINUE ELSE DO 20 K = 0, WL-1 DO 20 L = 0, WL-1 DO 20 J = 0, WR-1 V2(WL*K+L)= V2(WL*K+L) + V1(WR*K+J) * V1(WR*L+J) 20 CONTINUE END IF RETURN END CCC SUBROUTINE GETDMR( WL, WR, V2, V1 ) INTEGER WL, WR REAL*8 V2(0:*), V1(0:*) IF( WL .GE. WR ) THEN DO 10 J = 0, WR-1 DO 10 K = 0, WR-1 DO 10 L = 0, WL-1 V2(WR*J+K)= V2(WR*J+K) + V1(WR*L+K) * V1(WR*L+J) 10 CONTINUE ELSE DO 20 L = 0, WL-1 DO 20 J = 0, WR-1 DO 20 K = 0, WR-1 V2(WR*J+K)= V2(WR*J+K) + V1(WR*L+K) * V1(WR*L+J) 20 CONTINUE END IF RETURN END CCC SUBROUTINE DIAGON( LV, N, A, E ) INTEGER LV, N, IERR REAL*8 W(0:10000), A(0:*), E(0:*) C IF( N .EQ. 0 ) THEN RETURN END IF IF( N .EQ. 1 ) THEN E(0)= A(0) A(0)= 1.0D0 IF( LV .GE. 3 ) THEN Write(6,*) '-------------------- Eigen Value', E(0) END IF RETURN END IF C C +----------------------------------------------+ C | Diagonalization of the matrix A: | C | ASL routine: Please read NEC mannual | C +----------------------------------------------+ C -NEC- CALL DCSMAA( A, N, N, E, W, IERR ) IF( IERR .NE. 0 ) THEN Write(6,*) 'ERROR in DCSMAA, Flag= ', IERR CALL ERRSTP( 'DIAGON' ) END IF IF( LV .GE. 3 ) THEN DO 10 I = 0, N-1 10 Write(6,*) '-------------------- Eigen Value', E(I) END IF RETURN END CCC C *========================== CAUTION =========================* C | Density Matrix Cut-off Routine (Subrouitne SELECT) | C | is NOT perfect now: Do check your result when the | C | ground state is collection of local states!! | C *============================================================* SUBROUTINE SELECT( LV, FLG, MX, N, ER, V, IP, FMAX ) INTEGER LV, FLG, MX, N, IP REAL*8 TP,COR, ER, V(0:100000,0:3), FMAX C IF( N .LE. MX ) THEN FLG = 0 IF( LV .GE. 2 ) Write(6,*) 'All States Preserved' IP = N FMAX = 0.0D0 RETURN END IF C FLG = 1 CALL VMOVE( N, 0, 3, V ) DO 30 I = 0, N-2 IP = I FMAX = V(I,0) DO 20 J = I+1, N-1 IF(FMAX.LT.V(J,0) ) THEN FMAX = V(J,0) IP = J END IF 20 CONTINUE IF( IP .NE. I ) THEN TP = V( I,0) V( I,0)= V(IP,0) V(IP,0)= TP END IF 30 CONTINUE C ------- Cut Off Ratio, COR -------- COR = 1.05D0 C ------- (Dummy for Debug) -------- IP = MIN( MX-1,N-1 ) V(IP+1,0)= V(IP,0) 40 CONTINUE IF( IP .EQ. 0 ) THEN IF( LV .GE. 3 ) Write(6,*) 'Only One State Is Selected.' IF( LV .GE. 2 ) THEN WRITE(6,*) 'Your result might be incorrect,' WRITE(6,*) 'because of minor error in the ' WRITE(6,*) 'subroutine SELECT. Do check your result!' END IF GOTO 45 END IF IF( V(IP,0) .LT. ER ) THEN IP = IP-1 GOTO 40 END IF IF( V(IP,0)/V(IP+1,0) .LT. COR ) THEN IP = IP-1 GOTO 40 END IF 45 CONTINUE C FMAX = MAX( V(IP,0), ER ) DO 50 I = 0, N-1 IF( V(I,3) .LT. FMAX ) V(I,3)= -2.0D0 50 CONTINUE IP = IP + 1 IF( LV .GE. 2 ) WRITE(6,100) IP, FMAX 100 FORMAT(1H ,' ', & 'Kept;', I5,' Min[DM]', E15.7 ) C RETURN END CCC SUBROUTINE NEWIDX( LV, N, Y, X, V ) INTEGER LV, N, Y(-1:101,0:3), X(-1:101,0:3) REAL*8 CNT, V(0:100000) DO 20 I = -1, N+1 CNT = Y(I,0) DO 10 J = 0, Y(I,0)-1 IF( V(Y(I,1)+J) .LT. -1.0D0 ) CNT = CNT -1 10 CONTINUE X(I,0)= CNT 20 CONTINUE CALL POINTR( LV, N, X ) IF( LV .GE. 2 ) THEN WRITE(6,100) (X(I,0),I=0,N) 100 FORMAT(1H ,15I5) END IF RETURN END CCC SUBROUTINE ROTATD( FLG, WO, WN, GM, SM, OM, TM ) INTEGER FLG, WO, WN REAL*8 SM(0:*), GM(0:*), OM(0:*), TM(0:*) IF( FLG .EQ. 0 ) THEN DO 10 J = 0, WO-1 DO 10 I = 0, WO-1 SM(WO*J+I)= GM(WO*J+I) 10 CONTINUE ELSE C ------ TM(I,K) ++ GM(I,J)*OM(J,K) CALL VEC0( WO*WN, TM ) DO 20 K = 0, WN-1 DO 20 J = 0, WO-1 DO 20 I = 0, WO-1 TM(WO*K+I)= TM(WO*K+I) + GM(WO*J+I)*OM(WO*K+J) 20 CONTINUE C ------ SM(I,K) ++ OM(J,I)*TM(J,K) CALL VEC0( WN*WN, SM ) DO 30 K = 0, WN-1 DO 30 I = 0, WN-1 DO 30 J = 0, WO-1 SM(WN*K+I)= SM(WN*K+I) + OM(WO*I+J)*TM(WO*K+J) 30 CONTINUE END IF RETURN END CCC SUBROUTINE ROTATP( FLG, Y, X, GM, SM, OM, QM, TM ) INTEGER FLG, Y(0:1), X(0:1) REAL*8 GM(0:*),SM(0:*),OM(0:*),QM(0:*),TM(0:*) IF( FLG .EQ. 0 ) THEN DO 10 J = 0, Y(1)-1 DO 10 I = 0, Y(0)-1 SM(Y(0)*J+I)= GM(Y(0)*J+I) 10 CONTINUE ELSE C ------ TM(I,K) ++ GM(I,J)*QM(J,K) CALL VEC0( Y(0)*X(1), TM ) DO 20 K = 0, X(1)-1 DO 20 J = 0, Y(1)-1 DO 20 I = 0, Y(0)-1 TM(Y(0)*K+I)= TM(Y(0)*K+I) + GM(Y(0)*J+I)*QM(Y(1)*K+J) 20 CONTINUE C ------ SM(I,K) ++ OM(J,I)*TM(J,K) CALL VEC0( X(0)*X(1), SM ) DO 30 K = 0, X(1)-1 DO 30 I = 0, X(0)-1 DO 30 J = 0, Y(0)-1 SM(X(0)*K+I)= SM(X(0)*K+I) + OM(Y(0)*I+J)*TM(Y(0)*K+J) 30 CONTINUE END IF RETURN END CCC SUBROUTINE CRESZE( W, X, BH, S ) INTEGER MP, EP, W, X(0:1) REAL*8 BH(0:*), S(0:*) CALL VEC0( W*W, BH ) DO 10 K = 0, 1 EP = K*X(0)*X(0) MP = K*(W *X(0)+X(0)) DO 10 J = 0, X(K)-1 DO 10 I = 0, X(K)-1 BH(MP+J*W+I)= S(EP+J*X(K)+I) 10 CONTINUE RETURN END CCC SUBROUTINE CRESZJ( W, X, BH ) INTEGER W, X(0:1) REAL*8 BH(0:*) CALL VEC0( W*W, BH ) DO 10 I = 0, X(0)-1 BH((W+1)*I)=-1.0D0/2 10 CONTINUE DO 20 I = X(0), X(0)+X(1)-1 BH((W+1)*I)= 1.0D0/2 20 CONTINUE RETURN END CCC SUBROUTINE CRESPE( Y, X, BH, S ) INTEGER MP, EP, Y(0:1), X(0:2) REAL*8 BH(0:*), S(0:*) CALL VEC0( Y(0)*Y(1), BH ) DO 10 K = 0, 1 MP = K*(Y(0)*X(1)+X(0)) EP = K* X(1)*X(0) DO 10 J = 0, X(K+1)-1 DO 10 I = 0, X(K )-1 BH(MP+J*Y(0)+I)= S(EP+J*X(K)+I) 10 CONTINUE RETURN END CCC SUBROUTINE CRESPJ( Y, X, BH ) INTEGER Y(0:1), X(0:2) REAL*8 BH(0:*) CALL VEC0( Y(0)*Y(1), BH ) DO 10 I = 0, X(1)-1 BH(X(0)+(Y(0)+1)*I)= 1.0D0 10 CONTINUE RETURN END CCC CC C SUBROUTINE BLOCKH( LV,NSB, SZ, MP, XL, XR, JH,HL,HR, SL,SR, BH) INTEGER OL,DL,ZL, LV,NSB, SZ(0:100, 0:1), XL(-1:101,0:3) INTEGER OR,DR,ZR, MP(0:100, 0:5), XR(-1:101,0:3) REAL*8 JH(0:4),SL(0:10000,0:5), HL(0:10000) REAL*8 BH(0:40000,0:1),SR(0:10000,0:5), HR(0:10000) C IF( MAX( MP(NSB,4),MP(NSB,5)) .LT. 40000 ) THEN CALL VEC0( MP(NSB,4), BH(0,0) ) CALL VEC0( MP(NSB,5), BH(0,1) ) ELSE Write(6,*) 'Memory Over in Block Hamiltonian Stack' CALL ERRSTP( 'BLOCKH' ) END IF C DO 20 I = 0, NSB-1 ZL = SZ( I,0)-1 ZR = SZ( I,1)-1 DL = XL(ZL,2) OL = XL(ZL,3) DR = XR(ZR,2) OR = XR(ZR,3) IF( MP(I,0) .EQ. 0 ) GOTO 10 CALL BLKHH( XL(ZL,0), HL(DL ), BH(MP(I,4),0) ) CALL BLKSD( XL(ZL,0), SL(DL,0), SL(DL,3), JH, BH(MP(I,4),0) ) CALL BLKSP( XL(ZL,0), SL(OL,1), SL(OL,4), JH, BH(MP(I,4),0) ) CALL BLKSM( XL(ZL,0), SL(OL,2), SL(OL,5), JH, BH(MP(I,4),0) ) C IF( LV .GE. 4 ) CALL SMTDSP(I, MP(I,0), BH(MP(I,4),0) ) 10 CONTINUE IF( MP(I,1) .EQ. 0 ) GOTO 20 CALL BLKHH( XR(ZR,0), HR(DR ), BH(MP(I,5),1) ) CALL BLKSD( XR(ZR,0), SR(DR,0), SR(DR,3), JH, BH(MP(I,5),1) ) CALL BLKSP( XR(ZR,0), SR(OR,1), SR(OR,4), JH, BH(MP(I,5),1) ) CALL BLKSM( XR(ZR,0), SR(OR,2), SR(OR,5), JH, BH(MP(I,5),1) ) C IF( LV .GE. 4 ) CALL SMTDSP(I, MP(I,1), BH(MP(I,5),0) ) 20 CONTINUE C IF( LV .GE. 3 ) Write(6,*) 'Block Hamiltonian is created.' RETURN END CCC SUBROUTINE SMTDSP( L, N, A ) INTEGER L, N REAL*8 A(0:*) WRITE(6,*) '--- Width:', L IF( N .GT. 12 ) RETURN DO 10 I = 0, N-1 WRITE(6,100) ( A(I+N*J), J=0, N-1 ) 100 FORMAT(1H ,12F6.2) 10 CONTINUE RETURN END CCC SUBROUTINE BLKHH( X, H, BH ) INTEGER AP, X(0:1), W, EP REAL*8 H(0:*), BH(0:*) W = X(0)+X(1) DO 10 I = 0, 1 AP = I*(W*X(0)+X(0)) EP = I* X(0)*X(0) DO 10 J = 0, X(I)-1 DO 10 K = 0, X(I)-1 BH(AP+W*J+K)= BH(AP+W*J+K) + H(EP+X(I)*J+K) 10 CONTINUE RETURN END CCC SUBROUTINE BLKSD( X, S0, S1, JH, BH ) INTEGER X(0:1), EP, AP, W REAL*8 F0,F1,S0(0:*), S1(0:*), JH(0:3), BH(0:*) C ---------------- End Point:S0, Joint Point:S1 W = X(0)+X(1) DO 10 I = 0, 1 C ---------------- N.N.Int.:F0, N.N.N.Int.:F1 F0 = JH(0)*JH(2)*(I-0.50D0) F1 = JH(1)*JH(2)*(I-0.50D0) AP = I*(W*X(0)+X(0)) EP = I* X(0)*X(0) DO 10 J = 0, X(I)-1 DO 10 K = 0, X(I)-1 BH(AP+W*J+K)= BH(AP+W*J+K) +F0*S1(EP+X(I)*J+K) +F1*S0(EP+X(I)*J+K) 10 CONTINUE RETURN END CCC SUBROUTINE BLKSP(X, S0, S1, JH, BH ) INTEGER X(0:1), W, AP REAL*8 S0(0:*), S1(0:*), JH(0:4), BH(0:*), F0, F1 W = X(0)+X(1) F0 = JH(0)*JH(3)/2.0D0 F1 = JH(1)*JH(3)/2.0D0 AP = X(0)*W DO 10 J = 0, X(1)-1 DO 10 K = 0, X(0)-1 BH(AP+W*J+K)= BH(AP+W*J+K) +F0*S1(X(0)*J+K) +F1*S0(X(0)*J+K) 10 CONTINUE RETURN END CCC SUBROUTINE BLKSM(X, S0, S1, JH, BH ) INTEGER X(0:1), W, AP REAL*8 S0(0:*), S1(0:*), JH(0:4), BH(0:*), F0, F1 W = X(0)+X(1) F0 = JH(0)*JH(3)/2.0D0 F1 = JH(1)*JH(3)/2.0D0 AP = X(0) DO 10 J = 0, X(0)-1 DO 10 K = 0, X(1)-1 BH(AP+W*J+K)= BH(AP+W*J+K) +F0*S1(X(1)*J+K) +F1*S0(X(1)*J+K) 10 CONTINUE RETURN END CCC CC C SUBROUTINE MULT(LV, ID, IS, NSB, SZ, MP, XL,XR, JH,SL,SR,BH,V) INTEGER BA, LV, ID, IS, NSB, SZ(0:100,0:1), XL(-1:101,0:3) INTEGER YR, ZL, ZR, WL, WR, MP(0:100,0:5), XR(-1:101,0:3) REAL*8 JH(0:3), BH(0:40000,0:1),V(0:100000,0:3) REAL*8 SL(0:10000,0:5), SR(0:10000,0:5) C CALL VEC0( MP(NSB,3), V(0,ID) ) DO 20 I = 0, NSB-1 WL = MP(I,0) WR = MP(I,1) IF( MP(I,2) .EQ. 0 ) GOTO 10 BA = MP(I,3) ZL = SZ(I,0)-1 ZR = SZ(I,1)-1 CALL MULBLK(ID,IS,BA,WL,WR, BH(MP(I,4),0), BH(MP(I,5),1), V ) CALL MULMMD(ID,IS,BA,WR,XL(ZL,0),XR(ZR,0), JH,V ) CALL MULMRD(ID,IS,BA,WR,XL(ZL,0),XR(ZR,0), SR(XR(ZR,2),3), JH,V ) CALL MULLMD(ID,IS,BA,WR,XL(ZL,0),XR(ZR,0), SL(XL(ZL,2),3), JH,V ) IF( I .EQ. NSB-1 ) GOTO 10 YR = SZ(I,1)-2 CALL MULMMO(ID,IS,BA,WR,XL(ZL,0),XR(YR,0), JH,V ) CALL MULMRO(ID,IS,BA,WR,XL(ZL,0),XR(YR,0), SR(XR(YR,3),5), JH,V ) CALL MULLMO(ID,IS,BA,WR,XL(ZL,0),XR(YR,0), SL(XL(ZL,3),4), JH,V ) 10 CONTINUE 20 CONTINUE C IF( LV .GE. 5 ) Write(6,*) 'End of Matrix Multiplication.' RETURN END CCC SUBROUTINE MULBLK(ID,IS, BA, WL, WR, BHL, BHR, V ) INTEGER ID,IS, BA, WL, WR REAL*8 BHL(0:*), BHR(0:*), V(0:100000,0:3) C IF( WR .GT. WL ) THEN DO 10 I = 0, WL-1 DO 10 J = 0, WL-1 DO 10 K = 0, WR-1 V(BA+J*WR+K,ID)= V(BA+J*WR+K,ID) +BHL(I*WL+J)*V(BA+I*WR+K,IS) 10 CONTINUE DO 20 I = 0, WL-1 DO 20 J = 0, WR-1 DO 20 K = 0, WR-1 V(BA+I*WR+K,ID)= V(BA+I*WR+K,ID) +BHR(J*WR+K)*V(BA+I*WR+J,IS) 20 CONTINUE ELSE DO 30 K = 0, WR-1 DO 30 I = 0, WL-1 DO 30 J = 0, WL-1 V(BA+J*WR+K,ID)= V(BA+J*WR+K,ID) +BHL(I*WL+J)*V(BA+I*WR+K,IS) 30 CONTINUE DO 40 J = 0, WR-1 DO 40 K = 0, WR-1 DO 40 I = 0, WL-1 V(BA+I*WR+K,ID)= V(BA+I*WR+K,ID) +BHR(J*WR+K)*V(BA+I*WR+J,IS) 40 CONTINUE END IF RETURN END CCC SUBROUTINE MULMMD(ID,IS, BA, WR, XL, XR, JH, V ) INTEGER AP, ID,IS, BA, WR, XL(0:1), XR(0:1) REAL*8 FC, JH(0:3), V(0:100000,0:3) DO 10 I = 0, 1 DO 10 J = 0, 1 FC = JH(0)* JH(2)*(I-0.50D0)*(J-0.50D0) AP = I* WR* XL(0) + J*XR(0) + BA DO 10 K = 0, WR*(XL(I)-1), WR DO 10 L = 0, XR(J)-1 V(AP+K+L,ID)= V(AP+K+L,ID) + FC*V(AP+K+L,IS) 10 CONTINUE RETURN END CCC SUBROUTINE MULMRD(ID,IS, BA, WR, XL, XR, SR, JH, V ) INTEGER AP, EP,ID,IS, BA, WR, XL(0:1), XR(0:1) REAL*8 FC, JH(0:3), SR(0:*), V(0:100000,0:3) DO 10 I = 0, 1 DO 10 J = 0, 1 AP = I *WR*XL(0) + J*XR(0) + BA EP = J *XR(0)**2 DO 10 K = 0, XR(J)-1 DO 10 L = 0, XR(J)-1 FC = JH(1)*JH(2)*(I-0.50D0) *SR(EP+L+K*XR(J)) DO 10 M = 0, WR*(XL(I)-1), WR V(AP+M+L,ID)= V(AP+M+L,ID) + FC*V(AP+M+K,IS) 10 CONTINUE RETURN END CCC SUBROUTINE MULLMD(ID,IS, BA, WR, XL, XR, SL, JH, V ) INTEGER AP, EP,ID,IS, BA, WR, XL(0:1), XR(0:1) REAL*8 FC, JH(0:3), SL(0:*), V(0:100000,0:3) DO 10 I = 0, 1 DO 10 J = 0, 1 AP = I *WR*XL(0) + J*XR(0) + BA EP = I *XL(0)**2 DO 10 K = 0, XL(I)-1 DO 10 L = 0, XL(I)-1 FC = JH(1)*JH(2)*(J-0.50D0) *SL(EP+L+K*XL(I)) DO 10 M = 0, XR(J)-1 V(AP+WR*L+M,ID)= V(AP+WR*L+M,ID) + FC*V(AP+WR*K+M,IS) 10 CONTINUE RETURN END CCC SUBROUTINE MULMMO( ID, IS, BA, WV, XL, XR, JH, V ) INTEGER VP,HP,WH,ID, IS, BA, WV, XL(0:1), XR(-1:1) REAL*8 FC, JH(0:3), V(0:100000,0:3) VP = BA + WV* XL(0) HP = BA + WV*(XL(0)+XL(1)) + XR(-1) WH = XR(-1)+XR(0) FC = JH(0)*JH(3)/2.0D0 DO 10 I = 0, XL(1)-1 *VDIR NODEP DO 10 J = 0, XR(0)-1 V(VP+WV*I+J,ID)= V(VP+WV*I+J,ID) + FC*V(HP+WH*I+J,IS) V(HP+WH*I+J,ID)= V(HP+WH*I+J,ID) + FC*V(VP+WV*I+J,IS) 10 CONTINUE RETURN END CCC SUBROUTINE MULMRO( ID, IS, BA, WV, XL, XR, SR, JH, V ) INTEGER VP,HP,WH,ID, IS, BA, WV, XL(0:1), XR(-1:1) REAL*8 FC, JH(0:3), SR(0:*), V(0:100000,0:3) WH = XR(-1)+XR(0) DO 10 I = 0, 1 VP = BA + WV* XL(0) + I*XR( 0) HP = BA + WV*(XL(0)+XL(1)) + I*XR(-1) EP = I* XR(0)*XR(-1) DO 10 J = 0, XR(I-1)-1 DO 10 K = 0, XR(I)-1 FC = JH(1)*JH(3)*SR(EP+XR(I)*J+K)/2.0D0 *VDIR NODEP DO 10 L = 0, XL(1)-1 V(VP+WV*L+K,ID)= V(VP+WV*L+K,ID) + FC*V(HP+WH*L+J,IS) V(HP+WH*L+J,ID)= V(HP+WH*L+J,ID) + FC*V(VP+WV*L+K,IS) 10 CONTINUE RETURN END CCC SUBROUTINE MULLMO( ID, IS, BA, WV, XL, XR, SL, JH, V ) INTEGER VP,HP,WH,ID, IS, BA, WV, XL(0:2), XR(-1:1) REAL*8 FC, JH(0:3), SL(0:*), V(0:100000,0:3) WH = XR(-1)+XR(0) DO 10 I = 0, 1 VP = BA + I*WV*XL(0) HP = BA + WV*(XL(0)+XL(1)) + I*WH*XL(1) + XR(-1) EP = I* XL(0)*XL(1) DO 10 J = 0, XL(I+1)-1 DO 10 K = 0, XL(I)-1 FC = JH(1)*JH(3)*SL(EP+XL(I)*J+K)/2.0D0 *VDIR NODEP DO 10 L = 0, XR(0)-1 V(VP+WV*K+L,ID)= V(VP+WV*K+L,ID) + FC*V(HP+WH*J+L,IS) V(HP+WH*J+L,ID)= V(HP+WH*J+L,ID) + FC*V(VP+WV*K+L,IS) 10 CONTINUE RETURN END CCC CC------------------------------------------- C Subroutines For the Lanczos Diagonalization C ------------------------------------------- SUBROUTINE ERRSTP( A ) CHARACTER*6 A Write(6,*) 'Software Termination in SUBROUTINE ', A, '.' Write(6,*) 'Do check Program/Parameters.' STOP END CCC SUBROUTINE INPRO( NV, ID, IS, PRD, V) INTEGER NV, ID, IS REAL*8 PRD, V(0:100000,0:3) PRD = 0.0D0 DO 10 I = 0, NV-1 PRD = PRD + V(I,ID)*V(I,IS) 10 CONTINUE RETURN END CCC SUBROUTINE ADDER( NV, ID, IS1, IS2, A, B, V) INTEGER NV, ID, IS1, IS2 REAL*8 A, B, V(0:100000,0:3) DO 10 I = 0, NV-1 V(I,ID)= A*V(I,IS1) + B*V(I,IS2) 10 CONTINUE RETURN END CCC SUBROUTINE DIVER( NV, ID, IS, D, V) INTEGER NV, ID, IS REAL*8 D, V(0:100000,0:3) DO 10 I = 0, NV-1 V(I,ID)= V(I,IS) / D 10 CONTINUE RETURN END CCC SUBROUTINE VEC0( NV, V) INTEGER NV REAL*8 V(0:100000) DO 10 J = 0, NV-1 V(J)= 0.0D0 10 CONTINUE RETURN END CCC SUBROUTINE VMOVE( NV, ID, IS, V ) INTEGER NV, ID, IS REAL*8 V(0:100000,0:3) DO 10 I = 0, NV-1 V(I,ID)= V(I,IS) 10 CONTINUE RETURN END CCC SUBROUTINE VECINI( NV, V) INTEGER NV REAL*8 PRD, V(0:100000,0:3) CALL VEC0( NV, V(0,0) ) CALL VEC0( NV, V(0,1) ) CALL VEC0( NV, V(0,2) ) CALL VEC0( NV, V(0,3) ) DO 10 I = 0,NV-1 V(I,0)=( DBLE(IRAND()) + 00.1 )/SQRT( DBLE(NV) ) 10 CONTINUE CALL INPRO( NV, 0, 0, PRD, V ) PRD = SQRT(PRD) CALL DIVER( NV, 0, 0, PRD, V ) RETURN END CCC CC C SUBROUTINE LANZOS(LV, NSB,LLS,SZ, MP, XL,XR, JH, SL,SR, BH, E, V) INTEGER LS,MLS,LK,LV, NSB, NV, MP(0:100,0:5) INTEGER NQ,LLS,SZ(*), XL(*), XR(*) REAL*8 JH(*), BH(*), SL(*), SR(*), E(0:3) REAL*8 ENG, E2 , ALP(0:999), TE(0:999),V(0:100000,0:3) REAL*8 EPS, E2P, BET(0:999), ME(0:999,0:4), QE(0:2999) C IF( LV.GE.3 ) WRITE(6,*) '-vvvv- Sub. LANcZOS -vvvvv-' C +----------------------------------------+ C | LLS - Minimum Lanczos Step MLS | C | MLS - Maximum Lanczos Step 999 | C | EPS - Eigenvalue Convergence Parameter | C +----------------------------------------+ NV = MP(NSB,3) MLS = 190 EPS = 1.0E-10 C ++++++ Initialization Steps ++++++ E2 = 0.0D0 CALL VECINI( NV, V ) CALL VMOVE( NV,3,0,V ) CALL MULT(LV,1,0,NSB, SZ, MP, XL,XR, JH,SL,SR,BH,V) IF( LV .GE. 3 ) WRITE(6,*) IF( LV .GE. 3 ) WRITE(6,*) '-+- Eigenvalues -+-' C ++++++ LS = Lanczos Step, LK = 5*(LS/5) ++++++ LK = 0 LS = 0 10 CONTINUE E2P = E2 C ++++++ Five Laanczos Steps ++++++ DO 20 J = 0, 4 CALL MACRO1(LV,NSB,SZ,MP,XL,XR,JH,SL,SR,BH,V,TE,ALP, BET,LK,0 ) 20 CONTINUE C C ++++++ Diagonalization of the Tri-Diagonal matrix ++++++ CALL MALANZ( LK, LV, ALP,BET,ENG, E2, ME) LS = LS + 1 C ++++++ Convergence Check ++++++ IF( LS*5 .LE. LLS ) GOTO 10 IF( LS*5 .LE. MLS .AND. ABS(E2P-E2) .GE. EPS ) GOTO 10 LS = LS - 1 IF( LV .GE. 2 ) WRITE(6,100) ENG 100 FORMAT(1H ,' Energy Upper Bound by Lanczos = ',E30.20) C ++++++ Lanczos Vector Creation ++++++ DO 30 I = 0,999 TE(I) = ME(I,0) 30 CONTINUE C ++++++ Two-pass Method ++++++ CALL VMOVE( NV, 0, 3, V ) CALL ADDER( NV, 3, 3, 0, TE(1), 0.0D0, V ) CALL MULT(LV,1,0,NSB, SZ, MP, XL,XR, JH,SL,SR,BH,V) LK = 0 C +--------------------------------+ C | NQ - Stacked-Eigenvalue Number | C | QE - Stacked-Eigenvalue | C +--------------------------------+ NQ = 0 DO 50 I = 0, LS DO 40 J = 0, 4 CALL MACRO1(LV,NSB,SZ,MP,XL,XR,JH,SL,SR,BH,V,TE,ALP, BET,LK,1 ) 40 CONTINUE C C ++++++ Stack Low-Lying Eigenvalues ++++++ CALL MALAN2( LK, NQ, ALP, BET, QE(NQ) ) IF( LV .GE. 5 ) WRITE(6,*) LK, ' -th Second Pass' IF( NQ .GE. 2999 ) THEN WRITE(6,*) 'Eigenvalue Stack Over Flow.' CALL ERRSTP( 'LANCZOS' ) END IF 50 CONTINUE C C ++++++ Sort/Search of Low-Lying states ++++++ CALL EVSORT( NQ, QE ) CALL SEARCH( NQ, QE ) DO 60 I = 0, 3 60 E(I)= QE(I) C C ++++++ Inprove the Lanczos Vector via Inverse Iteration. ++++++ CALL INVERS(LV,NSB,SZ,MP,XL,XR,JH,SL,SR,BH,ENG,V) C ++++++ Inproved Ground State Energy ++++++ E(0)= ENG IF( LV .GE. 2 ) WRITE(6,200) ENG 200 FORMAT(1H ,' Ground State Energy by I.I. / C.G.=', F20.15) C RETURN END CCC SUBROUTINE MACRO1(LV,NSB,SZ,MP,XL,XR,JH,SL,SR, & BH,V,TE,ALP,BET,LK,ISW) INTEGER NV,LK,ISW,LV,NSB, MP(0:100,0:5) INTEGER XL(*), XR(*), SZ(*) REAL*8 PROD, ALP(0:999), BET(0:999), TE(0:999) REAL*8 SL(*), SR(*), BH(*), JH(*), V(0:100000,0:3) C NV = MP(NSB,3) C ++++++ A Lanczos Step ++++++ CALL INPRO( NV, 0, 1, PROD, V ) ALP(LK) = PROD PROD = -1.0 * PROD CALL ADDER( NV, 2, 1, 0, 1.0D0, PROD, V ) CALL INPRO( NV, 2, 2, PROD, V ) PROD = SQRT(PROD) BET(LK) = PROD CALL DIVER( NV, 2, 2, PROD, V ) CALL MULT(LV,1,2,NSB, SZ, MP, XL,XR, JH,SL,SR,BH,V) PROD = -1.0 * PROD CALL ADDER( NV, 1, 1, 0, 1.0D0, PROD, V ) CALL VMOVE( NV, 0, 2, V ) C ++++++ Second Path: Lanczos Vector Creation ++++++ IF( ISW .EQ. 1 ) CALL ADDER( NV, 3, 3, 0, 1.0D0, TE(LK), V) LK = LK + 1 RETURN END CCC SUBROUTINE MALANZ( LK, LV, ALP, BET, ENG, E2, VE ) INTEGER N, M, LK, LV, LNV, ISW, IW1(0:999),IERR REAL*8 ENG,E2, ALP(0:999), VE(0:999,0:4) REAL*8 EPS, BET(0:999), E(0:999), W1(0:5999) C N = LK EPS = -1.0 M = 5 LNV = 1000 ISW = -1 C C +----------------------------------------------------+ C | Diagonalization of the tri-diagonal matrix: | C | ASL routine: Please read NEC mannual. | C +----------------------------------------------------+ C -NEC- CALL DCSTSS( ALP,N,BET,EPS,E,M,VE,LNV,ISW,IW1,W1,IERR ) IF( IERR .NE. 0 ) THEN WRITE(6,100) IERR 100 FORMAT(1H ,' ERROR NO. in DCSTSS(ASL)',I10) END IF ENG = E(0) E2 = E(1) IF( LV .GE. 2 ) WRITE(6,200) (E(J),J = 0,2), LK 200 FORMAT(1H ,'***',3F20.15,'*** ',I3,'- Steps') RETURN END CCC SUBROUTINE MALAN2( LK, NQ, ALP, BET, QE ) INTEGER LK, NQ, IER REAL*8 TP, ALP(*), BET(*), QE(*), W1(0:5999) C -NEC-(Save Bet[n] because it is not conserved, actually.)- TP = BET(LK) C C +----------------------------------------------------+ C | Diagonalization of the tri-diagonal matrix: | C | ASL routine: Please read NEC mannual. | C +----------------------------------------------------+ C -NEC- CALL DCSTSN( ALP, LK, BET, -1.0D0, QE, LK/5 + 4, -1, W1, IER ) BET(LK) = TP NQ = NQ + LK/5 + 4 RETURN END CCC SUBROUTINE EVSORT( NQ, QE ) INTEGER IP, NQ REAL*8 MIN, TP, QE(0:*) DO 20 I = 0, NQ-2 IP = I MIN = QE(I) DO 10 J = I+1, NQ-1 IF( MIN .GT. QE(J) ) THEN MIN = QE(J) IP = J END IF 10 CONTINUE IF( I .NE. IP ) THEN TP = QE(I) QE(I) = QE(IP) QE(IP) = TP END IF 20 CONTINUE RETURN END CCC SUBROUTINE SEARCH( NQ, QE ) INTEGER MIP, NQ, NE, SEP REAL*8 EPS, DL, QE(0:*), E(0:100), DQ(0:100000) C C ................... Eigenvalue Pick-Up Levels ......... EPS = 1.0E-10 DL = 1.0E-5 C ................... Eigenvalue differential ......... DO 10 I = 0, NQ- 2 DQ(I)= QE(I+1)- QE(I) 10 CONTINUE C NE = 1 E(0) = QE(0) SEP = 1 20 CONTINUE IF( DQ(SEP) .LT. DL .AND. SEP .LT. NQ-2 ) THEN SEP = SEP + 1 GOTO 20 END IF 30 CONTINUE IF( DQ(SEP) .GT. EPS .AND. SEP .LT. NQ-2 ) THEN SEP = SEP + 1 GOTO 30 END IF IF( SEP .EQ. NQ-2 ) GOTO 50 MIP = SEP 40 IF( DQ(SEP) .LT. DL .AND. SEP .LT. NQ-2 ) THEN SEP = SEP + 1 IF( DQ(SEP) .LT. DQ(MIP) ) THEN MIP = SEP END IF GOTO 40 END IF IF( SEP .LT. NQ-2 ) THEN E(NE)= QE(MIP) NE = NE + 1 IF( NE .LE. 100) GOTO 30 END IF 50 CONTINUE C NQ = NE DO 60 I = 0, NE-1 QE(I)= E(I) 60 CONTINUE C RETURN END CCC CC C SUBROUTINE INVERS(LV,NSB,SZ,MP,XL,XR,JH,SL,SR,BH,LMD,V) INTEGER NV,MLR,LR,LV,NSB, MP(0:100,0:5), SZ(*) INTEGER XL(*), XR(*) REAL*8 BH(*), JH(*), SL(*), SR(*) REAL*8 DTR, LMD, DTP,EPS, V(0:100000,0:3) C IF( LV.GE.3 ) WRITE(6,*) '-vvvv- Sub. INVERS -vvvvv-' C +-----------------------------------------+ C | MLR - Maximum Inverse Iteration # | C | EPS - Eigenvector Convergence Parameter | C +-----------------------------------------+ LR = 0 MLR = 20 EPS = 0.5E-11 NV = MP(NSB,3) C CALL VMOVE( NV, 0, 3, V ) GOTO 20 10 CONTINUE CALL CGMETH( LV,NSB,SZ,MP,XL,XR,JH,SL,SR,BH,LMD,V ) 20 CALL VMOVE( NV, 1, 0, V ) CALL INPRO( NV, 1, 1, DTP, V ) DTP = SQRT(DTP) CALL DIVER( NV, 1, 1, DTP, V ) CALL MULT( LV, 2, 1,NSB, SZ, MP, XL,XR, JH,SL,SR,BH,V) CALL INPRO( NV, 2, 2, DTP, V ) DTP = SQRT(DTP) CALL DIVER( NV, 2, 2, DTP, V ) IF( LMD .LT. 0 ) THEN CALL ADDER( NV, 3, 2, 1, 1.0D0,1.0D0, V ) ELSE CALL ADDER( NV, 3, 2, 1,-1.0D0,1.0D0, V ) END IF CALL INPRO( NV, 3, 3, DTR, V ) IF( LV .GE. 3 ) WRITE(6,100) LMD, DTP, SQRT(DTR) 100 FORMAT(1H ,'Lmd=',E17.10,' Lmd2=',E17.10,' Difference=',E17.10) C ++++++ Ground State Energy Inprovement ++++++ IF( LMD .GT. 0 .AND. DTP .LT. LMD ) LMD = DTP IF( LMD .LT. 0 .AND. DTP .GT.-LMD ) LMD =-DTP C IF( ABS(LMD) .LE. EPS ) THEN WRITE(6,*)'Eigen value ',LMD,' might be Accidentaly Zero.' WRITE(6,*)'Please change chemical potential to get normal', & ' result.' IF( LV .GE. 3 ) THEN CALL ERRSTP( 'INVERS' ) END IF END IF C LR = LR + 1 IF( LR .GE. MLR ) THEN WRITE(6,*)'Eigen vector is not obtained in ',MLR,' steps.' WRITE(6,*)'Eigen Value = ', LMD, ' in INVERS' IF( LV .GE. 1 ) CALL ERRSTP( 'INVERS' ) END IF C ++++++ Convergence Check ++++++ IF( SQRT(DTR) .GE. EPS ) GOTO 10 C RETURN END CCC SUBROUTINE CGMETH(LV,NSB,SZ,MP,XL,XR,JH,SL,SR,BH,LMD,V) INTEGER NV, LR, LV,NSB, MP(0:100,0:5), MLR INTEGER XL(*), XR(*), SZ(*) REAL*8 BH(*), JH(*), SL(*), SR(*) REAL*8 DB, DR1, DR2, ALP, BET, EPS, DTP, LMD, V(0:100000,0:3) C C +-----------------------------------------+ C | MLR - Maximum Inverse Iteration # | C | EPS - Eigenvector Convergence Parameter | C +-----------------------------------------+ LR = 0 MLR = 20 EPS =1.0E-8 NV = MP(NSB,3) C (I) CALL INPRO( NV, 1, 1, DB, V) CALL VMOVE( NV, 2, 1, V) DO 10 I = 0,NV-1 V(I,0)= 0.0D0 10 CONTINUE C (II) CALL INPRO( NV, 1, 1, DR1, V) 20 CALL MULT( LV, 3, 2,NSB, SZ, MP, XL,XR, JH,SL,SR,BH,V) DO 30 I = 0, NV-1 V(I,3)= V(I,3) - LMD*V(I,2) 30 CONTINUE CALL INPRO( NV, 3, 2,DTP, V) ALP=DR1/DTP CALL ADDER( NV, 0, 0, 2, 1.0D0, ALP, V) CALL ADDER( NV, 1, 1, 3, 1.0D0,-ALP, V) CALL INPRO( NV, 1, 1, DR2, V) C IF( LV .GE. 5 ) WRITE(6,100) DR2 / DB, LR 100 FORMAT(1H ,' CG Conv. Parameter = ', E18.10,' -',I3,'-Step') IF( DR2/DB .LE. EPS ) GOTO 40 BET = DR2 / DR1 DR1 = DR2 LR = LR + 1 CALL ADDER( NV, 2, 1, 2,1.0D0, BET, V) IF( LR .GE. MLR ) THEN IF(LV.GE.4) WRITE(6,*)'Eigen vector is not obtained in ', & MLR,' CG steps.' GOTO 50 END IF GOTO 20 40 CONTINUE 50 IF( LV .GE. 5 ) WRITE(6,*) LR,' Rotations in CG' C RETURN END