路径积分蒙特卡洛模拟一维谐波势中原子的运动

Question

我正在尝试模拟原子在一维谐波势中的运动。我想计算原子在谐波势中的位置，然后绘制相对于（任意）时间单位的位置。我想使用路径积分蒙特卡洛方法 (PIMC) 进行模拟。

PIMC 是一种方法，其中一个量子粒子由 P 个“珠子”描述，其中 P 是所谓的 Trotter 数。通过这种方式，我们获得了同构经典系统的配分函数，即 P 珠的“项链/环聚合物”，所有这些都具有相同的势能。

Metropolis MC 的伪代码如下：

Put N*P particles at random positions;
REPEAT
    for i=1 to P DO
        select random P;
        select one of the N particles at P at random;
        generate random displacement of that particle;
        compute r = exp(-delta_tau*(H_new - H_old));
        accept/reject the displacement with normal Metropolis check;
    end for;
UNTIL FINISHED

到目前为止我写的代码只生成稀疏点，它似乎没有正确模拟位置。如果有人知道物理/PIMC 方法，你能帮我修复这段代码吗？我想它应该找到某种平衡点，但情节只显示了几个单独的点。

import math
import numpy as np
import random
import matplotlib.pyplot as plt
%matplotlib inline

def V(x, k=54.869):
    """
    1D potential for quantum harmonic oscillator.
    k = harmonic force constant
    x = position
    """
    return 1/2*k*x**2

def move(x, drmax):
    """
    Simple MC move.
    x = position (1D)
    drmax = maximum displacement parameter
    """
    ran = random.uniform(0, 1)
    pos = x + (ran - 0.5)*drmax
    return pos
    
def metropolis_check(r_old, r_new, k, kB, T):
    """
    Check wheter to accept or reject trial move.
    T = temperature
    kB = Boltzmann constant
    """
    r = math.nan
    if V(k, r_new) < V(k, r_old):
        r = r_new
    else:
        deltaV = V(k, r_new) - V(k, r_old)
        boltzmann = np.exp(-deltaV/(kB*T))
        ran = random.uniform(0, 1)
        if boltzmann > ran:
            r = r_new
        else:
            r = r_old
    return r
    
def PIMC(N, P, T):    
    """
    Metropolis Monte Carlo algorithm.
    N = number of particles
    P = number of beads per particle
    T = temperature
    """
    drmax = 1
    delta_time_slice = 10**(-6)
    k = 54.869 # harmonic force constant
    kB = 1.380649*10**(-23) #Boltzmann constant
    
    #put the NP particles at random positions
    positions = random.sample(range(-N*P, N*P), N*P)
    times = list(np.zeros(len(positions)))
    H = list(map(V, positions))
    
    for j in range(10): # loop so many times that position converges (?)
        for i in range(P):
            # first select time slice P at random
            # then select one of the N particles at time slice P at random
            time_slice = random.choice(range(P))
            times.append(time_slice)
            part = random.choice(range(N))
            bead = positions[time_slice]

            # generate random displacement of that bead (simple MC move)
            x = move(bead, drmax)
            H_old = H[i]
            H_new = V(k, x)
            H[i] = H_new
            r_old = positions[i - 1]
            r_new = np.exp(-delta_time_slice*(H_old - H_new))

            # accept or reject the displacement with metropolis check     
            r = metropolis_check(r_old, r_new, k, kB, T)
            positions.append(r)

    return positions, times

positions, times = PIMC(1, 5, 300)
plt.scatter(times, positions)
plt.ylabel("positions")
plt.xlabel("time")