Я пытаюсь реализовать свой код на Фортране в python, используя 'f2py'. Я извлек / извлек подпрограмму на Фортране из моего исходного кода модели Изинга на Фортране и провел вычисления, используя оба параметра с одинаковыми параметрами / условиями решетки и т. Д. c. Но мои результаты разные, что я не могу понять. Поскольку оба должны приводить к одинаковым кривым.
Выход:
- Размеры решетки = 64 x 64
- Диапазон температур 0,1 - 5,0
- Внешнее магнитное поле c Поле b = 0
- Константа связи j = 0
- Нет. MonteCarlo Simulations = 100 + 50
Оригинальный код Фортрана Результат:
На основе F2PY Результат:
Код:
Оригинальный код Фортрана:
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()