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

def dist(a1, a2):
    return np.linalg.norm(a1 - a2)

# Потенциал Леннарда-Джонса
def u(r):
    return 4 * eps * ((sigma / r)**12 - (sigma / r)**6)

def forces(coords, n_atoms):
    forces = np.zeros((n_atoms, 3))  # Трехмерные силы
    for i in range(n_atoms):
        for j in range(i + 1, n_atoms):
            r_ij = coords[i] - coords[j]
            r = np.linalg.norm(r_ij)
            if r > 0:
                force_magnitude = 24 * eps * (2 * (sigma / r)**12 - (sigma / r)**6) / r**2
                force_vector = force_magnitude * (r_ij / r)
                forces[i] += force_vector
                forces[j] -= force_vector
    return forces

def u_total(coords, n_atoms):
    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)
    return energy

mass = 1.0  # Масса атомов (в атомных единицах)
sigma = 1.0  # Параметр Леннарда-Джонса (в атомных единицах)
eps = 1.0  # Энергия (в атомных единицах)
t_max = 5000  # Максимальное количество итераций
time_step = 0.1  # Шаг по времени
energy_threshold = 1e-6  # Порог точности

coords = np.array([
    [-1.48299065476456, 1.76531716394081, -0.77268972225495],
    [-0.18040290485171, 2.22965591428082, -0.68325485345313],
    [0.87213948582557, 1.62148574017750, -1.39768218615107],
    [-1.73409785533136, 0.55308539835510, -1.46032742961606],
    [-0.62036597303262, -0.25804296469539, -1.98526599217501],
    [0.68791109855286, 0.29037837534202, -1.99953791196672],
    [-2.08660595759664, -1.67804437554509, -0.61980080272662],
    [-2.56397090983759, -0.26576127263212, -0.64111707410608],
    [2.01430494571546, 1.56508244846803, -0.53413195693634],
    [2.62058579236347, 0.31758429691267, -0.64281307312575],
    [1.79935952181604, -0.51512725294138, -1.54766425547717],
    [1.52563818697286, 1.74518052646804, 0.86856787583649],
    [0.21643576275457, 2.23451897453078, 0.77947013186982],
    [-0.81895028090721, 1.59624333112993, 1.46582749980378],
    [1.79269417591429, 0.55169217009261, 1.57269064998931],
    [0.68392374627209, -0.25912322438114, 2.02163795286504],
    [-0.62691898531055, 0.28303284007451, 2.00390882594822],
    [2.62028424421333, -0.27707490298611, 0.67002092800882],
    [-2.56699709370779, 0.28146380491787, 0.65464481137833],
    [-2.09632847490229, 1.69601546822453, 0.63457253851745],
    [-1.73536579134934, -0.53340689222586, 1.47603574077613],
    [0.87610268578258, -1.58932681808972, 1.42027561805953],
    [2.02025369584686, -1.52745035347449, 0.55974061524057],
    [-0.80752933592474, -1.57216182079846, -1.44768687236054],
    [0.22908419071840, -2.20547964856536, -0.75859445604328],
    [1.53615817527116, -1.70988038155071, -0.84423797020701],
    [-0.17162935410286, -2.20251827608069, 0.70307678250976],
    [-1.47665469432564, -1.74442130129554, 0.78907339604319],
])

# Количество атомов
n_atoms = len(coords)

# Начальные импульсы (нулевые для всех атомов)
p = np.zeros_like(coords)

# Массивы для графиков
kinetic_energy_history = []
potential_energy_history = []

# Сохранение координат для анимации
coords_history = []

prev_kinetic_energy = None
for n in range(t_max):
    # Сохранение текущих координат
    coords_history.append(coords.copy())

    # 3. Вычисляем полную потенциальную энергию и силы
    forces_vals = forces(coords, n_atoms)
    energy_potential = u_total(coords, n_atoms)

    # 4. Обновление импульсов по формуле leap-frog
    p_half_step = p + 0.5 * time_step * forces_vals  # p_i(t + 0.5)

    # 5. Обновление координат по формуле
    coords = coords + time_step * p_half_step / mass  # r_i(t + 1)

    # 6. Пересчёт сил на новых координатах
    forces_vals = forces(coords, n_atoms)

    # 7. Обновление импульсов на следующий шаг
    p = p_half_step + 0.5 * time_step * forces_vals  # p_i(t + 1)

    # 8. Вычисление кинетической энергии
    energy_kinetic = 0.5 * np.sum(p**2) / mass
    kinetic_energy_history.append(energy_kinetic)
    potential_energy_history.append(energy_potential)

    # 9. Проверка условия уменьшения кинетической энергии
    if prev_kinetic_energy is not None and energy_kinetic < prev_kinetic_energy:
        # 10. Проверка, если кинетическая энергия меньше порога
        if energy_kinetic < energy_threshold:
            print(f"Алгоритм завершен на итерации {n + 1} (Кинетическая энергия < {energy_threshold})")
            break
        else:
            # Если кинетическая энергия уменьшилась, но всё ещё больше порога, обнуляем импульсы
            p = np.zeros_like(coords)
            prev_kinetic_energy = None
            continue  # Переход к следующей итерации

    # Если кинетическая энергия не уменьшилась, продолжаем с шага 4
    prev_kinetic_energy = energy_kinetic

    # 10. Условие завершения, если кинетическая энергия меньше порога
    if energy_kinetic < energy_threshold:
        print(f"Алгоритм завершен на итерации {n + 1} (Кинетическая энергия < {energy_threshold})")
        break
else:
    print(f"Алгоритм завершен после максимального количества итераций: {t_max}")

print(f"Итого итераций: {n + 1}")
print("\nКонечные координаты атомов:")
print(coords)
print(f"Кинетическая энергия (K): {energy_kinetic:.6f}")
print(f"Потенциальная энергия (U): {energy_potential:.6f}")

# Анимация движения атомов
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='sandybrown')

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()

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

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

# График кинетической энергии
for i in range(1, len(kinetic_energy_history)):
    if kinetic_energy_history[i] < 0.02:
        ax[0].plot(time[i-1:i+1], kinetic_energy_history[i-1:i+1], linestyle='--', color='steelblue')
    else:
        ax[0].plot(time[i-1:i+1], kinetic_energy_history[i-1:i+1], color='steelblue')

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

plt.tight_layout()
plt.show()

здесь нужно задать параметры для атомов углерода в атомных единицах хартри где eps = 34.839 K, sigma = 3.805 ангстрема