如何增加 eqsteps 而不会在一维伊辛模型代码中出现太大波动？

Question

我已经为 100 个格点编写了 1-D Ising 模型的 Fortran 90 代码。我计算了磁化强度、能量、磁化率和比热，并使用 gnuplot 进行了绘图。我在 E~t 和 m~t 图上得到了几乎预期的结果，但我不确定 c~t 和 sus~t 图是否正确。

我还采取了 100000 eqsteps 来达到平衡。如果我将 eqsteps 减少到 100 或 200，我会得到太大的波动。请帮助我如何用较小的 eqsteps 获得结果。

program ising1d
        implicit none
        integer,allocatable,dimension(:)::lattice
    real,allocatable,dimension(:)::mag
        integer,dimension(100000)::seed=12345678
        integer::i,row,d,p,eqsteps
        integer::mcs,w
    real::j,t,rval1,rval2,average_m
    real::y,dE,sum_m,sum_E,average_E,E
    real::sum_m2,sum_E2,average_m2,average_E2,c,sus
    real,dimension(20000)::mm,EE
        
        j=1.0
        t=0.01
        d=100
        allocate(lattice(d))
        
        
        open (unit = 9, file = "ss.txt", action = "write", status = "replace")
    
    do while (t <= 10)
        
        sum_E=0
        sum_m=0
        sum_m2=0
        average_m2=0
        w=1
        
        do mcs=1,200
            do row=1,d
                
                call random_number(rval1)
                if (rval1 .ge. 0.5)then
                        lattice(row)=1.0
                else
                    lattice(row)=-1.0
                end if
            end do
!          This loop is for taking measurements 
        
!          This loop is for getting equilibrium     
                do eqsteps=1,100000
!                          print*,lattice
!                          choosing an random site to flip                          
                        call random_number(y)   
                        p=1+floor(y*d)
!                          print*,"the flipping site is  :",p

!                          boundary condition               
                        if(p==d) then
                                lattice(p+1)=lattice(1)
                        else if (p==1)then
                                lattice(p-1)=lattice(d)
                        else
                                lattice(p)=lattice(p)
                        end if
!                          energy change                
                        dE=2*(lattice(p))*(lattice(p-1)+lattice(p+1))
!                          print*,"the change of energy is  :",dE
!                          metropolis part      
                        if (dE .le. 0) then
                                lattice(p)=-lattice(p)
!                   print*, "flipped"
                        else if (dE .gt. 0) then
                                call random_number(rval2)
                                if (rval2 .lt. exp(-1.0*dE/t)) then
                                        lattice(p)=-lattice(p)
!                       print*, "flipped"
                                else
                                        lattice(p)=lattice(p)
!                       print*, " not flipped"
                                end if
                        else
                                lattice(p)=lattice(p)

                        end if
                end do
            mm(w)=abs(sum(lattice))
            sum_m = sum_m + mm(w)
            
            sum_m2=sum_m2 + mm(w)*mm(w)
            
            
            
            call calc_energy(d,row,E)
            EE(w)=E
            sum_E=sum_E+EE(w)
            
            sum_E2=sum_E2+EE(w)*EE(w) 


                w=w+1
        end do
        average_m=sum_m/(w*d)
        average_E=sum_E/(w*d)
        
        average_m2=sum_m2/(w*w*d*d)
        average_E2=sum_E2/(w*w*d*d)
        
        c=(average_E2-average_E*average_E)/(t*t)
        
        sus=(average_m2-average_m*average_m)/t
        
        write(9,*) t,average_m,average_E,c,sus
        print*,t,average_m,average_E,c,sus
        t = t + 0.01
        end do
        close(9)
        
    
    
contains
!energy calculation

subroutine calc_energy(d,row,E)
    integer,intent(in)::d,row
    real, intent(out)::E
    integer::k
    real::E00=0
    
    do k=2,d-1
        E00=E00+lattice(k)*(lattice(k-1)+lattice(k+1))
    end do
    E=-0.5*(E00+lattice(1)*(lattice(2)+lattice(d))+lattice(d)*(lattice(d-1)+lattice(1)))
    E00=0
end subroutine calc_energy
        
end program ising1d

我期望以较小的 eqstep 大小获得结果。