Phase Transition in 2D Ising Model in Fortran

Question

I am working on 2D Ising model using Metropolis-Montecarlo Algorithm. When I plot average magnetization vs temperature, the phase transition should be around 2.5K, but my phase transition is between 0.5 and 1.0. But a similar code written in python gives me correct results.

Output:

1. lattice dimensions = 25 x 25
2. No. of MonteCarlo simulations = 100
3. Temperature: 0.1 - 5.0
4. External Magnetic field B = 0

Fortran 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 = 25

    ! 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, id
    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

    if (row == 1) then
        top = dim
    else 
        top = row - 1
    end if 

    if (row == dim) then
        bottom = 1
    else 
        bottom = row - 1
    end if

    if (col == 1) then
        left = dim
    else 
        left = col - 1
    end if

    if (col == dim) then
        right = 1
    else 
        right = col - 1
    end if

    dE = 2 * j * array(row, col) * ((array(top, col) + array(bottom, col) + &
            array(row, left) + array(row, right)) + B * sum(array))
end subroutine EnergyCal

Are there any places, where can I reduce the number of lines of codes using any inbuilt function like periodic boundary conditions and improve the the speed of simulation?

Answer 1

Thank you @Ian Bush for hint about the error. I completely screwed my periodic boundary conditions, this code should be

! 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

or

if (row == 1) then
     top = dim
else 
    top = row - 1
end if 

if (row == dim) then
    bottom = 1
else 
    bottom = row + 1
end if

if (col == 1) then
    left = dim
else 
    left = col - 1
end if

if (col == dim) then
    right = 1
else 
    right = col + 1
end if