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