Участник:Ivan Shmakov/Параллельные вычисления в моделировании графена 2014

Хотя графен как объект теоретико-физического исследования не нов — по-видимому, первая рассматривающая одиночный слой графита работа была опубликована в 1947 г. — значительный интерес к проблеме возник только с появлением в 2004 г. практического способа получения «макроскопических» образцов материала. Среди предполагаемых применений можно найти и высокоэффективные фильтры для опреснения воды, и прозрачные электроды светоизлучающих электрохимических ячеек (как замена оксиду индия-олова.)

Однако, сколь угодно широкое применение «прямого» (ab initio) квантово-теоретического подхода в моделировании графена сдерживается большим количеством атомов (от сотен до десятков тысяч) для представляющих интерес моделей, что вынуждает прибегнуть к приближенным методам, одним из которых является моделирование межатомных связей потенциалом Морзе.

В данной работе рассматривается программное обеспечение для моделирования геометрии образцов графена и нахождения энергии связи в приближении потенциального взаимодействия по межатомным связям с использованием параллельных вычислительных систем, реализующих спецификацию MPI.

Формулировка задачи

$\$

Каждый из атомов (𝐫_𝑖, 𝑚_𝑖, 𝑖 = 1, …, 𝑁) считается находящимся в потенциальном поле 𝑉_𝑖, представленном как сумма полей 𝑉_𝑖𝑗, создаваемых связанными с ним (𝐷_e𝑖𝑗 ≠ 0) другими атомами. При использовании потенциала Морзе, выражения принимают следующий вид (где 𝐷_e𝑖𝑗, 𝜌_e𝑖𝑗, 𝑎_e𝑖𝑗 — параметры химической связи.)

${\begin{aligned}V_{i}&=\sum _{j}V_{ij},\quad {D_{\mathrm {e} ij}\neq 0};\\V_{ij}&=D_{\mathrm {e} ij}\left[1-\exp(-a_{ij}(R_{ij}-\rho _{\mathrm {e} ij}))\right]^{2};\\R_{ij}&={\sqrt {(\mathbf {r} _{i}-\mathbf {r} _{j})^{2}}}.\\\end{aligned}}$

Кроме того, может иметь смысл включить в потенциал следующие вращательные компоненты, а также компоненту дальнодействующего (независимого от химических связей) ван-дер-ваальсова взаимодействия.

${\begin{aligned}V_{ijk}^{\mathrm {a} }&={{\frac {1}{2}}K_{\mathrm {\theta } }\left[\cos \theta _{ijk}-\cos \theta _{0}\right]^{2}};\\V_{ijkl}^{\mathrm {t} }&={{\frac {1}{2}}K_{\mathrm {\varphi } }\left[1-\cos(2\varphi _{ijkl})\right]^{2}};\\V_{ij}^{\mathrm {vdW} }&={4\varepsilon \left[\left({\frac {\sigma }{R_{ij}}}\right)^{12}-\left({\frac {\sigma }{R_{ij}}}\right)^{6}\right]}.\\\end{aligned}}$

Найти (или, скорее, — уточнить) геометрию системы можно или решив задачу оптимизации (нахождения минимума потенциальной энергии системы; для чего применим, в частности, метод градиентного спуска), или же решив задачу прямого моделирования динамики системы в рамках классической механики (например, методами Рунге — Кутты.) Поскольку математически данные подходы в известной мере близки, ограничимся рассмотрением последнего из них.

Решение обыкновенных дифференциальных уравнений методами Рунге — Кутты требует приведения исходной системы к виду 𝐲̇ = 𝐟 (𝑦), для чего уравнения динамики записываются следующим образом.

$\left\{{\begin{aligned}{\dot {\mathbf {r} }}_{i}&={\mathbf {v} }_{i},\\{\dot {\mathbf {v} }}_{i}&=m_{i}\nabla V_{i}.\end{aligned}}\right.$

Для ∇𝑉_𝑖 получаем следующее выражение.

${\begin{aligned}\nabla V_{i}&=\sum _{j}\nabla V_{ij};\\\nabla V_{ij}&=2a_{ij}D_{\mathrm {e} ij}(1-\kappa _{ij})\kappa _{ij}{\frac {\mathbf {R} _{ij}}{R_{ij}}};\\\kappa _{ij}&=\exp(-a_{ij}(R_{ij}-\rho _{\mathrm {e} ij}));\\\mathbf {R} _{ij}&=\mathbf {r} _{i}-\mathbf {r} _{j}.\end{aligned}}$

Многопоточность

Пример топологии образца графена (10 колец, 34 атома.)

При выполнении вычислений «в несколько потоков», каждый поток моделирует динамику некоторой доли атомов образца, для чего вычисляются значения действующих на атом сил 𝑚_𝑖∇𝑉_𝑖 = 𝑚_𝑖 ∑_𝑗 ∇𝑉_𝑖𝑗. При этом, хотя равенство ∇𝑉_𝑖𝑗 = − ∇𝑉_𝑗𝑖 позволяет вычислять «напряженность» каждой связи однократно на каждой итерации, обмен получаемыми значениями между потоками усложняет код и предполагает дополнительные накладные расходы, связанные с использованием сети передачи данных. Как следствие, в текущем варианте кода «повторные» вычисления опускаются только в тех случаях, когда оба связанных атома «принадлежат» одному потоку. Это, в свою очередь, делает желательным такое распределение атомов между потоками, при котором количество связей между атомами различных потоков будет минимально.

Кроме того, представляется полезным использовать как можно более простой алгоритм обмена результатами каждой итерации (𝐫_𝑖, 𝐯_𝑖) между потоками. В данном случае, с каждым потоком помимо диапазона индексов моделируемых атомов [𝑚; 𝑛], связывается диапазон индексов атомов, данные о состоянии которых необходимы для такого моделирования, — [𝑚′; 𝑛′] ⊇ [𝑚; 𝑛] (где [𝑚′; 𝑛′] = [𝑚; 𝑛] возможно только в случае моделирования одним потоком, или же при моделировании образца с неодносвязным графом связей.)

В рамках топологии образцов, подобной представленной на рисунке, использование такого алгоритма предполагает моделирование каждым потоком целого количества строк. В данном примере, это распределение (для случая шести потоков) может быть следующим.

Поток	Моделирует	Требует
0	[0; 4]	[0; 9]
1	[5; 10]	[0; 16] ∖ {11}
2	[11; 16]	[6; 21]
3	[17; 22]	[11; 28] ∖ {23}
4	[23; 28]	[18; 33]
5	[29; 33]	[23; 33] ∖ {28}

Реализация

С учетом вышесказанного, был разработан программный код, реализующий следующий алгоритм.

Основной процесс	Вспомогательные процессы
Инициализация
Разбор аргументов командной строки. Загрузка данных. Формирование задач для всех процессов.
Отправка задач вспомогательным процессам.	Получение задачи для данного процесса, или запроса на завершение моделирования.
Вычисление вклада задачи основного процесса в начальные значения потенциальной и кинетической энергий. Инициирование «фиктивной» итерации (с нулевым шагом) для получения вклада задач вспомогательных процессов.
Главный цикл
Отправка параметров итерации (шаг по времени, флаг отмены последнего расчета.)	Получение параметров итерации.
	(При неактивном флаге отмены.) Ожидание и обработка запросов состояния от процессов с меньшим номером.
(При неактивном флаге отмены.) Отправка запросов состояния и получение ответов. Ожидание и обработка запросов состояния от процессов с бо́льшим номером.
Вычисление следующего состояния задачи для заданного шага. Вычисление вклада задачи в значения энергии системы.
Сбор вкладов в значения энергии системы.	Отправка вклада задачи в значения энергии системы.
Определение необходимости повтора итерации (изменение энергии более чем на 10⁻⁴ по сравнению с ожидаемым или изменение потенциальной энергии более чем на ⅛) и вычисление нового значения шага.
Завершение
Отправка запросов на завершение моделирования.
Сбор конечных состояний задач.	Отправка конечного состояния задачи.
Сохранение результирующего состояния системы.

В текущей редакции для решения системы дифференциальных уравнений используется метод «средней точки» (второй порядок точности.) В дальнейшем предполагается также изучить целесообразность применения классического метода Рунге — Кутты четвертого порядка.

Решение проверочных задач показало необходимость периодически «отнимать» от системы часть кинетической энергии частиц. В противном случае, рост кинетической энергии может существенно превысить падение потенциальной, вплоть до выхода частиц из соответствующих «потенциальных ям», что, по-видимому, объясняется накоплением погрешности вычислений. Для моделирования потери энергии на каждой итерации выполняется домножение скоростей частиц на величину 𝜇 (где 𝜏 — характерное время потери, ∆𝑡 — шаг моделирования):

$\mu ={\sqrt {1-\exp(-\tau /\Delta t)}}$

При решении проверочной задачи (22 атома, 6 колец, начальная энергия системы — 1.972 × 10⁻¹⁸ Дж, начальный шаг — 2 × 10⁻¹⁹ с), уменьшение полной энергии системы в порядка 450 раз (до 4.929 × 10⁻²¹ Дж) заняло 3.027 × 10⁻¹³ с. Расчет потребовал 64 278 итераций (средний шаг — 4.71 × 10⁻¹⁹ с) и занял порядка 6 с машинного времени.

«Черновая» реализация была опубликована 2014-02-26 и состоит из следующих компонент:

• http://am-1.org/~ivan/src/ns8xrqxiojtkiin4drpx4ne19y.c — моделирующий код;
• http://am-1.org/~ivan/src/gpqfyhge73tgcr4dkfp9xydju8.perl — вспомогательный код преобразования начального состояния системы из YAML в специфичный для моделирующего кода «двоичный» формат.

Также доступна документация.

См. также

• en:User:Ivan Shmakov/Potopt 2014
• b:en:Potopt Molecular Dynamics Suite