title: Different Ising Model results with Fortran and F2PY
body: <p>I am trying to implement my Fortran code in python using 'f2py'. I derived/extracted a Fortran subroutine from my original Fortran code of Ising model and did calculation using both with same lattice parameters/conditions, etc. But my results are different, which I am unable to understand. Since both should result in similar curves.</p>
  
  <p><strong>Output:</strong></p>
  
  <ol>
  <li>Lattice dimensions = 64 x 64</li>
  <li>Temperature range 0.1 - 5.0</li>
  <li>External Magnetic Field b = 0</li>
  <li>coupling constant j = 0</li>
  <li>No. of MonteCarlo Simulations = 100 + 50</li>
  </ol>
  
  <p>Original Fortran code Result:</p>
  
  <p><img src="https://i.stack.imgur.com/RRrqy.png" alt=""></p>
  
  <p>F2PY based Result:</p>
  
  <p><img src="https://i.stack.imgur.com/m4w2S.png" alt=""></p>
  
  <p><strong><em>Code:</em></strong></p>
  
  <p>Original Fortran Code:</p>
  
  <pre><code>program ising_model
      implicit none
      integer,  allocatable, dimension(:,:)   :: lattice
      real, allocatable, dimension(:)         :: mag
      integer                                 :: row, col, dim, mcs, sweeps, relax_sweeps
      real                                    :: j, dE, B, temp, rval
  
      ! specify the file name
      open(12,file='myoutput.txt')
  
      ! square - lattice specification
      dim = 64
  
      ! coupling constant
      j = 1.0
  
      ! magnetic field
      B = 0.0
  
      ! number of montecarlo simulations
      mcs = 100
  
      ! sparse averaging
      relax_sweeps = 50
  
      allocate(lattice(dim, dim))
      allocate(mag(mcs + relax_sweeps))
      call spin_array(dim, lattice)
      !call outputarray(dim, lattice)
  
      temp = 0.1
      do while (temp .le. 5.0)
          mag = 0
          do sweeps =  1, mcs + relax_sweeps
              ! One complete sweep
              do row = 1, dim
                  do col = 1, dim
                      call EnergyCal(row, col, j, B, lattice, dim, dE)
                      call random_number(rval)
                      if (dE .le. 0) then
                          lattice(row,col) = -1 * lattice(row,col)
                      elseif (exp(-1 * dE / temp) .ge. rval) then
                          lattice(row,col) = -1 * lattice(row,col)
                      else
                          lattice(row,col) = +1 * lattice(row,col)
                      end if
                  end do
              end do
              mag(sweeps) = abs(sum(lattice))/float((dim * dim))
          end do
          write(12,*) temp, sum(mag(relax_sweeps:))/float(mcs + relax_sweeps)
          temp = temp + 0.01
      end do
      !print*,''
      !call outputarray(dim, lattice)
  end program ising_model
  
  !--------------------------------------------------------------
  !   subroutine to print array
  !--------------------------------------------------------------
  subroutine outputarray(dim, array)
      implicit none
      integer                         :: col, row, dim
      integer, dimension(dim, dim)    :: array
      do row = 1, dim
          write(*,10) (array(row,col), col = 1, dim)
  10  format(100i3)
      end do
  end subroutine outputarray
  
  !--------------------------------------------------------------
  !   subroutine to fill the square lattice with spin randomly
  !--------------------------------------------------------------
  subroutine spin_array(dim, array)
      implicit none
      integer                         :: dim, row, col
      real                            :: rval
      integer, dimension(dim, dim)    :: array
      do row = 1, dim
          do col = 1, dim
              call random_number(rval)
              if (rval .ge. 0.5) then
                  array(row, col) = +1
              else
                  array(row, col) = -1
              end if
          end do
      end do
  end subroutine spin_array
  
  !--------------------------------------------------------------
  !   subroutine to calculate energy
  !--------------------------------------------------------------
  subroutine EnergyCal(row, col, j, B, array, dim, dE)
      implicit none
      integer, intent(in)                         :: row, col, dim
      real, intent(in)                            :: j, B
      integer                                     :: left, right, top, bottom
      integer, dimension(dim, dim), intent(in)    :: array
      real, intent(out)                           :: dE
  
      ! periodic boundry condtions
      left = col - 1
      if (col .eq. 1) left = dim
  
      right = col + 1
      if (col .eq. dim) right = 1
  
      top = row - 1
      if (row .eq. 1) top = dim
  
      bottom = row + 1
      if (row .eq. dim) bottom = 1
  
      dE = 2 * j * array(row, col) * ((array(top, col) + array(bottom, col) + &amp;
              array(row, left) + array(row, right)) + B)
  end subroutine EnergyCal
  </code></pre>
  
  <p>F2PY Extracted Code which is converted using F2PY:</p>
  
  <pre><code>subroutine mcmove(inlattice, dim, j, b, temp, finlattice)
      implicit none
      integer                                     :: row, col
      real, intent(in)                            :: j, B, temp
      integer, intent(in)                         :: dim
      integer, dimension(dim,dim), intent(in)     :: inlattice
      integer, dimension(dim,dim), intent(out)    :: finlattice
      real                                        :: rval, de
  
      finlattice = inlattice
  
      do row = 1, dim
          do col = 1, dim
              call EnergyCal(row, col, j, b, finlattice, dim, dE)
              call random_number(rval)
              if (dE .le. 0) then
                  finlattice(row,col) = -1 * finlattice(row,col)
              elseif (exp(-1 * dE / temp) .ge. rval) then
                  finlattice(row,col) = -1 * finlattice(row,col)
              else
                  finlattice(row,col) = +1 * finlattice(row,col)
              end if
          end do
      end do
  end subroutine mcmove
  
  !--------------------------------------------------------------
  !   subroutine to calculate energy
  !--------------------------------------------------------------
  subroutine EnergyCal(row, col, j, b, array, dim, dE)
      implicit none
      integer, intent(in)                         :: row, col, dim
      real, intent(in)                            :: j, b
      integer                                     :: left, right, top, bottom
      integer, dimension(dim, dim), intent(in)    :: array
      real, intent(out)                           :: dE
  
      ! periodic boundry condtions
      left = col - 1
      if (col .eq. 1) left = dim
  
      right = col + 1
      if (col .eq. dim) right = 1
  
      top = row - 1
      if (row .eq. 1) top = dim
  
      bottom = row + 1
      if (row .eq. dim) bottom = 1
  
      dE = 2 * j * array(row, col) * ((array(top, col) + array(bottom, col) + &amp;
              array(row, left) + array(row, right)) + b)
  end subroutine EnergyCal
  </code></pre>
  
  <p>Implemented F2PY + Python (if in case needed)</p>
  
  <pre><code># required packages
  import numpy as np
  import matplotlib.pyplot as plt
  
  # fortran module
  import ising_f90 as fmod
  print(fmod.__doc__)
  
  ## definitions ##
  def spin_field(rows, cols):
      ''' generates a configuration with spins -1 and +1'''
      return np.random.choice([-1, 1], size = (rows, cols))
  
  ## relavant details about the configuration ##
  
  # number of monte carlo simulations
  mcs = 500
  
  # sparse averaging
  relax_sweeps = 50
  
  # square lattice dimensions
  dim = 64
  
  # coupling constant
  j = 1.0
  
  # external magnetic field
  b = 0.0
  
  # temperature - heat bath
  lowTemp = 0.1
  upTemp = 5.0
  
  # initialisation of lattice with random spins
  inlattice = spin_field(dim, dim)
  
  avg_mg = []
  T      = []
  n = 0
  for temp in np.linspace(0.1,5.0,20):
      mag = np.zeros(mcs + relax_sweeps)
      for sweeps in range(0, mcs + relax_sweeps):
          finlattice = fmod.mcmove(inlattice,j,b,temp,dim)
          mag[sweeps] = abs(sum(sum(finlattice))) / (dim * dim)
          n += 1
      T.append(temp)
      avg_mg.append(sum(mag[relax_sweeps:]) / (sweeps+1))
  
  plt.plot(T,avg_mg,'k.--',label='{0} x {0}'.format(str(dim)))
  plt.legend()
  plt.xlabel('Temperature')
  plt.ylabel('Magnetization')
  plt.title('Ising Model 2D')
  plt.show()
  
  </code></pre>