(← 2004 年 5 月 6 日へのリンク)
え〜と、この私の Web Page も「更新されるのは日記だけ」となって久しいので、何か
他の所を Edit するよりも「新しい内容を盛り込むなら、まずココに掲示してから、後でリ
ンク集みたいなページを作る」ようになりつつある。パソコンのデータベース機能を信頼
して。
さて、という訳で、ウチの学生さんたちの為に、イジング模型用の CTMRG のソースを
ば〜んと掲示。このプログラムは、45 度回転した正方格子 Ising Model を 4 等分するよ
うな、vertex-type の CTMRG Program で、入力するのは反復回数、m、温度の 3 つ。
出力は、磁化とフリーエネルギー。他愛もないプログラムじゃ....中身の説明は、明日以降。
c123456789012345678901234567890123456789012345678901234567890123456789012
c ----- -----
c ----- Vertex CTMRG for Ising model: Chess Board -----
c ----- version. -----
c ----- -----
c ----- '2002/08/24, Ver.1.0 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 = 80 ! 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) 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) C(0:mXT) ! Corner Transfer Matrix
real(8) P(0:mXT) ! Half Row Transfer 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 ) ! Set Physical Parameters
call initbw( K, W ) ! 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)= 0.0d0
goto 20
c
10 continue
call rjoint( m, m, iq*m, C, P, G1 )
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 observ( iq, m, mg, Z, G1 )
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 ); use var_size
integer L, mX; real(8) K, T
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 T (J=1) '; read(5,*) T; K = 1.0d0 / T
write(6,*)' mX = ', mX, ' T = ',T
end
ccccc
subroutine initbw( fK, W ); use var_size
integer i,j,k,l; real(8) fK, W(0:*) ! i---j
do i = 0, iq-1; do j = 0, iq-1 ! | W |
do k = 0, iq-1; do l = 0, iq-1 ! l---k
W(((i*iq+j)*iq+k)*iq+l)= exp( 4*fK*(i-0.5d0)*(j-0.5d0)
& + 4*fK*(j-0.5d0)*(k-0.5d0)
& + 4*fK*(k-0.5d0)*(l-0.5d0)
& + 4*fK*(l-0.5d0)*(i-0.5d0) )
end do; end do; end do; end do; 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 observ( q, m, mg, Z, G )
integer i,j, q, m; real(8) mg, Z, G(0:*)
c
mg = 0.0d0 ; Z = 0.0d0
do i = 0, q-1
do j = 0, m-1
Z = Z + G(((i*m+j)*q+i)*m+j)
mg = mg + G(((i*m+j)*q+i)*m+j)*(2*i-1)
end do; end do; mg = mg / Z
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 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 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
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