Задача 2

Молекулярная динамика и гибрид с basin hopping

Лавров и Никонов (arXiv:2609.20777) ищут минимумы затухающей молекулярной динамикой (QMD) с той же граничной эвристикой, что у Amore–Zarate, и считают её для малых N лучше basin hopping. Для больших N предложен гибрид: выбрать «бассейн» (подобласть), спустить его молекулярной динамикой, принять или отвергнуть результат по Метрополису — и так по кругу. Мы сравнили эти методы при одинаковых затратах.

QMD по протоколу Лаврова–Никонова, N = 10…330

Всё как в статье: N_b по формуле Amore–Zarate, случайный старт, 150 000 шагов схемы Верле с Δt = 0,001 и трением b = 0,9, σ плавно поднимается до 1, затем финальная релаксация; 10 независимых запусков на каждое N. Единственное наше добавление — короткая релаксация FIRE перед стартом: без неё при N ≳ 200 близкие пары случайного старта разгоняют заряды до скоростей 10³–10⁵, и конфигурация разрушается.

3210запусков QMD (321 значение N × 10)

0энергий ниже рекорда Amore–Zarate

57 %запусков, где к концу изменилось число граничных зарядов

Чтобы сравнить с basin hopping честно, считаем ожидаемое время до рекорда: время одного запуска QMD, делённое на вероятность того, что запуск найдёт рекорд. Для BH берём медианное время из задачи 1 (тот же процессор, одно ядро).

Табл. 1. QMD и BH+TN (Amore–Zarate) при равных условиях
N	доля запусков QMD с рекордом	время запуска, с	QMD: ожид. время, с	BH+TN: медиана, с	QMD / BH
10–50	0,390	0,2	0,5	0,1	3,7
51–100	0,086	0,9	11	1,0	11
101–150	0,064	2,3	37	1,7	21
151–200	0,108	4,4	40	1,4	29
201–250	0,058	7,0	121	3,1	39
251–300	0,022	10	467	11	41
301–330	0,003	13	4005	101	40

Доля успешных запусков QMD и доля запусков с изменившимся числом граничных зарядов — **Рис. 1.** Сверху — доля из 10 запусков QMD, нашедших рекорд Amore–Zarate. Снизу — доля запусков, в которых к концу изменилось число граничных зарядов N_b: несколько внутренних зарядов «садятся» на границу, и запуск заканчивается уже в другом классе состояний (почему — см. ниже).

Основное состояние N = 300 и три типичных исхода QMD — **Рис. 2.** N = 300: основное состояние и три типичных исхода запуска QMD. У основного состояния все дефекты (красные — 5 соседей, синие — 7) лежат в пограничном слое, середина — чистая треугольная решётка. QMD застревает с дефектами в объёме (пары 5–7, кластеры) или с лишними зарядами на границе.

Откуда берётся дрейф N_b и как его убрать

Первая гипотеза — виноват плавный подъём σ до 1. Проверили напрямую: те же 20 запусков при десяти значениях N от 100 до 330 с подъёмом σ и с σ, зафиксированным до финальной релаксации. Результаты совпали почти запуск в запуск (различия — один запуск из 20 при двух N): подъём σ ни на что не влияет.

Настоящая причина — в старте. Если брать начальные радиусы по плотности Amore–Zarate (ρ ∝ 1/√(1−r²)), у края оказывается много зарядов. Уже после первой релаксации 2–3 из них прижаты к ограничению r = σ. Так они остаются всю динамику, а когда ограничение снимают, переходят на границу, и N_b вырастает. Если после старта к ограничению никто не прижат, дрейфа нет. Поэтому сравнили три старта: по плотности, по плотности, но поджатый к центру на 10 %, и равномерный по кругу (по 20 запусков на N).

Число запусков QMD с рекордом и с изменившимся N_b для трёх типов старта — **Рис. 3.** Запуски QMD (протокол Лаврова–Никонова) из трёх стартов. Поджатый или равномерный старт находит рекорд в 14 и 13 из 100 запусков против 4 из 100 у старта по плотности, а N_b меняется в 31–37 % запусков вместо 56 %.

Трудные N: все методы при одинаковом времени

Взяли значения N, на которых в задаче 1 basin hopping искал рекорд дольше всего или не нашёл (299, 300, 301, 304, 316, 318, 319, 324, 325, 327, 330), и дали каждому методу по 15 минут одного ядра, по два независимых прогона. Методы: BH с TN (Amore–Zarate), BH с FIRE, QMD с подъёмом σ и без него (повтор запусков со случайных стартов до конца времени), гибрид «бассейн → MD → Метрополис» (тепловой толчок зарядам подкруга, затухающая MD всех зарядов, спуск FIRE, Метрополис при T = 0,01).

Табл. 2. Время до рекорда Amore–Zarate, секунды (два прогона; «—» — не найден за 15 минут)
N	BH + TN (Amore–Zarate)	BH + FIRE	QMD (Лавров–Никонов)	QMD, σ фикс.	гибрид бассейн→MD→Метрополис
299	— / —	7 / —	— / —	— / —	— / —
300	— / 502	— / —	— / —	— / —	— / —
301	346 / 750	466 / —	— / —	— / —	— / —
304	40 / 30	— / 14	— / 804	— / 817	— / —
316	112 / 10	396 / —	— / —	— / —	— / 504
318	44 / 9	536 / 162	— / —	— / —	— / —
319	11 / 6	666 / 648	— / —	— / —	— / —
324	10 / 10	26 / 30	842 / 517	858 / 523	— / —
325	143 / 284	444 / 134	— / —	— / —	— / —
327	123 / 151	42 / —	— / —	— / —	— / —
330	— / —	— / —	— / —	— / —	— / —
найдено	17 из 22	13 из 22	3 из 22	3 из 22	1 из 22

Большие N: гибрид внутри divide & conquer

Схему научного руководителя — подкруг («бассейн») → спуск молекулярной динамикой → критерий Метрополиса — мы встроили в D&C Amore–Zarate и запускали со случайного старта при N = 1000 и 2000, по 40 минут (подробно — задача 1, табл. 2 и рис. 3).

MD с постоянным шагом (схема Лаврова–Никонова, Δt = 0,001) внутри подкруга сходится медленно. За 40 минут D&C успевает 700–900 итераций при N = 1000 и около 200 при N = 2000 против 1400–3100 с TN или FIRE, а итог — на 3–6 единиц выше рекорда.
MD с адаптивным шагом (FIRE) работает. При N = 2000 D&C + FIRE и D&C + FIRE + Метрополис опустились ниже рекорда Amore–Zarate (на 1,31 и 0,95). При N = 1000 лучший из всех прогонов — D&C + FIRE + Метрополис (+0,07).

Задача 2½. Можно ли не мельчить шаг схемы Верле

Почему молекулярная динамика «задыхается» при больших N. Самые жёсткие колебания — у границы, где заряды сидят плотнее всего. Мы посчитали спектр колебаний для основных состояний N = 100…2000: наибольшая частота растёт ровно как N (ω_max ≈ 0,85 N), а схема Верле устойчива только при Δt < 2/ω_max. Шаг Лаврова–Никонова Δt = 0,001 перестаёт быть устойчивым около N ≈ 2300; дальше шаг приходится уменьшать пропорционально 1/N, а число шагов — увеличивать.

Проверили два способа обойти это ограничение, не меняя самих минимумов:

Массы по жёсткости. Заряду даётся масса m_i ∝ Σ_j 1/r_ij³ — у края частицы тяжелее. Энергия от масс не зависит, а предельный шаг падает медленнее: как N^−0,73 вместо N^−1,02, при N = 2000 он в 2,7 раза больше. При N = 300 схема устойчива с Δt = 0,005 и сходится за 5 000 шагов вместо 25 000 при Δt = 0,001.
Адаптивный шаг (FIRE). Та же затухающая динамика, но шаг и трение подстраиваются по ходу. При N = 1000 она сходится за 500 шагов вместо 25 000 у схемы с постоянным шагом — в 50 раз быстрее. Поэтому в гибриде и в D&C мы используем FIRE.

Предельный шаг схемы Верле и число шагов до сходимости — **Рис. 4.** Слева: предельный шаг 2/ω_max для основных состояний N = 100…2000 при единичных массах (красный) и массах по жёсткости (синий); пунктир — шаг Лаврова–Никонова. Справа: N = 1000, число шагов затухающей динамики (b = 0,9) до сходимости в зависимости от шага; крестик — разлёт; чёрные квадраты — FIRE.

Статью Э. Г. Никонова о схеме с крупным шагом, упомянутую на встрече, мы в открытом доступе не нашли. Ближайшая — симметричная симплектическая схема 3-го порядка (Батгэрэл, Никонов, Пузынин, 2016); она требует вторых и третьих производных энергии и для кулоновской задачи с N ≳ 10³ дорога. Её проверку можно сделать отдельно, если пришлёте статью.

Выводы