Поверхность. Рентгеновские, синхротронные и нейтронные исследования, 2020, № 11, стр. 84-87

ВВЕДЕНИЕ

Гипердинамика [1] – один из методов ускоренного молекулярно-динамического моделирования. К этим методам относятся также температурно-ускоренная динамика [2] и метод параллельных реплик [3]. Перечисленные подходы были разработаны как способы динамического моделирования процессов, важную роль в которых играют так называемые редкие события. К таким событиям относятся переходы системы между локальными минимумами потенциальной энергии. Моделирование этих процессов методами классической молекулярной динамики (МД) представляет большую вычислительную сложность, так как продолжительность переходов мала по сравнению со временем, в течение которого система находится в окрестности потенциального минимума между переходами. Таким образом, большая часть ресурсов при подобном моделировании расходуется на моделирование поведения системы, находящейся около энергетического минимума, а не на моделирование переходов. Для решения данной проблемы были разработаны методы ускоренной молекулярной динамики. Они используются при моделировании различных систем начиная с кристаллов [4, 5] и заканчивая биомолекулами [6, 7].

Основная идея гипердинамики состоит в понижении барьеров между минимумами потенциальной энергии системы, т.е. в уменьшении энергии активации переходов. Это достигается путем введения, дополнительного слагаемого U_b в потенциал системы:

где U₀ – исходный потенциал межатомного взаимодействия в системе, U – смещенный потенциал, используемый при гипердинамическом моделировании. Существуют различные подходы к построению U_b, которые постоянно развиваются. Часть этих подходов представлена в [8–10]. Рассмотрим метод гипердинамики со слагаемым U_b, для построения которого не требуется большого количества вычислительных ресурсов. Такой способ построения потенциала U рассматривали ранее в [11, 12], его тестировали на двумерных атомных системах с применением парного потенциала. В настоящей работе проводится моделирование трехмерной атомной системы с использованием многочастичного потенциала для демонстрации применимости данного способа построения смещенного потенциала при моделировании реалистичных систем. Критерием корректности получаемых результатов является их соответствие данным, полученным методом МД.

МЕТОДИКА

При моделировании используется потенциал следующего вида:

где суммирование происходит по атомам системы, слагаемое U_b,i определяется выражением: где V₀ – параметр, на значение которого снижается высота потенциальных барьеров системы; ρ_i – расстояние от i-го атома до ближайшего минимума потенциальной энергии; R₁, R₂ – параметры потенциала; функция φ(ρ_i) обеспечивает непрерывность U_b,i, а также некоторого числа ее производных.

Особенности данного потенциала, а также ограничения, накладываемые на параметры R₁ и R₂, рассмотрены в [11]. Для всех вычислений, результаты которых описаны ниже, параметр R₁ принимал одно из следующих значений: 0.26, 0.27, 0.28, 0.3 Å, а параметр R₂ – 1.0, 1.1, 1.15 Å. Функция φ(ρ_i) – полиномиальный сплайн пятого порядка, состоящий из двух фрагментов.

Система представляла собой кристалл алюминия, включающий вакансию. Образец содержал 8191 атом и был искусственно сжат в направлениях [110] и [1$\overline 1 $0], так что средние расстояния между атомами в данных направлениях сократились на 1 и 0.5% соответственно. Таким образом, для вакансии в любой момент времени существовали конкурирующие переходы трех типов. Методика модельных сжатий кристалла позволяет проверить, будет ли гипердинамический подход давать такие же частоты для различных переходов, что и классический метод МД. В качестве исходного потенциала U₀ был применен многочастичный потенциал Акланда [13]. Для установления и поддержания постоянной температуры кристалла был использован метод “термической ванны” Берендсена [14], который уже применялся нами ранее [15, 16]. Интегрирование уравнений движения выполнялось согласно скоростному методу Верле. Моделирование проводили при температурах 650, 600 и 550 К. При каждой температуре было получено по 1000 переходов для каждого набора параметров R₁, R₂, V₀ методом гипердинамики и столько же методом МД. Для сравнения результатов, полученных различными методами, использовали метод однородности Смирнова [17]. Данный метод применяли для подтверждения гипотезы о подчинении полученных методами МД и гипердинамики выборок времен, проходящих между двумя последовательными переходами системы, одному закону распределения.

РЕЗУЛЬТАТЫ И ИХ ОБСУЖДЕНИЕ

Анализ результатов моделирования показал, что выборки времен переходов, которые дает метод гипердинамики, подчиняются тому же закону распределения, что и выборки, получаемые методом МД. Эта гипотеза была подтверждена на уровне значимости 0.1 при V₀, не превышающем 0.15 эВ, и на уровне значимости 0.05 при более высоких значениях V₀.

На рис. 1 представлена зависимость ускорения τ по времени вычисления, получаемого при использовании метода гипердинамики при сравнении с классическим методом МД, от значения параметра V₀. Здесь τ определено следующим образом:

где t_MD и t_HD – вычислительное время, затраченное на совершение системой 1000 переходов, при использовании методов МД или гипердинамики соответственно. В исследованиях при уменьшении моделируемой температуры использовали бóльшие значения V₀. Однако в [11] было показано, как при одинаковых значениях параметров смещенного потенциала при низкой моделируемой температуре могут быть получены менее точные значения, чем в случае высокой температуры. В представленном здесь исследовании ухудшения качества результатов удалось избежать путем уменьшения значения параметра R₁ при понижении температуры. Это стало возможным, поскольку при более низких температурах система колеблется в меньшей пространственной области около минимума энергии. С уменьшением параметра R₁ можно повысить значение параметра V₀ без потери качества результатов.

На рис. 1 видно, что при одинаковых значениях V₀ вычислительное ускорение увеличивается с ростом температуры. На графике с осью абсцисс V₀/k_BT (рис. 2), видно, что выполняется следующее соотношение:

где k_B – постоянная Больцмана, Т – моделируемая температура. Это хорошо согласуется с положениями метода гипердинамики, описанными в [1]. Коэффициент пропорциональности в данном случае будет тем меньше, чем больше времени будет расходоваться на нахождение расстояний ρ_i. Исследуемая в настоящей работе система содержит малый процент атомов, которые могут совершать перескоки в положение вакансии. В выражении (2) суммирование проводилось только по таким атомам и вычислялись только соответствующие расстояния ρ_i. В представленном подходе это возможно, так как параметр R₁ может быть различным для разных U_b,i. Полагалось, что значение R₁ достаточно велико для всех атомов, не соседствующих с вакансией, чтобы они всегда находились в области ρ_i < R₁. В этих случаях U_b,i равен V₀, и, таким образом, силы, действующие на атомы в смещенном потенциале, не отличаются от соответствующих сил в исходном потенциале. Однако подобные предположения не являются обязательными, особенно в случае низких температур и больших значений V₀. Хотя они позволяют получить дополнительное ускорение в несколько раз, данный выигрыш в ускорении не сопоставим со значением экспоненты в выражении (5). В частности, полученные результаты при Т = 300 К и V₀ = 0.3 эВ показывают, что величина отношения t_MD/t_HD порядка 10⁵, тогда как при Т = 650 К она составляет 10².

ЗАКЛЮЧЕНИЕ

Описанный в [11] способ построения смещенного потенциала был протестирован на трехмерной атомной системе с использованием многочастичного потенциала. Показано, что в этом случае метод гипердинамики позволяет получать результаты, сопоставимые с результатами, которые дает метод классической МД. Получено, что ускорение растет экспоненциально с увеличением параметров V₀ и 1/T. Таким образом, при низких температурах, когда метод МД не позволяет моделировать процессы переходов системы между состояниями, метод гипердинамики дает возможность такого моделирования.