Lanczos Diagonalization for S=1 Heisenberg Spin Ring
c123456789012345678901234567890123456789012345678901234567890123456789012
c ----- -----
c ----- Diagonalization of S = 1 Heisenberg Model -----
c ----- -----
c ----- H = (XX) + (YY) + DLT (ZZ) + AXI -----
c ----- + JNN(NNN) + .... -----
c ----- -----
c ----- '2003/04/04, Ver.0.1 by T.Nishino -----
c ----- -----
c
c ns: Config. # Reduced Config. #
c (Sz=0) (*About*)
c 3: 7 3
c 4: 19 7
c 5: 51 11
c 6: 141 23
c 7: 393 56
c 8: 1,107 138
c 9: 3,139 348
c 10: 8,953 895
c 11: 25,653 2,332
c 12: 73,789 6,149
c 13: 212,941 16,380
c 14: 616,227 44,016
c 15: 1,787,607 119,173 <--- vleh
c 16: 5,196,627 324,789
c 17: 15,134,931 890,290
c 18: 44,152,809 2,452,933
c 19: 128,996,853 6,789,308
c 20: 377,379,369 18,868,968
c 21: 1,105,350,729 52,635,749
c
c +---------------------------------------------------+
c | if NS > 15, increase vleh and vlen |
c | |
c | if NS > 32 modify nhx in "CRENUM" |
c +---------------------------------------------------+
c
module phys_vars ! Physical Parameters
integer :: TZ = 0 ! Total Sz
real(8) :: JNN = 1.0d0 ! NNN Interaction
real(8) :: DLT = 1.0d0 ! DLT*(Sz)*(Sz)
real(8) :: DSZ = 0.0D0 ! Onsite D-term
end module
c
c
module var_size ! Comon Variable Size
integer,parameter :: Lv = 1 ! Debug Level 0-4
integer,parameter :: ns = 10 ! Sys_Size: EVEN Number
integer,parameter :: nh = ns / 2 ! half size
integer,parameter :: nhx = 16 ! max hs (in CRENUM)
integer,parameter :: half = 3**nh
integer,parameter :: nsub = 9953 ! (in CRESTA)
integer,parameter :: ntp = 131441 ! (in CRESTA)
integer,parameter :: vleh = 119173 ! half vector length
integer,parameter :: vlen = vleh*2 ! vector length
integer,parameter :: nex = 999 ! Exceptional State #
integer,parameter :: ldiv = 141072 ! Loop Div. (in MULT)
end module
c
c
module com_vars ! Comon Variables
use var_size
integer NST, ST(0:vleh)
integer NK, NX, EXS(0:nex)
real(8) CS(0:ns), DIA(0:vleh), NRM(0:nex)
real(8) SN(0:ns), VEC(0:vlen,0:3)
end module
c
c
module lan_vars ! Used in Lanczos
integer LK, NV
real(8) ENG, ALP(0:999), MER(0:999,0:4)
real(8) BET(0:999), TER(0:999)
end module
c
c
program main; use com_vars; use phys_vars
c
integer BH(0:half,0:2) ; write(6,*); write(6,*)
c
write(6,*) '*** Ns, Tz, JNN, DLT', ns, TZ, JNN, DLT
c
do NK = 0, ns/2; write(6,*) 'Total Momentum = ',
& 2*NK, '/', ns, '*pi'
call CRENUM( BH )
call CRESTA( BH ); if( NST .eq. 1 ) then; write(6,*)
& 'tribial'; goto 10; end if
call ORTSET( )
call NORMAL( )
call CREDIA( )
call LANZOS( )
c call VECDMP( )
call SZDIST( BH )
10 continue
end do; end program main
ccc
subroutine CRENUM( BH ); use var_size
integer B0(0:2), T(0:nhx), A0, A1, A2, A3, A4, A5
integer B1(0:2),BH(0:half,0:2), A6, A7, A8, A9, AA, AB
integer B2(0:2), I, AC, AD, AE, AF
c
do I = 0, nhx; if( I .LT. nh ) then; T(I)= 1
else; T(I)= 0; end if; end do
B0(0)= 1; B0(1)= 0; B0(2)= 0
B1(0)= 0; B1(1)= 1; B1(2)= 0
B2(0)= 0; B2(1)= 0; B2(2)= 1
c
do I = 0, half - 1 ; A0 = mod(I ,3)
A1 = mod(I/3 ,3); A2 = mod(I/3** 2,3); A3 = mod(I/3** 3,3)
A4 = mod(I/3** 4,3); A5 = mod(I/3** 5,3); A6 = mod(I/3** 6,3)
A7 = mod(I/3** 7,3); A8 = mod(I/3** 8,3); A9 = mod(I/3** 9,3)
AA = mod(I/3**10,3); AB = mod(I/3**11,3); AC = mod(I/3**12,3)
AD = mod(I/3**13,3); AE = mod(I/3**14,3); AF = mod(I/3**15,3)
c
BH(I,0)= T( 0)*B0(A0) +T( 1)*B0(A1) +T( 2)*B0(A2) +T( 3)*B0(A3)
& + T( 4)*B0(A4) +T( 5)*B0(A5) +T( 6)*B0(A6) +T( 7)*B0(A7)
& + T( 8)*B0(A8) +T( 9)*B0(A9) +T(10)*B0(AA) +T(11)*B0(AB)
& + T(12)*B0(AC) +T(13)*B0(AD) +T(14)*B0(AE) +T(15)*B0(AF)
BH(I,1)= T( 0)*B1(A0) +T( 1)*B1(A1) +T( 2)*B1(A2) +T( 3)*B1(A3)
& + T( 4)*B1(A4) +T( 5)*B1(A5) +T( 6)*B1(A6) +T( 7)*B1(A7)
& + T( 8)*B1(A8) +T( 9)*B1(A9) +T(10)*B1(AA) +T(11)*B1(AB)
& + T(12)*B1(AC) +T(13)*B1(AD) +T(14)*B1(AE) +T(15)*B1(AF)
BH(I,2)= T( 0)*B2(A0) +T( 1)*B2(A1) +T( 2)*B2(A2) +T( 3)*B2(A3)
& + T( 4)*B2(A4) +T( 5)*B2(A5) +T( 6)*B2(A6) +T( 7)*B2(A7)
& + T( 8)*B2(A8) +T( 9)*B2(A9) +T(10)*B2(AA) +T(11)*B2(AB)
& + T(12)*B2(AC) +T(13)*B2(AD) +T(14)*B2(AE) +T(15)*B2(AF)
end do; end
ccc
subroutine CRESTA( BH ); use com_vars; use phys_vars
integer NL, I, PT(0:ntp), BH(0:half,0:2)
integer LZ, J, TP(0:ntp), NLS(-nh:nh)
integer K, FG(0:ntp), LS(0:nsub,-nh:nh)
c
do I =-nh, nh ; NL = 0 ! Lower Loop
do J = 0, half -1; LZ = BH(J,2) - BH(J,0) ! Polarization
c
if( LZ .eq. I ) then; LS(NL,I) = J; NL = NL + 1; end if
end do; NLS(I) = NL
end do
c ! Higher Loop
NST = 0 ! Lower Polarization
do I = 0, half -1; LZ = TZ -(BH(I,2)-BH(I,0))
c
if( abs(LZ) .le. nh ) then
do J = 0, NLS(LZ) -1; PT(J)= I*half + LS(J,LZ)
FG(J)= 0; end do
c
TP(0:NLS(LZ)-1)= PT(0:NLS(LZ)-1)
c
do K = 1, ns-1
do J = 0, NLS(LZ) -1; TP(J)= mod(TP(J)*3,3**ns)
& + TP(J)/3**(ns-1)
if( PT(J) .gt. TP(J) ) FG(J) = 1
end do; end do
c --- Select States ---
do J = 0, NLS(LZ) -1
if( FG(J) .eq. 0 ) then; ST(NST) = PT(J); NST = NST + 1
end if
end do
end if
end do
c
if( Lv .ge. 1 ) write(6,*) 'Total # of basis=', NST
if( Lv .ge. 3 ) then; do I = 0, NST -1
call STDISP( ST(I), 0.0D0 ); end do; end if
end
ccc
subroutine STDISP( stat, F ); use var_size
integer I, TP, stat, LW(0:ns); real(8) F
TP = stat
do I = ns -1, 0, -1; LW(I)= mod(TP,3)
TP = TP/3
end do; write(6,100) stat, F, (LW(I), I=0, ns-1)
100 format(' ',I14,' ---', F20.10,' ', 40I1)
end
ccc
subroutine ORTSET( ); use com_vars; integer I; real(8) PI
c
if( Lv .ge. 2 ) write(6,*) '--- Subroutine ORTSET ---'
c
PI = 4.0D0*atan(1.0D0)
c
do I = 0, ns-1; SN(I)= SIN( 2.0D0*I*NK*PI/ns )
CS(I)= COS( 2.0D0*I*NK*PI/ns )
if( Lv .ge. 2 ) write(6,100) SN(I), CS(I)
100 format(' ',' SIN=', F10.5, ' COS=', F10.5 )
end do; end
ccc
subroutine NORMAL( ); use com_vars
integer RT, TP(0:vleh); real(8) F
c
if( Lv .ge. 2 ) write(6,*) '--- Subroutine NORMAL ---'
c
VEC(0:NST-1,0)= 0.0D0; TP(0:NST-1)
& = ST(0:NST-1)
c
do I = 1, ns -1 ! --- Mark Self Similar States ---
do J = 0, NST -1
TP(J)= mod(TP(J)*3,3**ns) + TP(J)/3**(ns-1)
c
if( ST(J) .eq. TP(J) ) VEC(J,0) = 1.0D0
end do; end do
c --- Pick up Self Similar States ---
NX = 0
do I = 0, NST -1
if( VEC(I,0) .gt. 1.0E-8 ) then; EXS(NX)= I; NX = NX + 1
end if
end do
c --- Normalize ---
do I = 0, NX -1; RT = ST(EXS(I))
F = 1.0D0
do J = 1, ns -1; RT = RT/3**(ns-1) + mod(RT*3,3**ns)
if( RT .eq. ST(EXS(I)) ) F = F + CS(J)
end do
if( F .gt. 1.0E-8 ) then; NRM(I) = 1.0D0 / sqrt( F )
else; NRM(I) = 0.0D0
end if
end do
c
if( NX .gt. nex ) then
write(6,*)'Degenerate Basis # > ', nex
write(6,*)'...Increase 999 in the module var_size'
stop
end if
c
if( Lv .ge. 1 ) write(6,200) NX
200 format(' ',' # of Auto Correlated Basis=', I10 )
c
if( Lv .ge. 2 ) then
do I = 0, NX -1; call STDISP( ST(EXS(I)), NRM(I) ); end do
end if
end
ccc
subroutine CREDIA( ); use com_vars; use phys_vars; integer I,J
c
do J = 0, NST -1 ! *** Boundary Part ***
DIA(J)=
& (mod(ST(J)/3**(ns-2),3)-1)*(mod(ST(J)/3**(ns-1),3)-1)*DLT
& + (mod(ST(J)/3**(ns-1),3)-1)*(mod(ST(J) ,3)-1)*DLT
& + (mod(ST(J)/3**(ns-2),3)-1)*(mod(ST(J) ,3)-1)*DLT*JNN
& + (mod(ST(J)/3**(ns-1),3)-1)*(mod(ST(J)/3 ,3)-1)*DLT*JNN
& + (mod(ST(J)/3**(ns-2),3)-1)**2 *DSZ
& + (mod(ST(J)/3**(ns-1),3)-1)**2 *DSZ
end do
c
do I = 0, ns -3; do J = 0, NST -1
DIA(J)= DIA(J)
& + (mod(ST(J)/3**I,3)-1)*(mod(ST(J)/3**(I+1),3)-1)*DLT
& + (mod(ST(J)/3**I,3)-1)*(mod(ST(J)/3**(I+2),3)-1)*DLT*JNN
& + (mod(ST(J)/3**I,3)-1)**2 *DSZ
end do; end do
c
if( Lv .ge. 3 ) then
do I = 0, NST-1; write(6,100) I, DIA(I); end do
100 format(' ', I4, F15.7 )
end if
end
ccc
subroutine VECDMP( ); use com_vars; integer IP; real(8) FTP
if( Lv .eq. 0 ) return
if( NST .gt. 500 ) then; write(6,*) ' Too Many Basis '; return
end if
c
do I = 0, NST -1
if( NK.eq.0 .or. NK.eq.ns/2 ) then
VEC(I,3) = VEC(I,1)*VEC(I,1)
else; VEC(I,3) = VEC(I,1)*VEC(I,1) + VEC(NST+I,1)*VEC(NST+I,1)
end if; end do
C
do I = 0, MIN(25,NST) -1
FTP = 0.0D0
IP = 0
do J = 0, NST -1
if( FTP .LT. VEC(J,3) ) then; FTP = VEC(J,3); IP = J
end if; end do
c
if( FTP .gt. 1.0E-8 ) call STDISP( ST(IP), VEC(IP,3) )
VEC(IP,3) = 0.0D0
end do; end
ccc
subroutine SZDIST( BH ); use com_vars
integer TP, M, TH, TL, BH(0:half,0:2); real(8) FC, F2, F1, F0
c
M = ns /2
TP = 3** M
F0 = 0.0D0; F1 = 0.0D0; F2 = 0.0D0
c
if( NK.eq.0 .or. NK.eq.nh ) then
do I = 0, NST -1
TH = ST(I)/TP
TL = mod(ST(I),TP)
FC = VEC(I,1) * VEC(I,1)
F0 = F0 + FC*( BH(TH,0) + BH(TL,0) )
F1 = F1 + FC*( BH(TH,1) + BH(TL,1) )
F2 = F2 + FC*( BH(TH,2) + BH(TL,2) )
end do
else
do I = 0, NST -1
TH = ST(I)/TP
TL = mod(ST(I),TP)
FC = VEC(I,1) *VEC(I,1) + VEC(NST+I,1) *VEC(NST+I,1)
F0 = F0 + FC*( BH(TH,0) + BH(TL,0) )
F1 = F1 + FC*( BH(TH,1) + BH(TL,1) )
F2 = F2 + FC*( BH(TH,2) + BH(TL,2) )
end do
end if
c
write(6,100) F0/ns, F1/ns, F2/ns
100 format(' ','v 0 ^ Ratio = ', 3F15.7)
end
ccc
subroutine INPRO( NV, ID, IS, PROD); use com_vars
integer NV, ID, IS; real(8) PROD
PROD = 0.0D0; do I = 0, NV-1; PROD = PROD + VEC(I,ID)*VEC(I,IS)
end do; end
ccc
subroutine DISPLN( val ); use var_size
real(8) val(0:*); write(6,100) (val(I),I=0,ns-1)
100 format(' ',8F9.6)
end
ccc
subroutine VECINI( ); use com_vars; real(8) PROD
c
if( NK.eq.0 .or. ns.eq.NK*2 ) then; VEC(0:NST-1,0:3)= 0.0d0
c
do I = 0, NST-1
VEC(I,0)= (mod(I,13)+mod(I,19)+00.1) / sqrt(dble(NST))
end do
*VDIR NODEP
do I = 0, NX-1; VEC(EXS(I),0)
& = VEC(EXS(I),0)*NRM(I); end do
c
call INPRO( NST, 0, 0, PROD )
c
PROD = sqrt(PROD); VEC(0:NST-1,0) = VEC(0:NST-1,0)/PROD
c
else; VEC(0:NST*2-1,0:3)= 0.0d0
c
do I = 0, NST-1
VEC(I ,0)= (mod(I,13)+mod(I,19)+00.1) / sqrt(DBLE(NST))
VEC(I+NST,0)= (mod(I,17)+mod(I,23)+00.1) / sqrt(DBLE(NST))
end do
*VDIR NODEP
do I = 0, NX-1
VEC(EXS(I) ,0)= VEC(EXS(I) ,0)*NRM(I)
VEC(EXS(I)+NST,0)= VEC(EXS(I)+NST,0)*NRM(I); end do
c
call INPRO( NST*2, 0, 0, PROD )
c
PROD = sqrt(PROD); VEC(0:NST*2-1,0) = VEC(0:NST*2-1,0)/PROD
end if
end
ccc
subroutine SCALE( ID, IS ); use com_vars
integer ID, IS
c
if( Lv .ge. 3 ) write(6,*) '--- Subroutine SCALE ---'
c
if( NK.eq.0 .or. ns.eq.NK*2 ) then
VEC(0:NST-1,ID)= 0.0D0
*VDIR NODEP
do I = 0, NX -1
VEC(EXS(I),IS)= VEC(EXS(I),IS)*NRM(I); end do
else
VEC( 0:NST*2-1,ID)= 0.0D0
*VDIR NODEP
do I = 0, NX -1
VEC(EXS(I) ,IS)= VEC(EXS(I) ,IS)*NRM(I)
VEC(EXS(I)+NST,IS)= VEC(EXS(I)+NST,IS)*NRM(I); end do
end if
end
ccc
subroutine RESCAL( ID, IS ); use com_vars
integer ID, IS
c
if( Lv .ge. 3 ) write(6,*) '--- Subroutine RESCAL ---'
c
if( NK.eq.0 .or. ns.eq.NK*2 ) then
*VDIR NODEP
do I = 0, NX -1
VEC(EXS(I),ID)= VEC(EXS(I),ID)*NRM(I)
if( NRM(I) .gt. 1.0E-8 ) then
VEC(EXS(I),IS)= VEC(EXS(I),IS)/NRM(I)
end if
end do
else
*VDIR NODEP
do I = 0, NX -1
VEC(EXS(I) ,ID)= VEC(EXS(I) ,ID)*NRM(I)
VEC(EXS(I)+NST,ID)= VEC(EXS(I)+NST,ID)*NRM(I)
if( NRM(I) .gt. 1.0E-8 ) then
VEC(EXS(I) ,IS)= VEC(EXS(I) ,IS)/NRM(I)
VEC(EXS(I)+NST,IS)= VEC(EXS(I)+NST,IS)/NRM(I)
end if
end do
end if
end
ccc
subroutine MULT( ID, IS ); use com_vars; use phys_vars
integer IP, NW, ID, IS, VLE, IX
integer V0(0:ldiv), V2(0:ldiv), V4(0:ldiv), V5(0:ldiv)
real(8) F0(0:ldiv), F1(0:ldiv)
c
call SCALE( ID, IS )
c
do I = 0, NST-1, ldiv; IX = I; VLE = min(ldiv,NST)
c
V0(0:VLE-1)= ST(IX:IX+VLE-1) ! --- State Copy ---
c
do 80 N = 0, ns-1 ! --- Loop: (Shift)^N ---
c
if( Lv .GE. 3 ) write(6,*) 'Rotation Level =', N
c
do K = 0, ns-1; IP = K ! --- bond index ---
c ! --- NN_S+S- Part ---
call NNSPSM( IP, VLE, V0, V2, F0 )
call COMPRS( ST(NST-1), IX, VLE, NW, V2, F0, V4, V5, F1 )
call BISECT( NW, V2, V4 )
c
if( NK .eq. 0 .or. NK .eq. ns/2 ) then
*VDIR NODEP
do L = 0, NW -1
if( V4(L) .eq. ST(V2(L)) ) then
VEC(V5(L) ,ID) = VEC(V5(L) ,ID)
& + VEC(V2(L) ,IS)*F1(L)*CS(N)
end if; end do
else
*VDIR NODEP
do L = 0, NW -1
if( V4(L) .eq. ST(V2(L)) ) then
VEC(V5(L) ,ID) = VEC(V5(L) ,ID)
& + VEC(V2(L) ,IS)*F1(L)*CS(N)
& - VEC(V2(L)+NST,IS)*F1(L)*SN(N)
VEC(V5(L)+NST,ID) = VEC(V5(L)+NST,ID)
& + VEC(V2(L) ,IS)*F1(L)*SN(N)
& + VEC(V2(L)+NST,IS)*F1(L)*CS(N)
end if; end do
end if
c ! --- NNN_S+S- Part ---
call NNNPM( IP, VLE, V0, V2, F0 )
call COMPRS( ST(NST-1), IX, VLE, NW, V2, F0, V4, V5, F1 )
call BISECT( NW, V2, V4 )
c
if( NK .eq. 0 .or. NK .eq. ns/2 ) then
*VDIR NODEP
do L = 0, NW -1
if( V4(L) .eq. ST(V2(L)) ) then
VEC(V5(L) ,ID) = VEC(V5(L) ,ID)
& + VEC(V2(L) ,IS)*F1(L)*CS(N)*JNN
end if; end do
else
*VDIR NODEP
do L = 0, NW -1
if( V4(L) .eq. ST(V2(L)) ) then
VEC(V5(L) ,ID) = VEC(V5(L) ,ID)
& + VEC(V2(L) ,IS)*F1(L)*CS(N)*JNN
& - VEC(V2(L)+NST,IS)*F1(L)*SN(N)*JNN
VEC(V5(L)+NST,ID) = VEC(V5(L)+NST,ID)
& + VEC(V2(L) ,IS)*F1(L)*SN(N)*JNN
& + VEC(V2(L)+NST,IS)*F1(L)*CS(N)*JNN
end if; end do
end if
c --- NN_S-S+ Part ---
call NNSMSP( IP, VLE, V0, V2, F0 )
call COMPRS( ST(NST-1), IX, VLE, NW, V2, F0, V4, V5, F1 )
call BISECT( NW, V2, V4 )
C
if( NK .eq. 0 .or. NK .eq. ns/2 ) then
*VDIR NODEP
do L = 0, NW -1
if( V4(L) .eq. ST(V2(L)) ) then
VEC(V5(L) ,ID) = VEC(V5(L) ,ID)
& + VEC(V2(L) ,IS)*F1(L)*CS(N)
end if; end do
else
*VDIR NODEP
do L = 0, NW -1
if( V4(L) .eq. ST(V2(L)) ) then
VEC(V5(L) ,ID) = VEC(V5(L) ,ID)
& + VEC(V2(L) ,IS)*F1(L)*CS(N)
& - VEC(V2(L)+NST,IS)*F1(L)*SN(N)
VEC(V5(L)+NST,ID) = VEC(V5(L)+NST,ID)
& + VEC(V2(L) ,IS)*F1(L)*SN(N)
& + VEC(V2(L)+NST,IS)*F1(L)*CS(N)
end if; end do
end if
c --- NNN_S-S+ Part ---
call NNNMP( IP, VLE, V0, V2, F0 )
call COMPRS( ST(NST-1), IX, VLE, NW, V2, F0, V4, V5, F1 )
call BISECT( NW, V2, V4 )
C
if( NK .eq. 0 .or. NK .eq. ns/2 ) then
*VDIR NODEP
do L = 0, NW -1
if( V4(L) .eq. ST(V2(L)) ) then
VEC(V5(L) ,ID) = VEC(V5(L) ,ID)
& + VEC(V2(L) ,IS)*F1(L)*CS(N)*JNN
end if; end do
else
*VDIR NODEP
do L = 0, NW -1
if( V4(L) .eq. ST(V2(L)) ) then
VEC(V5(L) ,ID) = VEC(V5(L) ,ID)
& + VEC(V2(L) ,IS)*F1(L)*CS(N)*JNN
& - VEC(V2(L)+NST,IS)*F1(L)*SN(N)*JNN
VEC(V5(L)+NST,ID) = VEC(V5(L)+NST,ID)
& + VEC(V2(L) ,IS)*F1(L)*SN(N)*JNN
& + VEC(V2(L)+NST,IS)*F1(L)*CS(N)*JNN
end if; end do
end if
c
end do
c
if( N .eq. ns-1 ) GOTO 80 ! --- State Translation ---
c
do J = 0, VLE -1; V0(J)= mod(V0(J)*3,3**ns) + V0(J)/3**(ns-1)
end do
80 continue
end do
c
call RESCAL( ID, IS )
C --- Diagonal Part ---
if( NK.eq.0 .or. NK.eq.ns/2 ) then
VEC( 0: NST-1,ID)
& = VEC( 0: NST-1,ID) +DIA( 0: NST-1)*VEC( 0: NST-1,IS)
else
VEC( 0: NST-1,ID)
& = VEC( 0: NST-1,ID) +DIA( 0: NST-1)*VEC( 0: NST-1,IS)
VEC(NST:2*NST-1,ID)
& = VEC(NST:2*NST-1,ID) +DIA( 0: NST-1)*VEC(NST:2*NST-1,IS)
end if
C
if( Lv .ge. 3 ) then
write(6,*) 'Stopped by Debug Switch Lv=', Lv; stop
end if
end
ccc
subroutine NNSPSM( IP, VLE, V0, V2, F0 ); use var_size
integer IP, VLE, V0(0:*),V2(0:*)
integer P0, P1, T0, T1; real(8) F0(0:*)
c
c --- NN_S+[T1] S-[T0] Part ---
P0 = 3** IP; P1 = 3**(IP+1)
if( IP .eq. ns-1 ) then; P0 = 1 ; P1 = 3**(ns-1); end if
c
do I = 0, VLE -1 ; T0 = mod(V0(I)/P0,3)
T1 = mod(V0(I)/P1,3)
if( T1.NE.2 .and. T0.NE.0 ) then ; V2(I)= V0(I) + P1 - P0
F0(I)= 1.0D0
else ; F0(I)= 0.0D0
end if
end do; end
ccc
subroutine NNNPM( IP, VLE, V0, V2, F0 ); use var_size
integer IP, VLE, V0(0:*),V2(0:*)
integer P0, P1, T0, T1; real(8) F0(0:*)
c
c --- NNN_S+[T1] S-[T0] Part ---
P0 = 3** IP; P1 = 3**(IP+2)
if( IP .eq. ns-2 ) then; P0 = 1 ; P1 = 3**(ns-2); end if
if( IP .eq. ns-1 ) then; P0 = 3 ; P1 = 3**(ns-1); end if
c
do I = 0, VLE -1 ; T0 = mod(V0(I)/P0,3)
T1 = mod(V0(I)/P1,3)
if( T1.NE.2 .and. T0.NE.0 ) then ; V2(I)= V0(I) + P1 - P0
F0(I)= 1.0D0
else ; F0(I)= 0.0D0
end if
end do; end
ccc
subroutine NNSMSP( IP, VLE, V0, V2, F0 ); use var_size
integer IP, VLE, V0(0:*),V2(0:*)
integer P0, P1, T0, T1 ; real(8) F0(0:*)
c
c --- NN_S-[T1] S+[T0] Part ---
P0 = 3** IP; P1 = 3**(IP+1)
if( IP .eq. ns-1 ) then; P0 = 1 ; P1 = 3**(ns-1); end if
c
do I = 0, VLE -1 ; T0 = mod(V0(I)/P0,3)
T1 = mod(V0(I)/P1,3)
if( T1.ne.0 .and. T0.ne.2 ) then ; V2(I)= V0(I) - P1 + P0
F0(I)= 1.0D0
else ; F0(I)= 0.0D0
end if
end do; end
ccc
subroutine NNNMP( IP, VLE, V0, V2, F0 ); use var_size
integer IP, VLE, V0(0:*),V2(0:*)
integer P0, P1, T0, T1 ; real(8) F0(0:*)
c
c --- NNN_S-[T1] S+[T0] Part ---
P0 = 3** IP; P1 = 3**(IP+2)
if( IP .eq. ns-2 ) then; P0 = 1 ; P1 = 3**(ns-2); end if
if( IP .eq. ns-1 ) then; P0 = 3 ; P1 = 3**(ns-1); end if
c
do I = 0, VLE -1 ; T0 = mod(V0(I)/P0,3)
T1 = mod(V0(I)/P1,3)
if( T1.ne.0 .and. T0.ne.2 ) then ; V2(I)= V0(I) - P1 + P0
F0(I)= 1.0D0
else ; F0(I)= 0.0D0
end if
end do; end
ccc
subroutine COMPRS( MST, IX, VLE, NW, V2, F0, V4, V5, F1 )
integer MST, IX, VLE, NW, V2(0:*), V4(0:*), V5(0:*)
real(8) F0(0:*), F1(0:*)
C
NW = 0
do I = 0, VLE -1
if( abs(F0(I)) .gt. 1.0E-16 .and. V2(I) .le. MST ) then
V4(NW) = V2(I)
F1(NW) = F0(I)
V5(NW) = IX + I
NW = NW + 1
end if
end do; end
ccc
subroutine BISECT( NW, V2, V4 ); use com_vars
integer DW, NW, V2(0:*), V4(0:*)
C
DW = NST - NST/2
do I = 0, NW -1
if( V4(I) .ge. ST(DW) ) then; V2(I) = DW
else; V2(I) = 0
end if; end do
C
20 CONTINUE
DW = DW - DW/2
do I = 0, NW -1
if( V4(I) .ge. ST( MIN(NST-1,V2(I)+DW) ) )
& V2(I) = MIN(NST-1,V2(I)+DW)
end do
if( DW .NE. 1 ) GOTO 20
end
ccc
subroutine LANZOS( ); use com_vars; use lan_vars
integer LS; real(8) E2, E2P
c
if( NK.eq.0 .or.NK.eq.ns/2 ) then; NV = NST
else; NV = 2*NST
end if
VEC(0:NV-1,0:3)= 0.0d0
call VECINI( ); VEC(0:NV-1, 3)= VEC(0:NV-1,0)
c
call MULT( 1,0 )
c
LK = 0; LS = 0; E2 = 0.0D0
10 continue
E2P = E2
do J = 0, 4; call MACRO1( ); LK = LK + 1; end do
c
call MALANZ( E2 )
LS = LS + 1
if( LS.le.40 .and. ABS(E2P-E2).ge.1.0E-16) goto 10
LS = LS - 1
c write(6,100) ENG
100 format(' ','Eng= ', F20.15)
c
if( abs(ENG) .lt. 1.0E-6 ) then
write(6,*) 'Energy is accidentally Zero'
write(6,*) 'do not believe the following data'
return
end if
c
TER(0:LK) = MER(0:LK,0); VEC(0:NV-1,0)= VEC(0:NV-1,3)
VEC(0:NV-1,3)= TER( 0)*VEC(0:NV-1,3)
call MULT( 1,0 )
LK = 0
do I = 0, LS
do J = 0, 4 ; call MACRO1( ); LK = LK + 1
VEC(0:NV-1,3)
& = VEC(0:NV-1,3)+ TER(LK)*VEC(0:NV-1,0)
end do; end do
c
call INVERS( ); write(6,110) ENG
110 format(' ', 'Eng= ', F20.15)
end
ccc
subroutine MACRO1( ); use com_vars; use lan_vars; use phys_vars
real(8) PROD
c
call INPRO( NV, 0, 1, PROD ); ALP(LK) = PROD
PROD = -PROD
VEC(0:NV-1,2)= VEC(0:NV-1,1) + PROD*VEC(0:NV-1,0)
c
call INPRO( NV, 2, 2, PROD ); PROD = SQRT(PROD)
BET(LK) = PROD
VEC(0:NV-1,2) = VEC(0:NV-1,2) /PROD
c
call MULT( 1, 2 )
PROD = -PROD
VEC(0:NV-1,1)= VEC(0:NV-1,1) + PROD*VEC(0:NV-1,0)
VEC(0:NV-1,0)= VEC(0:NV-1,2)
end
ccc
subroutine MALANZ( E2 ); use lan_vars; use com_vars
integer N; real(8) E2, TP, E(0:999)
c
N = LK
TP = BET(N); call TRIDIG( ALP, N, BET, E, 5, MER ); BET(N)= TP
c
ENG = E(0) ; E2 = E(1)
if( Lv .GE. 0 ) write(6,110) (E(J),J=0,2)
110 format(' ','*****',3F20.15,'*****')
end
ccc
subroutine INVERS( ); use com_vars; use lan_vars
integer LROT; real(8) DTP, DTR
LROT = 0
VEC(0:NV-1,0)= VEC(0:NV-1,3)
goto 20
10 continue
call CGMETH( )
20 VEC(0:NV-1,1)= VEC(0:NV-1,0)
call INPRO( NV, 1, 1, DTP )
DTP = sqrt( DTP ); VEC(0:NV-1,1)= VEC(0:NV-1,1)/DTP
call MULT( 2, 1 )
call INPRO( NV, 2, 2, DTP )
DTP = sqrt( DTP ); VEC(0:NV-1,2)= VEC(0:NV-1,2)/DTP
c
if( ENG .lt. 0 ) then; VEC(0:NV-1,3)
& = VEC(0:NV-1,2)+ VEC(0:NV-1,1)
else ; VEC(0:NV-1,3)
& =-VEC(0:NV-1,2)+ VEC(0:NV-1,1)
end if
call INPRO( NV, 3, 3, DTR ); write(6,100) ENG, DTP, sqrt(DTR)
100 format(' ', 'Eng= ',F20.15, ' ? ', E20.10, ' D ',E20.10)
c
if( ENG .gt. 0 .and. DTP .lt. ENG ) ENG = DTP
if( ENG .lt. 0 .and. DTP .gt.-ENG ) ENG =-DTP
c
LROT = LROT + 1
if( LROT .ge. 20 ) then
write(6,*) 'Did not converge in INVERS'; stop
end if
if( sqrt(DTR) .ge. 0.5E-11 ) goto 10
end
ccc
subroutine CGMETH( ); use com_vars; use lan_vars
integer LR; real(8) EPS, DB, ALP_, BET_, DTP, DR1, DR2
c
call INPRO( NV, 1, 1, DB ); VEC(0:NV-1,2)= VEC(0:NV-1,1)
VEC(0:NV-1,0)= 0.0D0
EPS = 1.0E-8; LR = 0
c
call INPRO( NV, 1, 1, DR1 )
20 call MULT( 3, 2 ); VEC(0:NV-1,3)
& = VEC(0:NV-1,3) - ENG *VEC(0:NV-1,2)
call INPRO( NV, 3, 2, DTP )
ALP_ = DR1 / DTP ; VEC(0:NV-1,0)
& = VEC(0:NV-1,0) + ALP_*VEC(0:NV-1,2)
VEC(0:NV-1,1)
& = VEC(0:NV-1,1) - ALP_*VEC(0:NV-1,3)
call INPRO( NV, 1, 1, DR2 )
c
if( DR2/DB .LE. EPS ) goto 40
c
BET_ = DR2 / DR1; DR1 = DR2
LR = LR + 1 ; VEC(0:NV-1,2)
& = VEC(0:NV-1,1) + BET_*VEC(0:NV-1,2)
if( LR .GE. 20 ) then
c write(6,110)
110 format(' ','Did not converge in CG')
goto 50
end if
goto 20
40 continue
50 write(6,120) LR
120 format(' After ',I10,' Rotation in CG')
end
ccc
subroutine TRIDIG( A, N, B, E, M, VE )
integer MWN,NC,N,M,SP,BP,WN(0:999)
real(8) EPS,A(0:*),B(0:*),W1(0:5999),E(0:*)
real(8) MIN,MDL,BB(0:999),VE(0:999,0:2)
c
BP = 0; EPS = -1.0d0; MIN = abs(A(0)) + abs(B(0))
c
do I = 0, N-3
MIN = max(MIN,abs(B(I ))+abs(A(I+1))+abs(B(I+1)))
BB(I)= B(I)*B(I)
end do
MIN = max(MIN,abs(B(N-2))+abs(A(N-1)))
BB(N-2)= B(N-2)*B(N-2)
C
call COUNT(N,NC,-MIN,A,BB); MIN =-MIN; MWN = NC
call COUNT(N,NC,-MIN,A,BB); SP = 0; W1(SP)=-MIN; WN(SP)= NC
C
20 continue
MDL =(MIN+W1(SP)) / 2
if(W1(SP)-MIN.lt.EPS.or.W1(SP)-MDL.lt.abs(MDL)*1.0E-17.or.
& MDL-MIN.LT.abs(MDL)*1.0E-17) goto 25
C
call COUNT(N,NC,MDL,A,BB)
if(MWN-NC.eq.0)then; MIN = MDL; MWN = NC; goto 20; end if
C
SP = SP + 1; W1(SP)= MDL; WN(SP)= NC
goto 20
C
25 continue
C
do I = 0, MWN - WN(SP) - 1; E(BP+I)= MDL; end do
BP = BP + MWN - WN(SP)
C
if(BP.LT.M) then
40 continue
if(WN(SP).eq.WN(SP-1)) then; SP = SP-1; goto 40; end if
c
MIN = W1(SP); MWN = WN(SP); SP = SP - 1
goto 20
end if
C
call INVERT(N,E(0),A,B,W1,VE)
end
CCC
subroutine COUNT(N,NC,LAM,A,BB); integer I, N,NC
real(8) G(0:999),GRN,EPS,LAM,A(0:*),BB(0:*)
C
NC = 0; EPS = 1.0E-18; GRN = 1.0E+18; G(0) = LAM - A(0)
c
do I = 1,N- 1
if(abs(G(I-1)).gt.EPS) then; G(I) = LAM -A(I) -BB(I-1)/G(I-1)
else; G(I) = GRN
end if; end do
C
do I = 0,N-1; if(G(I).le.0) NC = NC + 1; end do
end
CCC
subroutine INVERT(N,E,A,B,W1,VE); integer N
real(8) PROD,MROD
real(8) E,A(0:*),B(0:*),W1(0:999,0:5),VE(0:999,0:2)
C
VE(0,0)= 1.0D0; VE(1:N-1,0)= 0.0D0
C
call CHRSKY(N,E,A,B,W1)
C
20 continue; VE(0:N-1,1)= VE(0:N-1,0)
c
call VSOLVE(N,W1,VE)
C
PROD = 0.0D0
do I = 0,N-1; PROD = PROD + VE(I,0)*VE(I,0); end do
C
PROD = sqrt(PROD); VE(0:N-1,0)= VE(0:N-1,0) / PROD
C
PROD = 0.0D0
MROD = 0.0D0
do I = 0,N-1
PROD = PROD +(VE(I,1)-VE(I,0))*(VE(I,1)-VE(I,0))
MROD = MROD +(VE(I,1)+VE(I,0))*(VE(I,1)+VE(I,0))
end do
C -WR- write(6, * )' MROD, PROD ', MROD, PROD
if(PROD.gt.1.0E-16.and.MROD.gt.1.0E-16) goto 20
end
CCC
subroutine CHRSKY(N,E,A,B,W1); integer N
real(8) E,A(0:*),B(0:*),W1(0:999,0:5)
C
W1(0,0)= A(0)- E; do I = 0,N - 2
W1(I ,1)= B(I)/W1(I,0)
W1(I+1,0)= A(I+1)- E - B(I)*W1(I,1)
end do; end
CCC
subroutine VSOLVE(N,W1,VE); integer N
real(8) VE(0:999,0:2),W1(0:999,0:5)
C
do I = 1,N-1; VE(I,0)= VE(I,0) - W1(I-1,1)*VE(I-1,0); end do
C
do I = 0,N-1; VE(I,0)= VE(I,0)/W1(I,0); end do
C
do I = N -2,0,-1; VE(I,0)= VE(I,0) - W1(I,1)*VE(I+1,0); end do
end