Каркас программы для структурно-атомарного моделирования процесса роста кристалла с помощью метода Монте-Карло Текст научной статьи по специальности «Компьютерные и информационные науки»

УДК 544.412.3+519.237.3

Г. Ю. Аверчук*, М. Н. Шишов, Э. М. Кольцова, Е. С. Куркина

Российский химико-технологический университет им. Д.И. Менделеева, Москва, Россия 125480, Москва, ул. Героев Панфиловцев, д. 20 * e-mail: altermn@gmail.com

КАРКАС ПРОГРАММЫ ДЛЯ СТРУКТУРНО-АТОМАРНОГО МОДЕЛИРОВАНИЯ ПРОЦЕССА РОСТА КРИСТАЛЛА С ПОМОЩЬЮ МЕТОДА МОНТЕ-КАРЛО

С использованием техники обобщённого программирования разработан каркас программы для выполнения структурно-атомарного моделирования процесса роста кристалла. На основании базовых сущностей, предоставляемых каркасом, можно реализовать программный код, в котором моделируемый процесс будет выражен в терминах кинетической схемы. Заложенная кинетическая схема рассчитывается с помощью стохастического метода Монте-Карло.

Ключевые слова: каркас программы, моделирование, метод Монте-Карло, рост кристалла, латеральные взаимодействия, алмаз.

Моделирование процесса роста кристалла заданной конфигурации ставит ряд проблем, обусловленных необходимостью более глубокого понимания процесса роста. Рост кристалла происходит благодаря взаимодействию атомов кристалла с окружающей средой, в результате чего атомы окружающей среды встраиваются в кристаллическую решётку. Встраивание осуществляется благодаря силам, возникающим на поверхности растущего кристалла [1]. При более детальном рассмотрении процесса взаимодействия атомов можно выделить стабильные атомарные структуры, наличие которых определяет скорости элементарных стадий роста кристалла. Локальные скорости элементарных процессов эволюционируют благодаря изменению локального окружения взаимодействующих атомарных структур. Множество всех возможных элементарных процессов, имеющих место на поверхности растущего кристалла, определяет состав и скорость роста моделируемого кристалла и может быть представлено в виде кинетической схемы [2].

Определение адекватной кинетической схемы представляется достаточно сложной задачей и требует возможности проверки составленной кинетической схемы с последующим сравнением результата расчёта схемы с результатами эксперимента. Проверочная математическая модель может быть построена с использованием системы ОДУ или с помощью метода клеточного автомата [3], или же с использованием стохастической имитационной модели, рассчитываемой с помощью метода Монте-Карло [4]. Чаще всего для каждого конкретного процесса строится самостоятельная модель, в которую заложены все необходимые ограничения. Настоящая работа направлена на создание программного каркаса, предоставляющего

базовые сущности, в терминах которых описывается моделируемая схема.

Каркас программы разработан с использованием техники обобщённого программирования [5] и предоставляет шаблонные структуры данных, комбинируя которые могут быть получены сущности, необходимые для моделирования кинетической схемы процесса роста кристалла. К таким сущностям относятся атомарные структуры различной конфигурации, а так же элементарные стадии, участвующие в процессе роста. Шаблонные структуры данных оснащены алгоритмами поиска самих себя в заданном множестве атомов поверхности кристалла. Множество атомов описывается моделью многокомпонентного решёточного газа, имеющего строгую периодику согласно используемой кристаллической решётке с возможностью выделения атомов, не входящих в состав кристалла, находящихся на границе решётки. Любой атом рассматриваемого множества обязательно имеет связь с другим атомом этого множества. Используемые атомы могут находиться в различных состояниях, обусловленных матрицей состояний, к которой должна быть применена операция транзитивного замыкания. Указанная матрица используется в алгоритмах поиска шаблонных структур данных, предоставляемых разработанным каркасом. Матрица состояний атомов может быть получена в автоматическом режиме с помощью анализатора процесса роста кристалла [2].

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

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

алгоритмической реализации которого выбран динамический метод [6]. Программный каркас включает этот алгоритм, группируя все возможные в системе события по типу, что обеспечивает возможность создания алгоритма, скорость работы которого не зависит от размера системы.

После того, как наиболее вероятное событие выбрано, алгоритм каркаса производит вызов метода, построенного на базе предоставляемых каркасом программных инструментов. Такой метод изменяет валентное состояние взаимодействующих атомов и от каждого измененного атома запускает алгоритм поиска новых образованных структур и элементарных реакций, с ними связанных, а также алгоритм удаления испорченных структур и реакций, от них зависящих. Испорченные атомарные структуры появляются в виду изменения состояния реагирующих атомов, поскольку атом в одном и том же валентном состоянии может принадлежать различным атомарным структурам. Поиск атомарных структур начинается с наименьших структур, определённых на этапе анализа кинетической схемы [2]. После того как атомарная структура найдена, она запускает поиск структур и реакций, зависимых от неё. Если найденная атомарная структура влияет на скорости элементарных процессов, имеющих место вокруг неё, то она изменяет эти процессы, производя соответствующее обновление множества всех возможных событий. Если найденная реакция зависит от близлежащих атомарных структур вокруг реагентов, то производится поиск таких структур. В случае успеха скорость реакции так же изменяется. Такие реакции называются реакциями с латеральными взаимодействиями.

Благодаря тому, что множества структур и реакций кинетической схемы процесса роста кристалла упорядочены в результате анализа [2], а транзитивная матрица состояний атомов позволяет определить принадлежность атома какому-либо состоянию за постоянное время, сложность алгоритма поиска выражена логарифмической зависимостью от количества структур и реакций. Моделирующее приложение,

построенное на базе разработанного программного каркаса [7], обладает высокой производительностью, что обусловлено использованием эффективных моделей и алгоритмов описанных выше.

Результаты моделирования могут быть представлены в виде графиков, отражающих скорость послойного роста кристалла (рис. 1), или же в виде снимков моделируемой поверхности в заданные интервалы времени (рис. 3). Снимки поверхности сохраняются в формате МоШ1е У3000, который может быть визуализирован любой популярной программой для просмотра молекул. После остановки моделирующего процесса предоставляется отчёт о частоте возникновения тех или иных событий за всё время моделирования. Формат полученных данных, их полнота, а также размер исходного кристалла определяются во время запуска моделирующего приложения. Исходная информация о температуре окружающей среды и концентраций газовых реагентов может быть задана посредством конфигурационного файла.

time (s)

Рис. 1. Кривые послойного роста кристалла алмаза размером 20х20 атомов в плоскости (100)

0.12 0.1 # 0.08 2 0.06 о 0.04 0.02 0

0 0 1 0.2 0.3 0 4 0 5 0.6 0.7 0 8 0 9 1 time (s)

Рис. 2. Кривые изменения концентраций атомов, принадлежащих димерам, на поверхности роста кристалла размером 20х20 атомов в плоскости (100)

Рис. 3. Снимки поверхности моделируемого кристалла алмаза размером 20х20 атомов в плоскости (100) в разные моменты времени

На основании разработанного программного каркаса был реализован программный код, соответствующий процессу роста алмазной плёнки. Процесс был представлен 32 элементарными стадиями кинетической схемы, которые использовали 30 атомарных структур, составленных из атомов, число возможных состояний которых равнялось 36. Значения скоростей элементарных процессов были взяты из литературы [7, 8]. В результате моделирования процесса роста кристалла алмаза были получены кривые скорости роста для каждого отдельно взятого слоя (рис. 1, 2), а также снимки поверхности моделируемого кристалла (рис. 3, 4). Снимки поверхности были визуализированы с помощью МагапБрасе [9].

Рис. 4. Термодинамически устойчивая макроструктура - димерный ряд, возникающая в плоскости (100) кристалла алмаза

В кинетической схеме, используемой для получения представленных результатов, имели место реакции с латеральными взаимодействиями, наличие которых обусловило возникновение на поверхности устойчивых макро-структур - так называемых димерных рядов [10]. Ряды димеров представлены на рисунке 4, а на рисунке 3 выражены более тёмными областями, обрамлёнными светлыми полосками атомов водорода. Как видно из рисунка 3,с, димерные ряды изменяют своё направление на каждом следующем растущем слое. Возникновение условий формирования димерных рядов является причиной образования более ровных краёв растущего слоя (рис. 3,ф, а так же, если верить графикам (рис. 1, 2), крайне необходимо для инициализации роста следующего слоя кристалла. Программная модель, построенная на базе каркаса, позволяет проводить исследования влияния скоростей реакций на морфологию и скорость роста кристалла алмаза, а так же выделять наиболее важные элементарные стадии, ответственные за формирование кристаллических структур определённой формы. Наличие программного каркаса позволяет моделировать гораздо более сложные системы, нежели это позволяют сделать таргетированные модели [4].

Разработанный программный каркас [7] позволяет производить нано-масштабное моделирование процесса роста кристалла заданной конфигурации. Исследования выполнены при финансовой поддержке РФФИ в рамках научного проекта 14-07-00960.

Аверчук Глеб Юрьевич, ведущий программист кафедры информационных компьютерных технологий РХТУ им. Д. И. Менделеева, Россия, Москва

Шишов Максим Николаевич, студент третьего курса факультета информационных технологий и управления РХТУ им. Д. И. Менделеева, Россия, Москва

Кольцова Элеонора Моисеевна, д.т.н., профессор, заведующая кафедрой информационных компьютерных технологий РХТУ им. Д. И. Менделеева, Россия, Москва

Куркина Елена Сергеевна, д.ф-м.н., профессор кафедры информационных компьютерных технологий РХТУ им. Д. И. Менделеева, Россия, Москва

Литература

1. Дерягин Б. В., Чураев Н. В., Муллер В. М. Поверхностные силы. - М.: Наука, 1985. - 398 с.

2. Аверчук Г. Ю., Куркина Е. С., Кольцова Э. М. Параллельный алгоритм поиска структур и реакций для моделирования процесса роста монокристаллической алмазной плёнки // Высокопроизводительные параллельные вычисления на кластерных системах. Материалы XIII Всероссийской конференции (Н.Новгород, 14-16 ноября 2013г.) / Под ред. проф. В.П. Гергеля. -Нижний Новгород: Изд-во Нижегородского госуниверситета, 2013, С. 9-13.

3. Аверчук Г. Ю., Кольцова Э. М. Моделирование процесса химического осаждения из газовой фазы монокристаллической алмазной плёнки / Успехи в химии и химической технологии: сб. науч. тр. Т. 25, № 1 (117). М.: РХТУ им. Д.И. Менделеева, 2011. С. 66-72.

4. Netto A., Frenklach M. Kinetic Monte Carlo simulations of CVD diamond growth—Interlay among growth, etching, and migration // Diamond Relat. Mater. 2005. 14. P. 1630-1646.

5. Александреску А. Современное проектирование на С++. - М.: Издательский дом "Виллиамс", 2002. - 336 с.

6. Аверчук Г. Ю., Куркина Е. С. Сравнительный анализ методов Монте-Карло на примере расчета сложной динамики решеточной модели химической реакции // Вестник Московского Университета, Вычислительная математика и кибернетика, 2013. №3. С. 6-17.

7. Репозиторий с исходным кодом программного комплекса для моделирования процесса роста кристалла. URL: https://github.com/newmen/versatile-diamond (дата обращения: 14.05.2014).

8. Skokov S., Weiner B., Frenklach M. Elementary Reaction Mechanism for Growth of Diamond (100) Surfaces from Methyl Radicals // J. Phys. Chem. 1994. 98. P. 7073-7082.

9. Copyright © 1998-2013 ChemAxon Ltd. URL: http://www.chemaxon.com (дата обращения 14.05.2014).

10. Frenklach M., Skokov S. Surface Migration in Diamond Growth // J. Phys. Chem. B. 1997. 101. P. 3025-3036.

Averchuk Gleb Yurievich*, Shishov Maksim Nicolaevich, Koltsova Eleonora Moiseevna, Kurkina Elena Sergeevna

D.I. Mendeleev University of Chemical Technology of Russia, Moscow, Russia.

* e-mail: altermn@gmail.com

FRAMEWORK FOR ATOMICALLY STRUCTURAL MODELING OF CRYSTAL GROWTH USING MONTE-CARLO METHOD

Abstract

Generic programming techniques used in developed framework for performing atomically structural modeling of crystal growth process. Framework provides elemental terms of kinetic scheme, you can implement the code based on them, which will simulate growth process. Implicit kinetic scheme is calculated using a stochastic Monte Carlo method.

Key words: simulation framework, crystal growth, lateral interactions, nano-structures, diamond.