VMC 算法中的意外行为

Question

我正在用 Fortran 语言编写二维 N 粒子硬球气体被困在谐波势中的模拟。我的代码应该最大化以下功能： 其中第一项是粒子之间的相互作用项，第二项是粒子与简谐势之间的相互作用项。这是通过以下函数在我的代码中实现的：

real*8 function trial_WF(R)
    ....
    trial_WF = 1
    !atoms interaction term
    do i_atom=1, Natoms
        do j_atom = 1, Natoms
            if (i_atom /= j_atom) then 
                dist = norm2(R(i_atom,:)-R(j_atom,:))
                trial_WF = trial_WF*twobody_corr(dist)
            end if 
        end do 
    end do
    !harmonic oscillator interaction term
    do i_atom = 1, Natoms
        dist = norm2(R(i_atom,:))
        trial_WF = trial_WF*harmonic_GS(dist)
    end do 

蒙特卡洛是通过以下方式实现的：

call init_random_seed()
allocate(walker(Natoms,DIM,2))
NEW = 2; OLD = 1  !init swappers
do MC_step = - NStabSteps, NMCsteps
    do step = 1, NThermSteps!thermalization steps to avoid correlated results
        call diffuse(walker(:,:,OLD),walker(:,:,NEW),sigma)
        TWF(OLD) = trial_WF(walker(:,:,OLD)) 
        TWF(NEW) = trial_WF(walker(:,:,NEW))

        !metropolis question
        call random_number(c1)
        if( (TWF(NEW)/TWF(OLD))**2 > c1  ) then 
            OLD = 3 - OLD; NEW = 3 - NEW !swap NEW <--> OLD
            Nacceptances = Nacceptances + 1
        end if 
        COUNTER = COUNTER + 1 !total cycles
    end do 
    !update accumulators
    if(MC_step > 0) then 
        E      = Elocal(walker(:,:,OLD))
        E_avg  = E_avg + (E)/NMCsteps
        E2_avg = E2_avg + (E**2)/NMCsteps
        call print_states(walker(:,:,OLD)) !energy acc/ density profile
    end if 
end do

问题： 所有粒子似乎都表现正常，被困在彼此排斥的谐波势的中心，但最后一个粒子在正 y 轴中经历了漂移，如下图所示：

当然，我考虑过数据泄漏，因为它基本上只影响我位置的最后一个数组单元，但由于某些原因，情况似乎并非如此：

• 如果我关闭粒子交互评论两个循环，这种奇怪的行为就会消失
• 如果我用一些空单元扩展粒子数组，则不会发生任何变化
• 如果我不将值固定为某个常数（如 0.6），而不是随机大都会问题，则这种原本不太可能出现的意外行为会再次消失
• 如果我更改代码中的粒子数，漂移仍然存在

还有：

• 没有理由漂移应该沿 y 轴不对称，因为相互作用项取决于粒子距离的模，因此无法区分 x、y 及其各自的符号。
• 我还考虑了随机数之间的一些奇怪的相关性，但我的 init_random_seed() 应该非常强大，并且如果我更改粒子数，也不会发生任何变化。

Answer 1

编辑，我发现了错误：我在代码中覆盖了一个变量，这导致了粒子的奇怪行为。无论如何，感谢您的评论，也是因为我在尝试写下“最小的、可重现的示例”时发现了它