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: !
  !###############################################################
  
  module dou_ex
  contains
  
      subroutine calcu_V( L, N, J_H, lambda, fai, theta, V )
          implicit none
          integer :: L, N, i, j, i_minus, i_plus, j_minus, j_plus, N_i_minus, N_i_plus, N_j_minus, N_j_plus, N_ij, N_i_plus_j_minus, N_i_minus_j_plus
          real*8 :: J_H, lambda, fai(N), theta(N), V
          V = 0
          do i = 0, L-1, 1
              do j = 0, L-1, 1
                  !period of boundary 
                  i_plus = mod(i+1, L)
                  j_plus = mod(j+1, L)
                  if (i==0) then
                      i_minus = L-1
                  else
                      i_minus = i-1
                  end if 
                  if (j==0) then
                      j_minus = L-1
                  else
                      j_minus = j-1
                  end if 
                  N_ij = i*L+j+1
                  N_i_minus = i_minus*L+j+1
                  N_i_plus = i_plus*L+j+1
                  N_j_minus = i*L+j_minus+1
                  N_j_plus = i*L+j_plus+1
                  N_i_minus_j_plus = i_minus*L+j_plus+1
                  N_i_plus_j_minus = i_plus*L+j_minus+1
  
                  V = V-1.0/2*J_H*sin(theta(N_ij))*cos(fai(N_ij))*sin(theta(N_i_minus))*cos(fai(N_i_minus))
                  V = V-1.0/2*J_H*sin(theta(N_ij))*cos(fai(N_ij))*sin(theta(N_i_plus))*cos(fai(N_i_plus))
                  V = V-1.0/2*J_H*sin(theta(N_ij))*cos(fai(N_ij))*sin(theta(N_j_minus))*cos(fai(N_j_minus))
                  V = V-1.0/2*J_H*sin(theta(N_ij))*cos(fai(N_ij))*sin(theta(N_j_plus))*cos(fai(N_j_plus))
                  V = V-1.0/2*J_H*sin(theta(N_ij))*cos(fai(N_ij))*sin(theta(N_i_plus_j_minus))*cos(fai(N_i_plus_j_minus))
                  V = V-1.0/2*J_H*sin(theta(N_ij))*cos(fai(N_ij))*sin(theta(N_i_minus_j_plus))*cos(fai(N_i_minus_j_plus))
  
                  V = V-1.0/2*J_H*sin(theta(N_ij))*sin(fai(N_ij))*sin(theta(N_i_minus))*sin(fai(N_i_minus))
                  V = V-1.0/2*J_H*sin(theta(N_ij))*sin(fai(N_ij))*sin(theta(N_i_plus))*sin(fai(N_i_plus))
                  V = V-1.0/2*J_H*sin(theta(N_ij))*sin(fai(N_ij))*sin(theta(N_j_minus))*sin(fai(N_j_minus))
                  V = V-1.0/2*J_H*sin(theta(N_ij))*sin(fai(N_ij))*sin(theta(N_j_plus))*sin(fai(N_j_plus))
                  V = V-1.0/2*J_H*sin(theta(N_ij))*sin(fai(N_ij))*sin(theta(N_i_plus_j_minus))*sin(fai(N_i_plus_j_minus))
                  V = V-1.0/2*J_H*sin(theta(N_ij))*sin(fai(N_ij))*sin(theta(N_i_minus_j_plus))*sin(fai(N_i_minus_j_plus))
  
                  V = V-1.0/2*J_H*cos(theta(N_ij))*cos(theta(N_i_minus))
                  V = V-1.0/2*J_H*cos(theta(N_ij))*cos(theta(N_i_plus))
                  V = V-1.0/2*J_H*cos(theta(N_ij))*cos(theta(N_j_minus))
                  V = V-1.0/2*J_H*cos(theta(N_ij))*cos(theta(N_j_plus))
                  V = V-1.0/2*J_H*cos(theta(N_ij))*cos(theta(N_i_plus_j_minus))
                  V = V-1.0/2*J_H*cos(theta(N_ij))*cos(theta(N_i_minus_j_plus))
  
                  V = V+1.0/2*lambda*(-1.0)*( sin(theta(N_i_minus))*sin(fai(N_i_minus))*cos(theta(N_ij))-cos(theta(N_i_minus))&
                      *sin(theta(N_ij))*sin(fai(N_ij)) )
                  V = V+1.0/2*lambda*( sin(theta(N_i_plus))*sin(fai(N_i_plus))*cos(theta(N_ij))-cos(theta(N_i_plus))&
                      *sin(theta(N_ij))*sin(fai(N_ij)) )
                  V = V+1.0/2*lambda*(-1.0)*( (1.0/2)*( sin(theta(N_j_minus))*sin(fai(N_j_minus))*cos(theta(N_ij))&
                      -cos(theta(N_j_minus))*sin(theta(N_ij))*sin(fai(N_ij)) ) + sqrt(3.0)/2*( cos(theta(N_j_minus))&
                      *sin(theta(N_ij))*cos(fai(N_ij))-sin(theta(N_j_minus))*cos(fai(N_j_minus))*cos(theta(N_ij)) ) )
                  V = V+1.0/2*lambda*( (1.0/2)*( sin(theta(N_j_plus))*sin(fai(N_j_plus))*cos(theta(N_ij))&
                      -cos(theta(N_j_plus))*sin(theta(N_ij))*sin(fai(N_ij)) ) &
                      + sqrt(3.0)/2*( cos(theta(N_j_plus))*sin(theta(N_ij))*cos(fai(N_ij))-sin(theta(N_j_plus))*cos(fai(N_j_plus))*cos(theta(N_ij)) ) )
                  V = V+1.0/2*lambda*( sin(theta(N_i_plus_j_minus))*sin(fai(N_i_plus_j_minus))*cos(theta(N_ij))&
                      -cos(theta(N_i_plus_j_minus))*sin(theta(N_ij))*sin(fai(N_ij)) ) &
                      +1.0/2*lambda*(-1.0)*( (1.0/2)*( sin(theta(N_i_plus_j_minus))*sin(fai(N_i_plus_j_minus))*cos(theta(N_ij))-cos(theta(N_i_plus_j_minus))*sin(theta(N_ij))*sin(fai(N_ij)) ) &
                      + sqrt(3.0)/2*( cos(theta(N_i_plus_j_minus))*sin(theta(N_ij))*cos(fai(N_ij))-sin(theta(N_i_plus_j_minus))*cos(fai(N_i_plus_j_minus))*cos(theta(N_ij)) ) )
                  V = V+1.0/2*lambda*(-1.0)*( sin(theta(N_i_minus_j_plus))*sin(fai(N_i_minus_j_plus))*cos(theta(N_ij))&
                      -cos(theta(N_i_minus_j_plus))*sin(theta(N_ij))*sin(fai(N_ij)) ) &
                      +1.0/2*lambda*( (1.0/2)*( sin(theta(N_i_minus_j_plus))*sin(fai(N_i_minus_j_plus))*cos(theta(N_ij))-cos(theta(N_i_minus_j_plus))&
                      *sin(theta(N_ij))*sin(fai(N_ij)) ) + sqrt(3.0)/2*( cos(theta(N_i_minus_j_plus))*sin(theta&
                      (N_ij))*cos(fai(N_ij))-sin(theta(N_i_minus_j_plus))*cos(fai(N_i_minus_j_plus))*cos(theta(N_ij)) ) )
  
              end do
          end do
  
  
      end subroutine calcu_V
  
  
negative.code: 
      end subroutine calcu_deltaV
      
  end module modu
  
  real(8) function ran()
  !----------------------------------------------!
  