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