anchor.diff:          x = ran()
           theta(d,i) = x*PI
       end do
  -    !ij = 0
  +    call calcu_V( L, J_H, lambda, fai(d,:), theta(d,:), H, V )
   do ii = 1, num_T_para
  +    energy_aver = 0.0
  +    energy_aver_B = 0.0
       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()
  
positive.code:     read(12,*) pa
      close(12)
      B_z = pa*t
      T_para = (/ 1.0, 0.95, 0.9, 0.85, 0.8, 0.75, 0.7, 0.65, 0.6, 0.55, 0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.15, 0.1, 0.05/)
      open(unit = 10, file = 'T.out')
      do i = 1, num_T_para
          write(10,*) T_para(i)
      end do
      close(10)
  
      call initran(1)
  
      do i = 1, N
  
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'
  