CC BY
72
і-
ФИЗИКО - МАТЕМАТИЧЕСКИЕ НАУКИ
УДК 519.688/544 Д. В. МЫШЛЯВЦЕВ
П. В. СТИШЕНКО
Омский государственный технический университет
МОДИФИКАЦИЯ АЛГОРИТМА МЕТРОПОЛИСА ДЛЯ МОДЕЛИРОВАНИЯ МЕТАЛЛИЧЕСКИХ НАНОЧАСТИЦ
В данной работе мы предлагаем ряд модификаций алгоритма Метрополиса, повышающих его эффективность применительно к задаче о поиске равновесной формы металлических наночастиц. Наша реализация предлагаемого алгоритма, позволила найти равновесную форму наночастицы из 4033 атомов за 2 часа вычислений на современном микропроцессоре при использовании неаддитивных потенциалов.
Ключевые слова: алгоритм Метрополиса, металлические наночастицы, неаддитивный потенциал.
Введение
Уникальные свойства наночастиц (далее НЧ) известны достаточно давно, но обнаруженная в 1989 году повышенная каталитическая активность НЧ золота [ 1 ] вызвала ещё больший интерес к исследованиям в данной области. Благодаря своим особым, по сравнению со сплошными поверхностями, свойствам, НЧ нашли применение в фотонике, электронике, медицине, разработке химических сенсоров [2]. В силу малого размера их исследование экспериментальным
путем весьма сложно и дорого, и не всегда возможно напрямую измерить представляющие интерес показатели и свойства процессов, протекающих с участием НЧ. Поэтому особенно важной частью исследований в этой области является компьютерное моделирование. Огромное количество работ, в которых с помощью компьютерного моделирования изучались НЧ различных типов, подтверждает наличие потребности в таком виде исследований [3]. Обзор этих работ показывает, что большой интерес представляет форма НЧ в состоянии термодинамического
ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 1 (107) 2012 ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 1 (107) 2012
равновесия, в значительной степени определяющая их особые свойства. При моделировании НЧ применяются методы молекулярной динамики, теории функционала плотности (далее ББТ), генетические алгоритмы, а также алгоритмы типа Метрополиса, основанные на использовании метода Монте-Карло, как для расчета равновесного состояния, так и кинетики исследуемых систем. Несмотря на потенциальную мощность методов молекулярной динамики, их применение ограничено системами, быстро достигающими состояния термодинамического равновесия, поскольку время моделирования напрямую зависит от времени протекания процесса в реальности. Поиск состояния с минимальной энергией для НЧ, состоящих из десятков или сотен атомов, выполняется в рамках ББТ. Но интерес представляют НЧ размером до тысяч атомов. В этом случае чаще всего применяется алгоритм Метрополиса и его модификации, иногда в комбинации с генетическими алгоритмами.
С момента опубликования оригинального алгоритма Метрополиса [4] появилось множество его модификаций для различных молекулярных систем. Описание многих из них можно найти в учебниках [5 — 7]. Применение той или иной модификации может в разы изменить скорость моделирования, сделав его доступным для интересующих исследователя систем. Несмотря на большую практику применения алгоритмов типа Метрополиса для моделирования металлических НЧ и других подобных объектов (например, нанопроволок), насколько нам известно, на сегодня нет общепринятого способа реализации алгоритма, который бы использовался для данного типа молекулярных систем. Также нельзя сказать, что существует какой-то программный код, который используется в большинстве таких работ. Обычно приводится лишь беглое описание примененного метода, указывается используемый потенциал межатомного взаимодействия и, иногда, ограничение предельного радиуса взаимодействий. Наблюдается явный недостаток работ с описанием алгоритма достаточно детальным, чтобы можно было оценить его производительность. Это затрудняет оценку практической применимости других модификаций, поскольку неясно, в комбинации с какими структурами данных и алгоритмическими приемами их предстоит использовать.
В данной работе мы приводим описание модификации алгоритма Метрополиса для НЧ с детальным описанием ключевых особенностей алгоритма, оказывающих наибольшее влияние на его производительность. В предлагаемом алгоритме используются как известные подходы к ускорению алгоритма, так и разработанные нами приёмы, которые, насколько нам известно, ранее не описывались в литературе.
Алгоритм Метрополиса для металлических наночастиц
Классический алгоритм Метрополиса состоит из следующих шагов:
1. Взять систему в некотором случайном состоянии.
2. Обозначить текущее состояние системы как г. В соответствии с базовой матрицей марковского процесса выбрать некоторое состояние /.
3. Вычислить разницу энергии системы ДЕ// в состояниях г и /.
4. С вероятностью ехр( — ДЕ//квТ), где кв — постоянная Больцмана, а Т — температура системы, перевести систему в состояние /. Перейти к пункту 2.
Каждый шаг этого алгоритма можно реализовать различными способами, в зависимости от требований, предъявляемых к алгоритму, и характеристик исследуемых систем и процессов. Очевидно, что на первом шаге не много возможностей по модификации алгоритма. Второй шаг алгоритма полностью определяется базовой матрицей марковского процесса. Далее попытку перехода из текущего состояния в новое будем называть пробным шагом.
В решеточной модели возможны три варианта пробных шагов для моделирования диффузии атомов: прыжок в случайную точку, прыжок в соседнюю ячейку (в ближайшую или находящуюся в пределах определенной координационной сферы) и прыжки между случайными элементами списков активных атомов и вакансий. В первом случае велика вероятность попытки выполнить переход в ячейку решетки не имеющую рядом с собой никаких атомов. Поскольку вероятность успеха такой попытки очень мала, это сильно снижает долю успешных шагов. В случае прыжков в пределах некоторой сферы шанс на заведомо неудачную попытку гораздо меньше, но ограничение максимального расстояния перемещения атома может осложнить обход фазового пространства из-за «застреваний» в областях локальных минимумов энергии или метастабильных состояний. Используя динамические списки активных атомов и вакансий можно избежать недостатков обоих подходов. При таком подходе нужно учитывать, что ведение этих списков требует дополнительных вычислительных затрат, которые тем чувствительнее, чем выше процент успешных шагов, поскольку актуализация списков должна производиться после каждого успешного шага. Кроме того, в этом случае важно следить за соблюдением условия детального баланса. Поскольку число активных атомов и вакансий меняется после каждого успешного шага, базовая матрица теряет свою симметричность. Чтобы компенсировать это нарушение детального баланса, на шаге 4 используется смещенная вероятность принятия пробного шага. Например, в работе [8] для этого применялось следующее выражение:
p = min
NaNb
N^N*
exp
- АЕ
(1)
где Nai(fj — количество поверхностных атомов до (i) и после (f) пробного шага, Nbjff) — соответствующее количество активных вакансий, DE — изменение энергии системы, вызываемое пробным шагом, kB — постоянная Больцмана, T — температура системы.
Насколько нам известно, сегодня нет общепринятого и описанного в литературе способа ведения списков активных атомов и вакансий. В тоже время при моделировании ведется активная работа с этими списками и выборка эффективного контейнера для этих целей значительно влияет на скорость моделирования. При моделировании выполняются следующие операции со списками активных атомов и вакансий:
— выбор случайного элемента на каждом пробном шаге;
— вставка и удаление элементов после каждого успешного шага.
Массив, дерево или список имеют сложность не меньше 0(log N) как минимум для одной из этих операций. Для решеточной модели возможно обеспечить константную сложность всех операций. В нашей реализации алгоритма этот контейнер состоит из двух массивов R и L и индекса TAIL. Массив R
2
kgT
имеет заведомо достаточную длину, чтобы вместить все активные атомы (вакансии). Массив L является развёрткой кристаллической решётки. Каждый элемент массива L содержит либо индекс элемента в массиве R, либо —1, означающее, что в данной ячейке решётки нет активного атома. Массив R содержит индексы элементов массива L, которые могут быть легко преобразованы в декартовы координаты. При этом соблюдается условие, что R[i]=j, если и только если L[j] = i. В массиве R элементы, ссылающиеся на активные атомы (вакансии) решётки, расположены непрерывно от начала массива и до элемента, на который указывает TAIL. Остальные элементы массива не используются. Когда в ходе моделирования число активных атомов меняется, TAIL сдвигается в ту или другую сторону.
В предлагаемой структуре данных выбор случайного элемента списка осуществляется генерацией случайного числа от 1 до TAIL, а удаление и вставка по заданным координатам выполняется сначала в массиве L, с последующей коррекцией массива R и индекса TAIL. Очевидно, что все требуемые операции в таком контейнере имеют сложность 0(1).
На следующем шаге алгоритма Метрополиса, вычисляется изменение энергии системы, вызываемое выбранным ранее пробным шагом. Расчет изменения потенциальной энергии в результате пробного шага это самый или, один из самых затратных с вычислительной точки зрения шагов алгоритма, поэтому ему необходимо уделять особое внимание.
Вычисление энергии системы может осуществляться с использованием различных эмпирических потенциалов межатомного взаимодействия, которые выбираются в зависимости от характеристик моделируемой системы, требуемой точности, от исследуемых явлений, и доступных вычислительных ресурсов. В некоторых случаях можно использовать простые парные потенциалы, такие как потенциал Лен-нарда-Джонса. Несмотря на свою простоту, такие потенциалы успешно используются для моделирования благородных газов, и даже некоторых свойств металлических НЧ [9].
С точки зрения компьютерного моделирования, важной особенностью таких потенциалов является их аддитивный характер. Это дает возможность легко получить изменение энергии системы, вызванное пробным шагом. Выражение полной энергии в общем виде для таких потенциалов выглядит так:
Ет = II Е(г/),
(2)
£Ч /єЧ;
где Y — множество всех атомов системы, Y — множество атомов в окрестности атома i, r. — расстояние между атомами i и j, E — функция энергии парного взаимодействия. Здесь и далее под окрестностью атома понимается пространство в пределах некоторого расстояния называемого радиусом отсечки. Очевидно, что изменение ET можно вычислить как разницу между E(rj}) для тех i и j, для которых изменилось расстояние rij.
Присущая аддитивным потенциалам модель парных взаимодействий плохо описывает такие свойства металлов как константы упругости, энергия формирования вакансии, реконструкция и релаксация поверхности. Для исследования в этих областях применяются многочастичные потенциалы, теряющие свойство аддитивности: метод погруженного атома (EAM/MEAM), tight-binding (TB) потенциалы (Гупты [10], TB-SMA [11], Саттона-Чена [12]). Они учиты-
вают непарный характер взаимодействий в металлах и в настоящее время широко применяются как при моделировании методом Монте-Карло, так и с методами молекулярной динамики. TB потенциалы хорошо зарекомендовали себя именно в моделировании наноразмерных структур, а именно: НЧ, нанопроволок, длинных атомных цепочек (LAC).
Общая энергия системы для TB модели вычисляется по формуле
E = I I Rj -II Bi
/ЄЧІ
/єЧі
(З)
Здесь и — вклад отталкивания и притяжения
в общую энергию взаимодействия между атомами I и]. Выражения для Я^ и Б^ у разных потенциалов отличаются. В случае ТБ-БМА они выглядят следующим образом:
Ri = Aie-pA/ го-1), b, = /Г'/(г // го-«,
(4)
(5)
где А, р, — параметры потенциала, определяемые типом атомов і и /. В решеточной модели удобно задавать потенциал в виде таблицы значений Я/ и В/ для различных типов атомов и расстояний т/, а не пользоваться формулами (4) и (5) напрямую.
Формулу (3) можно переписать в следующем виде
E = II(R(rt))-У IB/)
іЄЧ /єЧ;
(б)
ієЧ \ jєЧ.
Отсюда видно, что для составляющей притяжения при перемещении 1-го атома необходимо пересчитать
не только корень
Z B(rij), но и корни /X В(гц)
j-f, V j-fi
для всех i таких, что 1 е f. Если |f,| = M , то таких корней также будет M. Это означает, что если для парных потенциалов расчет изменения энергии имеет сложность O(M), то для TB потенциалов — 0(M2).
Насколько нам известно, ранее в литературе эта особенность многочастичные потенциалов явно не обсуждалась, хотя в целом считается, что их вычислительная стоимость выше, чем у парных [13].
Если при моделировании металлических НЧ не применяются специальные техники для повышения процента успешных шагов, а также если используется решеточная модель расположения атомов, то этот процент будет очень низким на протяжении большей части времени моделирования. Благодаря этому имеет смысл хранить результаты промежуточных вычислений энергии системы для ускорения последующих вычислений. Такой прием называется мемоизацией. Если для каждого i-го атома хранить
актуальное значение суммы X В(г^), то при переме-
j-f,
щении 1-го атома можно вычислить её изменение как разницу B(rJ—B(r) Очевидно, что в отличие от полного пересчета суммы, эта операции имеет сложность не O(M), а 0(1). Это делает общую сложность вычисления изменения энергии не 0(M2), а O(M), то есть такой же, как и у аддитивных потенциалов. Конечно, необходимо поддерживать актуальность
суммы X В(тц) для всех атомов частицы, что требует
j-f,
гєЧ
ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 1 (107) 2012 ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК №1 (107) 2012
*
Рис. 1. Нанесенная наночастица платины из 4033 атомов, полученная с помощью предлагаемой модификации алгоритма Метрополиса (светлые сферы обозначают атомы платины, тёмные — атомы подложки)
дополнительных затрат. Поэтому не следует ожидать ускорения вычислений в М раз. В разработанной нами реализации алгоритма Метрополиса при моделировании в рамках решеточной модели платиновой НЧ применение такой мемоизации позволило ускорить процесс моделирования примерно в 15 раз (для М = 78).
Реализация алгоритма
С применением описанных выше оптимизаций и техник реализация алгоритма Метрополиса для моделирования нанесенных НЧ выглядит следующим образом:
1. Заполнить модельную решётку атомами частицы и подложки; рассчитать и Б{ для г е ¥ ; составить списки активных атомов и вакансий.
2. Выбрать случайные элементы из списков активных атомов и вакансий, посчитать суммарную энергию Е1 атомов, находящихся в окрестностях этих двух точек.
3. Энергия взаимодействия выбранного атома вычитается из аддитивных составляющих энергии в ячейках его соседей.
4. Атом перемещается в выбранную вакансию, а значения аддитивных составляющих энергии его новых соседей увеличиваются в соответствии с формулами (4) и (5).
5. Вычисляется суммарная энергия Еа атомов, находящихся в окрестностях начальной и конечной точек расположения атома.
6. Вычисляется разница АЕ=Е—Е1, которая используется для принятия или отклонения данного перехода в соответствие с распределением Больцмана.
7. Если пробный переход был принят, то ближайшие соседи исходной и конечной ячеек данного перехода, при необходимости, добавляются или удаляются из списков активных атомов и вакансий. Иначе — положение атома и значения аддитивных составляющих энергии вокруг атома и вакансии возвращаются в исходные состояния.
8. Если не достигнуто состояние термодинамического равновесия, то алгоритм возвращается к пункту 2.
С использованием данной реализации мы определили равновесную форму свободных и нанесенных
частиц платины размером из 4033 атомов при температуре 1162K. На процессоре Intel Xeon X5472 3 GHz реализованный алгоритм выполнял около 90 000 пробных шагов в секунду, из них около 500 успешных. Полученная нанесённая НЧ изображена на рис. 1.
Для достижения состояния термодинамического равновесия частиц потребовалось до 2000 шагов Монте-Карло, на каждом из которых предпринималось 72000 попыток (определяется размером решетки) выполнить акт диффузии. Чтобы убедиться, что частица пришла в состояние термодинамического равновесия, моделирование продолжалось до 5000 шагов. Весь процесс моделирования занимал не более двух часов, что очевидно можно считать приемлемыми вычислительными затратами при проведении исследований. Приведенный алгоритм успешно использовался нами в работе [ 14] для вычисления равновесной структуры адсорбирующей поверхности.
Заключение
Предлагаемая реализация алгоритма Метрополиса для моделирования свободных и нанесенных НЧ расширяет границы применимости алгоритма, позволяя за несколько часов процессорного времени решить задачу о поиске равновесной формы частиц из нескольких тысяч атомов. Требуемая производительность алгоритма была достигнута за счёт разработанной структуры данных для ведения списков активных атомов и вакансий с константной сложностью всех требуемых операций, а также благодаря предложенному приёму мемоизации при вычислении изменения энергии системы с помощью неаддитивных потенциалов.
Библиографический список
1. Haruta M., Yamada N., Kobayashi T. and Iijima S. Gold catalysts prepared by coprecipitation for low-temperature oxidation of hydrogen and of carbon monoxide // Journal of Catalysis. 1989, Vol. 115, pp. 301-309.
2. Metal nanoparticles: synthesis, characterization, and applications / ed. by Feldheim D.L., Foss C.A. - New York: Marcel Dekker, 2001. — 338 p.
3. D. J. Wales, J. P. K. Doye, A. Dullweber, M. P. Hodges, F. Y. Naumkin F. Calvo, J. HernSndez-Rojas and T. F. Middleton. The Cambridge Cluster Database, URL: http://www-wales.ch. cam.ac.uk/CCD.html .
4. Metropolis N., Ulam S. The Monte Carlo method // Journal of the American Statistical Association. 1949, Vol. 44, pp. 335 — 341.
5. Frenkel D. Understanding molecular simulation: from algorithms to applications / Frenkel D., Smit B. - San Diego: Academic Press, 2002. — 638 p.
6. Landau D.P., and Binder K. A Guide to Monte Carlo Simulations in Statistical Physics. — 2nd ed. — Cambridge: Cambridge University Press, 2000. — 449 p.
7. Leach A.R. Molecular Modelling: Principles and Applications. — 2nd ed. — Prentice Hall, 2001. — 773 p.
8. McKenna K.P., Sushko P.V. and Shluge A.L. Transient Atomic Configurations of Supported Gold Nanocrystallites at Finite Temperature // Journal of Physical Chemistry C, 2007, Vol. 111, pp. 2823-2826.
9. Mbller M., Albe K. Lattice Monte Carlo simulations of FePt nanoparticles: Influence of size, composition, and surface segregation on order-disorder phenomena // Physical Review B, 2005, Vol. 72, pp. 094203.
10. Gupta R. P. Lattice relaxation at a metal surface // Physical Review B, 1981, Vol. 23, pp. 6265-6270.
11. Cleri F., Rosatto V. Tight-binding potentials for transition metals and alloys // Physical Review B, 1993, Vol. 48, pp. 22-33.
12. Sutton A.P., Chen J. Long-range Finnis-Sinclair potentials // Philosophical Magazine Letters, 1990, Vol. 61, pp. 139-164.
13. Germann T. C., Kadau K., Lomdahl P. S. 25 Tflop/s Mul-tibillion-atom Molecular Dynamics Simulations and Visualization/ Analysis on BlueGene/L // Proceedings of the 2005 ACM/IEEE Conference on Supercomputing. 2005.
14. Myshlyavtsev A., Stishenko P. Monte Carlo model of CO adsorption on supported Pt nanoparticle // Applied Surface Science. 2010. Vol. 256, no. 17. Pp. 5376-5380.
МЫШЛЯВЦЕВ Александр Владимирович, доктор химических наук, профессор кафедры «Химическая технология переработки углеводородов», проректор по учебной работе.
СТИШЕНКО Павел Викторович, младший научный сотрудник кафедры высшей математики.
Адрес для переписки: pvstishenko@omgtu.ru
Статья поступила в редакцию 30.08.2011 г.
© А. В. Мышлявцев, П. В. Стишенко
УДК 51987 М. А. БАРЫШНИКОВ
А. А. КОЛОКОЛОВ
Омский государственный технический университет
Омский филиал Института математики им. С. Л. Соболева СО РАН
РЕШЕНИЕ НЕКОТОРЫХ ЗАДАЧ РАЗВОЗКИ НЕФТЕПРОДУКТОВ С ИСПОЛЬЗОВАНИЕМ ДИСКРЕТНОЙ ОПТИМИЗАЦИИ_______________________________
В работе рассматриваются задачи оптимальной развозки нефтепродуктов. Для их решения построены модель дискретной оптимизации на графе и модель целочисленного линейного программирования, разработаны и реализованы алгоритмы муравьиной колонии, а также гибридный алгоритм с применением процедур локального поиска. Приведены результаты вычислительного эксперимента с реальными исходными данными, которые показали эффективность предложенного подхода к решению указанных задач.
Ключевые слова: задача оптимальной развозки, дискретная оптимизация, целочисленное программирование, дискретная оптимизация, алгоритм муравьиной колонии, локальный поиск.
Введение
Решение задач оптимальной развозки продукции [1, 2] необходимо для повышения эффективности работы транспортных структур крупных предприятий. Данная проблема особенно актуальна для организаций, в которых транспортные перевозки являются одним из основных видов деятельности [3]. Многие из указанных задач относятся к области развозки нефтепродуктов [4].
Рассматриваемые задачи, как правило, являются ИР-трудными и на практике имеют большую размерность, вследствие чего нахождение точного решения представляется весьма трудоемким [5]. Использование эвристических алгоритмов позволяет получать приближенные решения с относительно небольшой погрешностью за приемлемое время [6].
В данной работе построены и реализованы эвристические алгоритмы, основанные на алгоритме муравьиной колонии [7, 8]. Приводятся результаты вы-
ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 1 (107) 2012 ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
