anchor.diff:  
   
   program skyrmionlattice
  -    use dou_ex 
  +    use  modu
       implicit none
       !#########################################
       !L, size of lattice
       !N = L*L, num of lattice
       !t, J_k, the parameter if H
       !##########################################
  -    integer, parameter :: L = 40, N = L*L, num_MC = 10**6, num_T_para = 15
  +    integer, parameter :: L = 40, N = L*L, num_MC = 4*10**5, num_T_para = 20
       real*8, parameter :: t = 1.0, PI = 3.141592654
  -    integer :: i, j, k, ii, p, d, di, dd, ij, aver_num
  -    real*8 :: ran, J_k, beta, miu, Tc, x, energyV, energyD
  -    real*8 :: V, deltaV, J_H, lambda, B_z, pa, V_new, deltaV_new
  +    integer :: i, j, k, ii, r, d, di, dd, ij, aver_num
  +    real*8 :: ran, beta, Tc, x, Hp, B_z
  +    real*8 :: V, deltaV, J_H, lambda, pa, V_new
  +    integer, dimension(1) :: mloc
       real*8, dimension(2,N) :: fai, theta
       real*8, dimension(num_T_para) :: T_para
  -    open( unit = 12, file = 'B_z.in' )
  +    real*8, dimension(N) :: H
  +    character(len=2) :: cha
  +    d = 1
  +    lambda = 0.4*PI*t
  +    J_H = t
  +    aver_num = 100
  +    open(unit = 12, file = 'B_z.in')
       read(12,*) pa
  -    !pa = 55.0/89
  -    d = 1
  -    lambda = 0.3*PI
  -    J_H = lambda/11.7810
       B_z = pa*t
  -    deltaV = 0
  -    V = 0
  -    deltaV_new = 0
  -    V_new = 0
  -    energyV = 0
  -    energyD = 0
  -    aver_num = 100
  -
  -    T_para = (/1.0, 0.5, 0.1, 0.08, 0.05, 0.002, 0.01, 0.008, 0.005 ,0.001, 0.0009, 0.008, 0.0007, 0.0006, 0.0005/)
  +    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)
  -    open( unit = 13, file = 'lattice.out' )
   
       do i = 1, N
           x = ran()
  
positive.code: 
      end subroutine calcu_deltaV
      
  end module dou_ex
  
  real(8) function ran()
  !----------------------------------------------!
  
negative.code:         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()
          fai(dd, r) = fai(d,r)+(x-0.5)*2*2*PI/12
          if ( fai(dd,r) < 0 ) then
              fai(dd,r) = fai(dd,r)+2*PI
          end if
          if ( fai(dd,r) > 2*PI ) then
              fai(dd,r) = fai(dd,r)-2*PI
          end if
  
          x = ran()
          theta(dd, r) = theta(d,r)+(x-0.5)*2*PI/6
          if ( theta(dd,r) < 0 ) then
              theta(dd,r) = theta(dd,r)+PI
          end if
          if (theta(dd,r) > PI ) then
              theta(dd,r) = theta(dd,r)-PI
          end if
  
          call calcu_Hp( r, L, J_H, lambda, fai(dd,:), theta(dd,:), Hp)
          deltaV = B_z*( cos(theta(dd,r))-cos(theta(d,r)) )
  
          x = ran()
          !write(*,*) x, '<?', exp( -beta*(Hp-H(r)+deltav) )
          if ( x > exp( -beta*(Hp-H(r)+deltaV) ) ) then
              fai(dd,r) = fai(d,r)
              theta(dd, r) = theta(d, r)
          else
              V = V+Hp-H(r)
              H(r) = Hp
          end if
          d = dd
          !if (ii == num_T_para .AND. i > num_MC-aver_num*N ) then
          if ( i > num_MC-aver_num*N ) then
              if ( mod(i+aver_num*N-num_MC, N) == 0 ) then
                  do j = 1, N
                      write(14,*) theta(d,j), fai(d,j)
                  end do
              end if
          end if
          !if ( r == N ) then
          !    r = 0
          !end if
          write(13,*) V
      end do
      close(13)
      close(14)
  end do
  
      write(*,*) 'end'
  end program skyrmionlattice
  