Задание по молекулярной динамике

1_Задание.md

Задание по молекулярной динамике

Молекулярная динамика — слишком красивая и наглядная тема, чтобы обойтись скучными задачами на CodeForces. Поэтому мы решили сделать задачи с визуализацией, чтобы вы увидели, как те или иные формулы выглядят "вживую".

Мы подготовили файл-заглушку, в которую вам нужно будет дописать свой код для интегрирования Verlet и для конкретных сил, которые предложены в заданиях ниже.

Установка нужных библиотек

Откройте любым способом командную строку, где доступен Python (например, в VS Code) и выполните команду:

pip install "scipy < 1.6" ratcave six

Если у вас Linux, добавьте опцию --user.

Запустите файл md_task.py, приложенный ниже, и убедитесь, что вы видите медленно ползущий вправо шарик.

Замечание о координатах

Ось z виртуального мира направлена из экрана к зрителю, поэтому объектам надо давать отрицательную z-координату, чтобы они были видны на экране.

Предлагается распологать атомы в плоскости z = -3.

Задача 1 - орбита

Напишите интегрирование по схеме Verlet. Для его проверки создайте атом в точке -1, 0, -3 и начальной скоростью 0, 0.5 ** 0.5, 0. Ускорение будет направлено в "центр" мира в точке origin = 0, 0, -3:

Задача 2 - сила Лоренца

Поместите атом с зарядом 1 в точку -2, 0, -3 с начальной скоростью 1/6, 5/6, 0. Подействуйте на него силой Лоренца:

Где v — скорость атома, а B = [1, 0, 0] — вектор напряженности магнитного поля.

В этой задаче надо использовать Velocity Verlet.

Задача 3 - электростатика

Теперь создайте два атома разного заряда (1 и -1) в точках -1, 0, -3 и 1, 0, -3. Дайте им начальные скорости 0, 1/6, 0 и 0, -1/6, 0.

Ускорение атома i будет зависеть от положения атома j по закону Кулона:

Где q_i и q_j — заряды атомов.

Можете добавить ещё несколько атомов с разным зарядом и посмотреть, как они будут взаимодействовать.

Задача 4 - силы ван-дер-Ваальса

Теперь к закону Кулона добавим ещё силы ван-дер-Ваальса. Положения будут такие же, начальные скорости нулевые, а к ускорению добавится такое слагаемое:

Здесь σ = 0.5, а ε = 0.1.

md_task.py

import numpy as np
import pyglet
import ratcave as rc

obj_reader = rc.WavefrontReader(rc.resources.obj_primitives)


class Atom:
    """
    Класс атома. Имеет следующие поля:

    - pos: положение атома, массив NumPy
    - prev_pos: предудыщее положение, массив NumPy
    - init_speed: начальная скорость, массив NumPy
    - charge: заряд, по умолчанию равен 0

    С массивами NumPy можно делать операции как с числами: складывать, вычитать и умножать на число.

    Длину вектора можно посчитать так:

    np.linalg.norm(v)

    Пример работы с такими векторами -- ниже в функции verlet.
    
    Текущая скорость атома может быть вычислена как разность pos и prev_pos, разделенная на dt.
    """

    def __init__(
        self, position: tuple, init_speed: np.array, charge: float = 0.0
    ) -> None:
        self.mesh: rc.Mesh = obj_reader.get_mesh("Sphere", position=position, scale=0.3)

        tex = rc.Texture.from_image(rc.resources.img_uvgrid)
        self.mesh.textures.append(tex)

        self.prev_pos = None
        self.init_speed = init_speed
        self.charge = charge

    @property
    def pos(self) -> np.array:
        return self.mesh.position._array

    @pos.setter
    def pos(self, v: np.array) -> None:
        self.mesh.position = v


# Шаг интегрирования -- 20 мс.
dt = 0.02

# Список всех атомов с начальными скоростями и зарядами
atoms = [Atom((-1, 0, -3), np.array([1 / 30, 0.0, 0.0]), 0)]


def acceleration(atom: Atom) -> np.array:
    """
    Эта функция должна вычислять ускорение для каждого атома.
    """
    pass


def verlet(*args) -> None:
    """
    Эта функция должна менять положения атомов с помощью интегрирования Verlet.
    """

    for atom in atoms:
        # Присвойте atom.pos новое положение с учетом ускорения.
        atom.pos += atom.init_speed * dt


# Технические детали

window = pyglet.window.Window()

scene = rc.Scene(meshes=[atom.mesh for atom in atoms])
scene.bgColor = 1, 0, 0


@window.event
def on_draw():
    with rc.default_shader:
        scene.draw()


pyglet.clock.schedule_interval(verlet, dt)
pyglet.app.run()