我正在编写一个脚本,该脚本应该使用伦纳德-琼斯势来模拟气体分子的相互作用。我面临两个问题。
np.sqrt()
计算分子之间的距离时,我只得到由 nan
值组成的列表......这不是很有帮助。np.sqrt()
替换 ()**0.5
后,我解决了这个问题。但我仍然得到错误的结果。这里讨论的行是以下一行
dist12 = (r12[0]**2 + r12[1]**2)**0.5
。 然而,摆脱 nan 警告并没有完全解决我的问题,因为后来我检查 if dist12 == 0: return 0, 0
这导致位置列表(除了前四个粒子的列表之外)组成的列表主要包含零。
所以代码流程可能有问题。但我无法找到问题所在。因此,让我首先简要描述一下预期的代码流程。
粒子初始化,从
class Czastka
中初始化16个粒子及其相应的初始条件。
simulate()
函数会迭代 num_step = 1000
遍历粒子列表,并计算由于与其他粒子相互作用而作用在每个粒子上的力的总和。为此,使用 lennard_jones_forces(x1, y1, x2, y2)
函数。所应用的条件是,如果粒子相距超过 2.5*sigma
,则它们不会感受到任何吸引力。
使用函数
leapforg_single_step(particle, force_x_comp, force_y_comp):
计算新的位置和速度。
然后使用
p1.update_position(p1_x_next, p1_y_next, box_size)
将其附加到粒子位置列表中。 这里还检查周期性边界条件。其背后的想法是,粒子是否会离开模拟/包含气体的框架/盒子。从另一面看来,这是通过模除来实现的self.x = new_x % box_size
。
正如我所说,我不明白是什么导致了纳米错误,也不明白为什么我的位置列表中最终大部分为零?这是我的完整代码:
import numpy as np
#Inicjalizacja paramatrów
particle_num, box_size, eps, sigma = 16, 8.0, 1.0, 1.0
dt, temp, k_B, m, M = 0.0001, 2.5, 1, 1, 1
radius = sigma/2
class Czastka:
#Inicjalizacja danych cząstki
def __init__(self, radius, x_pos, y_pos, x_vel, y_vel):
#Składowe wartości wektorowych
self.x, self.y = x_pos, y_pos
self.vx, self.vy = x_vel, y_vel
"""
NOTE: The goal of the simulation in the end is to get only these 4 lists
"""
#Positions lists
self.x_positions_lst = []
self.y_positions_lst = []
#Velocity lists
self.x_velocity_lst = []
self.y_velocity_lst = []
def append_postions(self, x_pos, y_pos):
self.x_positions_lst.append(x_pos)
self.y_positions_lst.append(y_pos)
def update_velocities(self, x_vel, y_vel):
self.x_velocity_lst.append(x_vel)
self.y_velocity_lst.append(y_vel)
# Stosujemy periodyczne warunki brzegowe poprzez liczenie modułowe
def update_position(self, new_x, new_y, box_size):
self.x = new_x % box_size
self.y = new_y % box_size
self.append_postions(self.x, self.y)
# Inicjalizacja cząstek
initial_x, initial_y = [1, 3, 5, 7], [1, 3, 5, 7]
initial_vx, initial_vy = [1, 1, 2, 0.5], [1, 0.5, 0.25, 4]
particle_lst = []
for i in range(0, 4):
for j in range(0, 4):
particle_lst.append(Czastka(radius, initial_x[j], initial_y[j], initial_vx[j], initial_vy[j]))
#print(len(particle_lst))
#Siła jaka działa na skutek odziaływań w potencjale Lennard-Jonesa
def lennard_jones_forces(x1, y1, x2, y2):
global sigma, eps
#Obliczenia pomocnicze
r12 = [x2-x1, y2-y1]
dist12 = (r12[0]**2 + r12[1]**2)**0.5
print(r12[0], ' ', r12[1])
if dist12 == 0:
return 0, 0
# Calculate Lennard-Jones force
sigma_over_dist12 = sigma / dist12
sigma_over_dist12_14 = sigma_over_dist12 ** 14
sigma_over_dist12_8 = sigma_over_dist12 ** 8
Force21 = -(48 * eps / (sigma ** 2)) * (sigma_over_dist12_14 - 0.5 * sigma_over_dist12_8)
Force21_x = Force21 * r12[0]
Force21_y = Force21 * r12[1]
#W tym momecie nasza funkcja jest gotowa ALE musimy sprawdzić czy zachodzi warunek odcięcia
#Żeby zwiększyć wydajnośc obliczniową powinno się dokonać tego sprwdzenia, przed obliczniem sił
if np.isnan(Force21_x) or np.isnan(Force21_y):
print(Force21_x, ' ', Force21_y)
print("Nan detected in force calculation")
if dist12 > (2.5*sigma):
#print('Cut off')
Force21_x, Force21_y = 0, 0
return Force21_x, Force21_y
else:
#print('Normal operation')
return Force21_x, Force21_y
# Obliczanie poprzez wykorzystanie algorytmu żabki nowych współrzędnych
def leapforg_single_step(particle, force_x_comp, force_y_comp):
global dt, m
#Obliczanie pół-krokowych prędkości
vx_half = particle.vx + 0.5 * (force_x_comp / m) * dt
vy_half = particle.vy + 0.5 * (force_y_comp / m) * dt
#Obliczanie pół-krokowych położeń
x_next = particle.x + vx_half * dt
y_next = particle.y + vy_half * dt
#Obliczanie nowych prędkości
vx_next = vx_half + 0.5 * (force_x_comp / m) * dt
vy_next = vy_half + 0.5 * (force_y_comp / m) * dt
return x_next, y_next, vx_next, vy_next
def simulate():
global box_size
num_steps = 1000
#Cała symulacja składa się z num_steps kroków
for step in range(num_steps):
#Obliczmy sumę sił działającą na każdą cząstkę
for i, p1 in enumerate(particle_lst):
force_x, force_y = 0, 0
for j, p2 in enumerate(particle_lst):
if i != j:
# Obliczam sume sił działająca na cząsteki p1 (bez interakcji z samym sobą)
f_x, f_y = lennard_jones_forces(p1.x, p1.y, p2.x, p2.y)
force_x += f_x
force_y += f_y
p1_x_next, p1_y_next, p1_vx_next, p1_vy_next = leapforg_single_step(p1, force_x, force_y)
p1.update_position(p1_x_next, p1_y_next, box_size)
p1.update_velocities(p1_vx_next, p1_vy_next)
simulate()
def simulation_results():
for i, particle in enumerate(particle_lst):
print('Particle number: {}. List lenght: {} and x-positions list: {}'.format(i, len(particle.x_positions_lst), particle.x_positions_lst))
simulation_results()
def animation():
pass
任何帮助将不胜感激!干杯!
我进行了一些轻微的重构,但无法重现该问题。即使我通过引发异常来检测相同的噪声,也不会出现 NaN。这里还有其他需要改进的地方(你根本没有正确使用 Numpy)。
import numpy as np
# Inicjalizacja paramatrów
particle_num, box_size, eps, sigma = 16, 8.0, 1.0, 1.0
dt, temp, k_B, m, M = 0.0001, 2.5, 1, 1, 1
radius = sigma/2
class Particle:
# Inicjalizacja danych cząstki
def __init__(self, radius: float, x_pos: float, y_pos: float, x_vel: float, y_vel: float) -> None:
# Składowe wartości wektorowych
self.x, self.y = x_pos, y_pos
self.vx, self.vy = x_vel, y_vel
"""
NOTE: The goal of the simulation in the end is to get only these 4 lists
"""
# Positions lists
self.x_positions_lst: list[float] = []
self.y_positions_lst: list[float] = []
# Velocity lists
self.x_velocity_lst: list[float] = []
self.y_velocity_lst: list[float] = []
def append_postions(self, x_pos: float, y_pos: float) -> None:
self.x_positions_lst.append(x_pos)
self.y_positions_lst.append(y_pos)
def update_velocities(self, x_vel: float, y_vel: float) -> None:
self.x_velocity_lst.append(x_vel)
self.y_velocity_lst.append(y_vel)
# Stosujemy periodyczne warunki brzegowe poprzez liczenie modułowe
def update_position(self, new_x: float, new_y: float, box_size: float) -> None:
self.x = new_x % box_size
self.y = new_y % box_size
self.append_postions(self.x, self.y)
# Obliczanie poprzez wykorzystanie algorytmu żabki nowych współrzędnych
def leapfrog_single_step(self, force_x_comp: float, force_y_comp: float) -> tuple[
float, # x next
float, # y next
float, # vx next
float, # vy next
]:
# Obliczanie pół-krokowych prędkości
vx_half = self.vx + 0.5 * (force_x_comp / m) * dt
vy_half = self.vy + 0.5 * (force_y_comp / m) * dt
# Obliczanie pół-krokowych położeń
x_next = self.x + vx_half * dt
y_next = self.y + vy_half * dt
# Obliczanie nowych prędkości
vx_next = vx_half + 0.5 * (force_x_comp / m) * dt
vy_next = vy_half + 0.5 * (force_y_comp / m) * dt
return x_next, y_next, vx_next, vy_next
# Siła jaka działa na skutek odziaływań w potencjale Lennard-Jonesa
def lennard_jones_forces(
x1: float, y1: float, x2: float, y2: float,
) -> tuple[float, float]:
# Obliczenia pomocnicze
r12 = [x2-x1, y2-y1]
dist12 = (r12[0]**2 + r12[1]**2)**0.5
# print(r12[0], ' ', r12[1])
if dist12 == 0:
return 0, 0
# Calculate Lennard-Jones force
sigma_over_dist12 = sigma / dist12
sigma_over_dist12_14 = sigma_over_dist12 ** 14
sigma_over_dist12_8 = sigma_over_dist12 ** 8
Force21 = -(48 * eps / (sigma ** 2)) * (sigma_over_dist12_14 - 0.5 * sigma_over_dist12_8)
Force21_x = Force21 * r12[0]
Force21_y = Force21 * r12[1]
# W tym momecie nasza funkcja jest gotowa ALE musimy sprawdzić czy zachodzi warunek odcięcia
# Żeby zwiększyć wydajnośc obliczniową powinno się dokonać tego sprwdzenia, przed obliczniem sił
if np.isnan(Force21_x) or np.isnan(Force21_y):
raise ValueError(f"Nan detected in force calculation: {Force21_x} {Force21_y}")
if dist12 > (2.5*sigma):
# print('Cut off')
Force21_x, Force21_y = 0, 0
return Force21_x, Force21_y
else:
# print('Normal operation')
return Force21_x, Force21_y
def simulate(particle_lst: list[Particle]) -> None:
num_steps = 1000
# Cała symulacja składa się z num_steps kroków
for step in range(num_steps):
# Obliczmy sumę sił działającą na każdą cząstkę
for i, p1 in enumerate(particle_lst):
force_x, force_y = 0, 0
for j, p2 in enumerate(particle_lst):
if i != j:
# Obliczam sume sił działająca na cząsteki p1 (bez interakcji z samym sobą)
f_x, f_y = lennard_jones_forces(p1.x, p1.y, p2.x, p2.y)
force_x += f_x
force_y += f_y
p1_x_next, p1_y_next, p1_vx_next, p1_vy_next = p1.leapfrog_single_step(force_x, force_y)
p1.update_position(p1_x_next, p1_y_next, box_size)
p1.update_velocities(p1_vx_next, p1_vy_next)
def simulation_results(particle_lst: list[Particle]) -> None:
for i, particle in enumerate(particle_lst):
print(
f'Particle number: {i}. '
f'List length: {len(particle.x_positions_lst)} and '
f'x-positions list: {particle.x_positions_lst}'
)
def main() -> None:
# Inicjalizacja cząstek
initial_x, initial_y = [1, 3, 5, 7], [1, 3, 5, 7]
initial_vx, initial_vy = [1, 1, 2, 0.5], [1, 0.5, 0.25, 4]
particle_lst = []
for i in range(0, 4):
for j in range(0, 4):
particle_lst.append(Particle(radius, initial_x[j], initial_y[j], initial_vx[j], initial_vy[j]))
# print(len(particle_lst))
simulate(particle_lst)
simulation_results(particle_lst)
if __name__ == '__main__':
main()