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