import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

global t
t = 0

def u(r):
    return eps * ((sigma / r)**12 - 2*(sigma / r)**6)

def u_total(coords_flatten):
    # Преобразуем одномерный массив обратно в двумерный массив
    coords = coords_flatten.reshape((n_atoms, 3))
    global t
    t += 1
    if t % 58 == 0: coords_history.append(coords.copy())
    energy = 0.0
    for i in range(n_atoms):
        for j in range(i + 1, n_atoms):
            r = np.linalg.norm(coords[i] - coords[j])
            energy += u(r)
    potential_energy_history.append(energy)        
    return energy

coords = np.array([
    [1.99992272624515, -1.55394839253328, -0.54515336499829],
    [2.63668346641465, -0.31510893362187, -0.63469640100073],
    [2.65077799853344, 0.24711589083848, 0.70558024763462],
    [1.52722605028587, -1.76549882311611, 0.84939064982003],
    [1.81217410578072, -0.58406639378151, 1.55983167952112],
    [0.20763155899675, -2.23389038812251, 0.76039088898736],
    [0.84979137822815, -1.60653773470833, -1.40820126560287],
    [-0.19246482654185, -2.21844920965800, -0.69281138810739],
    [1.81004579987672, 0.53684135146920, -1.48007737825631],
    [0.66348262873344, -0.26502813547325, -1.98256865641318],
    [2.04428642916793, 1.50137798643177, 0.62289484468002],
    [1.56225208385088, 1.72510934816412, -0.76668369031050],
    [0.69140583701290, 0.24637939208619, 2.07449375339418],
    [0.90493937459836, 1.58260934225070, 1.49831085368104],
    [-0.81893681672243, -1.60610643960112, 1.48460559894418],
    [-0.59925876442654, -0.27789435722249, 2.07781201717052],
    [-1.49984196042733, -1.71605327994024, -0.77518957869575],
    [-1.96751690853535, -1.51248142516120, 0.62261394576123],
    [-0.61357889319083, 0.29137118601298, -1.97173474502841],
    [-1.75518413772831, -0.51796557767552, -1.46873904841851],
    [0.25582327825727, 2.22602808710782, -0.66326891570148],
    [-0.12918566902651, 2.22047078759246, 0.79384882238460],
    [-0.79385069614173, 1.62434020639910, -1.37670854833150],
    [-1.72535221760213, 0.55951183773771, 1.58695819504051],
    [-2.57276587823076, -0.25948658554608, 0.72961372876720],
    [-2.57256736101931, 0.32189901908437, -0.60263318803483],
    [-1.44780510967951, 1.75119436529804, 0.89033622781008],
    [-1.93513347670964, 1.55926687568858, -0.50221528469694]
])

n_atoms = len(coords)
eps = 6.1628 / 27.212
sigma = 1.2431 / 0.529177
xtol = 1.0e-5

coords_history = []
coords_history.append(coords.copy())
potential_energy_history = []

# Преобразуем начальные координаты в одномерный массив
coords_flatten = coords.ravel()

# Оптимизация
result_minimize = opt.minimize(u_total, coords_flatten, method='Nelder-Mead', options={'xatol': xtol, 'disp': True})

# Преобразуем оптимизированные координаты обратно в двумерный массив
optimized_coords = result_minimize.x.reshape((n_atoms, 3))

print("\nКонечные координаты атомов:")
print(optimized_coords)
print(f"Потенциальная энергия (U): {result_minimize.fun:.6f}")

# Построение графика потенциальной энергии
time = np.arange(len(potential_energy_history))

fig, ax = plt.subplots(1, 1, figsize=(8, 6))

# График потенциальной энергии
ax.plot(time, potential_energy_history, label="Потенциальная энергия", color='palevioletred')
ax.set_xlabel('Время')
ax.set_ylabel('Потенциальная энергия')
ax.set_title('Потенциальная энергия от времени')

plt.tight_layout()
plt.show()

# Анимация движения атомов
coords_history = np.array(coords_history)

fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
atoms, = ax.plot([], [], 'o', color='steelblue')

def update(frame):
    positions = coords_history[frame]
    atoms.set_data(positions[:, 0], positions[:, 1])
    return atoms,

ani = FuncAnimation(fig, update, frames=len(coords_history), interval=50, blit=True)

# Сохранение анимации в GIF
ani.save("lennard_jones_simulation.gif", writer="pillow", fps=20)

plt.show()

в эту строчку if t % 58 == 0: coords_history.append(coords.copy()) необходимо добавить сохранение координат в текстовые файлы, предположим имеется папка coords_history, туда сохраняем текстовые файлы с координатами, отдельный текстовый файл на каждую итерацию, файлы должны быть пронумерованы