我正在使用 MC 模拟来模拟 Lennard-Jones 流体。该代码始终使用缩减单位。我想找到系统的势能。实施了周期性边界条件。我已经用 2 个粒子进行了模拟，以确定模拟是否给出了正确的结果。对于 2 个粒子，平衡分离为$2^{1/6}$，所以能量是-1。（拿$\sigma =1$) . 代码如下。

import numpy as np
import random
import matplotlib.pyplot as plt


# reduced units:
#T(reduced) = kT/epsilon | r(reduced) = r/sigma | U(reduced) = U/epsilon
# General Parameters
DIM=2
npart=2
L=10
volume=L**DIM
density=npart/volume
print("volume = ", volume, " density = ", density,"Number of atoms =",npart)
print ("L is " ,L)
T = 2; Nsteps = 50000; maxdr =0.0001;printfeq=100;DIM=2
#System parameters
beta=1/T



def E(dr2):
    """Returns LJ 6-12 interaction energy for a particular distance between
       2 particles"""

    return  4*((dr2)**(-6) - (dr2)**(-3)) # r is given in unit of sigma. dr2 is distance^2

def P(x):
    """ gives boltzman factor for position at x"""
    return np.exp(-(beta*E(x))) 

def PBC(L,pos):       

    """PBC check for dim dimension system with equal length L in all dimension. 
    INPUT : position array,length,dimension
    OUTPUT: New position
    """
    for k in range(DIM):
            if (pos[k]>0.5):
                pos[k]=pos[k]-1
            if (pos[k]<-0.5):
                pos[k]=pos[k]+1


    return (pos)


def distance(current_position):
    """Takes the current position array of the configuration and finds out 
    distance between each pairs. Neglectd if distance > rcutoff
    INPUT: Array of current position of each particles
    OUTPUT: Array containing distances between each pair of LJ particles.
    """
    Distances=[]
    for i in range (npart):
        for j in range (i+1,npart):
            dr=(current_position[i]-current_position[j])*L
            dr2=np.dot(dr,dr)

            if (dr2!=0):
                Distances.append(dr2)

    return Distances



Energy=[0 for _ in range (Nsteps)]
Distances=[0 for _ in range (Nsteps)]
current_position=np.zeros([npart,DIM])



#------------------ Initialise the Setup with particles distributed uniformly ------------

ip=-1
x=0
y=0
lim=int(np.sqrt(npart))+1
for i in range(0,lim):
    for j in range(0,lim):
        if(ip<npart):
            x=i*(1/lim)
            y=j*(1/lim)
            current_position[ip]=np.array([x,y])
            ip=ip+1
        else:
            break
MassCentre = np.sum(current_position,axis=0)/npart
current_position=current_position-MassCentre


Distances[0]=distance(current_position)



for i in Distances[0]:

    Energy[0]+=E(i)
print(Energy[0])


# -------------------------MC Simulation ----------------------------
rejected=0
for step in range(1,Nsteps):
    if (step%printfeq==0):
        print ("Completed ",step,"steps")
    trial_position=np.zeros([npart,DIM])

    trial_energy=0

    for i in range (npart):

        displacex=(random.uniform(0,1)-0.5)*maxdr 
        displacey=(random.uniform(0,1)-0.5)*maxdr 

        pos=current_position[i]+np.array([displacex,displacey])
        trial_position[i]=PBC(L,pos)
    Distances[step]=distance(trial_position)
    #print(trial_position)

    for i in ((Distances[step])):
        trial_energy+= E(i)
    if (trial_energy<Energy[step-1]):
        current_position=trial_position
        Energy[step]=trial_energy 
    else:
        delta=trial_energy-Energy[step-1]
        if (random.random()<P(delta)):

            current_position=trial_position
            Energy[step]=trial_energy 

        else:
            rejected+=1

            Energy[step]=Energy[step-1]



print(Energy)           
print ("Rejected moves ",rejected,"out of",Nsteps)
steps=np.arange(0,Nsteps)            
plt.figure(1)    
plt.plot(steps, Energy,'o',label='simulation result')
#plt.plot(steps, Distances,'b-',label='simulation result')
plt.xlabel(' steps')
plt.ylabel('Energy of configuration')
plt.show()

我使用了 50000 步，没有单独的平衡步骤。我得到的输出就像

我想知道这种行为是否正常。似乎并非如此。势能突然跃升至-1。还是我做错了什么？在MC中，我没有关于动能的信息，所以我无法检查能量守恒。提前谢谢你。