Потенциал межмолекулярного взаимодействия

Код Приближение к равновесию идеального газа, но есть ошибка, когда атомы газа сближаются их сильно отталкивает. Не могу понять в чём проблема.

Вот весь код

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.patches import Circle

#ПАРАМЕТРЫ МОДЕЛИРОВАНИЯ
num_particles = 80         # Количество частиц
box_size = 60.0            # Размер квадратной области
particle_radius = 0.535 * 3.4  # Радиус частиц (3.4 - коэффициент для перевода в реальные единицы)
dt = 0.025                 # Шаг времени
total_time = 15.0          # Общее время моделирования
speed = 2.5                # Начальная скорость частиц


# Инициализация частиц в узлах сетки
np.random.seed(42)
positions = []
for IY in np.arange(particle_radius, box_size, particle_radius*2 + particle_radius):
    if len(positions) < num_particles:
        for IX in np.arange(particle_radius, box_size, particle_radius*2 + particle_radius):
            if len(positions) < num_particles:
                positions.append([IX, IY])
positions = np.array(positions)

# Начальные скорости (случайные направления)
angles = np.random.uniform(0, 2*np.pi, num_particles)
velocities = np.column_stack((np.cos(angles), np.sin(angles))) * speed

# Параметры потенциала Леннард-Джонса
R = 3.4          # Диаметр частиц
r_cut = 2.5 * R  # Радиус обрезания потенциала
m = 1.0          # Масса частицы (условные единицы)

# Настройка графики
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111)
ax.set_xlim(0, box_size)
ax.set_ylim(0, box_size)
ax.set_aspect('equal')
ax.set_facecolor('honeydew')
ax.set_title("Молекулярная динамика частиц", fontsize=14, pad=20)

# Создание графических объектов частиц
particles = [Circle((0,0), particle_radius, fc='deepskyblue', ec='navy', lw=0.8) 
             for _ in range(num_particles)]
for p in particles:
    ax.add_patch(p)

# Массивы для расчета сил
old_forces = np.zeros((num_particles, 2))

def update(frame):
    global positions, velocities, old_forces
    
    # Расчет межчастичных сил
    forces = np.zeros((num_particles, 2))
    for i in range(num_particles):
        for j in range(i+1, num_particles):
            r_ij = positions[j] - positions[i]
            r = np.linalg.norm(r_ij)
            
            if r < r_cut and r > 0:
                # Потенциал Леннард-Джонса (7-6)
                F_magnitude = 24 * ((R/r)**7) * (2*(R/r)**6 - 1)
                forces[i] += F_magnitude * r_ij / r
                forces[j] -= F_magnitude * r_ij / r
    
    # Интегрирование по алгоритму Velocity Verlet
    positions += velocities * dt + 0.5 * forces/m * dt**2
    velocities += 0.5 * (old_forces + forces)/m * dt
    old_forces = forces.copy()
    
    # Обработка столкновений со стенками
    for i in range(num_particles):
        for dim in [0, 1]:
            if positions[i, dim] < particle_radius:
                velocities[i, dim] = abs(velocities[i, dim])
                positions[i, dim] = particle_radius
            elif positions[i, dim] > box_size - particle_radius:
                velocities[i, dim] = -abs(velocities[i, dim])
                positions[i, dim] = box_size - particle_radius
    
    # Обновление позиций и цвета частиц
    for i, p in enumerate(particles):
        p.center = positions[i]
        speed = np.linalg.norm(velocities[i])
        p.set_facecolor(plt.cm.plasma(speed/4))  # Цвет зависит от скорости
    
    return particles

# Создание и запуск анимации
ani = FuncAnimation(
    fig,
    update,
    frames=int(total_time/dt),
    interval=30,
    blit=True,
    repeat=False
)

plt.show()