Python Grand Canonical Ensemble中LennardJones 2D相互作用对相关分布函数的计算
Python calculation of LennardJones 2D interaction pair correlation distribution function in Grand Canonical Ensemble
编辑
我认为直方图的归一化存在问题，因为必须除以每个元素的半径。
我正在尝试使用 python3 计算二维 LennardJones(LJ) 系统的粒子数波动和径向分布函数。虽然我相信粒子波动是正确的，但对相关性 g(r) 在小距离内是正确的，但随后会爆炸（计算使用 numpy 的直方图方法）。
问题是，我无法找出为什么会出现这种行为 - 也许是对某种方法的误解？就这样，我把相关代码贴在下面，如果需要，我也可以上传其他部分代码或整个脚本。
首先请注意，由于我们使用的是 Grand-Canonical Ensemble，因此随着粒子数量的变化，存储粒子的数组也会发生变化 - 也许这是另一个错误点实现可能存在。
粒子删除或插入
def mcex(L,npart,particles,beta,rho0,V,en):
factorin=(rho0*V)/(npart+1)
factorout=(npart)/(V*rho0)
print("factorin=",factorin)
print("factorout",factorout)
# Produce random number and check:
rand = random.uniform(0, 1)
if rand <= 0.5:
# Insert a particle at a random location
x_new_coord = random.uniform(0, L)
y_new_coord = random.uniform(0, L)
new_particle = [x_new_coord,y_new_coord]
new_E = particleEnergy(new_particle,particles, npart+1)
deltaE = new_E
print("dEin=",deltaE)
# Acceptance rule for inserting
if(deltaE>10):
P_in=0
else:
P_in = (factorin) *math.exp(-beta*deltaE)
print("pinacc=",P_in)
rand= random.uniform(0, 1)
if rand <= P_in :
particles.append(new_particle)
en += deltaE
npart += 1
print("accepted insertion")
else:
if npart != 0:
p = random.randint(0, npart-1)
this_particle = particles[p]
prev_E = particleEnergy(this_particle, particles, p)
deltaE = prev_E
print("dEout=",deltaE)
# Acceptance rule for removing
if(deltaE>10):
P_re=1
else:
P_re = (factorout)*math.exp(beta*deltaE)
print("poutacc=",P_re)
rand = random.uniform(0, 1)
if rand <= P_re :
particles.remove(this_particle)
en += deltaE
npart = npart - 1
print("accepted removal")
print()
return particles, en, npart
Monte Carlo 相关部分：对于 1/10 运行，检查插入或删除粒子的可能性
# MC
for step in range(0, runTimes):
print(step)
print()
rand = random.uniform(0,1)
if rand <= 0.9:
#----------- change energies-------------------------
#........
#........
else:
particles, en, N = mcex(L,N,particles,beta,rho0,V, en)
# stepList.append(step)
if((step+1)%1000==0):
for i, particle1 in enumerate(particles):
for j, particle2 in enumerate(particles):
if j!= i:
# print(particle1)
# print(particle2)
# print(i)
# print(j)
dist.append(distancesq(particle1, particle2))
NList.append(N)
我们调用函数mcex的地方可能粒子数组没有正确更新:
def mcex(L,npart,particles,beta,rho0,V,en):
factorin=(rho0*V)/(npart+1)
factorout=(npart)/(V*rho0)
print("factorin=",factorin)
print("factorout",factorout)
# Produce random number and check:
rand = random.uniform(0, 1)
if rand <= 0.5:
# Insert a particle at a random location
x_new_coord = random.uniform(0, L)
y_new_coord = random.uniform(0, L)
new_particle = [x_new_coord,y_new_coord]
new_E = particleEnergy(new_particle,particles, npart+1)
deltaE = new_E
print("dEin=",deltaE)
# Acceptance rule for inserting
if(deltaE>10):
P_in=0
else:
P_in = (factorin) *math.exp(-beta*deltaE)
print("pinacc=",P_in)
rand= random.uniform(0, 1)
if rand <= P_in :
particles.append(new_particle)
en += deltaE
npart += 1
print("accepted insertion")
else:
if npart != 0:
p = random.randint(0, npart-1)
this_particle = particles[p]
prev_E = particleEnergy(this_particle, particles, p)
deltaE = prev_E
print("dEout=",deltaE)
# Acceptance rule for removing
if(deltaE>10):
P_re=1
else:
P_re = (factorout)*math.exp(beta*deltaE)
print("poutacc=",P_re)
rand = random.uniform(0, 1)
if rand <= P_re :
particles.remove(this_particle)
en += deltaE
npart = npart - 1
print("accepted removal")
print()
return particles, en, npart
最后，我们创建 g(r) 直方图
可能归一化或直方图方法的使用不尽如人意
RDF(N,particles,L)
具有以下功能：
def RDF(N,particles, L):
minb=0
maxb=8
nbin=500
skata=np.asarray(dist).flatten()
rDf = np.histogram(skata, np.linspace(minb, maxb,nbin))
prefactor = (1/2/ np.pi)* (L**2/N **2) /len(dist) *( nbin /(maxb -minb) )
# prefactor = (1/(2* np.pi))*(L**2/N**2)/(len(dist)*num_increments/(rMax + 1.1 * dr ))
rDf = [prefactor*rDf[0], 0.5*(rDf[1][1:]+rDf[1][:-1])]
print('skata',len(rDf[0]))
print('incr',len(rDf[1]))
plt.figure()
plt.plot(rDf[1],rDf[0])
plt.xlabel("r")
plt.ylabel("g(r)")
plt.show()
结果是：
粒子N数波动
和
[
但我们想要
您的 g(r) 未正确归一化。您需要将每个箱中找到的对数除以系统的平均密度乘以与该箱关联的环空面积，其中后者仅为 2 pi r dr，其中 r 是箱的中点，dr 是箱尺寸。据我所知，您的前置因子不包含 "r" 位。还缺少其他一些东西，但如果不知道所有其他常量包含什么，就很难说清楚。
编辑：here 是一个 link，它将指导您执行计算 2D 和 3D 径向分布函数的例程
虽然我已经接受了一个答案，但我会在这里发布更多细节。
要正确规范化对相关性，必须除以每个 "number of particles found at a certain distance" 或数学上距离的 delta 函数之和，必须除以它自身的距离。
首先要了解 numpy.histogram 是一个包含两个元素的数组，第一个元素是所有计数事件的数组，第二个元素是 bin 的向量，必须取第一个数组的每个元素，假设 np.histogram[0] 并将其与第二个数组的 np.histogram[1] 成对相乘。
也就是说，必须做到以下几点：
def RDF(N,particles, L):
minb=0
maxb=25
nbin=200
width=(maxb-minb)/(nbin)
rings=np.linspace(minb, maxb,nbin)
skata=np.asarray(dist).flatten()
rDf = np.histogram(skata, rings ,density=True)
prefactor = (1/( np.pi*(L**2/N**2)))
rDf = [prefactor*rDf[0], 0.5*(rDf[1][1:]+rDf[1][:-1])]
rDf[0]=np.multiply(rDf[0],1/(rDf[1]*( width )))
在最后一条乘法线之前，我们将 bins 居中，以便它们的数量等于第一个数组的元素数量（你有五个手指，但它们之间有四个中间间隙）
编辑 我认为直方图的归一化存在问题，因为必须除以每个元素的半径。
我正在尝试使用 python3 计算二维 LennardJones(LJ) 系统的粒子数波动和径向分布函数。虽然我相信粒子波动是正确的，但对相关性 g(r) 在小距离内是正确的，但随后会爆炸（计算使用 numpy 的直方图方法）。
问题是，我无法找出为什么会出现这种行为 - 也许是对某种方法的误解？就这样，我把相关代码贴在下面，如果需要，我也可以上传其他部分代码或整个脚本。
首先请注意，由于我们使用的是 Grand-Canonical Ensemble，因此随着粒子数量的变化，存储粒子的数组也会发生变化 - 也许这是另一个错误点实现可能存在。
粒子删除或插入
def mcex(L,npart,particles,beta,rho0,V,en):
factorin=(rho0*V)/(npart+1)
factorout=(npart)/(V*rho0)
print("factorin=",factorin)
print("factorout",factorout)
# Produce random number and check:
rand = random.uniform(0, 1)
if rand <= 0.5:
# Insert a particle at a random location
x_new_coord = random.uniform(0, L)
y_new_coord = random.uniform(0, L)
new_particle = [x_new_coord,y_new_coord]
new_E = particleEnergy(new_particle,particles, npart+1)
deltaE = new_E
print("dEin=",deltaE)
# Acceptance rule for inserting
if(deltaE>10):
P_in=0
else:
P_in = (factorin) *math.exp(-beta*deltaE)
print("pinacc=",P_in)
rand= random.uniform(0, 1)
if rand <= P_in :
particles.append(new_particle)
en += deltaE
npart += 1
print("accepted insertion")
else:
if npart != 0:
p = random.randint(0, npart-1)
this_particle = particles[p]
prev_E = particleEnergy(this_particle, particles, p)
deltaE = prev_E
print("dEout=",deltaE)
# Acceptance rule for removing
if(deltaE>10):
P_re=1
else:
P_re = (factorout)*math.exp(beta*deltaE)
print("poutacc=",P_re)
rand = random.uniform(0, 1)
if rand <= P_re :
particles.remove(this_particle)
en += deltaE
npart = npart - 1
print("accepted removal")
print()
return particles, en, npart
Monte Carlo 相关部分：对于 1/10 运行，检查插入或删除粒子的可能性
# MC
for step in range(0, runTimes):
print(step)
print()
rand = random.uniform(0,1)
if rand <= 0.9:
#----------- change energies-------------------------
#........
#........
else:
particles, en, N = mcex(L,N,particles,beta,rho0,V, en)
# stepList.append(step)
if((step+1)%1000==0):
for i, particle1 in enumerate(particles):
for j, particle2 in enumerate(particles):
if j!= i:
# print(particle1)
# print(particle2)
# print(i)
# print(j)
dist.append(distancesq(particle1, particle2))
NList.append(N)
我们调用函数mcex的地方可能粒子数组没有正确更新:
def mcex(L,npart,particles,beta,rho0,V,en):
factorin=(rho0*V)/(npart+1)
factorout=(npart)/(V*rho0)
print("factorin=",factorin)
print("factorout",factorout)
# Produce random number and check:
rand = random.uniform(0, 1)
if rand <= 0.5:
# Insert a particle at a random location
x_new_coord = random.uniform(0, L)
y_new_coord = random.uniform(0, L)
new_particle = [x_new_coord,y_new_coord]
new_E = particleEnergy(new_particle,particles, npart+1)
deltaE = new_E
print("dEin=",deltaE)
# Acceptance rule for inserting
if(deltaE>10):
P_in=0
else:
P_in = (factorin) *math.exp(-beta*deltaE)
print("pinacc=",P_in)
rand= random.uniform(0, 1)
if rand <= P_in :
particles.append(new_particle)
en += deltaE
npart += 1
print("accepted insertion")
else:
if npart != 0:
p = random.randint(0, npart-1)
this_particle = particles[p]
prev_E = particleEnergy(this_particle, particles, p)
deltaE = prev_E
print("dEout=",deltaE)
# Acceptance rule for removing
if(deltaE>10):
P_re=1
else:
P_re = (factorout)*math.exp(beta*deltaE)
print("poutacc=",P_re)
rand = random.uniform(0, 1)
if rand <= P_re :
particles.remove(this_particle)
en += deltaE
npart = npart - 1
print("accepted removal")
print()
return particles, en, npart
最后，我们创建 g(r) 直方图
可能归一化或直方图方法的使用不尽如人意
RDF(N,particles,L)
具有以下功能：
def RDF(N,particles, L):
minb=0
maxb=8
nbin=500
skata=np.asarray(dist).flatten()
rDf = np.histogram(skata, np.linspace(minb, maxb,nbin))
prefactor = (1/2/ np.pi)* (L**2/N **2) /len(dist) *( nbin /(maxb -minb) )
# prefactor = (1/(2* np.pi))*(L**2/N**2)/(len(dist)*num_increments/(rMax + 1.1 * dr ))
rDf = [prefactor*rDf[0], 0.5*(rDf[1][1:]+rDf[1][:-1])]
print('skata',len(rDf[0]))
print('incr',len(rDf[1]))
plt.figure()
plt.plot(rDf[1],rDf[0])
plt.xlabel("r")
plt.ylabel("g(r)")
plt.show()
结果是：
粒子N数波动
和
[
但我们想要
您的 g(r) 未正确归一化。您需要将每个箱中找到的对数除以系统的平均密度乘以与该箱关联的环空面积，其中后者仅为 2 pi r dr，其中 r 是箱的中点，dr 是箱尺寸。据我所知，您的前置因子不包含 "r" 位。还缺少其他一些东西，但如果不知道所有其他常量包含什么，就很难说清楚。
编辑：here 是一个 link，它将指导您执行计算 2D 和 3D 径向分布函数的例程
虽然我已经接受了一个答案，但我会在这里发布更多细节。
要正确规范化对相关性，必须除以每个 "number of particles found at a certain distance" 或数学上距离的 delta 函数之和，必须除以它自身的距离。
首先要了解 numpy.histogram 是一个包含两个元素的数组，第一个元素是所有计数事件的数组，第二个元素是 bin 的向量，必须取第一个数组的每个元素，假设 np.histogram[0] 并将其与第二个数组的 np.histogram[1] 成对相乘。
也就是说，必须做到以下几点：
def RDF(N,particles, L):
minb=0
maxb=25
nbin=200
width=(maxb-minb)/(nbin)
rings=np.linspace(minb, maxb,nbin)
skata=np.asarray(dist).flatten()
rDf = np.histogram(skata, rings ,density=True)
prefactor = (1/( np.pi*(L**2/N**2)))
rDf = [prefactor*rDf[0], 0.5*(rDf[1][1:]+rDf[1][:-1])]
rDf[0]=np.multiply(rDf[0],1/(rDf[1]*( width )))
在最后一条乘法线之前，我们将 bins 居中，以便它们的数量等于第一个数组的元素数量（你有五个手指，但它们之间有四个中间间隙）