Стохастическая внерешеточная модель релаксации нанесенной наночастицы: линейная временная сложность Текст научной статьи по специальности «Физика»

УДК 544.023.55:519.245

М. П. Красильников

СТОХАСТИЧЕСКАЯ ВНЕРЕШЕТОЧНАЯ МОДЕЛЬ РЕЛАКСАЦИИ НАНЕСЕННОЙ НАНОЧАСТИЦЫ:

ЛИНЕЙНАЯ ВРЕМЕННАЯ СЛОЖНОСТЬ

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

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

Введение

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

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

Однако большинство таких моделей (если не все) - это решеточные модели [3]. Программная реализация стохастической модели для простой кубической решетки и даже для объемно-центрированной решетки не представляет особой сложности, но построение такой модели для более сложных решеток (не говоря уже о решетке произвольной геометрии) является очень сложной алгоритмической задачей.

Также необходимо отметить, что большинство таких моделей (опять же если не все) фактически являются двумерными [4, 5]. Если двухмерность для массивных кристаллов имеет физическое обоснование (наличие выраженной относительно однородной поверхности) [1], то для частиц наномет-рового размера подобные модели вряд ли можно признать реалистичными.

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

узлах решетки. Обычно для моделирования используются сферически-симметричные потенциалы межатомного взаимодействия, несмотря на то, что таким характером взаимодействия невозможно объяснить, например, неплотные упаковки атомов [6]. Применение несферически-симметрич-ных потенциалов является достаточно редким (напр.: [6, 7]).

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

Модель

Несмотря на то, что предложенная нами модель является теоретической, мы придерживаемся параметров, близких к реальным. Это касается и размеров атомов и энергий их взаимодействия. Именно поэтому для моделирования мы выбрали хорошо изученную реакцию окисления O на металлах платиновой группы [8].

Описание модели

В модели заложены четыре сорта атомов: атомы подложки, металла, CO, O. Для каждого сорта атомов задан, в общем случае, неизотропный парный потенциал межатомного взаимодействия.

Каждый атом модели определяется координатами и свободной энергией. Координаты атомов не целые, поскольку предлагаемая модель не является решеточной. В силу того, что потенциалы межатомного взаимодействия не сферически-симме-тричные, для каждого атома модели задана его ориентация. Ориентация атома определяется углами поворота его локальной системы координат от-

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

Рис. 1. Общий вид пространства модели

Атомы массивной инертной подложки расположены в узлах простой кубической решетки и неподвижны. Период решетки может варьироваться от эксперимента к эксперименту и не связан с размерами других атомов модели. Верхний слой атомов подложки расположен в плоскости ОХУ. Вдоль осей ОХ и ОУ наложены циклические граничные условия.

В модели реализованы следующие процессы:

1) диффузия всех атомов модели, кроме атомов подложки;

2) адсорбция молекул СО;

3) диссоциативная адсорбция О2;

4) десорбция любого из атомов модели, кроме атомов подложки;

5) бесконечно быстрая реакция окисления СО.

Алгоритм моделирования

Для моделирования мы используем хорошо известный алгоритм Метрополиса в следующем изложении.

Первоначально в центре подложки формируется аморфная частица так, как это делается в [1012]. Согласно [10], для подготовки исходного образца в заданной области пространства формируется образец со случайным распределением координат атомов. Первоначально на подложке в ее центре располагается один атом, затем к нему случайным образом добавляются атом за атомом до тех пор, пока их число не достигнет требуемого количества (которое является параметром модели). Алгоритм построен так, что он сохраняет только те атомы, которые оказываются на расстояния более двух радиусов от других атомов модели, но в поле их потенциала (с учетом радиуса обрезания).

Далее, в бесконечном цикле стандартным обра-

зом разыгрывается вероятность одного из процессов: диффузии, адсорбции 02, адсорбции СО. Вероятности процессов определяются их частотами, которые являются параметрами модели. Вероятности диффузии задаются известными соотношениями и могут быть определены, например, так:

Р

ш1 =

адс

к1

адс

ка.дс рі + ка.дс єхр

где к1адс - константа адсорбции вещества /, Р, - парциальное давление вещества I в газовой фазе, а АЕ - энергия окружения адсорбированной молекулы [13]. При этом адсорбция кислорода является диссоциативной, а в процессе диффузии может произойти и десорбция любого из атомов модели (кроме атомов подложки), если диффундирующий атом преодолевает барьер десорбции. Энергетический барьер десорбции определяется как минимальная энергия взаимодействия двух атомов, расположенных на расстоянии, равном радиусу обрезания плюс радиус наименьшего из атомов модели.

Реакция окисления СО является бесконечно быстрой и может произойти как в процессе диффузии, так и сразу же после акта адсорбции. В процессе моделирования учитываются все попытки изменения состояния системы. Шагом Монте-Карло считается число таких попыток, равное общему числу атомов системы. Отметим, что в процессе моделирования число атомов системы изменяется.

Опишем подробнее алгоритм каждого из процессов, происходящих в модели.

Диффузия:

- случайным образом выбирается атом модели;

- вычисляется энергия его взаимодействия с соседями, расположенными на расстоянии не более чем радиус обрезания;

- свободная энергия атома распределяется между поступательным движением и вращательным движением;

- случайным образом выбирается вектор перехода, зависящий от энергии поступательного движения;

- случайным образом выбираются углы поворота атома, зависящие от энергии вращательного движения;

- вычисляется энергия взаимодействия данного атома с теми атомами, которые окажутся на расстоянии не больше радиуса обрезания в его новом месте с его новой ориентацией;

- если новая энергия меньше или равна старой, то новое состояние принимается с вероятностью 1;

- если новая энергия больше старой, то это со-

- ЛЕ

стояние принимается с вероятностью е кТ , где АЕ

- разность энергий конечного и исходного состояния системы, Т - температура, k - постоянная Больцмана;

- если образовалась молекула С02 (О и СО оказались рядом), то она учитывается, а составляющие ее атомы удаляются из модели, энергия образования молекулы СО2 передается случайному соседнему атому;

- если атом преодолел барьер десорбции, то он удаляется из модели.

Направление вектора перехода атома выбирается случайным образом, с равномерным распределением вероятности, а его длина тоже является случайной величиной, но распределенной по Максвеллу

Г ™ Л3/2

f (r )dr = 4жг‘

m

2nkT

e 2kT dr,

к - постоянная Бо-

где г - дистанция перехода, льцмана, Т - температура.

На рис. 2 представлено распределение дистанции перехода для температур 300 К и 600 К. При этом параметры распределения подобраны так, что при температуре 300 К наиболее вероятная дистанция перехода равна радиусу атома. Аналогично распределены углы поворота атомов. Наиболее вероятный угол поворота при температуре 300 К -п/2.

Рис. 2. Распределение дистанции перехода атома.

Дистанция перехода задана в радиусах атома

Адсорбция СО

- В пространстве случайным образом выбирается точка и направление движения атома.

- Атом помещается в эту точку и начинает двигаться в выбранном направлении.

- если выбранный луч на расстоянии не больше дистанции свободного пробега (параметр модели)

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

Адсорбция O2

Адсорбция O2 происходит аналогично адсорбции CO за тем исключением, что сразу после акта адсорбции на расстоянии, равном размеру молекулы O2, в случайном направлении размещается еще один атом O. Если местоположение второго атома занято, то адсорбция O2 отвергается, но попытка засчитывается.

Программная реализация модели

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

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

Программа написана на языке Object Pascal, и в ней каждый сорт атомов реализован как класс. Вся совокупность атомов в нашей модели представлена, как линейный массив атомов. Важным является тот факт, что при таком подходе каждый атом модели получает свой уникальный номер, а поскольку атомы определены следующим образом:

Atom = class(TObject, MovingAtom, RotatingA-tom...

property X : real read GetX write SetX;

property Y : real read GetY write SetY;

property Z : real read GetZ write SetZ;

то для каждого атома определены (в том числе) и его координаты.

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

что в ней может разместиться не более одного атома модели. При движении атома в ячейку его старого местоположения записывается -1, а в ячейку нового - его номер в массиве атомов.

Общее представление о методе реализации модели можно получить из рис. 3. На нем изображена взаимосвязь между атомами в линейном массиве и их положением в пространстве. Получив номер атома из линейного массива атомов, мы получаем его координаты, а следовательно, и индексы ячейки пространства. Зная индексы ячейки, в которой расположен выбранный атом и радиус обрезания потенциала, мы можем найти номера атомов, попавших в окрестность выбранного атома. Таким образом, для известного радиуса обрезания мы получаем определенное конечное число атомов, с которыми взаимодействует выбранный атом и к которым мы можем обратиться как к элементам линейного массива по их индексам, полученным из окрестности данного атома. Следовательно, полученный алгоритм будет иметь временную сложность п • 0(1), где 0(1) - функция, определяющая число атомов в сферической окрестности радиуса I; п - общее число атомов модели.

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

Другие параметры модели (такие как размеры атомов, радиус обрезания потенциала, температура, частоты записи результатов эксперимента и др.) передаются в программу через ini-файл.

Программа реализована в виде консольного приложения на языке Object Pascal в среде Delphi 6 серийный номер h5ph-pwmy8-cswef и включена в Реестр программ для ЭВМ [14].

Модельные эксперименты

Релаксация

Переход частицы от неравновесного состояния к равновесному принято называть релаксацией.

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

Несмотря на то, что созданная нами программа позволяет моделировать многие физико-химические процессы, в данной работе мы изучаем только процесс релаксации наночастицы.

Для тестирования модели нами был проведен ряд модельных экспериментов, для чего был использован хорошо изученный сферически-симме-тричный парный потенциал Леннарда-Джонса [16]

Рис. 3. Программное представление модели

Циклические граничные условия, заданные в модели, реализованы за счет естественного переполнения типа ShortInt, что позволяет еще несколько увеличить скорость работы программы, поскольку такой подход не требует дополнительных проверок при переходе атома через границы модельного ящика (model box).

Потенциалы межатомного взаимодействия реализованы в виде внешних функций. Функции, описывающие потенциалы межатомного взаимодействия, имеют следующий прототип:

function P(x0, y0, z0, Rx, Ry, Rz, R, x, y, z : real) : real; stdcall;

где x0, y0, z0 - координаты атома, который создает потенциал; Rx, Ry, Rz - его углы поворота; x, y, z - точка, в которой вычисляется значение потенциала; R - радиус твердого ядра атома.

U (r )= 4s

>12

V

a

V r.

a

Vr

J

со следующим размером твердого ядра о = 0.138 х 10-9, одинаковым для всех атомов системы. Глубина потенциала в точке его минимума для атомов металла была выбрана равной е = 1.65 • 10-14, а для атомов подложки выбиралась равной 0.21 • 10-14,

0.41 • 10-14, 1.65 • 10-14, 6.60 • 10-14 для разных модельных экспериментов.

На рис. 4 представлена исходная аморфная форма наночастицы. Результат моделирования (равновесная форма наночастицы) на рис. 5. Общее число атомов наночастицы - 1000, температура моделирования 600 К.

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

Рис. 4. Начальная аморфная форма

Рис. 5. Равновесная форма. 100 000-й шаг Монте-Карло

выполнение определенного числа шагов моделирования [2] или достижение определенной заранее заданной формы частицы [15]. Мы придерживаемся иной методики, которую впервые использовали в [11]. Суть методики состоит в том, что параллельно проводятся два эксперимента с идентичными параметрами, за исключением исходной формы частиц. В первом эксперименте в качестве начальной формы используется аморфная форма наночастицы, сформированная так, как это описано выше. Во втором атомы исходного кристалла располагаются в узлах простой кубической решетки с параметром решетки, равным двум радиусам атома металла. Далее, визуально мы определяем номер шага Монте-Карло, при котором обе частицы приняли одинаковую форму. Затем, начиная с этого шага Монте-Карло, рассматривая общие энергии обеих наночастиц как две случайные величины, мы находим номер шага Монте-Карло, начиная с которого с 1 %-й достоверностью выполняют следующие три статистические гипотезы:

1. Гипотеза принадлежности обеих случайных величин нормальному распределению, для этого мы использовали следующий критерий [17]

\л\< зЩл) и \Е\< 5ЩЁ) , где Л - асимметрия, Е - эксцесс случайной величины, В(Л), Б(Е)

- их дисперсии соответственно.

2. Гипотеза равенства дисперсий двух случайных величин но Стьюденту [18].

3. Гипотеза равенства средних двух случайных величин с неизвестными, но равными дисперсиями [18].

Номер шага Монте-Карло, начиная с которого впервые выполняются эти три гипотезы, мы и считаем временем наступления равновесия (временем релаксации).

Время выполнения 10 ООО шагов Монте-Карло для наночастицы, состоящей из 1ООО атомов на компьютере с процессором Intel Pentium Dual Е21бО 1.8Ghz и общим объемом ОЗУ в 2Gb, составило на более 1б мин. Для модели с аналогичными параметрами, но реализованной традиционными методами, время выполнения такого же числа шагов Монте-Карло составляло несколько часов [11].

Зависимость равновесной формы наночастицы от энергии взаимодействия с подложкой

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

Результаты моделирования представлены на рис. б-9.

Рис. 6. Почти полное отсутствие смачивания. е = 0.21 • 10-14

Рис. 7. Несколько более сильное смачивание. е = 0.41 • 10-14

Рис. 8. Идентичные атомы подложки и наночастицы. е = 1.65 • 10-14

— 2З —

Рис. 9. Полное смачивание. £ = 6.60 • 10-14

На каждом из рисунков представлен 20 000-й шаг Монте-Карло. В «лупе» в более крупном масштабе изображен левый «край» наночастицы. Отчетливо прослеживается зависимость угла смачивания от глубины потенциальной ямы атомов подложки.

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

На рис. 6-9 б - глубина потенциала в точке минимума.

Выводы

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

Получена равновесная форма нанесенного на инертный носитель нанокристалла, состоящего из 1000 атомов, с потенциалами межатомного взаимодействия Леннарда-Джонса.

Выявлена зависимость равновесной формы нанокристалла от соотношения потенциалов взаимодействия атомов металла и атомов подложки.

Автор считает своей приятной обязанностью выразить свою искреннюю благодарность за ценные замечания по содержанию и стилю статьи Ольге Леонидовне Бандман.

Список литературы

1. Henry C. R. Surface studies of supported model catalysts // Surf. Sci. Rep. 1998. Vol. 31. P. 231-325.

2. Черезов Е. А., Ковалев Е. В., Елохин В. И., Мышлявцев А. В. Статистическая решеточная модель бимолекулярной реакции на динамически изменяющейся поверхности объемно-центрированного кристалла металла // Кинетика и катализ. 2006. Т. 47. № 4. С. 485-497.

3. Zhdanov V. P., Kasemo B. Simulations of the reaction kinetics on nanometer supported catalyst particles // Surf. Sci. Rep. 2000. Vol. 39. № 2-4. Р. 25-104.

4. Елохин В. И., Мышлявцев А. В., Латкин Е. И., Реснянский Е. Д., Шейнин Д. Э., Бальжинимаев Б. С. Статистические решеточные модели физико-химических процессов в каталитических реакциях: автоколебания, адсорбция на шероховатой поверхности, кристаллизация // Кинетика и катализ. 1998. Т. 39. № 2. C. 264-287.

5. Elokhin V. I., Myshlyavtsev A. V. Catalytic processes: nanoscale simulations // Dekker Encyclopedia of Nanoscience and Nanotechnology, 2nd Ed. V. II. C.I. Contescu and K. Putyera, Eds.; Taylor & Francis: New York, 2009. P. 782-793.

6. Баранов М. А., Дубов Е. А., Дятлова И. В., Черных Е. В. Атомно-дискретное описание влияния анизотропных межатомных взаимодействий на упругие свойства ГПУ металлов // ФТТ. 2004. Т. 46. Вып. 2. С. 212-217.

7. Килин В. А., Килин Р. Ю., Зеличенко В. М. Использование модельных потенциалов для расчета электронной структуры фуллерена // Вестн. Томского гос. пед. ун-та (Tomsk State Pedagogical University Bulletin). 2006. Вып. 6(57). С. 170-171.

8. Engle T., Ertl G. Elementary steps in the catalytic oxidation of carbon monoxide on platinum metals // Adv. Cat. 1979. Vol. 28. P. 1-77.

9. Экштайн В. Компьютерное моделирование взаимодействия частиц с поверхностью твердого тела / пер. с англ. М.: Мир, 1995. С. 321: ил.

10. Лагунов В. А., Синани А. Б. Образование биструктуры твердого тела в компьютерном эксперименте // ФТТ. 1998. Т. 40. № 10. С. 1919-1924.

11. Красильников М. П., Мышлявцев А. В. Имитационное моделирование равновесной формы нанесенных нанокристаллов: учет дальних взаимодействий и ковалентных связей // Состояние и освоение природных ресурсов Тувы и сопредельных регионов Центральной Азии. Геоэкология природной среды и общества: научн. тр. ТувИКОПР СО РАН / отв. ред. д.г.-м.н. В. И. Лебедев. Кызыл: ТувИКОПР СО РАН, 2002. С. 110-115.

12. Krasilnikov M. P. Two simple lattice models of the equilibrium shape and the surface morphology of supported 3D crystallites // Bull. Nov. Comp. Center, Comp. Scienc. Novosibirsk, 2008. Iss. 27. P. 63-69.

13. Ковалев Е. В. Статистические решеточные модели процессов на нанесенных катализаторах с учетом изменения морфологии поверхности и формы активных частиц: автореф. дис. ... канд. хим. наук. Новосибирск, 2007. 19 с.

14. Программа для ЭВМ «Программа линейной временной сложности для стохастического внерешеточного моделирования каталитических реакций на активных металлических наночастицах» // Свидетельство о государственной регистрации программы для ЭВМ № 2011612867 внесено в Реестре программ для ЭВМ 8 апреля 2011 г.

15. Combe N., Jensen P., Pimpinelli A. Changing shape in nanoworld // Phys. Rev. Lett. 2000. Vol. 85. № 1. P. 110-113.

16. Гулд Х., Тобочник Я. Компьютерное моделирование в физике: в 2 ч. Ч. 1 / пер. с англ. М.: Мир, 1990. 349 с.: ил.

17. Калиткин Н. Н. Численные методы. М.: Наука, 1978. С. 586.

18. Шметтерер Л. Введение в математическую статистику. М.: Наука, 1976. 520 с.: ил.

19. Resnyanskii E. D., Latkin E. I., Myshlyavtsev A. V., Elokhin V. I. Monomolecular adsorption on rough surfaces with dynamically changing morphology // Chem. Phys. Lett. 1996. Vol. 248. P. 136-140.

Красильников М. П., ведущий инженер.

Тывинский институт комплексного освоения природных ресурсов СО РАН.

Ул. Интернациональная, 117а, Кызыл, Республика Тыва, Россия, 667007.

E-mail: kmp000@gmail.com

Материал поступил в редакцию 08.06.2011.

M. P Krasilnikov

STATISTICAL MODEL OF RELAxATION OF OFF-LATTICE SUPPORTED NANOPARTICLES:

THE LINEAR TIME COMPLExITY

A new method for constructing stochastic models of off-lattice physical and chemical processes occurring at the active metal supported nanoparticles deposited on an inert carrier. The proposed algorithm has linear time complexity.

The created on this base program permits the use, in general, anisotropic interatomic interaction potential. For changes in the potentials of interatomic interaction program does not require recompilation. The equilibrium shape of the Lennard-Jones nanoparticles was obtained. It is shown that there is a dependency of the equilibrium shape of the nanoparticles on the interaction parameters of the substrate and metal atoms.

Key words: relaxation of the nanoparticle, off-lattice model, Monte Carlo, linear time complexity.

Tyva Institute for Exploration of Natural Resources SB RAS.

Ul. Internatsionalnaya, 117a, Kyzyl, Republic of Tyva, Russia, 667007.

E-mail: kmp000@gmail.com