BSS-QMC for Half-filled Hubbard model

This program runs with NEC's ASL library
(Please contact me if you wish to know about ASL.)


*#RUN: F=/BMC1/RN30BIT(08) F=/BMC1/NEGU3D(07)
C     -----                                                 -----
C     -----  Q.M.C. Simulation of n-dimensional Negative U  -----
C     -----  Hubbard Model: Zero Temperature Version        -----
C     -----                 Open Boundary Condition         -----
C     -----  NL >>> Time Slice  /  NG >>> Time Group        -----
C     -----   L >>> Lattice Size                            ----- 
C 
      INTEGER   NS,NE, NL,NG, TST, MCS,L(0:2)
      REAL*8     BETA,  U,SU(0:1), DU(0:1)
      REAL*8      FT(0:256,0:256), ET(0:256,0:256)
      REAL*8      PR(0:256,0:256), PL(0:256,0:256)
C
C                      --- Lattice Size (Lx,Ly,Lz) ---
       L(0) =     1
       L(1) =     1
       L(2) =     2
C                      --- # of Lattice Point ---
         NS = L(0) * L(1) * L(2)
C
      IF( MOD(NS,2) .EQ. 1 ) THEN
          WRITE(6,*)' NS should be an even number '
          STOP
        END IF
C                      --- # of up spin electron ---
         NE =    NS / 2
C                      --- Trotter Number ---
         NL =   128
C                      --- QR decomposition frequency ---
         NG =     8
C                      --- # of MCS without observation ---
        TST =  2000
C                      --- # of MCS  ---
        MCS = 12004
C
       BETA = 0.125D0 * NL
          U = 8.000D0
C
               IF( MOD(NL,NG).NE.0 .OR. NL/NG .GT. 32 ) THEN
                   WRITE(6,*) 'Invalid Pair Nl, Ng ', NL, NG, NL/NG
                   STOP
                 END IF
C
      CALL CRETRI( L,NS,NE,NL,  BETA,      FT,ET,PR,PL       )
      CALL CREDIF(         NL,  BETA,    U,            SU,DU )
      CALL SIMULA( L,NS,NE,NL,NG,MCS,TST,U,FT,ET,PR,PL,SU,DU )
C
      WRITE(6, 99)   NS,L(0),L(1),L(2),NE,     BETA,    U
      WRITE(6, * ) 'End of Simulation'
   99 FORMAT(1H ,'NS,L_123,NE =',5I3,' BETA=',F8.5,' U=',F8.5)
C
      STOP
      END
CCC
CC
C
      SUBROUTINE  DTIME( TIME1, TIME2, ASCII )
      CHARACTER*8                      ASCII
      REAL*8             TIME1, TIME2
      IF(  TIME1 .GT. TIME2 ) THEN
          CALL WATCH( TIME2 )
        ELSE
          CALL WATCH( TIME1 )
        END IF
       WRITE(6,100) ABS( TIME1 - TIME2 ), ASCII
  100 FORMAT(1H ,' - - - - - - -    ', F27.10,' Sec has passed in',A10)
      RETURN
      END
CCC
      SUBROUTINE  WATCH( T1 )
      REAL*8             T1
      CALL CLOCK( T1 )
      RETURN
      END
CCC
      SUBROUTINE MATCLR( NA, NB, A )
      INTEGER            NA, NB
      REAL*8        A(0:256, 0:256)
      DO 20 J = 0, NB -8, 8
      DO 10 I = 0, NA -1
      A(I,J  )= 0.0D0
      A(I,J+1)= 0.0D0
      A(I,J+2)= 0.0D0
      A(I,J+3)= 0.0D0
      A(I,J+4)= 0.0D0
      A(I,J+5)= 0.0D0
      A(I,J+6)= 0.0D0
      A(I,J+7)= 0.0D0
   10 CONTINUE
   20 CONTINUE
      DO 40 J = (NB/8)*8, NB -1
      DO 30 I =        0, NA -1
        A(I,J)= 0.0D0
   30 CONTINUE
   40 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE  MINUS( NA, NB, A )
      INTEGER            NA, NB
      REAL*8        A(0:256, 0:256)
      DO 20 J = 0, NB -8, 8
      DO 10 I = 0, NA -1
      A(I,J  )= -A(I,J  )
      A(I,J+1)= -A(I,J+1)
      A(I,J+2)= -A(I,J+2)
      A(I,J+3)= -A(I,J+3)
      A(I,J+4)= -A(I,J+4)
      A(I,J+5)= -A(I,J+5)
      A(I,J+6)= -A(I,J+6)
      A(I,J+7)= -A(I,J+7)
   10 CONTINUE
   20 CONTINUE
      DO 40 J = (NB/8)*8, NB -1
      DO 30 I =        0, NA -1
        A(I,J)= -A(I,J)
   30 CONTINUE
   40 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE MATCPY( NA, NB, A, B )
      INTEGER            NA, NB
      REAL*8        A(0:256,0:256), B(0:256,0:256)
      DO 20 J = 0, NB -8, 8
      DO 10 I = 0, NA -1
      B(I,J  )= A(I,J  )
      B(I,J+1)= A(I,J+1)
      B(I,J+2)= A(I,J+2)
      B(I,J+3)= A(I,J+3)
      B(I,J+4)= A(I,J+4)
      B(I,J+5)= A(I,J+5)
      B(I,J+6)= A(I,J+6)
      B(I,J+7)= A(I,J+7)
   10 CONTINUE
   20 CONTINUE
      DO 40 J = (NB/8)*8, NB -1
      DO 30 I =        0, NA -1
        B(I,J)= A(I,J)
   30 CONTINUE
   40 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE MATMUL( NA, NB, NC, A, B, C )
      INTEGER            NA, NB, NC
      REAL*8 A(0:256,0:256), B(0:256,0:256), C(0:256,0:256)
      CALL MATCLR( NA, NC, C )
C     DO 30 K = 0, NB -8, 8
C     DO 20 J = 0, NC -1
C     DO 10 I = 0, NA -1
C       C(I,J)= C(I,J) + A(I, K  )*B(K  , J) + A(I, K+1)*B(K+1, J)
C    &                 + A(I, K+2)*B(K+2, J) + A(I, K+3)*B(K+3, J)
C    &                 + A(I, K+4)*B(K+4, J) + A(I, K+5)*B(K+5, J)
C    &                 + A(I, K+6)*B(K+6, J) + A(I, K+7)*B(K+7, J)
C  10 CONTINUE
C  20 CONTINUE
C  30 CONTINUE
C     DO 60 K = (NB/8)*8, NB -1
      DO 60 K =        0, NB -1
      DO 50 J =        0, NC -1
      DO 40 I =        0, NA -1
        C(I,J)= C(I,J) + A(I,K)*B(K,J)
   40 CONTINUE
   50 CONTINUE
   60 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE  MULDR( NA, NB, A, D )
      INTEGER            NA, NB
      REAL*8        A(0:256,0:256), D(0:256)
      DO 20 J = 0, NB -8, 8
      DO 10 I = 0, NA -1
      A(I,J  )= A(I,J  )*D(J  )
      A(I,J+1)= A(I,J+1)*D(J+1)
      A(I,J+2)= A(I,J+2)*D(J+2)
      A(I,J+3)= A(I,J+3)*D(J+3)
      A(I,J+4)= A(I,J+4)*D(J+4)
      A(I,J+5)= A(I,J+5)*D(J+5)
      A(I,J+6)= A(I,J+6)*D(J+6)
      A(I,J+7)= A(I,J+7)*D(J+7)
   10 CONTINUE
   20 CONTINUE
      DO 40 J = (NB/8)*8, NB -1
      DO 30 I =        0, NA -1
        A(I,J)= A(I,J)*D(J)
   30 CONTINUE
   40 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE MULDRC( NA, NB, A, D, B )
      INTEGER           NA, NB
      REAL*8       A(0:256,0:256), D(0:256), B(0:256,0:256)
      DO 20 J = 0, NB -8, 8
      DO 10 I = 0, NA -1
      B(I,J  )= A(I,J  )*D(J  )
      B(I,J+1)= A(I,J+1)*D(J+1)
      B(I,J+2)= A(I,J+2)*D(J+2)
      B(I,J+3)= A(I,J+3)*D(J+3)
      B(I,J+4)= A(I,J+4)*D(J+4)
      B(I,J+5)= A(I,J+5)*D(J+5)
      B(I,J+6)= A(I,J+6)*D(J+6)
      B(I,J+7)= A(I,J+7)*D(J+7)
   10 CONTINUE
   20 CONTINUE
      DO 40 J = (NB/8)*8, NB -1
      DO 30 I =        0, NA -1
        B(I,J)= A(I,J)*D(J)
   30 CONTINUE
   40 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE  MULDL( NA, NB, D, A )
      INTEGER           NA, NB
      REAL*8       A(0:256,0:256), D(0:256)
      DO 20 J = 0, NB -8, 8
      DO 10 I = 0, NA -1
      A(I,J  )= D(I)*A(I,J  )
      A(I,J+1)= D(I)*A(I,J+1)
      A(I,J+2)= D(I)*A(I,J+2)
      A(I,J+3)= D(I)*A(I,J+3)
      A(I,J+4)= D(I)*A(I,J+4)
      A(I,J+5)= D(I)*A(I,J+5)
      A(I,J+6)= D(I)*A(I,J+6)
      A(I,J+7)= D(I)*A(I,J+7)
   10 CONTINUE
   20 CONTINUE
      DO 40 J = (NB/8)*8, NB -1
      DO 30 I =        0, NA -1
        A(I,J)= D(I)*A(I,J)
   30 CONTINUE
   40 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE MULDLC( NA, NB, D, A, B )
      INTEGER            NA, NB
      REAL*8       A(0:256,0:256), D(0:256), B(0:256,0:256)
      DO 20 J = 0, NB -8, 8
      DO 10 I = 0, NA -1
      B(I,J  )= D(I)*A(I,J  )
      B(I,J+1)= D(I)*A(I,J+1)
      B(I,J+2)= D(I)*A(I,J+2)
      B(I,J+3)= D(I)*A(I,J+3)
      B(I,J+4)= D(I)*A(I,J+4)
      B(I,J+5)= D(I)*A(I,J+5)
      B(I,J+6)= D(I)*A(I,J+6)
      B(I,J+7)= D(I)*A(I,J+7)
   10 CONTINUE
   20 CONTINUE
      DO 40 J = (NB/8)*8, NB -1
      DO 30 I =        0, NA -1
        B(I,J)= D(I)*A(I,J)
   30 CONTINUE
   40 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE  MATTP( NA, NB, A, B )
      INTEGER           NA, NB
      REAL*8       A(0:256,0:256), B(0:256,0:256)
      DO 20 I = 0, NA -1
      DO 10 J = 0, NB -1
        B(J,I)= A(I,J)
   10 CONTINUE
   20 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE    QDR( NA, NB,  Q,  A )
      INTEGER            NA, NB
      REAL*8 Q(0:256,0:256),  R(0:256,0:256), A(0:256,0:256)
C
      CALL MATCLR( NA, NB, Q )
      CALL MATCLR( NB, NB, R )
      DO 50 I = 0, NB-1
      DO 10 J = 0, NA-1
        R(I,I)= R(I,I) + A(J,I)*A(J,I)
   10 CONTINUE
        R(I,I)=    SQRT( R(I,I) )
      DO 20 J = 0, NA-1
        Q(J,I)=   A(J,I)/R(I,I)
   20 CONTINUE
      IF( I .NE. NB-1 ) THEN
          DO 40 J = I+1, NB-1
          DO 30 K =   0, NA-1
            R(I,J)= R(I,J) + Q(K,I)*A(K,J)
   30     CONTINUE
          DO 35 K =   0, NA-1
            A(K,J)= A(K,J) - Q(K,I)*R(I,J)
   35     CONTINUE
   40     CONTINUE
        END IF
   50 CONTINUE
C
      RETURN
      END
CCC
      SUBROUTINE    LDQ( NA, NB,  A,  Q )
      INTEGER            NA, NB
      REAL*8 A(0:256,0:256),  Q(0:256,0:256), T(0:256,0:256)
      CALL  MATTP( NA, NB, A, T )
      CALL    QDR( NB, NA, A, T )
      CALL  MATTP( NB, NA, A, Q )
      RETURN
      END
CCC
      SUBROUTINE MATINV( NA, A )
      INTEGER            NA
      REAL*8    A(0:256,0:256)
C
      IF( NA .EQ. 1 ) THEN
          A(0,0)= 1.0D0 / A(0,0)
        ELSE
          CALL MATLDU( NA, A )
          CALL LDUINV( NA, A )
        END IF
C
      RETURN
      END
CCC
      SUBROUTINE MATLDU( N, A )
      INTEGER            N
      REAL*8            TP, A(0:256,0:256)
      DO 60 I = 0, N-1
           TP = 0.0D0
      DO 10 K = 0, I-1
           TP = TP + A(I,K)*A(K,I)
   10 CONTINUE
        A(I,I)= A(I,I) - TP
      DO 50 J = I+1, N-1
           TP = 0.0D0
      DO 20 K = 0, I-1
           TP = TP + A(J,K)*A(K,I)
   20 CONTINUE
        A(J,I)= A(J,I) - TP
           TP = 0.0D0
      DO 30 K = 0, I-1
           TP = TP + A(I,K)*A(K,J)
   30 CONTINUE
        A(I,J)= ( A(I,J)-TP ) / A(I,I)
   50 CONTINUE
   60 CONTINUE
      RETURN
      END
CCCC
      SUBROUTINE LRMINV( N, A )
      INTEGER            N
      REAL*8     TL, TR, A(0:256,0:256)
      DO 30 J =  1, N-1
      DO 20 I =  0, J-1
           TL =  A(J,I)
           TR =  A(I,J)
        A(J,I)= -A(J,I)
        A(I,J)= -A(I,J)
      IF( J .NE. N-1 ) THEN
          DO 10 K = J+1, N-1
            A(I,K)= A(I,K) - TR*A(J,K)
            A(K,I)= A(K,I) - TL*A(K,J)
   10     CONTINUE
        END IF
   20 CONTINUE
   30 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE MULTLR( N, A, B )
      INTEGER            N
      REAL*8   A(0:256,0:256), B(0:256,0:256)
      DO 40 J =  0, N-1
      DO 30 I =  0, N-1
      IF( J .GT. I ) THEN
            B(I,J)= 0.0D0
          DO 10 K = J, N-1
            B(I,J)= B(I,J) + A(I,K)*A(K,J)
   10     CONTINUE
        ELSE
            B(I,J)= A(I,J)
          DO 20 K = I+1, N-1
            B(I,J)= B(I,J) + A(I,K)*A(K,J)
   20     CONTINUE
        END IF
   30 CONTINUE
   40 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE LDUINV( N, A )
      INTEGER            N
      REAL*8   A(0:256,0:256), B(0:256,0:256)
      DO 20 I =   0, N-1
        A(I,I)= 1.0D0  / A(I,I)
      DO 10 J = I+1, N-1
        A(J,I)= A(J,I) * A(I,I)
   10 CONTINUE
   20 CONTINUE
      CALL LRMINV( N, A )
      DO 40 I =   0, N-1
      DO 30 J = I+1, N-1
        A(J,I)= A(J,I) * A(J,J)
   30 CONTINUE
   40 CONTINUE
      CALL MULTLR( N,    A, B )
      CALL MATCPY( N, N, B, A )
      RETURN
      END
CCC
CC
C
      SUBROUTINE CRETRI(L,     NS,NE,NL, BETA, FT,ET,PR,PL)
      INTEGER        M, L(0:2),NS,NE,NL, ADD
      REAL*8  TP1,  PI, FK(0:256,0:2)  , T(0:256,0:256), E(0:256)
      REAL*8  TP2,  DT, FT(0:256,0:256),ET(0:256,0:256)
      REAL*8  TP ,BETA, PR(0:256,0:256),PL(0:256,0:256)
C
           PI = 4.0D0 * ATAN(1.0D0)
      DO 30 I = 0, L(0)-1
      DO 20 J = 0, L(1)-1
      DO 10 K = 0, L(2)-1
          M   =  I*L(1)*L(2) + J*L(2) + K
       FK(M,0)= (I+1)*PI/(L(0)+1)
       FK(M,1)= (J+1)*PI/(L(1)+1)
       FK(M,2)= (K+1)*PI/(L(2)+1)
        E(M)  = -2*( COS(FK(M,0))+COS(FK(M,1))+COS(FK(M,2)) )
   10 CONTINUE
   20 CONTINUE
   30 CONTINUE
C
C-DB- DO 11 I = 0, NS-1
C-DB- WRITE(6,100)E(I),FK(I,0),FK(I,1),FK(I,2)
C  11 CONTINUE
C
      DO 50 I =   0, NS-2
            K = I
           TP = E(I)
      DO 40 J = I+1, NS-1
      IF( E(J) .LT. TP ) THEN
           TP = E(J)
            K =   J
        END IF
   40 CONTINUE
       TP     =  E(I)
        E(I)  =  E(K)
        E(K)  = TP
       TP     = FK(I,0)
       FK(I,0)= FK(K,0)
       FK(K,0)= TP
       TP     = FK(I,1)
       FK(I,1)= FK(K,1)
       FK(K,1)= TP
       TP     = FK(I,2)
       FK(I,2)= FK(K,2)
       FK(K,2)= TP
   50 CONTINUE
C
C-DB- DO 55 I = 0, NS-1
C-DB- WRITE(6,100)E(I),FK(I,0),FK(I,1),FK(I,2)
C 100 FORMAT(1H , F15.7, ' | ', 3F10.6)
C  55 CONTINUE
C
           TP = SQRT(2.0D0/(L(0)+1)) * SQRT(2.0D0/(L(1)+1))
     &                               * SQRT(2.0D0/(L(2)+1))
      DO 66 I = 0, NS  -1
      DO 64 J = 0, L(0)-1
      DO 62 K = 0, L(1)-1
      DO 60 M = 0, L(2)-1
        ADD   =  J*L(1)*L(2) + K*L(2) + M
      T(ADD,I)= TP*SIN((J+1)*FK(I,0)) * SIN((K+1)*FK(I,1))
     &                                * SIN((M+1)*FK(I,2))
   60 CONTINUE
   62 CONTINUE
   64 CONTINUE
   66 CONTINUE
C
      CALL MATCPY( NS, NE,  T, PR )
      CALL  MATTP( NS, NE, PR, PL )
      CALL MATCLR( NS, NS, FT     )
      CALL MATCLR( NS, NS, ET     )
C
           DT = BETA / NL
      DO 80 M = 0, NS-1
          TP1 = EXP( -1.0D0*DT*E(M) )
          TP2 = EXP(  1.0D0*DT*E(M) )
      DO 77 J = 0, NS-1
      DO 75 I = 0, NS-1
       FT(I,J)= FT(I,J) + T(I,M)*T(J,M)*TP1
       ET(I,J)= ET(I,J) + T(I,M)*T(J,M)*TP2
   75 CONTINUE
   77 CONTINUE
   80 CONTINUE
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE CREDIF( NL, BETA, U, SU,      DU    )
      INTEGER            NL
      REAL*8             DA, BETA, U, SU(0:1), DU(0:1)
C
      DA   = U*BETA/NL
      SU(0)= EXP(DA/2) - SQRT( EXP(DA)-1.0D0 )
      SU(1)= 1.0D0 /SU(0)
      DU(0)= SU(0) /SU(1) - 1.0D0
      DU(1)= SU(1) /SU(0) - 1.0D0
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE RINIT( NRING, P, Q, R, SHIFT )
      INTEGER           NRING, P, Q, R, SHIFT(0:9700)
                 P  =   NRING - 9689
                 Q  =   NRING - 4187
                 R  =   NRING -    0
      DO 10  I = 0,     NRING
      READ(8,*) SHIFT(I)
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE RANDOM( P, Q, R, NRING, ITP, SHIFT )
      INTEGER            P, Q, R, NRING, ITP, SHIFT(0:9700)
                 P = MOD(P +1,    NRING)
                 Q = MOD(Q +1,    NRING)
                 R = MOD(R +1,    NRING)
          ITP = IEOR(SHIFT(P),  SHIFT(Q))
      SHIFT(R)= ITP
      RETURN
      END
CCC
      SUBROUTINE DIST( NS,     NRING, P, Q, R, SHIFT, ISING )
      INTEGER          NS,     NRING, P, Q, R, SHIFT(0:9700)
      INTEGER         RND,     ISING(0:256)
      DO 80 I   =  0, NS-1
      CALL RANDOM( P, Q, R, NRING, RND, SHIFT )
      ISING(I  )= MOD( RND, 2 )
   80 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE DDDIST( NS, NL, NRING, P, Q, R, SHIFT, ISING )
      INTEGER            NS, NL, NRING, P, Q, R, SHIFT(0:9700)
      INTEGER                 M, ISING(0:256,0:256)
      DO 10 I = 0, NL-1
            M = I
             CALL  DIST( NS, NRING, P, Q, R, SHIFT, ISING(0,M) )
   10 CONTINUE
      RETURN
      END
CCC
CC
C
      SUBROUTINE  CREKN( NS     , ISING,   SU, DNK )
      INTEGER            NS     , ISING(0:256)
      REAL*8             SU(0:1),   DNK(0:256)
      DO 10 I =   0,   NS-1
        DNK(I)=  SU(ISING(I))
   10 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE  CREKI( NS     , ISING,   SU,  DIK )
      INTEGER            NS     , ISING(0:256)
      REAL*8             SU(0:1),   DIK(0:256)
      DO 10 I =   0,   NS-1
        DIK(I)=  SU(MOD(ISING(I)+1,2))
   10 CONTINUE
      RETURN
      END
CCC
CC
C
      SUBROUTINE   GCHK( NS,  G, OG )
      INTEGER            NS
      REAL*8       GFN, DIF, G(0:256,0:256), OG(0:256,0:256)
C
          GFN =   0.0D0
          DIF =   0.0D0
C
      DO 20 I =   0,NS-1
      DO 10 J =   0,NS-1
          GFN = MAX( GFN, ABS(G(J,I)) )
   10 CONTINUE
      DO 12 J =   0,NS-1
          DIF = MAX( DIF, ABS( G(J,I)-OG(J,I) ) )
   12 CONTINUE
   20 CONTINUE
      IF( DIF/GFN .GT. 1.0E-2 ) THEN
          WRITE(6,100) DIF, GFN
  100     FORMAT(1H ,' STOP by  ll G-OG ll and ll G ll ', 2F20.10 )
          CALL MATDSP( NS, NS, OG, ' OG^ ' )
          CALL MATDSP( NS, NS,  G, '  G^ ' )
          STOP
        END IF
      IF( DIF/GFN .GT. 1.0E-4 ) THEN
          WRITE(6,110) DIF, GFN, DIF/GFN
  110     FORMAT(1H ,' DIF ',F10.5,' MAXG ',F10.5,' Ratio ',F20.10 )
        END IF
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE MATDSP( NA, NB, A, ASCII )
      CHARACTER*8                   ASCII
      INTEGER            NA, NB
      REAL*8        A(0:256,0:256)
      IF( NB .LE. 11  ) THEN
           WRITE(6, * ) ' element in the Matrix ', ASCII
          DO 10 I = 0, NA-1
           WRITE(6,100) ( A(I,J), J=0, NB-1 )
  100     FORMAT(1H ,11F7.4)
   10     CONTINUE
        END IF
      RETURN
      END
CCC
CC
C
      SUBROUTINE  EVALG( NS, NE, UL, UR, G )
      INTEGER            NS, NE
      REAL*8    UL(0:256,0:256), UR(0:256,0:256)
      REAL*8     G(0:256,0:256),  A(0:256,0:256)
C
      CALL MATMUL( NE, NS, NE, UL, UR, G )
      CALL MATINV( NE,                 G )
      CALL MATMUL( NS, NE, NE, UR,  G, A )
      CALL MATMUL( NS, NE, NS,  A, UL, G )
      CALL  MINUS( NS, NS,             G )
C
      DO 10 I = 0, NS-1
        G(I,I)= G(I,I) + 1.0D0
   10 CONTINUE
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE LEFT( NS, NE, NL, NG ,ISING,SU,FT,PL,UL)
      INTEGER          NS, NE, NL, NG ,ISING(0:256,0:256)
      REAL*8  SU(0:1), FT(0:256,0:256),   UL(0:256,0:256,0:64)
      REAL*8          DNK(0:256)      ,    A(0:256,0:256)
      REAL*8           PL(0:256,0:256),    B(0:256,0:256)
C
      CALL MATCPY( NE, NS, PL, UL(0,0,NL/NG) )
C
      DO 20 M = NL/NG,  1, -1
      CALL MATMUL( NE, NS, NS, UL(0,0,M),FT,  A      )
      CALL  CREKN( NS,      ISING(0,M*NG-1), SU, DNK )
      CALL MULDRC( NE, NS, A, DNK, B                 )
C
      DO 10 I = 2, NG
      CALL MATMUL( NE, NS, NS,  B, FT, A             )
      CALL  CREKN( NS,      ISING(0,M*NG-I), SU, DNK )
      CALL MULDRC( NE, NS, A, DNK, B                 )
   10 CONTINUE
C
      CALL    LDQ( NE, NS, B, UL(0,0,M-1) )
   20 CONTINUE
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE BOOSTL( NS,NE,NG,TM,ISING,SU,FT,ET,UL,UR,G)
      INTEGER    MP, MN, NS,NE,NG,TM,    ISING(0:256,0:256)
      REAL*8    SU(0:1)  , FT(0:256,0:256), UL(0:256,0:256,0:64)
      REAL*8               ET(0:256,0:256), UR(0:256,0:256,0:64)
      REAL*8   DNK(0:256),  A(0:256,0:256),  B(0:256,0:256)
      REAL*8   DIK(0:256),  G(0:256,0:256),  C(0:256,0:256)
      REAL*8               OG(0:256,0:256)
C-WR-
C-WR- WRITE(6,*)'-  - BOOSTL',TM
      CALL  CREKN( NS,   ISING(0,TM), SU, DNK )
      CALL  CREKI( NS,   ISING(0,TM), SU, DIK )
      CALL MULDRC( NS, NS,    FT,DNK,  B )
      CALL MULDLC( NS, NS,   DIK, ET,  A )
      CALL MATMUL( NS, NS, NS, B,  G,  C )
      CALL MATMUL( NS, NS, NS, C,  A,  G )
C
       MP = TM / NG
       MN = MP +  1
      IF( MOD(TM,NG) .EQ. 0 ) THEN
          CALL MATMUL( NS,NS,NE, B, UR(0,0,MP), UR(0,0,MN) )
        ELSE
          CALL MATMUL( NS,NS,NE, B, UR(0,0,MN),  C )
          CALL MATCPY( NS,   NE, C, UR(0,0,MN)     )
        END IF
C
      IF( MOD(TM+1,NG) .EQ. 0 ) THEN
          CALL MATCPY( NS,NS,  G, OG   )
          CALL MATCPY( NS,NE, UR(0,0,MN),  A )
          CALL    QDR( NS,NE, UR(0,0,MN),  A )
          CALL  EVALG( NS,NE, UL(0,0,MN), UR(0,0,MN), G )
          CALL   GCHK( NS,     G, OG   )
        END IF
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE BOOSTR( NS, ISING, SU, FT, ET, G )
      INTEGER            NS, ISING(0:256)
      REAL*8   SU(0:1)  , FT(0:256,0:256)
      REAL*8              ET(0:256,0:256)
      REAL*8  DNK(0:256),  G(0:256,0:256), A(0:256,0:256)
      REAL*8  DIK(0:256),  B(0:256,0:256), C(0:256,0:256)
C
      CALL  CREKN( NS,  ISING, SU, DNK )
      CALL  CREKI( NS,  ISING, SU, DIK )
      CALL MULDRC( NS, NS, FT,DNK,   B )
      CALL MULDLC( NS, NS,DIK, ET,   A )
      CALL MATMUL( NS, NS, NS,  A,G, C )
      CALL MATMUL( NS, NS, NS,  C,B, G )
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE   POST( NS,NE,NG, TM,ISING,SU,FT, UL,UR,G)
      INTEGER    MP, MN, NS,NE,NG, TM, ISING(0:256)
      REAL*8    SU(0:1), FT(0:256,0:256), UL(0:256,0:256,0:64)
      REAL*8              A(0:256,0:256), UR(0:256,0:256,0:64)
      REAL*8   DNK(0:256),G(0:256,0:256),  B(0:256,0:256)
      REAL*8             OG(0:256,0:256),  C(0:256,0:256)
C
C-WR- WRITE(6,*)'-  - POST ',TM
      CALL  CREKN( NS, ISING, SU, DNK    )
      CALL MULDRC( NS,    NS, FT, DNK, B )
C
      MP = TM / NG + 1
      MN =      MP - 1
C
      IF( MOD(TM+1,NG) .EQ. 0 ) THEN
          CALL MATMUL( NE,NS,NS, UL(0,0,MP), B, UL(0,0,MN))
        ELSE
          CALL MATMUL( NE,NS,NS, UL(0,0,MN), B, C )
          CALL MATCPY( NE,NS, C, UL(0,0,MN)       )
        END IF
C
      IF( MOD(TM,NG) .EQ. 0 ) THEN
          CALL MATCPY( NS,NS, G, OG )
          CALL MATCPY( NE,NS,    UL(0,0,MN),  A            )
          CALL    LDQ( NE,NS, A, UL(0,0,MN)                )
          CALL  EVALG( NS,NE,    UL(0,0,MN), UR(0,0,MN), G )
          CALL   GCHK( NS,    G, OG )
        END IF
C
      IF( TM .EQ. 0 ) THEN
        END IF
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE SIMULA(LS,NS,NE, NL,NG,MCS,TST,U,FT,ET,PR,PL,SU,DU)
      INTEGER      LS(0:2),NS,NE, NL,NG,MCS,TST, SHIFT(0:9700)
      INTEGER     TM,NRING, P, Q,  R,    ISING(0:256,0:256)
      REAL*8   TIME1,TIME2,TB,TU, TX,TY,U
      REAL*8 SU(0:1),FT(0:256,0:256), UL(0:256,0:256,0:64)
      REAL*8 DU(0:1),ET(0:256,0:256), UR(0:256,0:256,0:64)
      REAL*8          G(0:256,0:256), PR(0:256,0:256)
      REAL*8                          PL(0:256,0:256)
C
      CALL WATCH( TIME1 )
      CALL WATCH( TIME2 )
C
      CALL MATCPY( NS, NE, PR, UR )
C
      NRING = 9700
      CALL  RINIT(         NRING, P, Q, R, SHIFT        )
      CALL DDDIST( NS, NL, NRING, P, Q, R, SHIFT, ISING )
      CALL   LEFT( NS, NE, NL,NG,ISING,SU,    FT, PL,UL )
      CALL  EVALG( NS, NE, UL,UR, G )
      CALL  DTIME( TIME1, TIME2, ' LEFT ' )
C
           TU = 0.0D0
           TB = 0.0D0
          TST = TST / 2
          MCS = MCS / 2
C
      WRITE(6,98)
   98 FORMAT(1H ,'    L0   |    T   |   KU   |',
     &            '   NN   |   CW   |   SC   |   E0   |' )
C
      DO 90 I = 0, MCS -1
      DO 10 J = 0,  NL -1
           TM = J
C 
      IF( TM .EQ. NL/2 .AND. I .GT. TST ) THEN
          CALL OBSERV( NS,     LS, U, G)
        END IF
C 
      CALL  WATCH( TX )
      CALL UPDATE( NS, NRING,P,Q,R, ISING(0,TM),SHIFT, DU, G)
      CALL  WATCH( TY )
        TU = TU +  TY - TX
      CALL BOOSTL( NS,NE,NG, TM, ISING, SU, FT,ET, UL, UR, G)
      CALL  WATCH( TX )
        TB = TB +  TX - TY
   10 CONTINUE
C
      DO 20 J = NL  , 1, -1
           TM = J
C >>>
      IF( TM .EQ. NL/2  .AND. I .GT. TST ) THEN
          CALL OBSERV( NS,     LS, U, G)
        END IF
C >>>
      CALL  WATCH(TX)
      CALL BOOSTR( NS, ISING(0,TM-1), SU, FT, ET, G)
      CALL  WATCH(TY)
        TB = TB + TY - TX
      CALL UPDATE( NS, NRING,P,Q,R, ISING(0,TM-1), SHIFT, DU, G)
      CALL  WATCH(TX)
        TU = TU + TX - TY
      CALL   POST( NS,NE, NG, TM-1, ISING(0,TM-1), SU,FT, UL,UR,G)
      CALL  WATCH(TY)
        TB = TB + TY - TX
C-DB- CALL  DEBUG( NS, NE, NL, TM-1, ISING, SU, FT, PL, PR, G)
   20 CONTINUE
   90 CONTINUE
C
      WRITE(6,110) 'UPDATE',TU
      WRITE(6,110) ' BOOST',TB
  110 FORMAT(1H ,'PASS TIME IN ',A8,' IS ',F20.10 )
      CALL  DTIME( TIME1, TIME2, 'SM-END' )
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE UPDATE( NS,NRING,P,Q,R,ISING,SHIFT,DU, G )
      INTEGER RN,SPN,SPI,NS,NRING,P,Q,R,ISING(0:256),SHIFT(0:9700)
      REAL*8     RND, M0,R0,  DU(0:1),     VL(0:256),VR(0:256)
      REAL*8                          G(0:256,0:256)
C
      DO 60  I = 0, NS-1
C
      CALL RANDOM( P, Q, R, NRING, RN, SHIFT )
           RND =  DBLE(RN)/2**30
           SPN = ISING( I)
           SPI =   MOD(SPN+1,2)
            R0 = 1.0D0 + ( 1.0D0-G(I,I) )*DU(SPI)
C
      IF( ( R0*R0*( DU(SPN)+1.0D0 ) ) .GT. RND ) THEN
          ISING(I)= SPI
               M0 = DU(SPI) / R0
          DO 10 J =    0, NS-1
             VL(J)= -G(J,I)
             VR(J)=  G(I,J)*M0
   10 CONTINUE
             VL(I)= VL(I) + 1.0D0
C
          DO 30 M = 0, NS-4, 4
          DO 20 L = 0, NS-1
          G(L,M  )= G(L,M  ) - VL(L) * VR(M  )
          G(L,M+1)= G(L,M+1) - VL(L) * VR(M+1)
          G(L,M+2)= G(L,M+2) - VL(L) * VR(M+2)
          G(L,M+3)= G(L,M+3) - VL(L) * VR(M+3)
   20     CONTINUE
   30     CONTINUE
C
          DO 50 M = (NS/4)*4, NS-1
          DO 40 L = 0, NS-1
            G(L,M)= G(L,M) - VL(L) * VR(M)
   40     CONTINUE
   50     CONTINUE
C
        END IF
C
   60 CONTINUE
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE DEBUG( NS, NE, NL, TM, ISING, SU, FT, PL,PR, G)
      INTEGER       II, NS, NE, NL, TM,       ISING(0:256,0:256)
      REAL*8   SU(0:1)        , FT(0:256,0:256), PR(0:256,0:256)
      REAL*8  DNK(0:256)      ,  A(0:256,0:256), PL(0:256,0:256)
      REAL*8    B(0:256,0:256),  R(0:256,0:256),  Q(0:256,0:256)
      REAL*8                     L(0:256,0:256),  O(0:256,0:256)
      REAL*8                     G(0:256,0:256), NG(0:256,0:256)
C
      CALL MATCPY( NS, NE, PR, R )
C
      DO 10 I = 0, TM-1
           II = I
      CALL  CREKN( NS, ISING(0,II), SU, DNK )
      CALL MULDRC( NS, NS,  FT,    DNK,   B )
      CALL MATMUL( NS, NS,  NE,  B,  R,   A )
      CALL MATCPY( NS, NE,           A,   R )
   10 CONTINUE
      CALL    QDR( NS, NE,   Q,  R )
C
      CALL MATCPY( NE, NS,  PL,  L )
C
      DO 20 I = NL-1, TM, -1
           II = I
      CALL  CREKN( NS, ISING(0,II), SU, DNK )
      CALL MULDRC( NS, NS,  FT,    DNK,   B )
      CALL MATMUL( NE, NS,  NS,  L,  B,   A )
      CALL MATCPY( NE, NS,           A,   L )
   20 CONTINUE
      CALL    LDQ( NE, NS,   L,  O )
C
      WRITE(6,*) '-DB- IN DEBUG ',TM
      CALL  EVALG( NS, NE,  O, Q, NG )
      CALL MATDSP( NS, NS,  G, 'UPDAT^')
      CALL MATDSP( NS, NS, NG, 'REEVA^')
      CALL   GCHK( NS,  G, NG          )
C
      RETURN
      END
CCC
CC
C
      SUBROUTINE OBSERV( NS, LS, U, G )
      INTEGER            NS, LS(0:2)
      REAL*8      SC, KU, T, CW, U, L0, NN, E
      REAL*8     G(0:256,0:256)
C                                    --- Electron Density ---
C-DB- CALL NUMBER( NS,     FN, G )
C                                    --- Kinetic Energy ---
      CALL TRANSF( NS, LS,  T, G )
C                                    --- Local Moment ---
      CALL LOCALM( NS,     L0, G )
C                                    --- (n_iup-1/2)(n_idw-1/2)> ---
      CALL ONSITE( NS,     KU, G )
C                                    --- (ni-1)(nj-1)> ---
      CALL NNNEIG(     LS, NN, G )
C                                    --- CDW Order Parameter ---
      CALL CDWCRR(     LS, CW, G )
C                                    --- S.C. Order Parameter ---
      CALL SUPCRR(     LS, SC, G )
C
      E = T - U*KU
C
C      WRITE( 6,100)      L0, T, KU, NN, CW, SC, E
       WRITE( 7,110)      L0, T, KU, NN, CW, SC, E
  100 FORMAT(1H ,   7F9.5 )
  110 FORMAT(1H ,   9F8.4 )
C
      RETURN
      END
CCC
      SUBROUTINE NUMBER( NS, FN,  G )
      INTEGER            NS
      REAL*8                 FN,  G(0:256,0:256)
C
           FN = 0.0D0
      DO 10 I =   0 , NS -1
           FN = FN + 1.0D0 - G(I,I)
   10 CONTINUE
           FN = FN / NS
C
      RETURN
      END
CCC
      SUBROUTINE TRANSF(  NS, LS, T, G )
      INTEGER     IS, ID, NS, LS(0:2)
      REAL*8              T,  G(0:256,0:256)
C
            T = 0.0D0
C
      DO 30 I = 0, LS(0) -  2
      DO 20 J = 0, LS(1) -  1
      DO 10 K = 0, LS(2) -  1
           IS = (I+1)*LS(1)*LS(2) + J*LS(2) + K
           ID =  I   *LS(1)*LS(2) + J*LS(2) + K
            T =  T + 2 * G(IS,ID)
   10 CONTINUE
   20 CONTINUE
   30 CONTINUE
C
      DO 60 I = 0, LS(0) -  1
      DO 50 J = 0, LS(1) -  2
      DO 40 K = 0, LS(2) -  1
           IS = I*LS(1)*LS(2) +  J   *LS(2) + K
           ID = I*LS(1)*LS(2) + (J+1)*LS(2) + K
            T =  T + 2 * G(IS,ID)
   40 CONTINUE
   50 CONTINUE
   60 CONTINUE
C
      DO 90 I = 0, LS(0) -  1
      DO 80 J = 0, LS(1) -  1
      DO 70 K = 0, LS(2) -  2
           IS = I*LS(1)*LS(2) + J*LS(2) + K
           ID = I*LS(1)*LS(2) + J*LS(2) + K+1
            T =  T + 2 * G(IS,ID)
   70 CONTINUE
   80 CONTINUE
   90 CONTINUE
C
            T = T / NS
C
      RETURN
      END
CCC
      SUBROUTINE LOCALM( NS, L0, G )
      INTEGER            NS
      REAL*8                 L0, G(0:256,0:256)
C
           L0 = 0.0D0
      DO 10 I = 0, NS-1
           L0 = L0 + 1.50D0 * G(I,I) * ( 1.0D0 - G(I,I) )
   10 CONTINUE
           L0 = L0 / NS
C
      RETURN
      END
CCC
      SUBROUTINE ONSITE( NS, KU, G )
      INTEGER            NS
      REAL*8                 KU, G(0:256,0:256)
C
           KU =  0.0D0
      DO 10 I =  0, NS-1
           KU = KU + ( 0.50D0-G(I,I) )**2
   10 CONTINUE
           KU = KU / NS
C
      RETURN
      END
CCC
      SUBROUTINE NNNEIG( LS, NN, G )
      INTEGER            LS(0:2)
      REAL*8         NN,  G(0:256,0:256)
C
           NN = 0.0D0
C
      DO 30 I = 0, LS(0) -  2
      DO 20 J = 0, LS(1) -  1
      DO 10 K = 0, LS(2) -  1
           IS = (I+1)*LS(1)*LS(2) + J*LS(2) + K
           ID =  I   *LS(1)*LS(2) + J*LS(2) + K
           NN = NN + ( 1.0D0-2*G(IS,IS) ) * ( 1.0D0-2*G(ID,ID) )
     &             +         2*G(IS,ID)*G(ID,IS)
   10 CONTINUE
   20 CONTINUE
   30 CONTINUE
C
      DO 60 I = 0, LS(0) -  1
      DO 50 J = 0, LS(1) -  2
      DO 40 K = 0, LS(2) -  1
           IS = I*LS(1)*LS(2) +  J   *LS(2) + K
           ID = I*LS(1)*LS(2) + (J+1)*LS(2) + K
           NN = NN + ( 1.0D0-2*G(IS,IS) ) * ( 1.0D0-2*G(ID,ID) )
     &             +         2*G(IS,ID)*G(ID,IS)
   40 CONTINUE
   50 CONTINUE
   60 CONTINUE
C
      DO 90 I = 0, LS(0) -  1
      DO 80 J = 0, LS(1) -  1
      DO 70 K = 0, LS(2) -  2
           IS = I*LS(1)*LS(2) + J*LS(2) + K
           ID = I*LS(1)*LS(2) + J*LS(2) + K+1
           NN = NN + ( 1.0D0-2*G(IS,IS) ) * ( 1.0D0-2*G(ID,ID) )
     &             +         2*G(IS,ID)*G(ID,IS)
   70 CONTINUE
   80 CONTINUE
   90 CONTINUE
C
           NN = NN / ( 3*LS(0)*LS(1)*LS(2)
     &                  -LS(0)*LS(1)-LS(1)*LS(2)-LS(2)*LS(0) )
C
      RETURN
      END
CCC
      SUBROUTINE CDWCRR( LS, CW, G )
      INTEGER    IS, NS, LS(0:2), SN(0:256)
      REAL*8         CW,  G(0:256,0:256)
C
      DO 30 I = 0, LS(0) -  1
      DO 20 J = 0, LS(1) -  1
      DO 10 K = 0, LS(2) -  1
           IS =  I   *LS(1)*LS(2) + J*LS(2) + K
        SN(IS)=(-1)**(I+J+K)
   10 CONTINUE
   20 CONTINUE
   30 CONTINUE
C
           NS = LS(0)*LS(1)*LS(2)
           CW = 0.0D0
      DO 50 I = 0, NS-1
      DO 40 J = 0, NS-1
           CW = CW +SN(I)*SN(J)*( 4*(1.0D0-G(I,I))*(1.0D0-G(J,J))
     &                           -2*       G(I,J) *       G(J,I) )
   40 CONTINUE
   50 CONTINUE
C
      DO 60 I = 0, NS-1
           CW = CW + 2*G(I,I)
   60 CONTINUE
C
           CW = CW / (NS*NS)
      RETURN
      END
CCC
      SUBROUTINE SUPCRR( LS, SC, G )
      INTEGER        NS, LS(0:2)
      REAL*8         SC,  G(0:256,0:256)
C
           NS = LS(0)*LS(1)*LS(2)
           SC = 0.0D0
      DO 20 I = 0, NS-1
      DO 10 J = 0, NS-1
           SC = SC + G(J,I)*G(J,I)
   10 CONTINUE
   20 CONTINUE
C
      DO 30 I = 0, NS-1
           SC = SC + (1.0D0-G(I,I))*(1.0D0-G(I,I)) 
     &             -        G(I,I) *       G(I,I)
   30 CONTINUE
C
           SC = SC / (NS*NS)
      RETURN
      END