在数值积分过程中，两个物体之间的距离保持不变

Question

我正试图为我的研究项目运行一个数值积分代码。它是一个三原子系统，只受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)

Answer 1

你的代码有一些细节，它没有运行的主要原因是因为这一行return np.sum(xAcc), np.sum(yAcc)，你从所有粒子之间的相互作用计算了一个加速度矩阵，其中你添加了一个粒子的加速度和另一个粒子的加速度的倒数，因为所有粒子的质量都一样，那么加速度也一样，那么你把矩阵的所有元素求和，所有的项都抵消掉了，而不是返回每个粒子的加速度，你返回所有粒子的加速度之和，假设它们都有相同的质量，那么它们的加速度之和就是0，即使它不是0，也可能是错误的，因为你把它们都混合在一起，所以所有的粒子都会有相同的加速度，朝着相同的方向运动

关于代码的详细信息

• 一,

在同样的功能中，你也有一些可以改进的东西，比如

if condition:
    A_statement
    B_statement
else:
    B_statement

这与

if condition:
    A_statement
B_statement

因为B_statement总是会被执行

• 二,

您不必在使用它的变量上使用global，在使用它们的地方它们仍然在范围内，基本上在您将在其他函数中使用该变量时使用它，基本上所有具有相同数量表格或更多表格的对象都知道该变量，这在第一行结束，比那一行少列表，这导致了下一点

• 三,

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，该变量包含值1

• 四,

    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)

其中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)

我不知道分子动力学的物理原理是否正确，当我在物理系时，我只对动力学系统进行了模拟，也许我的方法对你的系统来说过于精确或经典，你可以用它来比较