Fortran二维Ising模型的相变

Question

我工作的二维伊辛模型使用大都会-Montecarlo算法。当我绘制平均磁化强度随温度变化时，相变应该在2.5K左右，而我的相变在0.5到1.0之间。但是用python编写的类似代码给出了正确的结果。

输出：

1. 晶格维数=25x25
2. No.100
3. Temperature: 0.1 -5.0
4. 外部磁场B= 0

​

​

Fortran代码：

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

有什么地方，哪里可以减少代码行数使用任何内置的函数，如周期边界条件，并提高模拟的速度？

Answer 1

“谢谢你，”伊恩·布什(Ian)暗示了这个错误。我完全搞砸了我的周期边界条件，这个代码应该是

! 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

或

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