【发布时间】：2020-03-18 13:38:11
【问题描述】：
我正在使用 Metropolis-Montecarlo 算法研究 2D Ising 模型。当我绘制平均磁化强度与温度的关系图时，相变应该在 2.5K 左右，但我的相变在 0.5 和 1.0 之间。但是用 python 编写的类似代码给了我正确的结果。
输出：
- 格子尺寸 = 25 x 25
- 没有。蒙特卡洛模拟 = 100
- 温度：0.1 - 5.0
- 外磁场 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
有什么地方可以使用周期性边界条件等内置函数减少代码行数并提高仿真速度？
【问题讨论】：
-
30 多年科学编码的第 1 条提示：如果您得到错误的答案，请不要担心速度。除此之外，您还必须更详细地告诉我们您想要实现的目标 - Ising 模型对这里的大多数人来说可能几乎没有任何意义。
-
30 多年科学编码的第二个提示：始终在打开所有运行时检查的情况下进行开发。如果我使用 gfortran 并使用 -fcheck=all 编译并运行，我会得到“在文件 Ising.f90 Fortran 运行时错误的第 126 行：数组 'array' 的第 1 维的索引 '0' 低于 1 的下限”
-
你搞砸了你的边界条件。因为这看起来像家庭作业，所以我只想说，除非你是我的学生，否则我会失败一个基本上很好的代码，因为在每个子程序的调用点没有范围内的接口，总是使用模块或包含的子程序出于与您使用 Implicit None 完全相同的原因 - 请参阅我对 stackoverflow.com/questions/31630535/… 的回答
-
@IanBush 感谢您的帮助。我不明白你最后的评论。任何来源将不胜感激。我还是 Fortran 新手
-
@147875 评论建议将所有子例程和函数放入模块中。
标签： fortran physics montecarlo