anchor.diff:  !
   !###############################################################
   
  -module dou_ex
  +module modu
   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)) ) )
  -
  +    subroutine calcu_p_nearest( p, p_i_pl, p_i_mi, p_j_pl, p_j_mi, p_i_pl_j_mi, p_i_mi_j_pl, L )
  +        implicit none
  +        integer :: p, p_i_pl, p_i_mi, p_j_pl, p_j_mi, p_i_pl_j_mi, p_i_mi_j_pl, L
  +        integer :: i, j, i_pl, j_pl, i_mi, j_mi
  +
  +        i = int((p-1)/L)
  +        j = mod(p-1,L)
  +        i_pl = mod(i+1,L)
  +        j_pl = mod(j+1,L)
  +        if (i==0) then
  +            i_mi = L-1
  +        else
  +            i_mi = i-1
  +        end if
  +        if (j==0) then
  +            j_mi = L-1
  +        else
  +            j_mi = j-1
  +        end if
  +
  +        p_i_pl = i_pl*L+j+1
  +        p_i_mi = i_mi*L+j+1
  +        p_j_pl = i*L+j_pl+1
  +        p_j_mi = i*L+j_mi+1
  +        p_i_pl_j_mi = i_pl*L+j_mi+1
  +        p_i_mi_j_pl = i_mi*L+j_pl+1
  +
  +    end subroutine calcu_p_nearest 
  +
  +    subroutine calcu_Hp( p, L, J_H, lambda, fai, theta, Hp)
  +        implicit none
  +        integer :: p, p_i_pl, p_i_mi, p_j_pl, p_j_mi, p_i_pl_j_mi, p_i_mi_j_pl, L
  +        real*8 :: J_H, lambda, fai(L*L), theta(L*L), Hp
  +        Hp = 0.0
  +        call calcu_p_nearest( p, p_i_pl, p_i_mi, p_j_pl, p_j_mi, p_i_pl_j_mi, p_i_mi_j_pl, L )
  +
  +        Hp = Hp-J_H*( sin(theta(p))*sin(theta(p_i_pl))*cos(fai(p)-fai(p_i_pl))+cos(theta(p))*cos(theta(p_i_pl)) )
  +        Hp = Hp-J_H*( sin(theta(p))*sin(theta(p_i_mi))*cos(fai(p)-fai(p_i_mi))+cos(theta(p))*cos(theta(p_i_mi)) )
  +        Hp = Hp-J_H*( sin(theta(p))*sin(theta(p_j_pl))*cos(fai(p)-fai(p_j_pl))+cos(theta(p))*cos(theta(p_j_pl)) )
  +        Hp = Hp-J_H*( sin(theta(p))*sin(theta(p_j_mi))*cos(fai(p)-fai(p_j_mi))+cos(theta(p))*cos(theta(p_j_mi)) )
  +        Hp = Hp-J_H*( sin(theta(p))*sin(theta(p_i_pl_j_mi))*cos(fai(p)-fai(p_i_pl_j_mi))+cos(theta(p))*cos(theta(p_i_pl_j_mi)) )
  +        Hp = Hp-J_H*( sin(theta(p))*sin(theta(p_i_mi_j_pl))*cos(fai(p)-fai(p_i_mi_j_pl))+cos(theta(p))*cos(theta(p_i_mi_j_pl)) )
  +
  +        Hp = Hp+lambda*( sin(theta(p))*sin(fai(p))*cos(theta(p_i_pl))-cos(theta(p))*sin(theta(p_i_pl))*sin(fai(p_i_pl)) )
  +        Hp = Hp-lambda*( sin(theta(p))*sin(fai(p))*cos(theta(p_i_mi))-cos(theta(p))*sin(theta(p_i_mi))*sin(fai(p_i_mi)) )
  +        Hp = Hp+(sqrt(3.0)/2)*lambda*( cos(theta(p))*sin(theta(p_j_pl))*cos(fai(p_j_pl))&
  +            -sin(theta(p))*cos(fai(p))*cos(theta(p_j_pl)) )&
  +            +(1.0/2)*lambda*( sin(theta(p))*sin(fai(p))*cos(theta(p_j_pl))-cos(theta(p))*sin(theta(p_j_pl))*sin(fai(p_j_pl)) )
  +        Hp = Hp-(sqrt(3.0)/2)*lambda*( cos(theta(p))*sin(theta(p_j_mi))*cos(fai(p_j_mi))&
  +            -sin(theta(p))*cos(fai(p))*cos(theta(p_j_mi)) )&
  +            -(1.0/2)*lambda*( sin(theta(p))*sin(fai(p))*cos(theta(p_j_mi))-cos(theta(p))*sin(theta(p_j_mi))*sin(fai(p_j_mi)) )
  +        Hp = Hp+lambda*(1.0/2)*( sin(theta(p))*sin(fai(p))*cos(theta(p_i_pl_j_mi))&
  +            -cos(theta(p))*sin(theta(p_i_pl_j_mi))*sin(fai(p_i_pl_j_mi)) )&
  +            -lambda*(sqrt(3.0)/2)*( cos(theta(p))*sin(theta(p_i_pl_j_mi))*cos(fai(p_i_pl_j_mi))&
  +            -sin(theta(p))*cos(fai(p))*cos(theta(p_i_pl_j_mi)) )
  +        Hp = Hp-lambda*(1.0/2)*( sin(theta(p))*sin(fai(p))*cos(theta(p_i_mi_j_pl))&
  +            -cos(theta(p))*sin(theta(p_i_mi_j_pl))*sin(fai(p_i_mi_j_pl)) )&
  +            +lambda*(sqrt(3.0)/2)*( cos(theta(p))*sin(theta(p_i_mi_j_pl))*cos(fai(p_i_mi_j_pl))&
  +            -sin(theta(p))*cos(fai(p))*cos(theta(p_i_mi_j_pl)) )
  +        !write(*,*) Hp
  +
  +    end subroutine calcu_Hp
  +
  +    subroutine calcu_V( L, J_H, lambda, fai, theta, H, V )
  +        implicit none
  +        integer :: p, L, i, j
  +        real*8 :: J_H, lambda, fai(L*L), theta(L*L), H(L*L), V
  +        V = 0.0
  +
  +        do i = 1,L
  +            do j = 1,L
  +                p = (i-1)*L+j
  +                call calcu_Hp( p, L, J_H, lambda, fai, theta, H(p))
  +                V = V+H(p)
               end do
  -        end do
  -
  +        end do!
  +        V = V/2
  +        !write(*,*) V
   
       end subroutine calcu_V
   
  
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
  