Реалізація алгоритму оптимізації параметрів молекулярно-динамічного потенціалу REAXFF
Публікаційна активність
Наступний номер на сайті
Реалізація алгоритму оптимізації параметрів молекулярно-динамічного потенціалу REAXFF
Чисельне моделювання хімічно реактивних молекулярних систем є актуальним для цілого ряду прикладних наукових завдань. Для невеликих хімічних сполук (до 100 атомів) це можна робити методами квантової механіки (ab initio), проте для моделювання більших структур ці методи непридатні через високу ресурсоємність. Традиційні молекулярно-динамічні методи, які зазвичай використовуються для великомасштабних систем, зазвичай засновані на класичному потенціалі силових полів і не дозволяють добре відтворювати багатокомпонентні хімічно реактивні системи. Тому у своїй роботі, пов'язаній із вивченням систем зберігання водню, для моделювання деструкції гідридів легких металів із виділенням водню ми використовуємо комбінований підхід. Він ґрунтується на потенціалі реактивних силових полів ReaxFF (Reactive Force Field) [1, 2], параметризація якого виконується на основі точних квантових розрахунків для малих кластерів. Зазначимо, що потенціал ReaxFF спочатку розроблявся для опису хімічної реактивності, дисоціації та формування хімічних зв'язків, дефектів, поверхневих ефектів та ін.
Потенціал ReaxFF. Форма МД-потенціалу ReaxFF подібна до багатьох нереактивних силових полів, де енергія системи складається з різних енергетичних вкладів. Принципова відмінність полягає у використанні не жорсткого методу зв'язування, а методупорядків зв'язків. Порядок зв'язку характеризує кратність зв'язку між конкретною парою атомів у поєднанні та виражається формулою
,
де індекс l дорівнює σ, π, ππ (типи ковалентного зв'язку); rij - Відстань між атомами; pl1, pl2, rl – параметри [2]. Порядки зв'язків, що обчислюються при чисельному моделюванні для всіх атомів безпосередньо з міжатомних відстаней та оновлюються на кожному кроці за часом, дозволяють описати освіту та дисоціацію зв'язків під час моделювання. Загальний вид потенціалу [2] наступний:
+Econj+Ehbond+Elp+EvdWaals+ECoulomb. (1)
Кожен доданок у правій частині формули відповідає за окремий тип взаємодії: ковалентне (характеризується порядком хімічних зв'язків, складові Ebond, Eover, Eunder), тричасткове (валентні кути, Eval, Epen, Ecoa), чотиричасткове (двогранні валентні кути, Etors, Eco кулонівське (ECoulomb), ван-дер-ваальсове (EvdWaals), водневі зв'язки (Ehbond) та енергія неподілених електронних пар (Elp).
Повна енергія є сумою функцій відносних положень пар атомів рій, трійок рійк і четвірок рійкл. Більше того, кожна з цих функцій залежить від певної кількості параметрів, які визначаються тим, які хімічні елементи взаємодіють. Таким чином, весь потенціал залежить від численних (залежно від складності з'єднання їх кількість може бути більшою за 100) параметрів: EReaxFF=f(p1, p2, p3, …, pN). Аналіз докладного вигляду всіх доданків у формулі (1) дозволяє зробити висновок, що EReaxFF та його перша похідна є безперервними функціями координат атомів навіть за хімічних реакцій [2].
Метод однопараметричної оптимізації. Реалізація алгоритму оптимізації параметрів ґрунтується на методі однопараметричного пошуку [5]. Дані, за якими проводитьсяпошук (оптимізуючий набір, training set), виходять з розрахунків простих моделей хімічних сполук (до 50 атомів) методами квантової хімії або беруться з експерименту. Даними оптимізації є геометричні характеристики молекул (відстань між атомами, валентні кути, двогранні кути), коливальні частоти, ефективні заряди атомів та енергії зв'язку з'єднань, а також будь-які інші характеристики, які можна обчислити методами молекулярної динаміки, а потім порівняти з даними оптимізуючого набір. Процедура однопараметричного пошуку полягає у зміні параметрів потенціалу для мінімізації суми квадратів:
(2)
Тут xiQC/Lit – дані набору, що оптимізує (конкретні характеристики моделі або молекули, отримані методами квантової хімії або з літератури (експерименту)); xiReaxFF – обчислені з використанням ReaxFF значення, що відповідають цим характеристикам; σi – критерій прийнятності, який має розмірність xi та характеризує вклад різниці xiQC/Lit–xiReaxFF у суму квадратів. Сума (2) є так званою функцією помилки – відхилення розрахунків ReaxFF від даних оптимізує набору. Кожне xiReaxFF залежить від параметрів потенціалу p1, . pN, тому Err є функцією цих параметрів.
Важливо відзначити, що для обчислення значення xiReaxFF для кожної моделі набору, що оптимізує, при конкретних значеннях параметрів необхідно провести пошук мінімуму потенційної енергії (1) цієї моделі як функції координат складових її атомів з використанням якого-небудь багатопараметричного методу, наприклад методу сполучених градієнтів. Цю процедуру не слід плутати із процедурою мінімізації суми (2) як функції параметрів потенціалу.
Основою методу однопараметричної оптимізаціїє параболічна екстраполяція, яка виконується по черзі для кожного параметра pk при фіксованих інших. На малих інтервалах залежність між помилкою Err та значенням кожного параметра pk можна вважати параболічною, тому для визначення оптимального значення параметра необхідно обчислити Err при трьох різних значеннях: pk, pk(1–δk), pk(1+δk), де δk
Одна ітерація програми полягає у виконанні описаної процедури один раз для кожного ПК зі списку. Оскільки більшість параметрів у силовому полі ReaxFF пов'язані між собою, оптимальне значення кожного параметра зсувається щоразу, коли змінюється інший параметр. Це означає, що процес оптимізації має бути повторений доти, доки Err не перестане зменшуватися, що визначається критерієм збіжності ε. Зазначимо, що для кожного параметра pk повинен бути заданий допустимий інтервал зміни [Dkmin, Dkmax], щоб уникнути збіжності у нереалістичну область.
Також слід зазначити, що метод однопараметричного пошуку гарантує збіжність мінімум, але цей мінімум може виявитися не абсолютним, а тільки локальним.
Опис програми оптимізації параметрів потенціалу ReaxFF. Блок-схема алгоритму розробленої програми показано малюнку. Загальна структура програми включає ініціалізацію, читання вхідних даних, два вкладені цикли, запис проміжних та вихідних результатів. Зовнішній цикл організує ітерації по повному набору параметрів, а внутрішній забезпечує процедуру екстраполяції параметром pk. Внутрішній цикл містить розрахунки з потенціалом ReaxFF всіх моделей набору, що оптимізує, і є дуже ресурсомістким. Він розпаралелен за допомогою технології MPI, що дозволяє запускати програму на обчислювальних вузлах кластера.
Зовнішній цикл організує ітерації набору параметрів. За один крок цього циклу забезпечується прохід по списку p1, . pN і для кожного pk виконується одна процедура екстраполяції. Результат однієї ітерації – новий набір, який буде використовуватись у наступній ітерації. Цикл завершується, якщо на черговому кроці t буде задоволена умова збіжності Errt-Errt-1
На виході програми виходять два текстові файли: файл із оптимізованими параметрами потенціалу та log-файл (журнал роботи програми). Журнал програми містить докладну інформацію про всі ітерації внутрішнього циклу. На кожній ітерації до нього записуються ім'я параметра; три його значення; помилки (2) цих значень; коефіцієнти екстраполяційної параболи; екстремум параболи; значення параметра, отримане екстраполяцією; екстраполіроване та обчислене значення помилки параметра; остаточне значення параметра, вибране як оптимальне; його помилка (2).
Журнал дозволяє контролювати хід процесу оптимізації, за необхідності зупинити програму та продовжити її виконання з довільної ітерації.
Структура пакету. Вихідний код програми розроблено мовою C і має модульну структуру. Програмний пакет зібраний та доступний у виглядіархіву з вихідними кодами, інструкцією та сценарієм компіляції, документацією для користувача та прикладами вхідних файлів.
Тестування та використання. Функціональність програми перевірялася з прикладу оптимізації параметрів алюмінію та її гідридів. За початкове наближення при оптимізації були взяті параметри, включені до дистрибутиву пакета LAMMPS, які слід вважати досить грубими, оскільки вони дають погані результати для решіток металу та його гідридів.
Для оптимізації параметрів ми використовували набір, до якого увійшли кластери чистого алюмінію (Al2 – Al8, Al13), кластери гідриду (AlnH3n, n=1, . 8), а також кристали α-β- та γ-модифікацій гідриду алюмінію [6, 7], породжені з осередків із періодичними умовами. Додавання набір нескінченних кристалічних структур дозволило врахувати рівняння стану, тобто залежність енергії кристалів від об'єму осередку. Загалом набір включав 56 структур.
Загалом у оптимізації брало участь 79 параметрів потенціалу. Одна ітерація процедури оптимізації (за одним параметром) на обчислювальному кластері займала у разі близько 4 годин. За одну ітерацію величина помилки зменшувалась на значення від 0,1 до 1% порівняно з попередньою ітерацією. Програма зробила близько 4 000 ітерацій, цей час величина помилки скоротилася приблизно 500 раз.
Для перевірки якості отриманих параметрів було виконано моделювання кристалічних решіток гідридів алюмінію при початкових і оптимізованих параметрах з метою порівняння стабільності, щільності і регулярності структури при кімнатній температурі. Всі вихідні решітки синтезовані за даними розрахунків методами квантової хімії для кожної з модифікацій, кожна система містить близько 4000 атомів і міститься у вакуумі. Чисельнийексперимент полягав у спостереженні за поведінкою кожної системи в часі за постійної температури (300 К), яка підтримувалася термостатом Нозе-Хувера.
При моделюванні з початковим набором параметрів кристалічна структура всіх гідридів розпалася, щільність речовини знизилася майже вдвічі. На оптимізованому наборі параметрів регулярна структура повністю збереглася, хоча постійні грати кристалів трохи змінилися (1–6 %). В кінці експерименту щільності кристалів при оптимізованих параметрах відрізнялися від реальних на 9, 14 і 4% для -, - і -модифікацій відповідно.
Таким чином, розроблена програма пошуку параметрів потенціалу ReaxFF є досить універсальною і може бути використана для систем атомів довільного типу. Оптимізуючий набір може включати будь-які з'єднання, хоча набір характеристик цих сполук поки обмежений: до нього входять довжини зв'язків, величини валентних кутів і двогранних валентних кутів, енергії атомізації, теплоти освіти та ефективні заряди атомів. Модульна структура коду дозволяє легко додавати розширення за необхідності врахування додаткових умов чи характеристик моделей. Застосування розроблених методів для оптимізації параметрів потенціалу ReaxFF продемонстровано для MД-моделювання гідридів алюмінію.