Задача 1

Алгоритм Amore–Zarate: воспроизведение и улучшение

Amore и Zarate (Phys. Rev. E 108, 055302, 2023) нашли рекорды для всех N ≤ 330 и для выборочных N до 40 886. Мы реализовали их алгоритм по тексту статьи, проверили, воспроизводит ли он их энергии, и попробовали его ускорить и улучшить.

Как работает алгоритм

Четыре шага алгоритма Amore–Zarate на конфигурации N = 300 — **Рис. 1.** Шаги алгоритма на реальной конфигурации N = 300. (1) N_b зарядов фиксируются на окружности, остальные удерживаются внутри круга σ. (2) Выбирается подкруг; заряды вне него заморожены. (3) Заряды внутри получают случайный толчок и спускаются в ближайший минимум; ход принимается, если энергия уменьшилась. (4) Следующий подкруг ставится на заряд, сместившийся сильнее всех («погоня» за дефектом). В конце ограничение σ снимается и делается полная релаксация.

Для N < 1000 авторы использовали спуск усечённым методом Ньютона (TN) из случайных стартов и basin hopping (BH); для N ≥ 1000 — divide & conquer (шаги 2–4). Часть параметров в статье не указана; наш выбор перечислен на странице «Методы».

Воспроизведение при N = 10…330

Для каждого N брали N_b по формуле Amore–Zarate (если их лучшее N_b отличалось, считали и его), 10 случайных стартов с TN, затем BH с TN до достижения их энергии или до лимита времени (от 30 с при малых N до 12 мин при N = 330). Один прогон на каждое N, одно ядро.

314 из 321значений N: найдена та же энергия, что у Amore–Zarate

1,9 смедианное время до рекорда

0значений N, где нашлась энергия ниже рекорда

Время до рекорда Amore–Zarate в зависимости от N — **Рис. 2.** Время, за которое BH+TN находит энергию Amore–Zarate (одно ядро M4). До N ≈ 250 рекорд находится за секунды; дальше время быстро растёт, и для N = 299, 300, 301, 325, 327, 330 лимита не хватило; при N = 127 BH остановился в другом минимуме, выше на 2,9·10⁻⁵. Совпадение проверено доводкой Ньютоном (разница энергий < 10⁻¹¹·E).

Табл. 1. Воспроизведение по диапазонам N
N	значений	рекорд найден	медиана, с	максимум, с
10–100	91	91	0,5	4
101–200	100	99	1,4	14
201–300	100	98	5,2	248
301–330	30	26	101,0	645

Большие N: divide & conquer со случайного старта

Повторили опыт рис. 2 статьи: случайный старт с N_b из формулы, радиус подкруга растёт от 0,1 до 1 за 3100 итераций, внутри — TN. Для сравнения те же прогоны с FIRE вместо TN, с критерием Метрополиса вместо «только вниз» и с молекулярной динамикой внутри подкруга (гибрид из задачи 2). Лимит — 40 минут на прогон, одно ядро. N = 1000 (3 прогона на вариант) и N = 2000 (3 прогона для варианта статьи, по 1 для остальных).

Табл. 2. Энергия в конце прогона минус рекорд Amore–Zarate (зелёным — ниже рекорда)
Вариант	N = 1000	N = 2000
D&C + TN (как в статье)	+0,62, +0,89, +1,47	−3,97, −1,42, −0,37
D&C + FIRE	+0,39, +0,63, +0,75	−1,31
D&C + FIRE + Метрополис	+0,07, +0,53, +0,68	−0,95
D&C + MD	+3,24, +4,32, +5,24	+5,87
D&C + MD + Метрополис	+3,24, +3,56, +5,09	+5,87

Энергия в ходе divide and conquer для N = 1000 и 2000 — **Рис. 3.** Ход D&C: наилучшая энергия минус рекорд Amore–Zarate в зависимости от номера итерации (ось сжата около нуля). Первые ~700 итераций с маленьким подкругом почти ничего не меняют; энергия падает, когда радиус доходит до 0,3–0,5. При N = 2000 все пять прогонов с TN или FIRE уходят ниже нуля — ниже опубликованного рекорда. Варианты с MD внутри подкруга (фиолетовый, красный) за то же время успевают сделать в 4–7 раз меньше итераций.

Состояния в ходе одного прогона D&C при N = 1000 — **Рис. 4.** Как меняется состояние в одном прогоне D&C (N = 1000, TN). Уже после первого спуска из случайного старта середина — треугольная решётка; дальше D&C переставляет «шрамы» из пар дефектов 5–7 (красные и синие ячейки) и уводит их к границе, энергия падает с +4,56 до +0,62 над рекордом.

Новый рекорд для N = 2000

5 из 7 прогонов D&C со случайного старта при N = 2000 закончились ниже опубликованного рекорда Amore–Zarate (3 003 830,44969). Лучший — E = 3 003 826,477153 при том же N_b = 442, на 3,97 ниже рекорда; это их собственный алгоритм в нашей реализации, 40 минут одного ядра. Независимо найдено и состояние с N_b = 443 («пересадка» заряда из их конфигурации и D&C, табл. 3): E = 3 003 827,089483, на 3,36 ниже рекорда. Проверка обоих: после доводки Ньютоном остаточная сила ≤ 2,1·10⁻⁵, все собственные значения гессиана положительны (кроме нулевой моды вращения), все граничные заряды прижаты наружу — это строгие локальные минимумы.

N = 2000: опубликованный рекорд и новое состояние — **Рис. 5.** N = 2000: опубликованное состояние Amore–Zarate и наше лучшее (оба с N_b = 442). Различаются расположение и длина «шрамов» — цепочек дефектов 5–7 в объёме (красные — 5 соседей, синие — 7). У нового состояния меньше дефектов: 307 против 332 (с пятью соседями — 52 против 66).

Можно ли опустить их рекорды при N = 2000…40 886

Для больших N Amore и Zarate пробовали 1–2 значения N_b и опубликовали по одной конфигурации. Мы взяли их конфигурации и применили два хода: «пересадку» граничного заряда (N_b ± 1, 2) и D&C с крупными подкругами (20–60 минут на 9 потоках). Затем и их, и наши конфигурации довели одинаково глубоко (метод Ньютона–CG, пока прогнозируемый выигрыш шага не станет меньше 10⁻⁵) и сравнили энергии.

Табл. 3. Большие N: опубликованная энергия, та же конфигурация после доводки и наша лучшая
N	опубликовано	их конф. после доводки	наша	N_b	разница	как найдена
2 000	3 003 830,450	3 003 830,449	3 003 826,477	442	−3,97	D&C с нуля
5 000	19 087 485,442	19 087 488,159	19 087 486,958	820	−1,20*	их конфигурация + D&C
10 000	76 987 197,734	76 987 197,734	76 987 194,960	1307	−2,77	их конфигурация + D&C + GPU
14 180	155 297 471,174	155 297 471,119	155 297 465,724	1647	−5,40	их конфигурация + D&C
40 886	1 300 043 776,140	1 300 043 775,700	1 300 043 761,671	3356	−14,03	их конфигурация + D&C + GPU

Разница — наша энергия минус энергия их конфигурации после той же доводки. * Для N = 5000 опубликованное число (19 087 485,4425) на 2,7 ниже энергии координат из их же архива; ниже архивных координат мы опустились, ниже опубликованного числа — нет. Для остальных N опубликованные числа совпадают с энергией их координат.

Проверка найденных состояний. Для каждого файла координат независимым кодом (numpy + точная сумма math.fsum, без нашего C-ядра) пересчитаны: число зарядов, |r| ≤ 1 для всех, число зарядов на окружности и энергия — совпадение с таблицей до 5·10⁻¹⁵ относительной точности. Для N = 2000 и 5000 посчитан полный спектр гессиана: отрицательных собственных значений нет (кроме нулевой моды вращения), то есть это строгие локальные минимумы. Для больших N полный гессиан не помещается в память; там утверждение «энергия ниже рекорда» опирается на сам факт существования конфигурации с такой энергией.

Координаты найденных состояний (текст «x y», первые N_b строк — на окружности): N = 2000 N_b = 442 и N_b = 443, N = 5000, N = 10 000, N = 14 180, N = 40 886.

N = 40 886: дефекты по всему диску и увеличенный фрагмент у границы — **Рис. 6.** Состояние N = 40 886 (Amore–Zarate после доводки). Слева показаны только заряды с числом соседей ≠ 6: зёрна треугольной решётки разделены цепочками дислокаций, у края — плотный пограничный слой. Справа — увеличенный фрагмент у границы.

Продолжение на GPU (vast.ai, одна A100)

Цель — проверить, не остановились ли мы при больших N слишком рано. Шаг D&C перенесён на GPU (code/gpu_dc.py, CuPy, двойная точность). Один ход: подкруг радиуса 0,04–0,4 (от 20 до 3700 зарядов) получает толчок 0,3–0,4 местного шага решётки и спускается в минимум (FIRE), остальные заряды заморожены; изменение энергии считается точно. В раундах 2–3 ход может быть принят и с небольшим ростом энергии (критерий Метрополиса, T = 0,05); лучшее состояние запоминается. Старт — наши состояния из табл. 3.

Проверка ядра: энергия совпадает с точной суммой до 1,5·10⁻¹⁶ (отн.), силы — до 10⁻¹⁵. Скорость: полная энергия при N = 40 886 — 13 мс (одно ядро M4 — около 0,5 с), силы на 1000 зарядов подкруга — 0,37 мс (на ядре M4 — около 24 мс).

Табл. 4. Прогоны на GPU: 2 370 ходов за 111 минут
раунд	N	старт	радиус подкруга	приём хода	мин	ходов	улучшений	ΔE на GPU
1	40 886	наш рекорд	0,05–0,3	только вниз	29	120	0	0
1	40 886	наш рекорд	0,08–0,4	только вниз	31	80	0	0
1	14 180	наш рекорд	0,07–0,4	только вниз	31	220	0	0
2	40 886	наш рекорд	0,04–0,2	T = 0,05	40	262	2	−2,362
2	14 180	наш рекорд	0,06–0,3	T = 0,05	40	435	0	0
3	40 886	итог раунда 2	0,04–0,2	T = 0,05	40	223	5	−0,006
3	10 000	наш рекорд	0,06–0,3	T = 0,05	40	465	24	−0,492
3	5 000	наш рекорд	0,08–0,4	T = 0,05	40	565	0	0

Раунд 1 остановлен через 30 минут: ни одного улучшения. ΔE на GPU — до общей доводки.

Изменение энергии каждого хода D&C на GPU в зависимости от размера подкруга — **Рис. 7.** Каждый ход D&C на GPU: изменение энергии (заряды вне подкруга заморожены) в зависимости от числа зарядов в подкруге. Почти все перестройки поднимают энергию на 1–100 или возвращают систему в тот же минимум. При N = 5000 и 14 180 нет ни одного хода вниз, при N = 10 000 и 40 886 — есть.

Главное наблюдение. Пока заряды вне подкруга заморожены, вокруг перестроенного подкруга остаётся напряжение (силы у края до 10⁴), поэтому ΔE хода — только верхняя оценка выигрыша. После общей доводки методом Ньютона–CG новое состояние оказывается заметно ниже: при N = 40 886 — на 13,18 ниже нашего прежнего (на GPU было видно 2,36), при N = 10 000 — на 0,86 (на GPU 0,49). Контроль: та же доводка, повторённая для наших прежних состояний, не меняет их энергию (изменение −2,1·10⁻⁷ при N = 10 000 и −9,8·10⁻⁵ при N = 40 886) — прежние состояния были настоящими минимумами, новые — другие, более низкие минимумы. Новые координаты проверены независимым пересчётом (numpy + math.fsum, |r| ≤ 1, N_b тот же).

N = 40 886: где изменилось состояние после ходов на GPU — **Рис. 8.** Что изменилось при N = 40 886. Слева — дефекты нового состояния, рамки — места перестроек. Справа — эти места до и после (цвет — число соседей: красный 5, синий 7; кольца — заряды, сместившиеся больше чем на четверть шага решётки). Перестроены два участка дислокационных цепочек у края (r ≈ 0,94–0,96): заметно сдвинулись 33 заряда из 40 886, число дефектов почти не изменилось (3194 → 3197). Выигрыш даёт расположение дефектов, а не их число; остальные заряды лишь слегка подстраиваются при общей доводке.

Выводы

Алгоритм Amore–Zarate восстановлен по тексту статьи и воспроизводит их рекорды: 314 из 321 значений N ≤ 330 за секунды–минуты на одном ядре; трудные N (≥ 300) находятся при более длинном прогоне (задача 2, табл. 2).
Их же алгоритм в нашей реализации даёт энергии ниже их рекордов при больших N: N = 2000 со случайного старта — на 3,97; из их конфигураций — ниже их же конфигураций при всех пяти проверенных N (до 14,0 при N = 40 886) и ниже опубликованных чисел при четырёх из пяти (кроме N = 5000, см. табл. 3). Опубликованные конфигурации при N ≥ 2000 — не основные состояния.
Что помогает: крупные подкруги D&C (радиус 0,3–0,6), «пересадка» заряда между границей и объёмом (при N = 2000 найдено состояние с N_b = 443 вместо их 442), доводка до сходимости. FIRE вместо TN — сопоставимо по качеству и проще.
GPU (одна A100, 111 мин): D&C с Метрополисом и общей доводкой опустил наши состояния ещё на 0,86 (N = 10 000) и 13,18 (N = 40 886); при N = 5000 и 14 180 улучшений нет. Замороженное окружение занижает выигрыш хода в 2–6 раз: принимать ходы по энергии после общей релаксации — следующий шаг.