В вашем коде есть некоторые детали, основная причина того, что он не работает, заключается в том, что из-за этой строки 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
Нет необходимости делать эти два идентичные функции, единственная разница - это имя, они делают одно и то же, имя параметров функции go выходит за рамки, когда вы завершаете тело функции, поэтому вы можете повторно использовать единственную функцию для x и y .
Некоторые примеры
def fun(a):
return a
x = a # NameError: name 'a' is not defined
a не входит в область видимости, когда мы используем его для x, он существовал только в теле fun, поэтому вы можете повторно использовать это имя для другая переменная или параметр.
* 103 7 *
Это происходит потому, что a затеняется параметром функции с тем же именем, а
a = 1
def fun():
return a
print( fun() ) # = 1
потому что, когда fun() ищет a, он не находит его в своей области, поэтому смотрит в более высокую область и находит переменную 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)
Я не знаю, точно ли физика должна быть такой, какой она должна быть для молекулярной динамики, когда Я занимался физикой, я моделировал только динамические системы, может быть, мой подход слишком точный или классический для вашей системы, в любом случае вы можете использовать его для сравнения.