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:

lattice dimensions = 25 x 25
No. of MonteCarlo simulations = 100
Temperature: 0.1 - 5.0
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?

Number 1 tip from 30 years+ of scientific coding: Don't worry about speed if you are getting the wrong answer. Beyond that you are going to have to tell us in a lot more detail what you are trying to achieve - the Ising model probably means almost nothing to most people here. — Ian Bush, Mar 18 '20 at 14:04
Number 2 tip from 30 years+ of scientific coding: Always develop with all the run time checks turned on. If I use gfortran and compile with -fcheck=all and run I get "At line 126 of file Ising.f90 Fortran runtime error: Index '0' of dimension 1 of array 'array' below lower bound of 1" — Ian Bush, Mar 18 '20 at 14:08
You've screwed up your boundary conditions. As this looks like homework I'll say no more than that except if you were a student of mine I would fail an otherwise basically nice code for not having an interface in scope at the calling point for each subprogram, ALWAYS use modules or contained subprograms for exactly the same reasons that you use Implicit None - see my answer to https://stackoverflow.com/questions/31630535/why-should-i-use-interfaces/31634745#31634745 — Ian Bush, Mar 18 '20 at 14:25
@IanBush thank you for your help. I didn't understand your last comment. Any source would be grateful. I am still novice in Fortran — 147875, Mar 18 '20 at 17:27
@147875 Tho comment was suggesting to place all your subroutines and functions into modules. — Vladimir F Героям слава, Mar 18 '20 at 19:20

score 0 · Accepted Answer · answered Mar 18 '20 at 17:09

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

Phase Transition in 2D Ising Model in Fortran

1 Answers1