Fortran和F2PY的不同伊辛模型结果

Question

我正试图用“f2py”在python中实现我的Fortran代码。我从伊辛模型的原始Fortran代码中导出/提取了一个Fortran子程序，并使用相同的晶格参数/条件等进行了计算。但我的结果不同，我无法理解。因为两者都会产生相似的曲线

输出：

1. 晶格尺寸=64 x 64
2. 温度范围0.1-5.0
3. 外磁场b=0
4. 耦合常数j=0
5. 蒙特卡洛模拟的数量=100+50

原始Fortran代码结果：

基于F2PY的结果：

代码：

原始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 = 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) + &
            array(row, left) + array(row, right)) + B)
end subroutine EnergyCal

使用F2PY转换的F2PY提取代码：

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) + &
            array(row, left) + array(row, right)) + b)
end subroutine EnergyCal

实现了F2PY+Python（如果需要）

# 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()