anchor.diff:          theta(d,i) = x*PI
       end do
       !ij = 0
  -do ii = 1, num_T_para 
  +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
  -    p = 0
  -    call calcu_V( L, N, J_H, lambda, fai(d,:), theta(d,:), V )
  -    call calcu_deltaV(L, N, B_z, theta(d,:), deltaV)
  -    !ln_Weight = -beta*(V+deltaV)
  +    !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()
  -        !p = int(x*N)+1 
  -        p = p+1
  -        !Write(*,*) p
  +        x = ran()
  +        r = int(x*N)+1 
  +        !r = r+1
  +        !Write(*,*) r
           fai(dd,:) = fai(d,:)
           theta(dd,:) = theta(d,:)
           x = ran()
  -        fai(dd, p) = fai(d,p)+(x-0.5)*2*2*PI/12
  -        if ( fai(dd,p) < 0 ) then
  -            fai(dd,p) = fai(dd,p)+2*PI
  -        end if
  -        if ( fai(dd,p) > 2*PI ) then
  -            fai(dd,p) = fai(dd,p)-2*PI
  -        end if
  -
  -        x = ran()
  -        theta(dd, p) = theta(d,p)+(x-0.5)*2*PI/6
  -        if ( theta(dd,p) < 0 ) then
  -            theta(dd,p) = theta(dd,p)+PI
  -        end if
  -        if (theta(dd,p) > PI ) then
  -            theta(dd,p) = theta(dd,p)-PI
  -        end if
  -        call calcu_V( L, N, J_H, lambda, fai(dd,:), theta(dd,:), V_new )
  -        call calcu_deltaV(L, N, B_z, theta(dd,:), deltaV_new)
  -        !ln_Weight_new = -beta*(V+deltaV)
  -
  -        x = ran()
  -        !write(*,*) x, '<?', exp( ln_Weight_new-ln_Weight )
  -        !write(*,*) ln_Weight, ln_Weight_new
  -        !write(*,*) x, '<?', exp( -beta*(V_new+deltaV_new)+beta*(V+deltaV) )
  -        !write(*,*) -beta*(V+deltaV), -beta*(V_new+deltaV_new)
  -        if ( x > exp( -beta*(V_new+deltaV_new)+beta*(V+deltaV) ) ) then
  -            fai(dd,p) = fai(d,p)
  -            theta(dd, p) = theta(d, p)
  +        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_new
  -            deltaV = deltaV_new
  -            !ln_Weight = ln_Weight_new
  +            V = V+Hp-H(r)
  +            H(r) = Hp
           end if
           d = dd
  -        !write(*,*) i/(num_MC/100)
  -        if (ii == num_T_para .AND. i > num_MC-aver_num*N ) then
  +        !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
  -                energyV = energyV+V
  -                energyD = energyD+deltaV
  -                !ij = ij+1
  -                !write(13,*) ij
                   do j = 1, N
  -                    write(13,*) j, theta(d,j), fai(d,j)
  +                    write(14,*) theta(d,j), fai(d,j)
                   end do
               end if
           end if
  -
  -        if ( p == N ) then
  -            p = 0
  -        end if
  -
  +        !if ( r == N ) then
  +        !    r = 0
  +        !end if
  +        write(13,*) V
       end do
  +    close(13)
  +    close(14)
   end do
  -energyV = energyV/aver_num
  -energyD = energyD/aver_num
  -open(unit = 18, file = 'enenrgy.out')
  -write(18,*) energyV, energyD
  -
   
       write(*,*) 'end'
   end program skyrmionlattice
  
positive.code: 
  
  program skyrmionlattice
      use dou_ex 
      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
      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
      real*8, dimension(2,N) :: fai, theta
      real*8, dimension(num_T_para) :: T_para
      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/)
      call initran(1)
      open( unit = 13, file = 'lattice.out' )
  
      do i = 1, N
          x = ran()
  
negative.code: 
      end subroutine calcu_deltaV
      
  end module modu
  
  real(8) function ran()
  !----------------------------------------------!
  