Научная статья на тему 'Реализация алгоритма оптимизации параметров молекулярно-динамического потенциала ReaxFF'

Реализация алгоритма оптимизации параметров молекулярно-динамического потенциала ReaxFF Текст научной статьи по специальности «Математика»

CC BY
115
219
Поделиться
Область наук
Ключевые слова
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ / МОЛЕКУЛЯРНАЯ ДИНАМИКА / ПОТЕНЦИАЛ ВЗАИМОДЕЙСТВИЯ / ХИМИЧЕСКИ РЕАКТИВНЫЕ СИСТЕМЫ / РЕАКТИВНОЕ СИЛОВОЕ ПОЛЕ / ОПТИМИЗАЦИЯ ПАРАМЕТРОВ / ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ

Аннотация научной статьи по математике, автор научной работы — Шефов К. С., Степанова М. М.

Молекулярно-динамические методы с потенциалом ReaxFF (Reactive Force Field) позволяют получать достаточно хорошие результаты при моделировании больших многокомпонентных химически реактивных систем. Реализация потенциала ReaxFF имеется в программном пакете LAMMPS, но для корректных расчетов требуется предварительная оптимизация многочисленных параметров этого потенциала для конкретной системы. В данной работе представлены результаты разработки параллельной программы, реализующей алгоритм однопараметрической оптимизации параметров молекулярно-динамического потенциала ReaxFF. Для оптимизации используются данные точных квантовых расчетов набора простых моделей химических соединений (training set). Оптимизация основана на минимизации функции ошибки, представляющей собой сумму квадратичных отклонений расчетов ReaxFF от квантовых расчетов по всему набору. В статье описан алгоритм поиска параметров и представлена блок-схема разработанной программы. Для процедуры поиска минимума потенциальной энергии моделей оптимизирующего набора с потенциалом ReaxFF используется молекулярно-динамический симулятор LAMMPS, собранный в качестве библиотеки. Программа позволяет подбирать параметры потенциала для произвольных типов химически реактивных систем, что необходимо для создания корректных численных моделей химических реакций методами молекулярной динамики c помощью пакета

Похожие темы научных работ по математике , автор научной работы — Шефов К.С., Степанова М.М.,

Текст научной работы на тему «Реализация алгоритма оптимизации параметров молекулярно-динамического потенциала ReaxFF»

УДК 519.6 Дата подачи статьи: 28.02.2014

РЕАЛИЗАЦИЯ АЛГОРИТМА ОПТИМИЗАЦИИ ПАРАМЕТРОВ МОЛЕКУЛЯРНО-ДИНАМИЧЕСКОГО ПОТЕНЦИАЛА ReaxFF

К.С. Шефов, аспирант; М.М. Степанова, к.ф.-м.н., доцент

(Санкт-Петербургский государственный университет, Университетская наб., 7-9, г. Санкт-Петербург, 199034, Россия, k. s. shefov@gmail. com, mstep@mms. nw. ru)

Молекулярно-динамические методы с потенциалом ReaxFF (Reactive Force Field) позволяют получать достаточно хорошие результаты при моделировании больших многокомпонентных химически реактивных систем. Реализация потенциала ReaxFF имеется в программном пакете LAMMPS, но для корректных расчетов требуется предварительная оптимизация многочисленных параметров этого потенциала для конкретной системы. В данной работе представлены результаты разработки параллельной программы, реализующей алгоритм однопараметрической оптимизации параметров молекулярно-динамического потенциала ReaxFF. Для оптимизации используются данные точных квантовых расчетов набора простых моделей химических соединений (training set). Оптимизация основана на минимизации функции ошибки, представляющей собой сумму квадратичных отклонений расчетов ReaxFF от квантовых расчетов по всему набору. В статье описан алгоритм поиска параметров и представлена блок-схема разработанной программы. Для процедуры поиска минимума потенциальной энергии моделей оптимизирующего набора с потенциалом ReaxFF используется молекулярно-динамический симулятор LAMMPS, собранный в качестве библиотеки. Программа позволяет подбирать параметры потенциала для произвольных типов химически реактивных систем, что необходимо для создания корректных численных моделей химических реакций методами молекулярной динамики c помощью пакета LAMMPS на высокопроизводительных параллельных кластерах.

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

IMPLEMENTATION OF THE PARAMETER OPTIMIZATION ALGORITHM FOR MOLECULAR DYNAMIC FORCE FIELD ReaxFF

Shefov K. S., Postgraduate Student; Stepanova M. M., Ph.D. (Physics and Mathematics), Associate Professor (St. Petersburg State University, Quay Universitetskaya 7-9, St. Petersburg, 199034, Russian Federation,

k. s. shefov@gmail. com, mstep@mms. nw. ru)

Received 28.02.2014

Аbstract. Molecular dynamic methods that use ReaxFF (Reactive Force Field) allow obtaining sufficiently good results in simulating large multicomponent chemically reactive systems. ReaxFF is implemented in the standard software package LAMMPS. However, pre-optimization of numerous parameters of this force field for a specific system is required for correct computations. This paper presents a parallel program that implements one-parameter search algorithm for the molecular dynamic force field ReaxFF. The program uses data obtained by accurate quantum chemical computations of characteristics (training set) of chemical compounds simple models for the optimization process. The parameter search is based on the minimization of the error function. This is a sum of squared deviations of training set entries computed by ReaxFF from ones computed by quantum chemical methods. The paper describes the parameter search algorithm and shows the program block-scheme. Molecular dynamic simulator LAMMPS is used to search the minimum of ReaxFF potential energy (optimal geometry) of the training set models. The package LAMMPS is compiled as a library and linked to the program. The program allows searching ReaxFF parameters for arbitrary types of chemically reactive systems. It is necessary for specific numerical simulations of chemical reactions using the molecular dynamic package LAMMPS on high-performance parallel computers.

Keywords: numerical simulation, molecular dynamics, reactive force field, chemically reactive systems, parameter optimization, parallel algorithm.

Численное моделирование химически реактивных молекулярных систем актуально для целого ряда прикладных научных задач. Для небольших химических соединений (до 100 атомов) это можно делать методами квантовой механики (ab initio), однако для моделирования более крупных структур данные методы непригодны из-за высокой ресурсоемкости. Традиционные молекулярно-динамические методы, которые обычно используются для крупномасштабных систем, как правило, основаны на классическом потенциале силовых полей и не позволяют хорошо воспроизводить многокомпонентные химически реактивные системы. Поэтому в своей работе, связанной с изучением систем хранения водорода, для моделирования деструкции гидридов легких металлов

с выделением водорода мы используем комбинированный подход. Он базируется на потенциале реактивных силовых полей ReaxFF (Reactive Force Field) [1, 2], параметризация которого выполняется на основе точных квантовых расчетов для малых кластеров. Отметим, что потенциал ReaxFF изначально разрабатывался для описания химической реактивности, диссоциации и формирования химических связей, дефектов, поверхностных эффектов и др. Такой подход позволяет на современных высокопроизводительных параллельных вычислительных кластерах моделировать развивающиеся во времени системы с числом атомов порядка 106-108.

Программная реализация потенциала имеется в ПО ADF/BAND (лицензионный платный па-

кет) [3] и LAMMPS (свободно распространяемый пакет) [4]. Для адекватных расчетов с использованием ReaxFF применительно к системе конкретного типа (например углеводороды или гидриды металлов) предварительно необходимо осуществить достаточно трудоемкий процесс поиска/оптимизации многочисленных параметров этого потенциала, при которых он будет корректно описывать процессы, происходящие в этой системе. Однако в отличие от самого потенциала ReaxFF программ оптимизации его параметров в открытом доступе нет. Поэтому, опираясь на описание метода однопараметрического поиска [5] и реализацию потенциала ReaxFF в LAMMPS, авторами был разработан свой параллельный алгоритм оптимизации параметров, описание которого представлено в данной работе.

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

ВО, = ^ exp

Pi i

г r V2 j

Л

где индекс I равен о, п, пп (типы ковалентной связи); Гу - расстояние между атомами; рь pl2, rl - параметры [2]. Порядки связей, вычисляемые при численном моделировании для всех атомов непосредственно из межатомных расстояний и обновляемые на каждом шаге по времени, позволяют описать образование и диссоциацию связей во время моделирования. Общий вид потенциала [2] следующий:

EReaxFFi{rij), {>>}, {j}, {q,}, {BO,j})=

Ebond+Eover +Eunder +Eval+Epen+Ecoa+Etors+ +Econj+Ehbond+Elp+ EvdWaals+ECoulomb. (1)

Каждое слагаемое в правой части формулы отвечает за отдельный тип взаимодействия: кова-лентное (характеризуется порядком химических связей, слагаемые Ebond, Eover, Eunder), трехчастичное (валентные углы, Eval, Epen, Ecoa), четырехчастич-ное (двугранные валентные углы, Etors, Econj), ку-лоновское (Ecouiomb), ван-дер-ваальсовое (EvdWaais), водородные связи (Ehbond) и энергия неподеленных электронных пар (Elp).

Полная энергия является суммой функций относительных положений пар атомов r,j, троек rijk и четверок rijkl. Более того, каждая из этих функций зависит от некоторого числа параметров, которые определяются тем, какие химические элементы взаимодействуют. Таким образом, весь потенциал зависит от многочисленных (в зависимости от сложности соединения их количество может быть более 100) параметров: EReaxFF=f(p\, Р2, Рз, ■■■, Pn).

Анализ подробного вида всех слагаемых в формуле (1) позволяет сделать вывод, что EReaxFF и его первая производная являются непрерывными функциями координат атомов даже при химических реакциях [2].

Метод однопараметрической оптимизации. Реализация алгоритма оптимизации параметров основана на методе однопараметрического поиска [5]. Данные, по которым производится поиск (оптимизирующий набор, training set), получаются из расчетов простых моделей химических соединений (до 50 атомов) методами квантовой химии или же берутся из эксперимента. Данными оптимизации являются геометрические характеристики молекул (расстояния между атомами, валентные углы, двугранные углы), колебательные частоты, эффективные заряды атомов и энергии связи соединений, а также какие-либо другие характеристики, которые можно вычислить методами молекулярной динамики, а затем сравнить с данными оптимизирующего набора. Процедура однопараметрического поиска состоит в изменении параметров потенциала для минимизации суммы квадратов:

Ёгг (р,р,р, ...,pN) =

•± [( х

'(Pi, Р, Рз >•••> Pn )) 1 ]

(2)

Здесь x,QC,L't - данные оптимизирующего набора (конкретные характеристики модели или молекулы, полученные методами квантовой химии или из литературы (эксперимента)); x,ReaxFF - вычисленные с использованием ReaxFF значения, соответствующие этим характеристикам; о, - критерий приемлемости, который имеет размерность x, и характеризует вклад разности x1QC,L't-xReaxFF в сумму квадратов. Сумма (2) представляет собой так называемую функцию ошибки - отклонение расчетов ReaxFF от данных оптимизирующего набора. Каждое x,ReaxFF зависит от параметров потенциала р, . ., pN, поэтому и Err является функцией этих параметров.

Важно отметить, что для вычисления значения

ReaxFF

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

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

мального значения параметра необходимо вычислить Err при трех различных значениях: pk, Pt(1-Sk), pt(l+5k), где Sk<<1 - множитель смещения. Для построенного по этим точкам полинома второй степени методом экстраполяции производится поиск значения pkopt, соответствующего минимальному значению полинома на интервале [pk(l-Ak), Pt(1+Afc)], где 5k<Ak<<1 - множитель экстраполяции.

Одна итерация программы состоит в выполнении описанной процедуры один раз для каждого pk из списка. Так как большинство параметров в силовом поле ReaxFF связаны между собой, оптимальное значение каждого параметра сдвигается каждый раз, когда другой параметр изменяется. Это означает, что процесс оптимизации должен быть повторен до тех пор, пока Err не перестанет уменьшаться, что определяется критерием сходимости е. Отметим, что для каждого параметра pk должен быть задан допустимый интервал изменения [Dkmm, Dkmax], чтобы избежать сходимости в нереалистичную область.

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

Описание программы оптимизации параметров потенциала ReaxFF. Блок-схема алгоритма разработанной программы показана на рисунке. Общая структура программы включает инициализацию, чтение входных данных, два вложенных цикла, запись промежуточных и выходных результатов. Внешний цикл организует итерации по полному набору параметров, а внутренний обеспечивает процедуру экстраполяции по параметру pk. Внутренний цикл содержит расчеты с потенциалом ReaxFF всех моделей оптимизирующего набора и является весьма ресурсоемким. Он распараллелен с помощью технологии MPI, что позволяет запускать программу на вычислительных узлах кластера.

На первой стадии программа инициализирует три процесса и производит чтение из файлов входных параметров: начальных значений пара-метровpi0, ..., pN0, величин множителей смещения Si, ..., SN, величин множителей экстраполяции Ai, ..., An, допустимых границ изменения параметров потенциала Dkmm и Dkmax, k=1, ..., N, максимального числа итераций M и критерия сходимости е. Также обрабатываются входной файл моле-кулярно-динамического симулятора LAMMPS и файлы моделей с их характеристиками x,QC/Lli и начальными (взятыми из квантово-химических расчетов) конфигурациями атомов.

Внешний цикл организует итерации по набору параметров. За один шаг этого цикла обеспечивается проход по всему списку pi, ..., pN, и для каждого pk выполняется одна процедура экстраполяции. Результат одной итерации - новый набор

{pk }, который будет использоваться в следующей итерации. Цикл завершается, если на очередном шаге t будет удовлетворено условие сходимости \Err-Errt-1 |<е или достигнут лимит количества итераций M.

Внутренний цикл, предназначенный для самой трудоемкой процедуры экстраполяции, выполняется параллельно. На этом этапе каждый из трех процессов получает свое значение параметра (pk, pt(1-5k), pt(1+5k)) и вычисляет для него значение функции ошибки (2) (на рисунке Err-, Err0 и Err+). Для этого используется вызов симулятора

LAMMPS с реализацией потенциала ReaxFF, собранного в качестве библиотеки. С помощью библиотеки LAMMPS происходит поиск минимума потенциальной энергии (1) моделей оптимизирующего набора при заданных значениях параметров как функции координат составляющих их атомов. Отметим, что алгоритмы, входящие в состав библиотеки LAMMPS, сходятся зачастую не к минимуму, а к седловой точке некоторого порядка. Поэтому к библиотечным методам минимизации пакета LAMMPS авторы добавили свой модуль для вычисления матрицы Гессе (второй производной потенциала по координатам атомов). Его использование помогает избежать сходимости в седловые точки и гарантирует попадание в минимум. После этого вычисляются характеристики

ReaxFF

моделей оптимизирующего набора x,ReaxFF и значения функции (2) Err-, Err0 и Err+. На следующем шаге алгоритма процесс 0 получает от процессов 1 и 2 вычисленные ими значения Err' и Err+ и выполняет процедуру параболической экстраполяции (Fextr на рисунке). Полученное оптимальное значение параметра pk записывается в журнал и выходной файл параметров потенциала. Последовательность повторяется для всех N параметров в списке, после чего последнее оптимальное значение функции ошибки записывается в память как

Errt.

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

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

Структура пакета. Исходный код программы разработан на языке C и имеет модульную структуру. Программный пакет собран и доступен в виде архива с исходными кодами, инструкцией и сценарием компиляции, документацией для пользователя и примерами входных файлов.

Тестирование и использование. Функциональность программы проверялась на примере оптимизации параметров алюминия и его гидридов. За начальное приближение при оптимизации были взяты параметры, включенные в дистрибутив пакета LAMMPS, которые следует считать довольно грубыми, поскольку они дают плохие результаты для решеток металла и его гидридов.

Для оптимизации параметров мы использовали набор, в который вошли кластеры чистого алюминия (Al2 - Al8, AlX3), кластеры гидрида (AlnH3n, n=l, ..., 8), а также кристаллы а- ß- и у-модификаций гидрида алюминия [6, 7], порожденные из ячеек с периодическими условиями. Добавление в набор бесконечных кристаллических структур позволило учесть уравнение состояния, то есть зависимость энергии кристаллов от объема ячейки. В общей сложности набор включал 56 структур.

Всего в оптимизации участвовало 79 параметров потенциала. Одна итерация процедуры оптимизации (по одному параметру) на вычислительном кластере занимала в данном случае около 4 часов. За одну итерацию величина ошибки уменьшалась на значение от 0,1 до l % по сравнению с предыдущей итерацией. Программа совершила около 4 000 итераций, за это время величина ошибки сократилась приблизительно в 500 раз.

Для проверки качества полученных параметров было выполнено моделирование кристаллических решеток гидридов алюминия при начальных и оптимизированных параметрах с целью сравнения их стабильности, плотности и регулярности структуры при комнатной температуре. Все исходные решетки синтезированы по данным расчетов методами квантовой химии для каждой из модификаций, каждая система содержит около 4 000 атомов и помещается в вакуум. Численный эксперимент состоял в наблюдении за поведением каждой системы во времени при постоянной температуре (300 К), которая поддерживалась термостатом Нозе-Хувера.

При моделировании с начальным набором параметров кристаллическая структура всех гидридов распалась, плотность вещества снизилась почти в два раза. На оптимизированном наборе параметров регулярная структура полностью сохранилась, хотя постоянные решеток кристаллов немного изменились (1-6 %). В конце эксперимента плотности кристаллов при оптимизированных параметрах отличались от реальных на 9, l4 и 4 % для а-, ß- и у-модификаций соответственно.

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

Литература

1. van Duin A.C.T., Dasgupta S., Lorant F., Goddard W.A. III ReaxFF: A Reactive Force Field for Hydrocarbons. Journ. of Physical Chemistry A., 2001, vol. 105, pp. 9396-9409.

2. Nomura K., Kalia R.K., Nakano A., Vashishta P. A scalable parallel algorithm for large-scale reactive force-field molecular dynamics simulations. Computer Physics Communications, 2008, vol. 178, pp. 73-87.

3. ADF.URL: http://www.scm.com (дата обращения: 22 февраля 2014).

4. Программный пакет LAMMPS. URL: http://lammps.san-dia.gov (дата обращения: 22 февраля 2012).

5. van Duin A.C.T., Baas J.M.A., van de Graaf B. Delft Molecular Mechanics: A New Approach to Hydrocarbon Force Fields. Journ. of the Chemical Society, Faraday Transactions, 1994, vol. 90, pp. 2881-2895.

6. Turley J.W., Rinn H.W. Inorg. Chem., 1969, vol. 8, p. 18.

7. Yartys V.A., Denys R.V., Maehlen J.P., Frommen C., Fichtner M., Bulychev B.M. // Inorganic Chemistry, 2007, vol. 46, pp. 1051-1055.

References

1. van Duin A.C.T., Dasgupta S., Lorant F., Goddard W.A. III ReaxFF: a reactive force field for hydrocarbons. Journ. of Physical Chemistry A. 2001, vol. 105, pp. 9396-9409.

2. Nomura K., Kalia R.K., Nakano A., Vashishta P. A scalable parallel algorithm for large-scale reactive force-field molecular dynamics simulations. Computer Physics Communications. 2008, vol. 178, pp. 73-87.

3. ADF. Available at: http://www.scm.com (accessed Feb. 22, 2014).

4. LAMMPS Molecular Dynamics Simulator. Available at: http://lammps.sandia.gov (accessed Feb. 22, 2014).

5. van Duin A.C.T., Baas J.M.A., van de Graaf B. Delft Molecular Mechanics: a new approach to hydrocarbon force fields. Journ. of the Chemical Society, Faraday Transactions. 1994, vol. 90, pp. 2881-2895.

6. Turley J.W., Rinn H.W. Inorg. Chem. 1969, vol. 8, p. 18.

7. Yartys V.A., Denys R.V., Maehlen J.P., Frommen C., Fichtner M., Bulychev B.M. Inorganic Chemistry. 2007, vol. 46, pp. 1051-1055.