anchor.diff:      !N = L*L, num of lattice
       !t, J_k, the parameter if H
       !##########################################
  -    integer, parameter :: L = 40, N = L*L, num_MC = 4*10**5, num_T_para = 20
  +    integer, parameter :: L = 40, N = L*L, num_MC = 4*10**5, num_T_para = 30
       real*8, parameter :: t = 1.0, PI = 3.141592654
       integer :: i, j, k, ii, r, d, di, dd, ij, aver_num
       integer :: p, p_i_pl, p_i_mi, p_j_pl, p_j_mi, p_i_pl_j_mi, p_i_mi_j_pl
       real*8 :: ran, beta, Tc, x, Hp, B_z
  -    real*8 :: V, deltaV, J_H, lambda, pa, V_new
  +    real*8 :: V, deltaV, J_H, lambda, pa, energy_aver, energy_aver_B
       integer, dimension(1) :: mloc
       real*8, dimension(2,N) :: fai, theta
       real*8, dimension(num_T_para) :: T_para
  
positive.code:         x = ran()
          theta(d,i) = x*PI
      end do
      !ij = 0
  do ii = 1, num_T_para
      write(cha,'(i2)') ii
      open(unit = 13, file = 'energy'//adjustl(trim(cha))//'.out' )
      open( unit = 14, file = 'lattice'//adjustl(trim(cha))//'.out' )
      write(*,*) num_T_para-ii
      Tc = T_para(ii)*J_H
      beta = 1/Tc
      !r= 0
      call calcu_V( L, J_H, lambda, fai(d,:), theta(d,:), H, V )
  
      do i = 1, num_MC 
  
          dd = di(d)
          x = ran()
          r = int(x*N)+1 
          !r = r+1
          !Write(*,*) r
          fai(dd,:) = fai(d,:)
          theta(dd,:) = theta(d,:)
          x = ran()
  
negative.code:             call calcu_Hp( p_i_mi_j_pl, L, J_H, lambda, fai(dd,:), theta(dd,:), H(p_i_mi_j_pl))
          end if
          d = dd
          if ( i > num_MC-aver_num*N ) then
              if ( mod(i+aver_num*N-num_MC, N) == 0 ) then
                  energy_aver = energy_aver+V
                  call calcu_deltaV(L, N, B_z, theta(d,:), deltaV)
                  energy_aver_B = energy_aver_B+V+deltav
                  do j = 1, N
                      write(14,*) theta(d,j), fai(d,j)
                  end do
              end if
          end if
          write(13,*) V
      end do
      close(13)
      close(14)
      write(10,*) T_para(ii), energy_aver/aver_num, energy_aver_B/aver_num
  end do
  
      write(*,*) 'end'
  