Tensor Product Variational Appriximation for 3D Ising Model
@
c123456789012345678901234567890123456789012345678901234567890123456789012 c ----- ----- c ----- IRF-variational approximation for 3D Ising model ----- c ----- assisted by CTMRG. ----- c ----- ----- c ----- '99/12/22, Ver.1.1 by T.Nishino ----- c ----- ----- module var_size integer,parameter :: iq = 2, mVL = iq**8, mXX=20 integer,parameter :: mIT =10000, mXT = mVL*mXX**2 end module c program main use var_size integer mX, L, LL, m(0:1), mN(0:1), q(0:1), i, j, N, info real(8) K, fP(0:1), G1(0:mXT,0:1) real(8) x(0:3), fC(0:1), W(0:mVL,0:3), G2(0:mXT,0:1) real(8) y(0:3), eP(0:1), P(0:mXT,0:1), G3(0:mXT,0:1) real(8) u, eC(0:1), C(0:mXT,0:1), G4(0:mXT,0:1) real(8) v, Z(0:1,0:mIT), F(0:mIT,0:1), G5(0:mXT,0:1) real(8) eps, mg(0:1,0:mIT) c call getprm( L, mX, K, eps ) ! --- Parameter Imput --- call initw( q, m, K, W, P, C ) ! --- Initialize TM --- c LL = 2; eP(0:1)= 0.0d0; eC(0:1)= 0.0d0 goto 20 c 10 continue call joint1( 0,0,1, q, m, C, P, G1 ) ! --- Conserve G1 --- call joint1( 1,0,1, q, m, P, G1, G2 ) call joint1( 0,0,2, q, m, C, G2, G3 ) call jointX( q, m, G3, G2 ) ! --- Conserve G2 --- call obsrv1( x,y, q, W, G2 ) c c -------------- Improve IRF weight ----------------------- c G3(0:q(1)**4-1,1) = G2(0:q(1)**4-1,1)*W(0:q(1)**4-1,3) c call atorame( q(0), G3(0,1), G4(0,1) ) N = q(0)**4 G3(0:N-1,1) = 0.0d0 do i = 0, N-1 G3(0:N-1,1)= G3(0:N-1,1) + G4(i*N:i*N+N-1,1)*W(i,2) end do c G3(0:N-1,1) = G3(0:N-1,1) / G2(0:N-1,0) c call norm( N, G3(0,1), u ) c call vratio( q, u, G2, W ) c W(0:N-1,2)= W(0:N-1,2) + eps * G3(0:N-1,1) c call norm( N, W(0,2), v ) call bigmac( q, W ) call vratio( q, v, G2, W ) c c ----------------------------------------- Create DMat --- call extP( q, m, W, P, G2 ) ! --- Conserve G2 --- call extC( q, m, G1, G2, G3 ) ! --- Conserve G3 --- call joint3( q, m, G3, G4 ) ! --- Expire G1 --- call joint3( q, m, G4, G5 ) call obsrv2( q, m, mg(0,LL), G5 ) call credm( q, m, G5, G4 ) ! --- Dens. Mat. --- call crepf( q, m, Z(0,LL), G4 ) c F(LL,0:1)= eP(0:1)*8 + eC(0:1)*4 + log(Z(0:1,LL)) c c ------------------------------------- Display Results --- if(LL .gt. 8) & write(6,100) LL, v/u-1.0d0, & mg(0,LL )/Z(0,LL ), mg(1,LL )/Z(1,LL ), y(0), & ( F(LL,0)+F(LL-4,0)-2*F(LL-2,0) & -F(LL,1)-F(LL-4,1)+2*F(LL-2,1) )/8 100 format('#',I4,E11.4, F15.10, F15.10, F12.7, F15.10) c call diag( q, m, G4, G1, G5 ) c mN(0:1)=min( mX, q(0:1)*m(0:1) ) c call renop( q, m, mN, G2, P, G4 ) ! --- Expire G2 --- call renoc( q, m, mN, G3, C, G4 ) ! --- Expire G3 --- m(0:1) = mN(0:1) c call symchp( q(0), m(0), P(0,0), ' -a0- ' ) ! # call symchp( q(1), m(1), P(0,1), ' -a1- ' ) ! Debug call symchc( q(0), m(0), C(0,0), ' -a0- ' ) ! Routines call symchc( q(1), m(1), C(0,1), ' -a1- ' ) ! # c 20 continue call norm( q(0)**2 * m(0)**2, P(0,0), fP(0) ) call norm( q(1)**2 * m(1)**2, P(0,1), fP(1) ) call norm( q(0) * m(0)**2, C(0,0), fC(0) ) call norm( q(1) * m(1)**2, C(0,1), fC(1) ) c eC(0:1)= eC(0:1) + 2*eP(0:1) + log(fC(0:1)) eP(0:1)= eP(0:1) + log(fP(0:1)) LL = LL + 2 if( LL .LT. L ) goto 10 end program ccccc subroutine errstp( ch ) character(6) ch write(6,*) '*** Error stop in subroutine -',ch,'- ***' stop "Error_stop" end ccccc subroutine getprm( L, mX, K, eps ) use var_size integer L, mX; real(8) K, eps c write(6,*) 'Input Max. Iteration # L'; read(5,*) L if( L .GT. mIT ) call errstp('L >mIT') c write(6,*) 'Input Maximum Dim. mX>= m'; read(5,*) mX if( mX .GT. mXX ) call errstp('mX>mXX') c write(6,*) 'Input improvement fac.eps'; read(5,*) eps c write(6,*) 'Input J/T K '; read(5,*) K write(6,*) 'L = ',L,' mX = ',mX, ' K = ',K write(6,*) 'eps= ', eps end ccccc subroutine initw( q, m, K, W, P, CC ) use var_size integer q(0:*), m(0:*), a, b, c, d real(8) K, f(0:1,0:1), W(0:mVL,0:*), P(0:mXT,0:*), CC(0:mXT,0:*) c integer :: s1(0:3)=(/0,1,0,1/), s2(0:3)=(/0,0,1,1/) c f(0,0)= exp( K/4 ); f(0,1)= exp( -K/4 ) ! --- Bond Factors --- f(1,1)= exp( K/4 ); f(1,0)= exp( -K/4 ) c q(0)= iq do a = 0, q(0) - 1 ! d----c do b = 0, q(0) - 1 ! | W2 | = Initial Face Factor do c = 0, q(0) - 1 ! a----b do d = 0, q(0) - 1 W(((a*q(0)+b)*q(0)+c)*q(0)+d,2)= f(a,b)*f(b,c)*f(c,d)*f(d,a) end do; end do; end do; end do c q(1)= iq**2 do a = 0, q(1) - 1 ! d----c do b = 0, q(1) - 1 ! | W3 | = Cube Factor do c = 0, q(1) - 1 ! a----b do d = 0, q(1) - 1 W((( a *q(1)+ b) *q(1)+ c) *q(1)+ d, 3)= & W(((s1(a)*q(0)+s1(b))*q(0)+s1(c))*q(0)+s1(d),2)* & W(((s2(a)*q(0)+s2(b))*q(0)+s2(c))*q(0)+s2(d),2)* & f( s1(a), s2(a))* f(s1(b), s2(b) )* & f( s1(c), s2(c))* f(s1(d), s2(d) ) end do; end do; end do; end do c call bigmac( q, W ) c do a = 0, 1 m(a)=1 ! --- Fixed Boundary Condition --- do b = 0, q(a)-1 do c = 0, q(a)-1 ! 0------c 0------c P(b*q(a)+c,a) ! | P0/1 | = | W0/1 | & =W(b*q(a)+c,a) ! 0------b 0------b end do; end do ! (use rotated W) do b = 0, q(a)-1 ! +------0 0------0 CC(b,a)= W(b,a) ! | C0/1 | = | W0/1 | end do ! 0------b 0------b end do ! (use rotated W) end ccccc subroutine bigmac( q, W ) use var_size integer i, j, q(0:*); real(8) W(0:mVL,0:*) real(8),allocatable :: tp(:) c ! d----c d----c d----c W(0:q(0)**4-1,0) = ! | W0 | = | W2 | * | W2 | & W(0:q(0)**4-1,2)**2 ! a----b a----b a----b c allocate( tp(0:q(0)**8-1) ); call atorame( q(0), W(0,3), tp ) c do i = 0, q(0)**4-1 ! d1---c1 d12--c12 d2---c2 d----c do j = 0, q(0)**4-1 ! | W2 | * | W3 | * | W2 | = | W1 | tp(i*q(0)**4+j)= W(i,2)! a1---b1 a12--b12 a2---b2 a----b & * tp(i*q(0)**4+j)* W(j,2) end do; end do c call torame( q(0), tp, W(0,1) ); deallocate( tp ) end ccccc subroutine torame( q, W, Y ) integer a,b,c,d,e,f,g,h, q; real(8) W(0:*), Y(0:*) c do a = 0, q-1; do b = 0, q-1; do c = 0, q-1; do d = 0, q-1 do e = 0, q-1; do f = 0, q-1; do g = 0, q-1; do h = 0, q-1 c Y(((((((a*q+e)*q+b)*q+f)*q+c)*q+g)*q+d)*q+h) & = W(((((((a*q+b)*q+c)*q+d)*q+e)*q+f)*q+g)*q+h) end do; end do; end do; end do; end do; end do; end do; end do end ccccc subroutine atorame( q, W, Y ) integer a,b,c,d,e,f,g,h, q; real(8) W(0:*), Y(0:*) c do a = 0, q-1; do b = 0, q-1; do c = 0, q-1; do d = 0, q-1 do e = 0, q-1; do f = 0, q-1; do g = 0, q-1; do h = 0, q-1 c Y(((((((a*q+b)*q+c)*q+d)*q+e)*q+f)*q+g)*q+h) & = W(((((((a*q+e)*q+b)*q+f)*q+c)*q+g)*q+d)*q+h) end do; end do; end do; end do; end do; end do; end do; end do end ccccc subroutine joint1( p1, p2, p3, q, m, A, B, C ) use var_size integer i,j,k,l,n, p1, p2, p3, q(0:*), m(0:*) real(8) A(0:mXT,0:*), B(0:mXT,0:*), C(0:mXT,0:*) c do i = 0, 1 C(0:m(i)**2 *q(i)**(p1+p3+1)-1,i)= 0.0d0 do j = 0, m(i)*q(i)** p1-1 do k = 0, q(i) -1 ! +-----L-----+ do l = 0, m(i)*q(i)** p2-1 ! | 2l | do n = 0, m(i)*q(i)** p3-1 ! | 1 | 3 | C((j*q(i)+k)*q(i)** p3*m(i)+n,i) = ! J-j---k---n-N & C((j*q(i)+k)*q(i)** p3*m(i)+n,i) + & A((j*q(i)+k)*q(i)** p2*m(i)+l,i) * & B((l*q(i)+k)*q(i)** p3*m(i)+n,i) end do; end do; end do; end do end do end ccccc subroutine jointX( q, m, A, B ) use var_size integer i,j,k,l,n, p1, p2, q(0:*), m(0:*) real(8) A(0:mXT,0:*), B(0:mXT,0:*) do i = 0, 1 B(0:q(i)**4-1,i)= 0.0d0 ! +--p1 p1-----+ do j = 0, q(i)-1 ! | | | | do k = 0, q(i)-1 ! | l l--n | do l = 0, q(i)-1 ! | | | | do n = 0, q(i)-1 ! | k--j j | do p1 = 0, m(i)-1 ! | | | | do p2 = 0, m(i)-1 ! +-----p2 p2--+ c B( (( j *q(i)+k)*q(i)+l)*q(i)+ n,i) & = B( (( j *q(i)+k)*q(i)+l)*q(i)+ n,i) & + A((((p2*q(i)+j)*q(i)+k)*q(i)+l)*m(i)+p1,i) & * A((((p1*q(i)+l)*q(i)+n)*q(i)+j)*m(i)+p2,i) c end do; end do; end do; end do; end do; end do end do end ccccc subroutine joint3( q, m, A, B ) use var_size integer q(0:*), m(0:*), i, c, d, e, f real(8) A(0:mXT,0:*), B(0:mXT,0:*) c do i = 0, 1 ! +---f B(0:q(i)**3*m(i)**2-1,i)=0 ! | A | do c = 0, q(i)* m(i)-1 ! e---d do d = 0, q(i) -1 ! | A | do e = 0, q(i)* m(i)-1 ! +---c do f = 0, q(i)* m(i)-1 B((c*q(i)+d)*q(i)*m(i)+f,i) = B((c*q(i)+d)*q(i)*m(i)+f,i) & + A((c*q(i)+d)*q(i)*m(i)+e,i) & * A((e*q(i)+d)*q(i)*m(i)+f,i) end do; end do; end do; end do end do end ccccc subroutine obsrv1( x, y, q, W, G ) use var_size integer :: i, q(0:*), dg(0:4)=(/1,2,4,8,1/) real(8) z, x(0:*), y(0:*), W(0:mVL,0:*), G(0:mXT,0:*) c x(0:3) = 0.0d0; y(0:3) = 0.0d0; z = 0.0d0 c do i = 0, q(0)**4-1 z = z + G(i,0)*W(i,0) x(0:3)= x(0:3) + G(i,0)*W(i,0) * (mod(i/dg(0:3),2)*2-1) y(0:3)= y(0:3) + G(i,0)*W(i,0) * (mod(i/dg(0:3),2)*2-1) & * (mod(i/dg(1:4),2)*2-1) end do x(0:3) = x(0:3) / z y(0:3) = y(0:3) / z end ccccc subroutine obsrv2( q, m, mg, G ) use var_size integer i,j,k, q(0:*), m(0:*) real(8) mg(0:*), G(0:mXT,0:*) c mg(0:1)= 0.0d0 do i = 0, 1 do j = 0, q(i) -1 do k = 0, q(i)*m(i)-1 mg(i)= mg(i) + G((k*q(i)+j)*q(i)*m(i)+k,i) * (mod(j,2)*2-1) end do; end do; end do end ccccc subroutine vratio( q, u, G2, W ) use var_size integer i, j, q(0:*) real(8) a(0:1), u, G2(0:mXT,0:*), W(0:mVL,0:*) c do i = 0, 1 a(i)= 0.0d0 do j = 0, q(i)**4-1 a(i)= a(i) + W(j,i)*G2(j,i) end do end do u = a(1)/a(0) end ccccc subroutine extP( q, m, W, P, G ) use var_size integer q(0:*), m(0:*), a, b, c, d, e, f, i real(8) W(0:mVL,0:*), P(0:mXT,0:*), G(0:mXT,0:*) c do i = 0, 1 do a = 0, m(i)-1 ! f----e---d do b = 0, q(i)-1 ! | P | W | do c = 0, q(i)-1 ! a----b---c do d = 0, q(i)-1 ! do e = 0, q(i)-1 ! careful to the a-b order do f = 0, m(i)-1 G(((((b*m(i)+a)*q(i)+c)*q(i)+d)*q(i)+e)*m(i)+f,i) = & P(((a*q(i)+b)*q(i) +e)*m(i)+f,i) * & W(((b *q(i)+c)*q(i)+d)*q(i)+e, i) end do; end do; end do; end do; end do; end do end do end ccccc subroutine extC( q, m, G1, G2, G3 ) use var_size integer q(0:*), m(0:*), a, b, c, d, e, f, i real(8) G1(0:mXT,0:*), G2(0:mXT,0:*), G3(0:mXT,0:*) c do i = 0, 1 G3(0:q(i)**3*m(i)**2-1,i) = 0 do a = 0, q(i) *m(i)-1 ! +-------f do b = 0, q(i) -1 ! | G1 | do c = 0, q(i) -1 ! e---d---c do d = 0, q(i) -1 ! | G2 | do e = 0, m(i)-1 ! A---a---b do f = 0, m(i)-1 G3( ((a*q(i)+b)*q(i)+c) *m(i)+f,i) = & G3( ((a*q(i)+b)*q(i)+c) *m(i)+f,i) + & G2((((a*q(i)+b)*q(i)+c)*q(i)+d)*m(i)+e,i) * & G1( ((e*q(i)+d)*q(i)+c) *m(i)+f,i) end do; end do; end do; end do; end do; end do end do end ccccc subroutine credm( q, m, A, B ) use var_size integer q(0:*), m(0:*), i, c, d, e real(8) A(0:mXT,0:*), B(0:mXT,0:*) c do i = 0, 1 B(0:q(i)**2*m(i)**2-1,i)=0 do c = 0, q(i)* m(i)-1 do d = 0, q(i) -1 do e = 0, q(i)* m(i)-1 B(c*q(i)* m(i)+e,i)= B( c *q(i)*m(i)+e,i) & + A((c*q(i)+d)*q(i)*m(i)+e,i) end do; end do; end do end do end ccccc subroutine crepf( q, m, Z, A ) use var_size integer q(0:*), m(0:*), i, j real(8) Z(0:*), A(0:mXT,0:*) do i = 0, 1 Z(i)= 0.0d0 do j = 0, q(i)*m(i)-1 Z(i)= A(j*q(i)*m(i)+j,i) + Z(i) end do; end do end ccccc subroutine diag( q, m, A, E, tp ) use var_size integer q(0:*), m(0:*), N, info real(8) A(0:mXT,0:*), E(0:mXT,0:*), tp(0:*) c do i = 0, 1 N = q(i)*m(i) call DSYEV('V','U', N, A(0,i), N, E(0,i), tp, (N+2)*N, info ) c call dcsmaa( A(0,i), N, N, E(0,i), tp, info ) c if( info .ne. 0 ) then write(6,*) 'dsyev_info = ', info call errstp(' diag ') end if c write(6,'(f30.10)') (E(j,i),j=N-1,N-min(N-1,8)-1,-1) end do end cccccc subroutine renop( q, m, mN, G, P, R ) use var_size integer i, a, b, c, q(0:*), m(0:*), N, mN(0:*) real(8) G(0:mXT,0:*), P(0:mXT,0:*), R(0:mXT,0:*) c do i = 0, 1 N = q(i) *m(i) ! c---c c---c P(0:q(i)**2*N*mN(i)-1,i)= 0.0d0 ! | | | | c ! c P | = c G b do a = 0, mN(i)-1 ! | | | |R\ do b = 0, N -1 ! C---a C---B--a do c = 0, q(i)**2*N-1 P(a*q(i)**2*N+c,i) = P(a*q(i)**2*N+c,i) & + G(b*q(i)**2*N+c,i) & * R( (N-a-1)*N+b,i) end do; end do; end do c ! a---a a---a G(0:q(i)**2*mN(i)**2-1,i) = 0.0d0 ! | | | | c ! | G | = c P | do a = 0, q(i)**2*mN(i)-1 ! | | /R| | do b = 0, mN(i)-1 ! b---a b--C---a do c = 0, N -1 G(a*mN(i)+b,i) = G( a* mN(i)+b,i) & + P( a* N +c,i) & * R((N-b-1)*N+c,i) end do; end do; end do c P(0:q(i)**2*mN(i)**2-1,i) = G(0:q(i)**2*mN(i)**2-1,i) end do end ccccc subroutine renoc( q, m, mN, G, C, R ) use var_size integer i, a,b, d, q(0:*), m(0:*), N, mN(0:*) real(8) G(0:mXT,0:*), C(0:mXT,0:*), R(0:mXT,0:*) c do i = 0, 1 ! d--d--d d--d--d N = q(i) *m(i) ! | | | | C(0:q(i)*N*mN(i)-1,i) = 0.0d0 ! | C | = | G b do a = 0, mN(i)-1 ! | | | |R\ do b = 0, N -1 ! +-----a +-----b--a do d = 0, q(i)*N -1 C(a*q(i)*N+d,i)= C( a*q(i)*N+d,i) & + G( b*q(i)*N+d,i) & * R((N-a-1)*N+b,i) end do; end do; end do c ! b G(0:q(i)*mN(i)**2-1,i) = 0.0d0 ! |R\ do a = 0, q(i)*mN(i)-1 ! b-----a d--d--a do b = 0, mN(i)-1 ! | G | = | C | do d = 0, N -1 ! +-----a +-----a G(a*mN(i)+b,i)= G( a* mN(i)+b,i) & + C( a* N +d,i) & * R((N-b-1)*N +d,i) end do; end do; end do c C(0:q(i)*mN(i)**2-1,i) = G(0:q(i)*mN(i)**2-1,i) end do end ccccc subroutine norm( 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 V(0:N-1)= V(0:N-1) / W end ccccc subroutine symchp( q, m, P, ch ) integer i,j,k,l, q, m; real(8) P(0:*), S; character(6) ch S = 0.0d0 do i = 0, m-1 do j = 0, q-1 do k = 0, q-1 do l = 0, m-1 S = S + ABS( P(((i*q+j)*q+k)*m+l)-P(((l*q+k)*q+j)*m+i) ) end do; end do; end do; end do if( S .gt. 1.0E-8 ) write(6,*) 'sym P = ', S, ch end ccccc subroutine symchc( q, m, C, ch ) integer i,j,k, q, m; real(8) C(0:*), S; character(6) ch S = 0.0d0 do i = 0, m-1 do j = 0, q-1 do k = 0, m-1 S = S + ABS( C((i*q+j)*m+k)-C((k*q+j)*m+i) ) end do; end do; end do; if( S .gt. 1.0E-8 ) write(6,*) 'sym C = ', S, ch end