在数值积分过程中两个物体之间的距离保持不变
Distance between two bodies remains constant during numerical integration process
我正在尝试 运行 我的研究项目的数值积分代码。它是一个 3 原子系统，只受到 Lennard-Jones 力。但是 r_x 变量在此过程中保持为 0。不幸的是我不知道为什么。这是代码的输出：
[[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]]
当我检查所有变量的值时，我看到 r_x 只有一个值，并且在此过程中为零。
import numpy as np
np.seterr(invalid = "ignore")
m = 1
x = np.array([1, 9, 15])
y = np.array([16, 22, 26])
def GetLJForce(r, epsilon, sigma):
return 48 * epsilon * np.power(sigma, 12) / np.power(r, 13) - 24 * epsilon * np.power(sigma, 6) / np.power(r, 7)
def GetAcc(xPositions, yPositions):
global xAcc
global yAcc
xAcc = np.zeros((xPositions.size, xPositions.size), dtype=object)
yAcc = np.zeros((xPositions.size, xPositions.size), dtype=object)
for i in range(0, xPositions.shape[0]-1):
for j in range(i+1, xPositions.shape[0]-1):
global r_x
r_x = xPositions[j] - xPositions[i]
r_y = yPositions[j] - yPositions[i]
global rmag
rmag = np.sqrt(r_x*r_x + r_y*r_y)
if(rmag[0]==0 or rmag[1]==0 or rmag[2]==0):
rmag += 1
force_scalar = GetLJForce(rmag, 0.84, 2.56)
force_x = force_scalar * r_x / rmag
force_y = force_scalar * r_y / rmag
xAcc[i,j] = force_x / m
xAcc[j,i] = - force_x / m
yAcc[i,j] = force_y / m
yAcc[j,i] = - force_y / m
else:
force_scalar = GetLJForce(rmag, 0.84, 2.56)
force_x = force_scalar * r_x / rmag
force_y = force_scalar * r_y / rmag
xAcc[i,j] = force_x / m
xAcc[j,i] = - force_x / m
yAcc[i,j] = force_y / m
yAcc[j,i] = - force_y / m
return np.sum(xAcc), np.sum(yAcc)
def UpdatexPos(x, v_x, a_x, dt):
return x + v_x*dt + 0.5*a_x*dt*dt
def UpdateyPos(y, v_y, a_y, dt):
return y + v_y*dt + 0.5*a_y*dt*dt
def UpdatexVel(v_x, a_x, a1_x, dt):
return v_x + 0.5*(a_x + a1_x)*dt
def UpdateyVel(v_y, a_y, a1_y, dt):
return v_y + 0.5*(a_y + a1_y)*dt
def RunMD(dt, number_of_steps, x, y):
global xPositions
global yPositions
xPositions = np.zeros((number_of_steps, 3))
yPositions = np.zeros((number_of_steps, 3))
v_x = 0
v_y = 0
a_x = GetAcc(xPositions, yPositions)[0]
a_y = GetAcc(xPositions, yPositions)[1]
for i in range(number_of_steps):
x = UpdatexPos(x, v_x, a_x, dt)
y = UpdateyPos(y, v_y, a_y, dt)
a1_x = GetAcc(xPositions, yPositions)[0]
a1_y = GetAcc(xPositions, yPositions)[1]
v_x = UpdatexVel(v_x, a_x, a1_x, dt)
v_y = UpdateyVel(v_y, a_y, a1_y, dt)
a_x = np.array(a1_x)
a_y = np.array(a1_y)
xPositions[i, :] = x
yPositions[i, :] = y
return xPositions, yPositions
sim_xpos = RunMD(0.1, 10, x, y)[0]
sim_ypos = RunMD(0.1, 10, x, y)[1]
print(sim_xpos)
您的代码有一些细节，它不是 运行 的主要原因是因为这一行 return np.sum(xAcc), np.sum(yAcc)，您根据所有粒子之间的相互作用计算了加速度矩阵将一个粒子的加速度和该加速度的倒数加到另一个粒子上，因为所有粒子的质量都相同，所以加速度也相同，然后对矩阵的所有元素求和，所以所有项都抵消了，而不是 return 每个粒子的加速度你 return 所有粒子的加速度的总和，并且假设它们都具有相同的质量，即使它不是0 它会是错误的，因为你将它们全部混合在一起，所以所有粒子都将具有相同的加速度并沿相同的方向移动。
纯粹关于代码的细节
- 1
在同一个函数中，您还有一些可以改进的地方，例如
if condition:
A_statement
B_statement
else:
B_statement
与
相同
if condition:
A_statement
B_statement
因为 B_statement 总是会被执行。
- 2
你不必在你使用它的变量上使用global，你使用它们的地方它们仍然在范围内，基本上当你要在其他函数中使用该变量时使用它，基本上所有具有相同数量或更多列表的东西都知道该变量，并且在第一行结束的列表比那个少，这导致下一个点。
- 3
def UpdatexVel(v_x, a_x, a1_x, dt):
return v_x + 0.5*(a_x + a1_x)*dt
def UpdateyVel(v_y, a_y, a1_y, dt):
return v_y + 0.5*(a_y + a1_y)*dt
这两个完全相同的函数不需要做，唯一不同的是名字，它们做的是一样的，当你结束函数体时，函数的参数名称超出范围，所以你可以为 x 和 y 重用一个函数。
一些例子
def fun(a):
return a
x = a # NameError: name 'a' is not defined
a 在我们将它用于 x 时不在范围内，它仅存在于 fun 的正文中，因此您可以将该名称重用于其他变量或参数。
a = 1
def fun(a):
return 2 + a
print( fun(7) ) # = 9
发生这种情况是因为 a 被同名函数参数遮蔽，而
a = 1
def fun():
return a
print( fun() ) # = 1
因为当 fun() 查找 a 时它没有在其作用域内找到它，所以它查找更高的作用域，并找到一个变量 a 保存值 1
- 4
return xPositions, yPositions
sim_xpos = RunMD(0.1, 10, x, y)[0]
sim_ypos = RunMD(0.1, 10, x, y)[1]
这个比较小，但是 python 有 tuples 所以这个重复计算 RunMD （你不想要）的代码可以简化为
return xPositions, yPositions
sim_xpos, sim_ypos = RunMD(0.1, 10, x, y)
其中 return 由 RunMD 编辑的对被分配给对（python 的两个元素的元组）sim_xpos 和 sim_ypos.
如果我在哪里重写部分代码，我会删除 numpy 来测试算法，然后使用 numpy 来向量化操作并使代码更高效，它看起来像这样
import math
def cart_to_polar(x, y):
rho = math.sqrt(x**2 + y**2)
phi = math.atan2(y, x)
return rho, phi
def polar_to_cart(rho, theta):
x = rho * math.cos(theta)
y = rho * math.sin(theta)
return x, y
def GetLJForce(r, epsilon, sigma):
return 48 * r * epsilon * ( ( ( sigma / r ) ** 14 ) - 0.5 * ( ( sigma / r ) ** 8 ) ) / ( sigma ** 2 )
def GetAcc(x_positions, y_positions, m):
xAcc = [0]*len(x_positions)
yAcc = [0]*len(x_positions)
for i in range(0, len(x_positions)-1):
for j in range(i+1, len(x_positions)):
# calculations are made from the point of view of i, and then flipped for j
delta_x = x_positions[j] - x_positions[i]
delta_y = y_positions[j] - y_positions[i]
radius, theta = cart_to_polar(delta_x, delta_y)
# in case two particles are at the same place
# this is some times called a cutoff radius
if radius==0: radius = 1e-10
force_mag = GetLJForce(radius, 0.84, 2.56)
# since the polar coordinates are centered in the i particle the result is readilly usable
# particle j interaction with particle i
force_x, force_y = polar_to_cart(force_mag, theta)
xAcc[i] += force_x / m[i]
yAcc[i] += force_y / m[i]
# flip the sing of the force to represent the
# particle i interaction with particle j
force_x, force_y = polar_to_cart(-force_mag, theta)
xAcc[j] += force_x / m[j]
yAcc[j] += force_y / m[j]
return xAcc, yAcc
def update_pos(x, v, a, dt):
return x + v*dt + 0.5 * a * dt ** 2
def update_vel(v, a, dt):
return v + a * dt
def runMD(dt, x, y, v_x, v_y, m):
num_particles = len(x)
a_x, a_y = GetAcc(x, y, m)
for i in range(num_particles):
v_x[i] = update_vel(v_x[i], a_x[i], dt)
v_y[i] = update_vel(v_y[i], a_y[i], dt)
for i in range(num_particles):
x[i] = update_pos(x[i], v_x[i], a_x[i], dt)
y[i] = update_pos(y[i], v_y[i], a_y[i], dt)
return x, y
# number of particles in the system
num_particles = 3
# mass of the particles
m = [1] * num_particles
# starting positions
x = [1, 9, 15]
y = [16, 22, 26]
# starting velocities
v_x = [0] * num_particles
v_y = [0] * num_particles
# number of steps of the simulation
number_of_steps = 10
# the simulation
for i in range(number_of_steps):
x, y = runMD(0.1, x, y, v_x, v_y, m)
print(x, y)
我不知道物理学是否正是分子动力学应该如何，当我在物理学的时候我只模拟了动力系统，也许我的方法对于你的系统来说太精确或经典了，无论哪种方式你可以用它来比较。
我正在尝试 运行 我的研究项目的数值积分代码。它是一个 3 原子系统，只受到 Lennard-Jones 力。但是 r_x 变量在此过程中保持为 0。不幸的是我不知道为什么。这是代码的输出：
[[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]
[ 1. 9. 15.]]
当我检查所有变量的值时，我看到 r_x 只有一个值，并且在此过程中为零。
import numpy as np
np.seterr(invalid = "ignore")
m = 1
x = np.array([1, 9, 15])
y = np.array([16, 22, 26])
def GetLJForce(r, epsilon, sigma):
return 48 * epsilon * np.power(sigma, 12) / np.power(r, 13) - 24 * epsilon * np.power(sigma, 6) / np.power(r, 7)
def GetAcc(xPositions, yPositions):
global xAcc
global yAcc
xAcc = np.zeros((xPositions.size, xPositions.size), dtype=object)
yAcc = np.zeros((xPositions.size, xPositions.size), dtype=object)
for i in range(0, xPositions.shape[0]-1):
for j in range(i+1, xPositions.shape[0]-1):
global r_x
r_x = xPositions[j] - xPositions[i]
r_y = yPositions[j] - yPositions[i]
global rmag
rmag = np.sqrt(r_x*r_x + r_y*r_y)
if(rmag[0]==0 or rmag[1]==0 or rmag[2]==0):
rmag += 1
force_scalar = GetLJForce(rmag, 0.84, 2.56)
force_x = force_scalar * r_x / rmag
force_y = force_scalar * r_y / rmag
xAcc[i,j] = force_x / m
xAcc[j,i] = - force_x / m
yAcc[i,j] = force_y / m
yAcc[j,i] = - force_y / m
else:
force_scalar = GetLJForce(rmag, 0.84, 2.56)
force_x = force_scalar * r_x / rmag
force_y = force_scalar * r_y / rmag
xAcc[i,j] = force_x / m
xAcc[j,i] = - force_x / m
yAcc[i,j] = force_y / m
yAcc[j,i] = - force_y / m
return np.sum(xAcc), np.sum(yAcc)
def UpdatexPos(x, v_x, a_x, dt):
return x + v_x*dt + 0.5*a_x*dt*dt
def UpdateyPos(y, v_y, a_y, dt):
return y + v_y*dt + 0.5*a_y*dt*dt
def UpdatexVel(v_x, a_x, a1_x, dt):
return v_x + 0.5*(a_x + a1_x)*dt
def UpdateyVel(v_y, a_y, a1_y, dt):
return v_y + 0.5*(a_y + a1_y)*dt
def RunMD(dt, number_of_steps, x, y):
global xPositions
global yPositions
xPositions = np.zeros((number_of_steps, 3))
yPositions = np.zeros((number_of_steps, 3))
v_x = 0
v_y = 0
a_x = GetAcc(xPositions, yPositions)[0]
a_y = GetAcc(xPositions, yPositions)[1]
for i in range(number_of_steps):
x = UpdatexPos(x, v_x, a_x, dt)
y = UpdateyPos(y, v_y, a_y, dt)
a1_x = GetAcc(xPositions, yPositions)[0]
a1_y = GetAcc(xPositions, yPositions)[1]
v_x = UpdatexVel(v_x, a_x, a1_x, dt)
v_y = UpdateyVel(v_y, a_y, a1_y, dt)
a_x = np.array(a1_x)
a_y = np.array(a1_y)
xPositions[i, :] = x
yPositions[i, :] = y
return xPositions, yPositions
sim_xpos = RunMD(0.1, 10, x, y)[0]
sim_ypos = RunMD(0.1, 10, x, y)[1]
print(sim_xpos)
您的代码有一些细节，它不是 运行 的主要原因是因为这一行 return np.sum(xAcc), np.sum(yAcc)，您根据所有粒子之间的相互作用计算了加速度矩阵将一个粒子的加速度和该加速度的倒数加到另一个粒子上，因为所有粒子的质量都相同，所以加速度也相同，然后对矩阵的所有元素求和，所以所有项都抵消了，而不是 return 每个粒子的加速度你 return 所有粒子的加速度的总和，并且假设它们都具有相同的质量，即使它不是0 它会是错误的，因为你将它们全部混合在一起，所以所有粒子都将具有相同的加速度并沿相同的方向移动。
- 1
在同一个函数中，您还有一些可以改进的地方，例如
if condition:
A_statement
B_statement
else:
B_statement
与
相同if condition:
A_statement
B_statement
因为 B_statement 总是会被执行。
- 2
你不必在你使用它的变量上使用global，你使用它们的地方它们仍然在范围内，基本上当你要在其他函数中使用该变量时使用它，基本上所有具有相同数量或更多列表的东西都知道该变量，并且在第一行结束的列表比那个少，这导致下一个点。
- 3
def UpdatexVel(v_x, a_x, a1_x, dt):
return v_x + 0.5*(a_x + a1_x)*dt
def UpdateyVel(v_y, a_y, a1_y, dt):
return v_y + 0.5*(a_y + a1_y)*dt
这两个完全相同的函数不需要做，唯一不同的是名字，它们做的是一样的，当你结束函数体时，函数的参数名称超出范围，所以你可以为 x 和 y 重用一个函数。
一些例子
def fun(a):
return a
x = a # NameError: name 'a' is not defined
a 在我们将它用于 x 时不在范围内，它仅存在于 fun 的正文中，因此您可以将该名称重用于其他变量或参数。
a = 1
def fun(a):
return 2 + a
print( fun(7) ) # = 9
发生这种情况是因为 a 被同名函数参数遮蔽，而
a = 1
def fun():
return a
print( fun() ) # = 1
因为当 fun() 查找 a 时它没有在其作用域内找到它，所以它查找更高的作用域，并找到一个变量 a 保存值 1
- 4
return xPositions, yPositions
sim_xpos = RunMD(0.1, 10, x, y)[0]
sim_ypos = RunMD(0.1, 10, x, y)[1]
这个比较小，但是 python 有 tuples 所以这个重复计算 RunMD （你不想要）的代码可以简化为
return xPositions, yPositions
sim_xpos, sim_ypos = RunMD(0.1, 10, x, y)
其中 return 由 RunMD 编辑的对被分配给对（python 的两个元素的元组）sim_xpos 和 sim_ypos.
如果我在哪里重写部分代码，我会删除 numpy 来测试算法，然后使用 numpy 来向量化操作并使代码更高效，它看起来像这样
import math
def cart_to_polar(x, y):
rho = math.sqrt(x**2 + y**2)
phi = math.atan2(y, x)
return rho, phi
def polar_to_cart(rho, theta):
x = rho * math.cos(theta)
y = rho * math.sin(theta)
return x, y
def GetLJForce(r, epsilon, sigma):
return 48 * r * epsilon * ( ( ( sigma / r ) ** 14 ) - 0.5 * ( ( sigma / r ) ** 8 ) ) / ( sigma ** 2 )
def GetAcc(x_positions, y_positions, m):
xAcc = [0]*len(x_positions)
yAcc = [0]*len(x_positions)
for i in range(0, len(x_positions)-1):
for j in range(i+1, len(x_positions)):
# calculations are made from the point of view of i, and then flipped for j
delta_x = x_positions[j] - x_positions[i]
delta_y = y_positions[j] - y_positions[i]
radius, theta = cart_to_polar(delta_x, delta_y)
# in case two particles are at the same place
# this is some times called a cutoff radius
if radius==0: radius = 1e-10
force_mag = GetLJForce(radius, 0.84, 2.56)
# since the polar coordinates are centered in the i particle the result is readilly usable
# particle j interaction with particle i
force_x, force_y = polar_to_cart(force_mag, theta)
xAcc[i] += force_x / m[i]
yAcc[i] += force_y / m[i]
# flip the sing of the force to represent the
# particle i interaction with particle j
force_x, force_y = polar_to_cart(-force_mag, theta)
xAcc[j] += force_x / m[j]
yAcc[j] += force_y / m[j]
return xAcc, yAcc
def update_pos(x, v, a, dt):
return x + v*dt + 0.5 * a * dt ** 2
def update_vel(v, a, dt):
return v + a * dt
def runMD(dt, x, y, v_x, v_y, m):
num_particles = len(x)
a_x, a_y = GetAcc(x, y, m)
for i in range(num_particles):
v_x[i] = update_vel(v_x[i], a_x[i], dt)
v_y[i] = update_vel(v_y[i], a_y[i], dt)
for i in range(num_particles):
x[i] = update_pos(x[i], v_x[i], a_x[i], dt)
y[i] = update_pos(y[i], v_y[i], a_y[i], dt)
return x, y
# number of particles in the system
num_particles = 3
# mass of the particles
m = [1] * num_particles
# starting positions
x = [1, 9, 15]
y = [16, 22, 26]
# starting velocities
v_x = [0] * num_particles
v_y = [0] * num_particles
# number of steps of the simulation
number_of_steps = 10
# the simulation
for i in range(number_of_steps):
x, y = runMD(0.1, x, y, v_x, v_y, m)
print(x, y)
我不知道物理学是否正是分子动力学应该如何，当我在物理学的时候我只模拟了动力系统，也许我的方法对于你的系统来说太精确或经典了，无论哪种方式你可以用它来比较。