python - Metropolis monte carlo 理想气体模拟量归零
问题描述
我正在尝试使用 Metropolis Monte Carlo 算法根据 Clapeyron 方程“pv=nkbT”对理想气体进行简单的模拟。这是一个非常简单的例子,我考虑二维分子,彼此之间没有相互作用,能量等于 E = pV 其中 V 是包含所有分子的圆的面积。
我的问题是,经过很少的蒙特卡罗步骤,我的气体体积总是几乎为零,无论我放了多少分子或压力。我不知道我的代码中是否有错误,或者是因为我没有没有任何分子相互作用。所有的帮助都会很重要,这是我的代码
我的结果显示在下面的图中,x 轴是蒙特卡罗步骤,y 轴是音量,我预计结果是在更多的步骤之后音量的一些非零恒定值。
import numpy as np
import random
import matplotlib.pyplot as plt
def centroid(arr):
length = arr.shape[0]
sum_x = sum([i.x for i in arr])
sum_y = sum([i.y for i in arr])
return sum_x/length, sum_y/length
class Molecule:
def __init__(self, xpos, ypos):
self.pos = (xpos, ypos)
self.x = xpos
self.y = ypos
class IdealGas:
# CONSTANTS
def __init__(self, n,full_radius, pressure, T, number_of_runs):
gas = []
for i in range(0, n):
gas.append(Molecule(random.random() * full_radius,
random.random() * full_radius))
self.gas = np.array(gas)
self.center = centroid(self.gas)
self.small_radius = full_radius/4.
self.pressure = pressure
self.kbT = 9.36E-18 * T
self.n = n
self.number_of_runs = number_of_runs
def update_pos(self):
random_molecule = np.random.choice(self.gas)
old_state_x = random_molecule.x
old_state_y = random_molecule.y
old_radius = np.linalg.norm(np.array([old_state_x,old_state_y])-np.array([self.center[0],self.center[1]]))
energy_old = np.pi * self.pressure * old_radius**2
random_molecule.x = old_state_x + (random.random() * self.small_radius) * np.random.choice([-1, 1])
random_molecule.y = old_state_y + (random.random() * self.small_radius) * np.random.choice([-1, 1])
new_radius = np.linalg.norm(np.array([random_molecule.x,random_molecule.y])-np.array([self.center[0],self.center[1]]))
energy_new = np.pi * self.pressure * new_radius**2
if energy_new - energy_old <= 0.0:
return random_molecule
elif np.exp((-1.0 * (energy_new - energy_old)) / self.kbT) > random.random():
return random_molecule
else:
random_molecule.x = old_state_x
random_molecule.y = old_state_y
return random_molecule
def monte_carlo_step(self):
gas = []
for molecule in range(0, self.n):
gas.append(self.update_pos())
self.gas = np.array(gas)
#self.center = centroid(self.gas)
return self.gas
def simulation(self):
volume = []
for run in range(self.number_of_runs):
step_gas = self.monte_carlo_step()
step_centroid = centroid(step_gas)
step_radius = max([np.linalg.norm(np.array([gas.x,gas.y])-np.array([step_centroid[0],step_centroid[1]]))
for gas in step_gas])
step_volume = np.pi * step_radius**2
volume.append(step_volume)
return volume
Gas = IdealGas(500, 50, 1000, 300, 100)
vol = Gas.simulation()
plt.plot(vol)
plt.show()
解决方案
如果新半径低于旧半径,您只允许您的分子移动:
if energy_new - energy_old <= 0.0:
相当于:
np.pi * self.pressure * new_radius**2 <= np.pi * self.pressure * old_radius**2
那是:
abs(new_radius) <= abs(old_radius)
所以所有的分子都去中心。
对我来说,你的假设太强了:你固定了温度、压力和分子数量。根据理想气体方程,体积v=nRT/p也是常数。如果体积可以改变,那么压力或温度就必须改变。在您的模拟中,允许压力变化意味着压力和体积的乘积是恒定的,因此能量是恒定的,因此分子可以在任意大的体积中自由移动。
顺便说一句,我认为分子应该初始化为:
(random.random() - 0.5) * full_radius
这样就占据了零附近的所有平面。
推荐阅读
- openid-connect - ASP.Net Core OIDC 如何识别新登录
- python-3.x - 在 Python 中将特定的中文标点符号替换为对应的英文标点符号
- ruby - 在视图中使用“助手”方法
- docker - Redis 哨兵故障转移
- vue.js - 我应该如何将 VeeValidate 3 用于手机上的浏览器?
- javascript - 如何正确处理移动设备的链接标签?
- python - sendVideo 给出错误“错误请求:无法获取 HTTP URL 内容”
- node.js - Bcrypt 没有将散列密码分配给 body.password 变量
- java - 数据验证失败时请求映射方法不运行(Spring MVC)
- python - 将序列化程序错误写入记录器。Django 休息