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

Материал из Викиверситета
Перейти к навигации Перейти к поиску

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

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

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

Формулировка задачи[править]

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

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

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

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

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

Многопоточность[править]

Пример топологии образца графена (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⁻⁴ по сравнению с ожидаемым или изменение потенциальной энергии более чем на ⅛) и вычисление нового значения шага.
Завершение
Отправка запросов на завершение моделирования.
Сбор конечных состояний задач. Отправка конечного состояния задачи.
Сохранение результирующего состояния системы.

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

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

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

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

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

См. также[править]

Примечания[править]

  1. Wallace, P. R. (1947). «The Band Theory of Graphite». Physical Review 71 (9). DOI:10.1103/PhysRev.71.622.
  2. (2009) «This Month in Physics History: October 22, 2004: Discovery of Graphene». APS News 18 (9).
  3. Cohen-Tanugi, David; Grossman, Jeffrey C. (2012). «Water Desalination across Nanoporous Graphene». Nano Letters 12 (7): 3602–3608. DOI:10.1021/nl3012853. ISSN 1530-6984. PMID 22668008.
  4. Matyba, P.; et al. (2010). «Graphene and Mobile Ions: The Key to All-Plastic, Solution-Processed Light-Emitting Devices». ACS Nano 4 (2): 637–42. DOI:10.1021/nn9018569. PMID 20131906.
  5. Guo, Y. N.; Karasawa N., Goddard W. A. (1991). «Prediction of fullerene packing in C60 and C70 crystals». Nature 351: 464–467.
  6. Le, M. Q.; Batra, R. C. (2013). «Single-edge crack growth in graphene sheets under tension». Computational Materials Science 69: 381–388.