Verlet Алгоритм в Python

Question

Я следую инструкциям, которые я нашел по этой ссылке ( http://gdrcorelec.ups-tlse.fr/files/python_verlet.pdf). Все работает отлично (потенциал Леннарда-Джонса, молекулярная геометрия, вычисление матрицы расстояний, потенциал для нескольких атомов, полная энергия), но я не могу заставить части кода ускорения и алгоритма верлета работать должным образом. Как я могу решить эту проблему? Благодарю.

Это параметры, которые я печатаю, когда спрашивается в коде.

Введите межатомное расстояние: 0,3 Количество атомов: 3 Масса атома: 10 Масса атома: 20 Масса атома: 15 Положение x: 0 Положение y: 0 Положение z: 0 Положение x: 0 Положение y: 0 Положение z: 0.3 Положение x: 0.1 Положение y: 0.2 Положение z: -0,3 Скорость атома: 0,1 Скорость атома: 0,2 Скорость атома: 0,3

КОД

import numpy as np

#Parâmetros iniciais
mass = 39.948
E = 0.0661
S = 0.3345

r = float(input("Enter inter-atomic distance: "))

#Potencial de Lennard-Jones
V = 4*E*((S/r)**12-(S/r)**6)
print(V)

#Número de átomos
Natoms = int(input("Number of atoms: "))

#Massa dos átomos
massa = []
t1 = 0
for i in range(Natoms):
    while t1 < Natoms:
            x1 = float(input("Mass of atom:"))
            massa = np.append(massa,x1)
            t1 += 1

print("Massa =", massa)

#Coordenadas iniciais dos átomos
pos = []
for i in range(Natoms):
    x1 = float(input("Position x:"))
    x2 = float(input("Position y:"))
    x3 = float(input("Position z:"))
    pos1 = np.array([x1, x2, x3])
    pos = np.append(pos, pos1, axis=0)

pos = np.array(pos).reshape((Natoms,3))
print("Posição =", pos)

dist1 = []
for mylist in pos:
    for item in mylist:
        dist1 = np.append(dist1, item)

print("D1 =", dist1)


dist = [[0] * Natoms for i in range(Natoms)]
for i in range(Natoms):
    for j in range(Natoms):
        if i < j:
            dist[i][j] = np.linalg.norm(pos[i]-pos[j])
        elif i > j:
            dist[i][j] = np.linalg.norm(pos[i]-pos[j])
        else:
            dist[i][j] = np.linalg.norm(pos[i]-pos[j])

print("Matriz de distâncias entre as partículas = ", dist)

dist2 = []
for mylist in dist:
    for item in mylist:
        dist2 = np.append(dist2, item)

print("D2 =", dist2)

dist3 = []
for num in dist2:
    if num not in dist3:
        dist3.append(num)

print("D3 = ", dist3)

#Remover as distâncias calculadas para as mesmas partículas
for x in dist3:
    if x==0.0:
        dist3.remove(x)

print("Distâncias = ", dist3)

#Cálculo do potencial total do sistema
V01 = []
for i in dist3:
    V02 = 4*E*((S/i)**12-(S/i)**6)
    V01 = np.append(V01,V02)

#print(V01)

Vr = sum(V01)
print("Vr = ", Vr)

#Cálculo da energia total
veloc = []
t2 = 0
for i in range(Natoms):
    while t2 < Natoms:
        y1 = float(input("Veloc of atom:"))
        veloc = np.append(veloc,y1)
        t2 += 1

print("Velocidades = ", veloc)

T01 =[]
for i in veloc:
    for j in massa:
        T02 = (1/2)*(j*i**2)
        T01 = np.append(T01, T02)

print(T01)

Ttot = sum(T01)
print("Ttot = ", Ttot)

E = Ttot + Vr

print("Energia = ", E)

#Cálculo da aceleração
dpos = 1
dist4 = dist2 + dpos
dist5 = dist2 - dpos
print("D4 = ", dist4)
print("D5 = ", dist5)

V04 = []
for i in dist4:
    V041 = 4*E*((S/i)**12-(S/i)**6)
    V04 = np.append(V04,V041)    

print("V04 = ", V04)

V05 = []        
for j in dist5:
    V051 = 4*E*((S/j)**12-(S/j)**6)
    V05 = np.append(V05,V051)

print("V05 = ", V05)

z = [x-y for x,y in list(zip(V04,V05)) ]
print("z=",z)

#dvdpos=[]
s = 2*dpos
dvdpos =[x/s for x in z]
#print(dvdpos)

dvdpos = np.array(dvdpos).reshape((Natoms,3))
print("dvdpos = ", dvdpos)


for i in massa:
    a = [-1/i*x for x in dvdpos]
    a = np.array(a).reshape((Natoms,3))
print("a=",a)

a2 = []
for mylist in a:
    for item in mylist:
        a2 = np.append(a2, item)
 print("a2 =", a2)

#Algoritmo de Verlet
Nsteps = 10
dt = 0.1
r_new = []
v_new = []
veloc2 = np.tile(veloc,Natoms)
print("V2=", veloc2)

t3 = 0
while t3 < Nsteps:
    r_new2 = dist2 + veloc2*dt + a2*dt**2/2
    r_new = np.append(r_new,r_new2)
    t3 += 1
    for j in massa:
    a3 = [-1/j*x for x in r_new]
    v_new2 = veloc2 + (a2 + a3)/2 * dt
    v_new = np.append(v_new,v_new2)