CTMRG for 2D Ising Model
c123456789012345678901234567890123456789012345678901234567890123456789012
c ----- -----
c ----- Vertex CTMRG for Ising model: Tutorial for -----
c ----- MC Students. -----
c ----- -----
c ----- '2005/04/29, Ver.0.1 by T.Nishino -----
c ----- -----
module var_size
integer,parameter :: mIT = 10000 ! Maximum Iteration Number
integer,parameter :: iq = 2 ! Ising (2-state) spins
integer,parameter :: mXX = 40 ! Maximum value of 'm'
integer,parameter :: mBW = iq**4 ! Size of the Boltzmann Weight
integer,parameter :: mXT = iq**3
& *mXX**2 ! Variable Size
end module
c
program main; use var_size
integer L, LL ! System Size L, Max. System Size LL
integer m, mX ! Block Spin State m, Max. Block Spin State mX
integer mN ! The value of m in the NeXT iteration.
integer tp ! Temporary Variable
c
real(8) K ! Inverse Temperature 1/T
real(8) eps ! Decouple element (***)
real(8) mg ! Spontaneous Magnetization
real(8) u ! Bond Energy
real(8) Z ! Normalized Partition Function
real(8) eP, fP ! Normalization Const. for HRTM
real(8) eC, fC ! Normalization Const. for CTM
real(8) F(0:mIT) ! Free Energy
c
real(8) W(0:mBW) ! Boltzmann Weight
real(8) S(0:mBW) ! Ising Spin Weight
real(8) Y(0:mBW) ! Decouple Weight (***)
real(8) C(0:mXT) ! Corner Transfer Matrix
real(8) D(0:mXT) ! Corner Decouple Matrix (***)
real(8) P(0:mXT) ! Half Row Transfer Matrix
real(8) H(0:mXT) ! Half Row Decouple Matrix (***)
real(8) G1(0:mXT) ! Work Area
real(8) G2(0:mXT) ! Work Area
real(8) G3(0:mXT) ! Work Area
real(8) G4(0:mXT) ! Work Area
c
call getprm( L, mX, K, eps ) ! Set Physical Parameters
call initbw( K, W, S ) ! Initialize Boltzmann Weights
c ! Set CTM, HRTM, etc.
LL = 2; eC = 0.0d0; C(0)= 1.0d0
m = 1; eP = 0.0d0; P(0)= 1.0d0; P(1)= 1.001d0
goto 20
c
10 continue
call rjoint( m, m, iq *m, C, P, G1 )
call rjoint( m*iq, m, iq *m, P, G1, G2 )
call rjoint( m, m, iq**2*m, C, G2, G3 )
call jointX( iq, m, G3, G2, C )
call symchW( iq, G2, C )
c
call observ( mg, Z, W, S, G2 )
c
call extP( iq, m, W, P, G2 )
call extC( iq, m, G1, G2, G3 )
call sjoint( iq *m, iq*m, G3, C )
call sjoint( iq *m, iq*m, C, G1 )
call trace( iq *m, G1, Z )
c
call diag( iq*m, G1, C, P ) ! C/P_temp
mN = min( iq*m, mX )
tp = iq*m*(iq*m-mN )
call renoP( iq,iq*m, mN, G2, P, G1(tp) )
call renoC( iq*m, mN, G3, C, G1(tp) ); m = mN
F(LL)= eP*8 + eC*4 + log(Z)
c
if(LL.gt.8) then; write(6,100) LL,mg,(F(LL)+F(LL-4)-2*F(LL-2))/8
else; write(6,100) LL,mg,0.0d0
end if
100 format('#',I4, F16.10, F16.10 )
c
20 continue
call symchP( iq, m, P, G1 ) ! Check the Symmetry of HRTM
call symchC( m, C, G1 ) ! Check the Symmetry of CTM
call norm( iq *m**2, P, fP ) ! Normalize HRTM
call norm( m**2, C, fC ) ! Normalize CTM
eC= eC + 2*eP + log(fC) ! Normalization Const for CTM
eP= eP + log(fP) ! Normalization Const for HRTM
LL = LL + 2; if( LL .LT. L ) goto 10
end program main
ccccc
subroutine errstp( ch ) ; character(8) ch
write(6,*) 'errstp <',ch,'>'; stop "Input_Error"; end
ccccc
subroutine getprm( L, mX, K eps); use var_size
integer L, mX; real(8) K, eps
write(6,*)'Max. Iteration # L'; read(5,*) L ; if( L.GT.mIT )
& call errstp( 'L > mIT')
write(6,*)'Maximum Dim. mX>= m'; read(5,*) mX; if( mX.GT.mXX )
& call errstp( 'mX > mXX')
write(6,*)'Parameter K = J/T '; read(5,*) K
write(6,*)'Parameter eps '; read(5,*) eps
write(6,*)' mX = ', mX,' K = ',K,' eps = ',eps
end
ccccc
subroutine initbw( K, W, S); use var_size
integer i, ibt; real(8) Xi, K, W(0:*), S(0:*)
c
Xi = exp(2*K) + sqrt(exp(2*K)**2-1)
do i = 0, iq**4-1
W(i)= Xi**(ibt(i)-2)/2**3 + Xi**(-ibt(i)+2)/2**3
S(i)= Xi**(ibt(i)-2)/2**3 - Xi**(-ibt(i)+2)/2**3
end do
end
ccccc
integer function ibt( D ); integer D
ibt = mod(D,2) + mod(D/2,2) + mod(D/4,2) + mod(D/8,2); end
cccccc
subroutine exchg2( m1, m2, A, B )
integer i, j, m1, m2; real(8) A(0:*), B(0:*)
do i = 0, m1-1
do j = 0, m2-1; B(j*m1+i) = A(i*m2+j)
end do; end do; end
cccccc
subroutine exchg3( m1, m2, m3, A, B )
integer i, j, m1, m2, m3; real(8) A(0:*), B(0:*)
do i = 0, m1-1
do j = 0, m2-1; B((j*m1+i)*m3:(j*m1+i)*m3 + m3-1)
& = A((i*m2+j)*m3:(i*m2+j)*m3 + m3-1)
end do; end do; end
cccccc
subroutine trace( m, A, Tr )
integer i, m; real(8) A(0:*), Tr; Tr = 0.0d0
do i = 0, m-1; Tr = Tr + A(i*m+i)
end do; end
cccccc
subroutine utrace( m1, m2, A, B )
integer i, m1, m2; real(8) A(0:*), B(0:*)
B(0:m2-1)= 0.0d0
do i = 0, m1-1; B(0:m2-1)= B(0:m2-1) + A(i*m2:i*m2 + m2-1)
end do; end
cccccc
subroutine sjoint( m1, m2, A, B )
integer i,j, m1, m2; real(8) A(0:*), B(0:*)
B(0:m2*m2-1)= 0.0d0
do i = 0, m1-1
do j = 0, m2-1; B(j*m2:j*m2 + m2-1) = ! *m1* m2
& B(j*m2:j*m2 + m2-1) + A(i*m2+j) ! *m1* m2
& * A(i*m2:i*m2 + m2-1) ! * * m2 m2
end do; end do; end
cccccc
subroutine ujoint( m1, m2, m3, A, B, C )
integer i,j, m1, m2, m3; real(8) A(0:*), B(0:*), C(0:*)
C(0:m2*m3-1)= 0.0d0
do i = 0, m1-1
do j = 0, m2-1; C(j*m3:j*m3 + m3-1) = ! *m1* m2
& C(j*m3:j*m3 + m3-1) + A(i*m2+j) ! *m1* m3
& * B(i*m3:i*m3 + m3-1) ! * * m2 m3
end do; end do; end
cccccc - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
subroutine rjoint( m1, m2, m3, A, B, C )
integer m1, m2, m3
real(8) A(m2,m1), B(m3,m2), C(m3,m1)
C = MATMUL(B,A)
end
cccccc - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
subroutine jointX( q, m, E, F, tp )
integer q, m; real(8) E(0:*), F(0:*), tp(0:*)
c ! iabj +-i-+---+
call exchg2( m *q*q, m, E,tp ) ! jiab | d |
call exchg3( m,m, q*q, tp, E ) ! jicd + a c +
call ujoint( m*m, q*q, q*q,tp, E, F ) ! abcd | b |
end ! +---+-j-+
cccccc
subroutine inpro( N, f, A, B )
integer i, N; real(8) f, A(0:*), B(0:*); f=0.0d0
do i = 0, N-1; f = f + A( i ) *B( i )
end do; end
cccccc
subroutine observ( mg, Z, W, S, G ); use var_size
real(8) mg, Z, W(0:*), S(0:*), G(0:*)
call inpro( iq**4, Z, W, G )
call inpro( iq**4, mg, S, G ); mg = mg / Z
end
cccccc
subroutine extP( q, m, W, P, G )
integer q, m; real(8) W(0,*), P(0:*), G(0:*)
c
call exchg2( m, q*m, P, G ) ! iaj -> [a]ji
call ujoint( q, m*m, q**3, G, W, P ) ! [a] bcd
call exchg2( m,m *q**3, P, G ) ! ibcdj <---- jibcd
call exchg3( m, q,m *q**2, G, P ) ! bicdj
c ! i b
G(0:m**2*q**3-1) = P(0:m**2*q**3-1) ! +-a-+-c
end ! j d
cccccc
subroutine extC( q, m, G1, G2, G3 )
integer q, m; real(8) G1(0:*), G2(0:*), G3(0:*)
c
call exchg3( m, q, m, G1, G3 ) ! iaj -> [ai]j
call ujoint( q*m, m, q**2*m, G3, G2, G1 ) ! [ai] bck
call exchg3( m, q, q *m, G1, G3 ) ! bjck <--- jbck
end
ccccc
subroutine extD( q, m, Y, H, D, G )
integer i,j,k,l, q, m; real(8) Y(0:*), H(0:*), D(0:*), G(0:*)
c
call rjoint( m, m, m, H, D, G )
call rjoint( m, m, m, G, H, D )
c
do i = 0, q-1; do j = 0, m-1
do k = 0, q-1; do l = 0, m-1
G(((i*m+j)*q + k)*m+l) = Y(i*q+k)*D(j*m+l)
end do; end do; end do; end do
D(0:q*q*m*m-1) = G(0:q*q*m*m-1)
end
ccccc
subroutine extH( q, m, Y, H, G )
integer i,j,k,l, q, m; real(8) Y(0:*), H(0:*), G(0:*)
c
do i = 0, q-1; do j = 0, m-1
do k = 0, q-1; do l = 0, m-1
G(((i*m+j)*q + k)*m+l) = Y(i*q+k)*H(j*m+l)
end do; end do; end do; end do
H(0:q*q*m*m-1) = G(0:q*q*m*m-1)
end
ccccc
subroutine diag( N, A, E, tp )
integer N, info; real(8) A(0:*), E(0:*), tp(0:*)
c
call DSYEV('V','U', N, A, N, E, tp, (N+2)*N, info )
c-SX- call dcsmaa( A, N, N, E, tp, info )
if( info .ne. 0 ) then
write(6,*) 'dsyev_info = ', info; call errstp(' -diag- ')
end if
end
cccccc
subroutine renoP( q, qm, mN, G, P, R )
integer q, qm, mN; real(8) G(0:*), P(0:*), R(0:*)
call rjoint( mN, qm, q*qm, R, G, P ) ! x[ai]
call exchg3( mN, q, qm, P, G ) ! [ai]bcj -> xbcj
call exchg2( q*mN, qm, G, P ) ! bxcj
call rjoint( mN, qm, q*mN, R, P, G ) ! y[cj] <- cjbx
c ! [cj]bx
P(0:q*mN**2-1) = G(0:q*mN**2-1) ! y bx
end
ccccc
subroutine renoC( qm, mN, G, C, R )
integer qm, mN; real(8) G(0:*), C(0:*), R(0:*)
call rjoint( mN, qm, qm, R, G, C ) ! x[ai] y[bj]
call exchg2( mN, qm, C, G ) ! [ai]bj -> xbj -> [bj]x
call rjoint( mN, qm, mN, R, G, C ) ! y x
end
cccccc
subroutine norm( N, V, w )
integer N; real(8) V(0:*), w
call maxel( N, V, w ); V(0:N-1)= V(0:N-1)/w
end
ccccc
subroutine maxel( N, V, w )
integer i, N; real(8) V(0:*), w; w = 0.0d0
do i = 0, N-1; if( abs(V(i)) .GT. w ) w = abs(V(i))
end do; end
ccccc
ccccc --------------------------------------
ccccc ---------- Debug Routines ------------
ccccc --------------------------------------
ccccc
subroutine symchP( q, m, P, tp )
integer i,j,k, q, m; real(8) P(0:*), tp(0:*), T, S; S = 0.0d0
do i = 0, m-1
do j = 0, q-1
do k = 0, m-1 ; tp((i*q+j)*m+k)
& = ( P((i*q+j)*m+k) + P((k*q+j)*m+i) )/2
if( S .lt. ABS( P((i*q+j)*m+k) - P((k*q+j)*m+i) ) )
& S = ABS( P((i*q+j)*m+k) - P((k*q+j)*m+i) )
end do; end do; end do
c
P(0:m*q*m-1) = tp(0:m*q*m-1); call maxel( m*q*m, P, T)
c
if( S/T .gt. 1.0E-12 ) write(6,*) 'sym P = ', S/T, q, m
end
cccccc
subroutine symchC( m, C, tp )
integer i,j, m; real(8) C(0:*), tp(0:*), T, S; S = 0.0d0
do i = 0, m-1
do j = 0, m-1 ; tp(i*m+j) =( C(i*m+j) + C(j*m+i) )/2
if( S .lt. ABS( C(i*m+j) - C(j*m+i) ) )
& S = ABS( C(i*m+j) - C(j*m+i) )
end do; end do
c
call maxel( m*m, C, T); C(0:m*m-1)= tp(0:m*m-1)
c
if( S/T .gt. 1.0E-12 ) write(6,*) 'sym C = ', S/T, m
end
cccccc
subroutine symchW( q, W, tp )
integer i,j,k,l, q; real(8) W(0:*), tp(0:*), T, S; S = 0.0d0
do i = 0, q-1
do j = 0, q-1
do k = 0, q-1
do l = 0, q-1 ; tp(((i*q+j)*q+k)*q+l) =
& ( W(((i*q+j)*q+k)*q+l)+ W(((j*q+k)*q+l)*q+i)
& + W(((k*q+l)*q+i)*q+j)+ W(((l*q+i)*q+j)*q+k)
& + W(((i*q+l)*q+k)*q+j)+ W(((l*q+k)*q+j)*q+i)
& + W(((k*q+j)*q+i)*q+l)+ W(((j*q+i)*q+l)*q+k) )/8
if( S .lt. ABS( W(((i*q+j)*q+k)*q+l)- W(((j*q+k)*q+l)*q+i) ) )
& S = ABS( W(((i*q+j)*q+k)*q+l)- W(((j*q+k)*q+l)*q+i) )
end do; end do; end do; end do
c
call maxel( q**4, W, T ); W(0:q**4-1) = tp(0:q**4-1)
if( S/T .gt. 1.0E-12 )
& write(6,*) ' sym W = ', S/T, q
end