This is the source code of CTMRG for q=5 potts
model
represented as a 5**4-vertex model.
This program runs with NEC's ASL library
(Please contact me if you wish to know about ASL.)
C123456789012345678901234567890123456789012345678901234567890123456789012
C ----- -----
C ----- Corner Transfer Matrix Renormalization Group Method -----
C ----- for the Square Lattice q-state Potts Model at Tc -----
C ----- -----
C ----- '96/ 3/18, Ver.0.9 by T.Nishino -----
C ----- -----
C ----- W---@---P---@---P---@---W @ --- Renormalized -----
C ----- | | | | Link Variable -----
C ----- @ o o | -----
C ----- | | | | o --- Bare Link -----
C ----- P---o---V---o---V---o---P Variable -----
C ----- | | | | -----
C ----- @ o o | W --- Corner Weight -----
C ----- | | | | -----
C ----- P---o---V---o---V---o---P P --- Vertex -----
C ----- | | | | Operator -----
C ----- @ o o | -----
C ----- | | | | V --- Vertex Weight -----
C ----- W---@---P---@---P---@---W -----
C
PARAMETER( iq=5, MXX=100 )
INTEGER M, N, MX, L
REAL*8 EKC, ER, W(0:iq*MXX), VW(0:iq**4), SW(0:2*iq**4)
REAL*8 C(0:iq*MXX,0:iq*MXX), P(0:iq*MXX,0:iq*MXX,0:iq-1)
REAL*8 R(0:iq*MXX,0:iq*MXX), O(0:iq*MXX,0:iq*MXX,0:iq-1)
C
CALL GETPRM( L, MX, EKC )
CALL INITTM( M, VW, SW, C, O, EKC )
K = 1
GOTO 20
10 CONTINUE
CALL EXP1( M, W, P, O, R )
CALL OBSERV( M, W, O, C, R, VW, SW )
C --- EXPCTM in Observ ---
CALL EXPVTO( M, P, O, VW )
20 CONTINUE
CALL MATCPY( iq*M, C, R )
CALL DIAGAL( iq*M, W, R )
CALL SELECT( M, N, W, MX, ER )
CALL DIAGPA( iq*M, N, W, C, R )
CALL RENORM( M, N, P, O, R )
CALL NORMW( N, W )
CALL NORMP( N, P )
M = N
K = K + 1
IF( K .LT. L ) GOTO 10
STOP
END
CCC
CC
C
SUBROUTINE GETPRM( L, MX, EKC )
PARAMETER( iq=5, MXX=100 )
INTEGER L, MX
REAL*8 EKC
Write(6,*) 'Input Max. Iteration # L'
Read(5,*) L
C L = 10
C Write(6,*) 'Input Boltsmann Fact. Exp(K)'
C Read(5,*) EKC
EKC = SQRT(1.0D0*iq) + 1.0D0
Write(6,*) 'Input Maximum Dim. mx >= m, '
Read(5,*) MX
C MX = 10
Write(6,100) L, MX, EKC
100 FORMAT(1H , 'L =',I4,' MX =',I3,' Exp(K) =',F15.7 )
IF( MX .GT. MXX ) THEN
Write(6,*) 'MX is larger than MXX'
STOP
END IF
RETURN
END
CCC
CC
C
SUBROUTINE INITTM( M, VW, SW, C, O, EKC )
PARAMETER( iq=5, MXX=100 )
INTEGER M, A, B, BC
REAL*8 EKC, VW(0:*), C(0:iq*MXX,0:*)
REAL*8 SW(0:*), O(0:iq*MXX,0:iq*MXX,0:*), U(0:iq-1,0:iq-1)
M = 1
DO 10 I = 0, iq-1
DO 10 J = 0, iq-1
IF( J .EQ. I ) THEN
U(J,I)= EKC + SQRT( (EKC+iq-1.0D0) * (EKC-1.0D0) )
ELSE
U(J,I)= 1.0D0
END IF
10 CONTINUE
C
C Fixed Boundary Condition / Free Boundary Condition
BC = 1
C BC = iq
C i i i i
C | | | |
C C---j O---k L---V---k L---S---k
C | | |
C j j j
DO 50 I = 0, iq-1
DO 50 J = 0, iq-1
C(I,J)= 0.0D0
DO 20 L = 0, BC-1
C(I,J)= C(I,J)+ U(I,L)*U(J,L)
20 CONTINUE
DO 50 K = 0, iq-1
O(I,J,K)= 0.0D0
DO 30 L = 0, BC-1
O(I,J,K)= O(I,J,K)+U(I,L)*U(J,L)*U(K,L)
30 CONTINUE
DO 50 L = 0, iq-1
A = iq*iq*iq*I + iq*iq*J + iq*K + L
SW(A)= U(I,0)*U(J,0)*U(K,0)*U(L,0)
SW(A+iq**4)= U(I,1)*U(J,1)*U(K,1)*U(L,1)
VW(A)= 0.0D0
DO 40 B = 0, iq-1
VW(A)= VW(A) + U(I,B)*U(J,B)*U(K,B)*U(L,B)
40 CONTINUE
50 CONTINUE
RETURN
END
CCC
CC
C
SUBROUTINE EXP1( M, W, P, O, R )
PARAMETER( iq=5, MXX=100 )
INTEGER M
REAL*8 W(0:*), P(0:iq*MXX,0:iq*MXX,0:*)
REAL*8 R(0:iq*MXX,0:*), O(0:iq*MXX,0:iq*MXX,0:*)
C
C k k
C | |
C i---O---j = i---P---j * W(j)
DO 10 K = 0, iq-1
DO 10 J = 0, M-1
DO 10 I = 0, M-1
O(I,J,K)= P(I,J,K) * W(J)
10 CONTINUE
C L n
C L n | |
C i---R---k = i---O---j---P---k
C
CALL MATCLR( iq*M, R )
DO 20 N = 0, iq-1
DO 20 K = 0, M-1
DO 20 J = 0, M-1
DO 20 L = 0, iq-1
DO 20 I = 0, M-1
R(L*M+I,N*M+K)= R(L*M+I,N*M+K) + O(I,J,L)*P(J,K,N)
20 CONTINUE
RETURN
END
CCC
CC
C
SUBROUTINE OBSERV( M, W, O, C, R, VW, SW )
PARAMETER( iq=5, MXX=100 )
INTEGER M
REAL*8 MG, Y, Z, Z1, SW(0:*), R(0:iq*MXX,0:*), C(0:iq*MXX,0:*)
REAL*8 E, Y0, Z0, VW(0:*), W(0:*), O(0:iq*MXX,0:iq*MXX,0:*)
C
CALL EXPCTM( M, C, R, SW )
CALL TRACE1( M, W, C, R, Y0 )
CALL QUADRA( M, C, O(0,0,0) )
CALL EXPCTM( M, C, R, SW(iq**4) )
CALL QUADRA( M, C, O(0,0,1) )
CALL EXPCTM( M, C, R, VW )
CALL TRACE1( M, W, C, R, Y )
MG = Y0 / Y
CALL QUADRA( M, C, R )
CALL TRACE2( M, R, O(0,0,0), Z0 )
CALL TRACE2( M, R, O(0,0,1), Z1 )
CALL TRACE3( M, R, Z )
E = ( Z0 + (iq-1)*Z1 ) / Z
C
C Write(6,100) MG, E
Write(6,100) (iq*MG-1.0D0) / (iq-1.0D0),
& E - (1.0D0+1.0D0 / SQRT(1.0D0*iq) ) / 2.0D0
100 FORMAT(1H ,F12.10,' ',F13.10)
RETURN
END
CCC
SUBROUTINE TRACE1( M, W, C, R, Z )
PARAMETER( iq=5, MXX=100 )
INTEGER M
REAL*8 Z, W(0:*), C(0:iq*MXX,0:*), R(0:iq*MXX,0:*)
C
C C-----+---W(k)
C | | |
C | | |
C +-----+--n--+
C | | |
C | L |
C W(i)----+-----R
Z = 0.0D0
DO 10 N = 0, iq-1
DO 10 L = 0, iq-1
DO 10 K = 0, M-1
DO 10 I = 0, M-1
Z = Z + C(L*M+I,N*M+K) * R(L*M+I,N*M+K) * W(I) * W(K)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE TRACE2( M, R, O, Z )
PARAMETER( iq=5, MXX=100 )
INTEGER M
REAL*8 Z, R(0:iq*MXX,0:*), O(0:iq*MXX,0:*)
Z = 0.0D0
DO 10 J = 0, iq*M-1
DO 10 I = 0, iq*M-1
Z = Z + R(I,J) * O(I,J)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE TRACE3( M, R, Z )
PARAMETER( iq=5, MXX=100 )
INTEGER M
REAL*8 Z, R(0:iq*MXX,0:*)
Z = 0.0D0
DO 10 J = 0, iq*M-1
DO 10 I = 0, iq*M-1
Z = Z + R(I,J) * R(I,J)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE EXPCTM( M, C, R, VW )
PARAMETER( iq=5, MXX=100 )
INTEGER M, S, T, A
REAL*8 VW(0:*), R(0:iq*MXX,0:*), C(0:iq*MXX,0:*)
C
C t k
C | |
C s---V--n--+ tk
C | | --> s |
C L | i---C
C i---+-----R
CALL MATCLR( iq*M, C )
DO 10 S = 0, iq-1
DO 10 N = 0, iq-1
DO 10 T = 0, iq-1
DO 10 L = 0, iq-1
A = iq*iq*iq*S + iq*iq*N + iq*T + L
DO 10 K = 0, M-1
DO 10 I = 0, M-1
C(S*M+I,T*M+K) = C(S*M+I,T*M+K) + VW(A) * R(L*M+I,N*M+K)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE QUADRA( M, C, A )
PARAMETER( iq=5, MXX=100 )
INTEGER M
REAL*8 C(0:iq*MXX,0:*), A(0:iq*MXX,0:*)
CALL MATCLR( iq*M, A )
DO 10 K = 0, iq*M-1
DO 10 J = 0, iq*M-1
DO 10 I = 0, iq*M-1
A(I,K)= A(I,K) + C(I,J) * C(J,K)
10 CONTINUE
RETURN
END
CCC
CC
C
SUBROUTINE EXPVTO( M, P, O, VW )
PARAMETER( iq=5, MXX=100 )
INTEGER M, S, A
REAL*8 VW(0:*),P(0:iq*MXX,0:iq*MXX,0:*),O(0:iq*MXX,0:iq*MXX,0:*)
C
C i L iL
C | | |
C P---k---V---s ---> O---s
C | | |
C j n jn
DO 10 I = 0, iq-1
CALL MATCLR( iq*M, O(0,0,I) )
10 CONTINUE
DO 20 K = 0, iq-1
DO 20 S = 0, iq-1
DO 20 L = 0, iq-1
DO 20 N = 0, iq-1
A = iq*iq*iq*K + iq*iq*S + iq*L + N
DO 20 J = 0, M-1
DO 20 I = 0, M-1
O(L*M+I,N*M+J,S)= O(L*M+I,N*M+J,S) + P(I,J,K) * VW(A)
20 CONTINUE
RETURN
END
CCC
CC
C
SUBROUTINE SELECT( M, N, W, MX, ER )
PARAMETER( iq=5, MXX=100 )
INTEGER M, N, MX
REAL*8 EPS, ER, R, W(0:*)
C Write(6,*)
C Write(6,100) (W(I),I=iq*M-1,0,-1)
C 100 FORMAT(1H ,' Eigen Values ',E15.7)
EPS = 1.0E-40
R = 1.001
N = iq*M
IF( N .LE. MX .AND. W(0) .GT. EPS ) THEN
ER = W(0)
RETURN
END IF
10 CONTINUE
N = N - 1
IF( N .EQ. 1 ) RETURN
IF( W(iq*M-N-1) .LT. EPS .AND. N. LE. MX .AND.
& W(iq*M-N ) .GT. EPS ) RETURN
IF( W(iq*M-N-1) .LT. EPS ) GOTO 10
IF( N .LE. MX .AND. W(iq*M-N)/W(iq*M-N-1) .GT. R ) THEN
ER = W(iq*M-N)
RETURN
END IF
GOTO 10
END
CCC
CC
C
SUBROUTINE RENORM( M, N, P, O, R )
PARAMETER( iq=5, MXX=100 )
INTEGER M, N
REAL*8 P(0:iq*MXX,0:iq*MXX,0:*), R(0:iq*MXX,0:*)
REAL*8 O(0:iq*MXX,0:iq*MXX,0:*)
C k k
C | |
C i---+---j ---> i---|
C | |R |
C i---O---j--L i---P---L
DO 10 I = 0, iq-1
CALL MATCLR( iq*M, P(0,0,I) )
10 CONTINUE
DO 20 K = 0, iq-1
DO 20 L = 0, N-1
DO 20 J = 0, iq*M-1
DO 20 I = 0, iq*M-1
P(I,L,K)= P(I,L,K) + O(I,J,K) * R(J,L)
20 CONTINUE
C k k
C | |
C i---+ ---> |
C R| | |
C L--i---P---j L---O---j
DO 30 I = 0, iq-1
CALL MATCLR( N, O(0,0,I) )
30 CONTINUE
DO 40 K = 0, iq-1
DO 40 J = 0, N-1
DO 40 L = 0, N-1
DO 40 I = 0, iq*M-1
O(L,J,K)= O(L,J,K) + P(I,J,K) * R(I,L)
40 CONTINUE
DO 50 I = 0, iq-1
CALL MATCPY( N, O(0,0,I), P(0,0,I) )
50 CONTINUE
RETURN
END
CCC
CC
C
SUBROUTINE NORMW( N, W )
INTEGER N
REAL*8 W(0:*)
DO 10 I = 1, N-1
W(I)= W(I) / W(0)
10 CONTINUE
W(0)= 1.0D0
RETURN
END
CCC
SUBROUTINE NORMP( N, P )
PARAMETER( iq=5, MXX=100 )
INTEGER N
REAL*8 FMX, P(0:iq*MXX,0:iq*MXX,0:*)
FMX = 0.0D0
DO 10 K = 0, iq-1
DO 10 J = 0, N-1
DO 10 I = 0, N-1
IF( FMX .LT. ABS( P(I,J,K) ) ) FMX = ABS( P(I,J,K) )
10 CONTINUE
DO 20 K = 0, iq-1
DO 20 J = 0, N-1
DO 20 I = 0, N-1
P(I,J,K)= P(I,J,K) / FMX
20 CONTINUE
RETURN
END
CCC
SUBROUTINE MATCPY( N, A, B )
PARAMETER( iq=5, MXX=100 )
INTEGER N
REAL*8 A(0:iq*MXX,0:*), B(0:iq*MXX,0:*)
DO 10 J = 0, N-1
DO 10 I = 0, N-1
B(I,J)= A(I,J)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE MATCLR( N, A )
PARAMETER( iq=5, MXX=100 )
INTEGER N
REAL*8 A(0:iq*MXX,0:*)
DO 10 J = 0, N-1
DO 10 I = 0, N-1
A(I,J)= 0.0D0
10 CONTINUE
RETURN
END
CCC
SUBROUTINE DIAGAL( N, W, R )
PARAMETER( iq=5, MXX=100 )
INTEGER N, IERR
REAL*8 R(0:iq*MXX,0:*), W(0:*), W1(0:5*iq*MXX)
CALL DCSMSN( R, iq*MXX+1, N, -1.0D0, W, N, -1, W1, IERR )
IF( IERR .NE. 0 ) THEN
WRITE(6,*) 'IERR = ', IERR,' in DIAGAL'
STOP
END IF
RETURN
END
CCC
SUBROUTINE DIAGPA( M2, N, W, C, R )
PARAMETER( iq=5, MXX=100 )
INTEGER M2, N, IERR, IW(0:iq*MXX)
REAL*8 C(0:iq*MXX,0:*), W(0:*), R(0:iq*MXX,0:*),W1(0:iq*8*MXX)
CALL DCSMSS(C,iq*MXX+1,M2,-1.0D0,W,N,R,iq*MXX+1,1,IW,W1,IERR)
IF( IERR .NE. 0 ) THEN
WRITE(6,*) 'IERR = ', IERR,' in DIAGPT'
STOP
END IF
RETURN
END