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