Поверхность. Рентгеновские, синхротронные и нейтронные исследования, 2019, № 7, стр. 109-112

ВВЕДЕНИЕ

Методы ускоренной молекулярной динамики были разработаны как более эффективное, чем молекулярная динамика (МД), средство динамического моделирования процессов, связанных с так называемыми редкими событиями. Такими событиями являются термически активируемые переходы, исследование которых важно для понимания процессов, определяющих, например, диффузию в атомных системах. Промежутки времени между такими переходами существенно больше продолжительности отдельно взятого перехода. Основа методов ускоренной молекулярной динамики – сокращение времени между двумя последовательными редкими событиями. Обязательным условием при этом является сохранение отношений вероятностей осуществления различных событий. Также критически важна возможность оценки промежутков времени между последовательными переходами. Это необходимо для понимания не только динамики системы, но и характерных масштабов времени, в которых она наблюдается, а также возможности оценки таких характеристик, как, например, коэффициент диффузии.

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

ОБЪЕДИНЕНИЕ МЕТОДОВ ТУД И ГД

Пусть U(r₁, r₂, …, r_N) – потенциал системы, состоящей из N атомов. Будем использовать следующий смещенный потенциал U *(r₁, r₂, …, r_N):

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

Для того чтобы не происходило тепловое расширение исследуемой системы, а также для улучшения оценок изменения частот нормальных мод системы при нагревании, используем метод ТУД с измененным взаимодействием между атомами системы. В такой модификации метода на i-й атом системы будет действовать сила:

где U(r₁, r₂, …, r_n) – U_i,min соответствует изменению энергии системы в случае перемещения i-го атома в ближайший энергетический минимум. Функция А(ρ_i) имеет следующий вид: где T_low – температура, при которой находится исследуемая система; T_high – температура, до которой система нагревается при ТУД-моделировании.

С учетом (1) сила, действующая на i-й атом, запишется в виде:

Моделирование описанной выше системы происходит в соответствии с методом ТУД. Вычисления выполняются с измененным взаимодействием между атомами системы согласно (5). При попытке перехода системы из одного состояния в другое данная попытка фиксируется, а система возвращается в исходное состояние. Принятие некоторой попытки перехода и перевод системы в новое состояние осуществляется по алгоритму, описанному в [2, 3]. Однако из-за изменения взаимодействия между атомами согласно (5) время t_low, приписываемое переходу после его принятия, рассчитывается следующим образом:

где E – высота энергетического барьера в системе с первоначальным потенциалом U(r₁, r₂, …, r_N).

Чтобы дополнительно ускорить вычислительный процесс, силу F_i для атомов, которые не участвуют в переходах, определяют следующим образом:

ОПИСАНИЕ МОДЕЛИ И ПОЛУЧЕННЫЕ РЕЗУЛЬТАТЫ

Для проверки предложенного метода моделирования ТУГД были проведены вычисления, в которых степень достоверности полученных результатов определяли путем их сравнения с результатами, полученными методом классической молекулярной динамики. Для получения достаточного количества событий термически активированных переходов моделирование атомной системы проводили в двумерном пространстве. Система представляла собой двумерный кристалл меди, содержащий вакансию. При моделировании был использован потенциал Морзе [7], который на расстояниях меньше 1 Å сшивался посредством полинома пятой степени с потенциалом Зиглера–Бирзака–Литтмарка [8]. Размер кристалла составил 16 × 16 параметров решетки. В горизонтальном и вертикальном направлениях были установлены периодические граничные условия. Для поддержания постоянной температуры кристалла использовали метод термической ванны Берендсена [9]. Мгновенную температуру Т определяли следующим образом:

где ${{\bar {E}}_{k}}$ – мгновенное значение средней кинетической энергии атомов системы. Выражение (8) отражает двумерный характер моделирования. Для определения высот энергетических барьеров использовался метод эластичной упругой ленты (NEB – Nudged Elastic Band) [10].

Вдоль определенного направления постоянная решетки кристалла была искусственно уменьшена на 0.25%. Таким образом, доступные в любой момент времени шесть переходов вакансии разделились на два типа: два перехода вдоль направления, в котором было произведено сжатие, и четыре в других направлениях. В выражении (6) показатель степени 3/2 был заменен единицей вследствие перехода к двумерному пространству. Параметры R₁ и R₂ были приняты равными, соответственно, 0.3 и 0.85 Å.

В каждом эксперименте было получено 500 переходов на расстояниях не меньше одной постоянной решетки от атомов диссипативных слоев до начального и конечного положения атома, совершившего переход. В результате моделирования сравнивали функции распределения времени, затраченного на атомный переход, полученные с помощью объединенного метода ТУГД, и соответствующие функции распределения, полученные методом классической МД. Для проверки нулевой гипотезы о соответствии эмпирических распределений, полученных обоими методами, одному и тому же закону использовали критерий однородности Смирнова [11].

Результаты моделирования представлены в табл. 1. Из таблицы видно, что результаты, полученные методом ТУГД, лучше соответствуют данным МД-моделирования, чем результаты моделирования методами ТУД и ГД, используемыми по отдельности.

Следует также отметить, что моделирование методом ТУД дает менее точные результаты, чем метод ГД. Это объясняется тем, что увеличение температуры системы согласно (5) приводит к изменению частот нормальных мод системы. Это также видно из данных, полученных объединенным методом ТУГД – повышение параметра T_high значительно ухудшает точность по сравнению с повышением параметра V₀. При этом параметры для вычислений подбирали таким образом, чтобы во всех случаях, кроме последнего в таблице, были близкие ускорения вычислений по сравнению с методом МД.

ЗАКЛЮЧЕНИЕ

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

Особенно эффективным по сравнению с оригинальными методами ускоренной молекулярной динамики предложенный метод может быть в случаях, когда метод температурно-ускоренной динамики может давать большую скорость вычислений, чем метод гипердинамики, но не позволяет получать удовлетворительную точность результатов при высоких значениях отношения T_high/T_low. В этом случае снижение скорости вычислений из-за выбора меньшего значения T_high может быть частично компенсировано увеличением параметра смещения потенциала V₀.

БЛАГОДАРНОСТИ

Работа выполнена с использованием оборудования Центра коллективного пользования сверхвысокопроизводительными вычислительными ресурсами МГУ имени М.В. Ломоносова [12].