Физика и химия стекла, 2023, T. 49, № 4, стр. 401-416

Набор самосогласованных параметров потенциала Леннард–Джонса для молекулярно-динамического моделирования боросиликатных стекол

Г. И. Макаров^1, *, К. С. Шилкова¹, А. В. Шунайлов¹, П. В. Павлов¹, Т. М. Макарова¹

¹ Южно-Уральский государственный университет
454080 Челябинск, пр. Ленина, 76, Россия

Поступила в редакцию 12.08.2022
После доработки 11.10.2022
Принята к публикации 07.12.2022

EDN: SGHZFG
DOI: 10.31857/S0132665122600315

Аннотация

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

Ключевые слова: молекулярная динамика, силовое поле, стекло, силикаты, оксиды, потенциал Леннард–Джонса

ВВЕДЕНИЕ

Современную технологию невозможно представить без разнообразных оксидных неорганических систем – алюмосиликатов, ферритов, разнообразных стекол – применяемых в качестве конструкционных, оптических, магнитных, электродных материалов. Такие материалы интенсивно исследуются и разрабатываются, причем в настоящее время экспериментальный подбор новых материалов все в большей степени дополняется рациональным дизайном [1], когда экспериментальному получению материала с заданными свойствами предшествует создание его достоверной модели, позволяющей оценить его характеристики. Ключом к объяснению и предсказанию свойств материала является его структура; и если для бездефектных кристаллических материалов она может быть с полной достоверностью установлена дифракционными методами, то, например, структура стекол с трудом поддается экспериментальному исследованию в силу отсутствия в них дальнего порядка.

Однако подобные слабо упорядоченные или несущие множественные дефекты системы могут быть смоделированы. Одним из методов, позволяющим моделировать структуру и свойства материалов, является метод классической молекулярной динамики (МД), симулирующий движение частиц атомно-молекулярной системы во времени, опираясь на принципы классической механики. Его применение требует задания потенциальной энергии системы с помощью совокупности уравнений потенциалов межатомных взаимодействий и их параметров – так называемых силовых полей или потенциальных наборов. К настоящему времени разработано несколько таких потенциальных наборов, успешно применяемых для моделирования стекол и кристаллических фаз алюмосиликатного состава. Во всех них межатомное взаимодействие представляется как комбинация электростатического взаимодействия и некоторого дисперсионного потенциала, задающего притяжение–отталкивание ионов. Таковы потенциалы Биста–Крамера–ван Сантена (BKS) [2], пересекающиеся между собой потенциалы Ху [3] и Гудмана [4], использующие потенциал Бэкингема в качестве дисперсионного потенциала, набор потенциалов Педоне [5], использующий в качестве дисперсионного потенциала комбинацию потенциала Морзе и отталкивающего потенциала $\frac{c}{{{{r}^{{12}}}}}.$ Однако, применяющиеся в них дисперсионные потенциалы на основе потенциалов Морзе и Бэкингема обеспечивают не самую высокую производительность расчетов МД с их применением и затрудняют использование GPGPU-ускорителей.

Одной из задач, связанной с моделированием оксидных материалов, является задача моделирования широко применяемых композитных материалов [6], использующих их в качестве наполнителя – главным образом с целью изучения взаимодействия органического полимерного связующего с неорганической подложкой. Описанные выше потенциалы и их аналоги, также основанные на потенциалах Морзе, Бэкингема и потенциалах более сложного вида, несовместимы с силовыми полями для моделирования органических соединений, такими как GAFF [7], GROMOS [8] или CHARMM [9], которые используют более производительный потенциал Леннард–Джонса 6–12. Кроме того, только для потенциала 6–12 решена задача обработки его дальнодействующей части с помощью метода сети частиц Эвальда [10].

Взаимодействие органических полимеров и низкомолекулярных соединений с неорганическими подложками может быть смоделировано с помощью Interface Force Field [11], полностью совместимого с силовыми полями для органических соединений. Однако это эффективное силовое поле основано на явном указании связей между атомами неорганического материала и охватывает только некоторые кристаллические тела, не позволяя моделировать in silico структуры стекол, что требует самопроизвольных структурных перестроек и перераспределения связей в моделируемой фазе.

Таким образом, создание потенциального набора, который, с одной стороны, успешно воспроизводил бы структуру и механические свойства оксидных материалов, подобно, например, потенциальному набору Педоне, и, с другой стороны, использовал бы в качестве дисперсионного потенциала потенциал Леннард–Джонса 6–12, обеспечивая этим совместимость с силовыми полями для органических соединений и высокую производительность расчетов МД – было бы весьма полезным. В настоящей работе мы приступили к решению этой задачи, построив набор потенциалов, воспроизводящих структуру и объемный модуль всестороннего сжатия некоторых бинарных оксидов элементов, содержащихся в распространенных стеклах и некоторых других материалах, с помощью комбинации электростатического потенциала и потенциала Леннард–Джонса 6–12. При этом, однако, мы не подбирали потенциалы Леннард–Джонса, воспроизводящие взаимодействия органических соединений с оксидными материалами.

МЕТОДЫ

Подбор параметров потенциала

Потенциальная энергия моделируемых неорганических кристаллов задавалась как сумма попарных взаимодействий, каждое из которых представляло собой сумму потенциала Леннард–Джонса и кулоновского взаимодействия:

(1)

$U\left( {{{r}_{{ij}}}} \right) = \frac{{{{z}_{i}}{{z}_{j}}{{e}^{2}}}}{{{{r}_{{ij}}}}} + {{\varepsilon }_{{ij}}}\left\{ {\frac{{\sigma _{{ij}}^{{12}}}}{{r_{{ij}}^{{12}}}} - \frac{{\sigma _{{ij}}^{6}}}{{r_{{ij}}^{6}}}} \right\},$

где ${{z}_{i}}$ – заряд иона, а $e = 0.6$ – заряд однозарядного катиона, которому кратны заряды остальных ионов с учетом знака. Эта величина заряда однозарядного иона используется в других силовых полях такого рода и восходит к работе Цунеюки [12], в которой были рассчитаны заряды $ - 1.2e$ атомов кислорода и $2.4e$ атомов кремния для моделирования α-кварца. При этом считалось, что между катионами металлов, кремния, бора действует только кулоновское отталкивание, и соответствующий ${{\varepsilon }_{{ij}}}$ принимается равным нулю; это решение основано на том, что в структурах кристаллов алюмосиликатов, оксидов и тому подобных фаз эти катионы связаны только кислородными мостиками.

Параметры ${{\sigma }_{{ij}}}$ и ${{\varepsilon }_{{ij}}}$ подбирались так, чтобы воспроизводить параметры элементарной ячейки, координаты и объемный модуль всестороннего сжатия (далее объемный модуль) соответствующего бинарного оксида (или силиката, в случае щелочных металлов). Источником кристаллических структур (табл. 1) служила Crystallography Open Database (COD) [13], значения объемного модуля для кристаллов извлекались из обзора [14], если не указано иное. Поиск параметров потенциала Леннард–Джонса выполнялся с помощью программы GULP [15] нижеследующим образом. Начальное приближение ${{\sigma }_{{ij}}}$ рассчитывалось из соответствующего межатомного расстояния ${{r}_{{ij}}},$ находимого из соответствующей пространственной структуры кристалла,

${{\sigma }_{{ij}}} = \frac{{{{r}_{{ij}}}}}{{\sqrt[6]{2}}},$

а начальное приближение ${{\varepsilon }_{{ij}}}$ произвольно принималось равным 100 кДж/моль. Затем эти значения попеременно уточнялись с помощью программы GULP [15], начиная с ${{\varepsilon }_{{ij}}},$ до тех пор, пока изменения этих величин на очередном шаге не становились малым. Полученные ${{\sigma }_{{ij}}}$ и ${{\varepsilon }_{{ij}}}$ затем последовательно совместно уточнялись до тех пор, пока не достигалось минимальное расхождение экспериментально измеренного и рассчитанного объемных модулей при сохранении не превышающего 1 процента расхождения геометрических параметров элементарной ячейки. Расчет геометрических характеристик и механических свойств кристалла при каждых рассматриваемых ${{\sigma }_{{ij}}}$ и ${{\varepsilon }_{{ij}}}$ выполнялся после минимизации свободной энергии Гиббса кристалла при температуре 298,15 К и давлении 100 кПа.

Таблица 1.

Сравнение экспериментальных и расчетных структурных параметров и объемных модулей использованных при построении потенциального набора кристаллических фаз

Кристалл	COD ID	Пространственная группа		a, Å	b, Å	c, Å	α, град.	β, град.	γ, град.	V, Å³	Объемный модуль, ГПа
SiO₂ α-кварц	9 000 775	P3₂21	Эксперимент	4.92	4.92	5.41	90.00	90.00	120.00	113.13	37.80
			Расчет	4.92	4.92	5.37	90.00	89.99	120.02	112.35	37.80
			Расхождение, %	0.00	0.01	–0.68	0.00	–0.01	0.02	–0.69	–0.00
GeO₂ аргутит	9 006 849	P4₂/mnm	Эксперимент	4.40	4.40	2.86	90.00	90.00	90.00	55.33	257.60
			Расчет	4.41	4.41	2.79	89.83	89.83	90.14	54.23	263.07
			Расхождение, %	0.34	0.34	–2.65	–0.19	–0.19	0.16	–1.99	2.12
TiO₂ рутил	9 015 662	P4₂/mnm	Эксперимент	4.59	4.59	5.92	90.00	90.00	90.00	124.87	215.00
			Расчет	4.61	4.61	5.93	90.00	90.00	89.99	126.01	355.75
			Расхождение, %	0.36	0.36	0.19	–0.00	–0.00	–0.01	0.91	65.46
ZrO₂ бадделеит	9007 448	P12₁/c₁	Эксперимент	5.15	5.21	5.32	90.00	99.23	90.00	140.88	159.00
			Расчет	5.25	5.31	5.34	90.00	98.21	90.00	147.38	159.81
			Расхождение, %	1.85	1.94	0.48	–0.00	–1.03	0.00	4.61	0.51
NaBSiO₄ малинкоит	9 004 616	P6₃	Эксперимент	13.90	13.90	7.70	90.00	90.00	120.00	1287.75
			Расчет	14.01	14.01	7.82	90.01	89.99	120.02	1329.37
			Расхождение, %	0.83	0.83	1.56	0.01	–0.01	0.02	3.23
Al₂O₃ корунд	2 300 375	R$\bar {3}$c	Эксперимент	4.76	4.76	12.99	90.00	90.00	120.00	254.42	253.50
			Расчет	4.70	4.70	13.25	89.96	90.02	119.98	253.96	255.24
			Расхождение, %	–1.09	–1.09	2.02	–0.04	0.02	–0.02	–0.18	0.69

Кристалл	COD ID	Пространственная группа		a, Å	b, Å	c, Å	α, град.	β, град.	γ, град.	V, Å³	Объемный модуль, ГПа
Ga₂O₃	2004987	C2/m	Эксперимент	12.21	3.04	5.80	90.00	103.83	90.00	208.85	199.00	^[23]
			Расчет	12.12	3.00	5.93	89.97	104.78	90.01	208.70	223.94
			Расхождение, %	–0.76	–1.11	2.26	–0.04	0.91	0.01	–0.07	12.53
Cr₂O₃ эсколаит	9016327	R$\bar {3}$c	Эксперимент	4.95	4.95	13.58	90.00	90.00	120.00	288.47	223.70
			Расчет	4.95	4.95	13.66	90.00	90.00	120.00	289.89	293.36
			Расхождение,	–0.06	–0.06	0.62	–0.00	0.00	0.00	0.49	31.14
Fe₂O₃ гематит	9000139	R$\bar {3}$c	Эксперимент	5.04	5.04	13.77	90.00	90.00	120.00	302.72	206.60
			Расчет	5.05	5.05	13.84	90.00	90.00	120.00	305.20	268.78
			Расхождение, %	0.16	0.16	0.50	–0.00	0.00	0.00	0.82	30.10
Sc₂O₃	4326667	Ia$\bar {3}$	Эксперимент	9.85	9.85	9.85	90.00	90.00	90.00	956.98	174.50	^[24]
			Расчет	9.88	9.88	9.88	90.00	90.00	90.00	963.01	208.31
			Расхождение, %	0.21	0.21	0.21	–0.00	–0.00	–0.00	0.63	19.38
Y₂O₃	1537882	Ia$\bar {3}$	Эксперимент	10.60	10.60	10.60	90.00	90.00	90.00	1191.35	149.50	^[25]
			Расчет	10.60	10.60	10.60	90.00	90.00	90.00	1191.68	155.89
			Расхождение, %o	0.01	0.01	0.01	–0.00	–0.00	–0.00	0.03	4.28
MnO манганозит	1514105	Fm$\bar {3}$m	Эксперимент	4.44	4.44	4.44	90.00	90.00	90.00	87.77	153.00
			Расчет	4.45	4.45	4.45	89.95	89.95	89.95	88.29	162.34
			Расхождение, %	0.20	0.20	0.20	–0.06	–0.06	–0.06	0.59	6.10

Кристалл	COD ID	Пространственная группа		a, Å	b, Å	c, Å	α, град.	β, град.	γ, град.	V, Å³	Объемный модуль, ГПа
FeO вюстит	1011169	Fm$\bar {3}$m	Эксперимент	4.40	4.40	4.40	90.00	90.00	90.00	85.12	223.70
			Расчет	4.40	4.40	4.40	89.96	89.96	89.96	85.32	223.70
			Расхождение, %	0.08	0.08	0.08	–0.05	–0.05	–0.05	0.24	0.00
CoO	1533087	Fm$\bar {3}$m	Эксперимент	4.26	4.26	4.26	90.00	90.00	90.00	77.47	183.30
			Расчет	4.28	4.28	4.28	90.00	90.00	90.00	78.67	187.07
			Расхождение, %	0.51	0.51	0.51	0.00	0.00	0.00	1.54	2.06
NiO бунзенит	4320496	Fm$\bar {3}$m	Эксперимент	4.19	4.19	4.19	90.00	90.00	90.00	73.77	205.00
			Расчет	4.20	4.20	4.20	89.94	89.94	89.94	73.85	214.06
			Расхождение, %	0.04	0.04	0.04	–0.06	–0.06	–0.06	0.11	4.42
BeO бромеллит	9011659	P6₃mc	Эксперимент	2.70	2.70	4.38	90.00	90.00	120.00	27.61	251.00
			Расчет	2.71	2.72	4.40	90.00	90.00	119.86	28.09	300.14
			Расхождение, %	0.33	0.76	0.49	0.00	0.00	–0.12	1.73	19.58
MgO периклаз	1011116	Fm$\bar {3}$m	Эксперимент	4.22	4.22	4.22	90.00	90.00	90.00	74.94	160.00
			Расчет	4.21	4.21	4.21	89.94	89.94	89.94	74.72	204.37
			Расхождение, %	–0.10	–0.10	–0.10	–0.07	–0.07	–0.07	–0.29	27.73
CaO известь	1011095	Fm$\bar {3}$m	Эксперимент	4.80	4.80	4.80	90.00	90.00	90.00	110.94	111.00
			Расчет	4.82	4.82	4.82	89.94	89.94	89.94	112.05	111.87
			Расхождение, %	0.33	0.33	0.33	–0.06	–0.06	–0.06	1.00	0.78

Кристалл	COD ID	Пространственная группа		a, Å	b, Å	c, Å	α, град.	β, град.	γ, град.	V, Å³	Объемный модуль, ГПа
SrO	7200689	Fm$\bar {3}$m	Эксперимент	5.16	5.16	5.16	90.00	90.00	90.00	137.49	87.30
			Расчет	5.16	5.16	5.16	89.95	89.95	89.95	137.52	87.27
			Расхождение, %	0.01	0.01	0.01	–0.06	–0.06	–0.06	0.02	–0.03
BaO	9008599	Fm$\bar {3}$m	Эксперимент	5.52	5.52	5.52	90.00	90.00	90.00	168.47	70.70
			Расчет	5.52	5.52	5.52	89.95	89.95	89.95	168.52	70.69
			Расхождение, %	0.01	0.01	0.01	–0.06	–0.06	–0.06	0.03	–0.02
Na₂SiO₃	2310858	Cmc2₁	Эксперимент	10.63	6.16	4.88	90.00	90.00	90.00	319.16	62.00	^[26]
			Расчет	10.81	6.11	4.83	90.04	89.99	90.03	319.14	61.81
			Расхождение, %	1.74	–0.79	–0.93	0.04	–0.01	0.04	–0.01	–0.30
KAlSiO₄ кальсилит	9001818	P31c	Эксперимент	5.16	5.16	8.71	90.00	90.00	120.00	200.51	59.63	^[19]
			Расчет	5.14	5.14	8.77	90.02	90.02	120.02	200.51	59.63
			Расхождение, %	–0.33	–0.33	0.70	0.02	0.02	0.02	0.00	0.00
AlPO4 берлинит	9006404	P3₂21	Эксперимент	4.94	4.94	10.95	90.00	90.00	120.00	231.77	29.30	^[27]
			Расчет	4.96	4.96	10.86	89.94	90.00	120.04	231.22	27.90
			Расхождение, %	0.38	0.28	–0.85	–0.07	–0.00	0.04	–0.24	–4.78

Необходимо также отметить, что построение набора параметров происходило в определенной последовательности. Сначала были найдены ${{\sigma }_{{ij}}}$ и ${{\varepsilon }_{{ij}}}$ для взаимодействий Si–O и O–O, воспроизводящих структуру, объемный модуль и модуль Юнга SiO₂ в форме α-кварца. Оксид кремния был выбран в качестве первоначальной точки потому, что силикатные кристаллы и стекла чрезвычайно распространены, так что разрабатываемый набор потенциалов должен хорошо воспроизводить структуру и механические свойства этих практически значимых материалов. Затем отыскивались ${{\sigma }_{{ij}}}$ и ${{\varepsilon }_{{ij}}}$ для взаимодействий X–O при неизменных параметрах взаимодействия O–O, где под X понимается соответствующий катион. При этом использовались структуры и экспериментально измеренные объемные модули таких оксидных фаз, которые были бы наиболее стабильны в нормальных условиях и содержали бы равнозначные контакты ионов X с ионами кислорода, так как потенциал Леннард–Джонса для взаимодействия X–O зависит только от расстояния между ионами. В наибольшей степени этим требованиям удовлетворяли бинарные оксиды, такие, как MgO (периклаз) или Al₂O₃ (корунд). Однако в некоторых случаях приходилось опираться на оксидные фазы более сложного состава, содержащие рассматриваемый ион X, для которого требуется подобрать параметры взаимодействия X–O, и другие ионы, для которых уже найдены параметры их взаимодействия с ионами кислорода. В случае щелочных металлов параметры взаимодействия X–O отыскивались для других кристаллов: метасиликата Na₂SiO₃ для натрия и кальсилита KAlSiO₄ при фиксированных параметрах взаимодействий Si–O, Al–O и O–O. Аналогичным способом параметры взаимодействия B–O искали, используя кристалл минерала малинкоита при неизменных взаимодействиях Na–O, Si–O и O–O, или параметры взаимодействия P–O, опираясь на структуру и свойства берлинита. Такой же прием применялся в работе [5]^.

Молекулярно-динамическое моделирование кристаллов

Молекулярно-динамические расчеты выполнялись в пакете GROMACS [16] версии 2019.4. Температура моделируемой системы управлялась термостатом масштабирования скоростей с добавочным стохастическим членом [17]. Постоянное давление было анизотропным и поддерживалось баростатом Берендсена [18] с периодом привязки 2 пс и изотермической сжимаемостью 4.65 ТПа^–1 (что соответствует объемному модулю кальсилита 59.6 ГПа [19]) при периодических граничных условиях по всем трем направлениям. Для обработки электростатических и Ван-дер-Ваальсовых взаимодействий применялся метод суммирования по Эвальду (PME) [20, 21].

Источником структур выбранных для моделирования кристаллов также служила Crystallography Open Database [13]. Для каждого из кристаллов (табл. 2) элементарная ячейка размножалась так, чтобы в полученной суперячейке a, b, c оказались примерно равны и в ней насчитывалось несколько тысяч атомов. Эта суперячейка последовательно подвергалась оптимизации потенциальной энергии методом Бройдена–Флетчера–Гольдфарба–Шанно и молекулярной динамике при постоянных объеме и температуре 298.15 К продолжительностью 100 пс при шаге интегрирования 1 фс. Полученное состояние вовлекалось в расчет молекулярной динамики при постоянных температуре 298.15 К и постоянном анизотропном давлении 1 бар продолжительностью 10 нс при шаге интегрирования 2 фс, в котором релаксировались длины ребер суперячейки. Конечное состояние уравновешивалось молекулярно-динамическим моделированием продолжительностью 6 нс при тех же шаге интегрирования и температуре, тогда как постоянное анизотропное давление в 1 бар, в отличие от предыдущего расчета, поддерживалось при ненулевых недиагональных элементах тензора сжимаемости, равных диагональным, что допускало изменения углов суперячейки. Для полученной таким способом траектории анализировались изменения параметров суперячейки a, b, c и ее объем V. Также рассчитывалась величина RMSD, понимаемая как среднее квадратичное смещение в данный момент времени t всех N атомов суперячейки:

(2)

${\text{RMSD}}\left( t \right) = \sqrt {\frac{{\sum\limits_i^N {{{{\left( {{{r}_{i}}\left( t \right) - {{r}_{{i,ref}}}} \right)}}^{2}}} }}{N}} .$

Таблица 2.

Сравнение экспериментальных структурных параметров пробных кристаллических фаз с расчетными, полученными при молекулярно-динамическом моделировании при постоянном давлении. Для МД расчетов показаны средние величины на участках с 4 по 10 нс траектории. Дана средняя величина RMSD на рассматриваемом участке траектории

Кристалл	COD ID	Симметрия	Суперячейка		a, нм	b, нм	c, нм	α	β	γ	V, нм³	RMSD, нм
Fe₃O₄	7 228 110	Fd$\bar {3}$m	3 × 3 × 3	Эксперимент	2.5080	2.5080	2.5080	90.00	90.00	90.00	15.8	0.0272
				Расчет	2.5616	2.5393	2.5626	90.00	89.79	90.13	16.7
				Расхождение, %	2.14	1.25	2.18	0.00	–0.24	0.15	5.7
Fe₂TiO₄	9 005 222	Fd$\bar {3}$m	3 × 3 × 3	Эксперимент	2.5632	2.5632	2.5632	90.00	90.00	90.00	16.8	0.0256
				Расчет	2.6032	2.6032	2.6032	90.00	90.00	90.00	17.6
				Расхождение, %	1.56	1.56	1.56	0.00	0.00	0.00	4.8
MgCr₂O₄	9 006 181	Fd$\bar {3}$m	3 × 3 × 3	Эксперимент	2.5003	2.5003	2.5003	90.00	90.00	90.00	15.6	0.0096
				Расчет	2.5021	2.5021	2.5021	90.00	90.00	90.00	15.7
				Расхождение, %	0.07	0.07	0.07	0.00	–0.00	0.00	0.2
Ca₃Cr₂Si₃O₁₂	9 007 149	Ia$\bar {3}$d	3 × 3 × 3	Эксперимент	3.5991	3.5991	3.5991	90.00	90.00	90.00	46.6	0.0139
				Расчет	3.6184	3.6184	3.6184	90.00	90.00	90.00	47.4
				Расхождение, %	0.54	0.54	0.54	0.00	0.00	0.00	1.6
Ca₃Ga₂Ge₃O₁₂	1 529 480	Ia$\bar {3}$d	3 × 3 × 3	Эксперимент	3.6768	3.6768	3.6768	90.00	90.00	90.00	49.7	0.0240
				Расчет	3.6351	3.6351	3.6351	90.00	90.00	90.00	48.0
				Расхождение, %	–1.14	–1.14	–1.14	0.00	0.00	0.00	–3.4
Ca₃Sc₂Ge₃O₁₂	1 530 425	Ia$\bar {3}$d	3 × 3 × 3	Эксперимент	3.7536	3.7536	3.7536	90.00	90.00	90.00	52.9	0.0263
				Расчет	3.7073	3.7073	3.7073	90.00	90.00	90.00	51.0
				Расхождение, %	–1.23	–1.23	–1.23	0.00	0.00	0.00	–3.7
Co₃Al₂Si₃O₁₂	9 010 332	Ia$\bar {3}$d	3 × 3 × 3	Эксперимент	3.4377	3.4377	3.4377	90.00	90.00	90.00	40.6	0.0261
				Расчет	3.4820	3.4820	3.4820	90.00	90.00	90.00	42.2
				Расхождение, %	1.29	1.29	1.29	0.00	0.00	0.00	3.9
Mn₃Al₂Si₃O₁₂	9 002 692	Ia$\bar {3}$d	3 × 3 × 3	Эксперимент	3.4845	3.4845	3.4845	90.00	90.00	90.00	42.3	0.0132
				Расчет	3.4987	3.4987	3.4987	90.00	90.00	90.00	42.8
				Расхождение, %	0.41	0.41	0.41	0.00	0.00	0.00	1.2
Sr₃Al₂Si₃O₁₂	1 538 398	Ia$\bar {3}$d	3 × 3 × 3	Эксперимент	3.6180	3.6180	3.6180	90.00	90.00	90.00	47.4	0.1924
				Расчет	3.9180	3.8646	3.8340	90.00	90.41	89.59	58.0
				Расхождение, %	8.29	6.82	5.97	0.00	0.46	–0.46	22.6
Y₃Al₅O₁₂	4 312 142	Ia$\bar {3}$d	3 × 3 × 3	Эксперимент	3.5970	3.5970	3.5970	90.00	90.00	90.00	46.5	0.0258
				Расчет	3.5495	3.5495	3.5495	90.00	90.00	90.00	44.7
				Расхождение, %	–1.32	–1.32	–1.32	0.00	0.00	0.00	–3.9

Кристалл	COD ID	Симметрия	Суперячейка		a, нм	b, нм	c, нм	α	β	γ	V, нм³	RMSD, нм
K₂Ca₃Si₃O₁₀	1 561 885	P$\bar {1}$	8 × 6 × 4	Эксперимент	4.5360	4.4250	4.7324	89.79	80.62	88.70	93.5	0.0467
				Расчет	4.6140	4.4365	4.8426	89.63	80.10	87.82	97.4
				Расхождение, %	1.72	0.26	2.33	–0.18	–0.64	–0.99	4.1
CaNiSi₂O₆	9 001 303	C2/c	4 × 5 × 8	Эксперимент	3.8944	4.4465	4.1824	90.00	105.83	90.00	69.7	0.0266
				Расчет	3.9569	4.4638	4.2133	90.00	106.43	90.00	71.4
				Расхождение, %	1.60	0.39	0.74	–0.00	0.57	0.00	2.4
CaZrSi₂O₇	9 004 208	C2	6 × 6 × 6	Эксперимент	4.1112	5.1954	2.8116	90.00	101.69	90.00	58.8	0.0343
				Расчет	4.1986	5.1607	2.8662	90.00	101.48	90.00	60.9
				Расхождение, %	2.13	–0.67	1.94	0.00	–0.20	–0.00	3.5
Be₂Y₂SiO₇	1 531 636	P$\bar {4}$2₁m	6 × 6 × 10	Эксперимент	4.3686	4.3686	4.7550	90.00	90.00	90.00	90.7	0.0490
				Расчет	4.4603	4.4623	4.7522	90.00	90.03	90.05	94.6
				Расхождение, %	2.10	2.15	–0.06	–0.00	0.03	0.05	4.2
KAlSiO₄	9 001 818	P31c	6 × 6 × 4	Эксперимент	3.0942	3.0942	3.4824	90.00	90.00	120.00	28.9	0.0186
				Расчет	3.0765	3.0758	3.4951	90.00	90.00	120.01	28.6
				Расхождение, %	–0.57	–0.60	0.37	0.00	–0.00	0.01	–0.8
NaBSiO₄	9 004 616	P6₃	4 × 4 × 6	Эксперимент	5.5584	5.5584	4.6200	90.00	90.00	120.00	123.6	0.0322
				Расчет	5.5881	5.5881	4.7057	90.00	90.00	120.00	127.3
				Расхождение, %	0.53	0.53	1.85	–0.00	0.00	0.00	2.9
K₂TiSi₃O₉	1 001 111	P6₃/m	6 × 6 × 4	Эксперимент	4.0644	4.0644	3.9688	90.00	90.00	120.00	56.8	0.0296
				Расчет	4.0944	4.0928	3.9824	90.21	90.60	87.19	57.8
				Расхождение, %	0.74	0.70	0.34	0.23	0.67	–27.35	1.8
Na₃YSi₂O₇	1 530 410	P6₃/m	4 × 4 × 3	Эксперимент	3.7688	3.7688	4.1370	90.00	90.00	120.00	50.9	0.0215
				Расчет	3.7867	3.7867	4.1370	90.00	90.00	120.00	51.4
				Расхождение, %	0.48	0.48	–0.00	–0.00	0.00	0.00	1.0
SrFe₂P₄O₁₄	2 003 248	P$\bar {1}$	8 × 6 × 5	Эксперимент	3.8360	4.2648	3.9150	89.30	87.59	73.11	61.2	0.0421
				Расчет	3.8587	4.3534	3.9443	89.73	89.11	72.51	63.2
				Расхождение, %	0.59	2.08	0.75	0.48	1.73	–0.81	3.2
CoFePO₅	1 525 963	Pnma	5 × 5 × 5	Эксперимент	3.6275	3.2235	3.7285	90.00	90.00	90.00	43.6	0.0352
				Расчет	3.7335	3.2493	3.7405	90.00	90.00	90.00	45.4
				Расхождение, %	2.92	0.80	0.32	0.00	0.00	0.00	4.1
KNi₄P₃O₁₂	1 520 998	Pmnn	7 × 2 × 4	Эксперимент	4.3064	3.2428	3.7936	90.00	90.00	90.00	53.0	0.0652
				Расчет	4.2852	3.4386	3.8577	90.00	90.00	90.00	56.8
				Расхождение, %	–0.49	6.04	1.69	0.00	0.00	0.00	7.3
NaFeP₂O₇	1 001 846	P2₁/c	5 × 5 × 4	Эксперимент	3.6620	3.9520	3.8300	90.00	111.86	90.00	51.4	0.0271
				Расчет	3.7112	3.9206	3.8635	90.00	111.46	90.00	52.3
				Расхождение, %	1.34	–0.79	0.87	–0.00	–0.36	0.00	1.7
Cr₃BO₆	7 209 462	Pnma	3 × 3 × 6	Эксперимент	2.9565	2.5203	2.6484	90.00	90.00	90.00	19.7	0.0233
				Расчет	2.9657	2.5458	2.6930	90.00	90.00	89.99	20.3
				Расхождение, %	0.31	1.01	1.68	–0.00	–0.00	–0.01	3.0
BaMnB₂O₅	2v241 822	P2₁/c	5 × 5 × 6	Эксперимент	4.1435	4.3285	3.9156	90.00	92.87	90.00	70.1	0.2724
				Расчет	4.3137	4.5108	3.4562	90.00	80.86	90.00	66.3
				Расхождение, %	4.11	4.21	–11.73	0.00	–12.93	0.00	–5.5
K₃Be₆B₉O₂₁	4 102 584	P2₁	6 × 3 × 5	Эксперимент	4.4862	3.9810	4.6630	90.00	111.27	90.00	77.6	0.2072
				Расчет	4.4266	3.9787	4.6076	90.16	110.81	89.56	75.9
				Расхождение, %	–1.33	–0.06	–1.19	0.17	–0.42	–0.49	–2.3
BaAl₂Ge₂O₈	2 106 270	I2/c	5 × 3 × 3	Эксперимент	4.3995	4.0113	4.4181	90.00	114.93	90.00	70.7	0.0519
				Расчет	4.2797	3.9805	4.3634	90.00	116.72	90.00	66.4
				Расхождение, %	–2.72	–0.77	–1.24	–0.00	1.56	0.00	–6.1
Na₄Zr₂Ge₃O₁₂	1 530 786	R$\bar {3}$	5 × 5 × 2	Эксперимент	4.7145	4.7145	4.5058	90.00	90.00	120.00	86.7	0.0359
				Расчет	4.6725	4.6725	4.4546	90.00	90.00	120.00	84.2
				Расхождение, %	–0.89	–0.89	–1.14	0.00	–0.00	0.00	–2.9
Mg₂TiO₄	9 001 694	P4₁22	4 × 4 × 3	Эксперимент	2.3908	2.3908	2.5248	90.00	90.00	90.00	14.4	0.0172
				Расчет	2.4220	2.4220	2.5076	90.00	90.00	90.00	14.7
				Расхождение, %	1.31	1.31	–0.68	0.00	0.00	0.00	1.9
Ba₂YFe₃O₈	4 309 913	P4/mmm	10 × 10 × 3	Эксперимент	3.9180	3.9180	3.5472	90.00	90.00	90.00	54.5	0.0626
				Расчет	3.9021	3.8900	3.7471	90.00	90.00	90.00	56.9
				Расхождение, %	–0.41	–0.72	5.64	0.00	0.00	0.00	4.5
BaCaGa₄O₈	2 008 101	Imm2	5 × 4 × 8	Эксперимент	4.1270	3.6832	4.1288	90.00	90.00	90.00	62.8	0.0307
				Расчет	4.1369	3.7606	4.0828	90.00	90.00	90.00	63.5
				Расхождение, %	0.24	2.10	–1.11	0.00	0.00	0.00	1.2
FeScO₃	4 122 825	R3c	8 × 8 × 3	Эксперимент	4.1576	4.1576	4.1808	90.00	90.00	120.00	62.6	0.0134
				Расчет	4.1371	4.1373	4.1911	90.00	90.00	120.00	62.1
				Расхождение, %	–0.49	–0.49	0.25	–0.00	0.00	–0.00	–0.7

Значения RMSD усреднялись с 4 по 10 нс траектории, в качестве опорных координат использовались исходные координаты атомов ${{r}_{{i,ref}}}$ в суперячейке. Средняя величина ${\text{RMSD}}$ использовалась как мера смещения атомов от их исходных положений в суперячейке.

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

Мы ставили своей целью найти параметры потенциалов Леннард–Джонса 6–12, позволяющих моделировать оксидные неорганические кристаллы и, главным образом, стекла того же состава совместно с органическими соединениями. При этом мы следовали широко используемому принципу построения таких потенциальных наборов, предполагающему, что в моделируемой системе катионы металлов, бора, кремния и т.п. элементов не взаимодействуют непосредственно, а только через посредство кислородных мостиков. При этом также предполагается, что анионы кислорода создают вокруг каждого из катионов координационную сферу, которая корректируется дисперсионным взаимодействием O–O. Эти соображения подкреплялись употребительными представлениями о структуре стекол, согласно которым стекла являются неупорядоченной пространственной сеткой, образованной атомами кислорода и кремния (или других элементов-стеклообразователей – бора или фосфора), причем связи между атомами кислорода и элемента-сеткообразователя имеют ковалентный характер. Внутри этой сетки по ионному механизму связываются ионы металлов, расположенные друг от друга на заметном, 0.3–0.5 нм, расстоянии и потому незначительно взаимодействующие друг с другом. Это косвенно подтверждается миграцией ионов из стекла в раствор, что является физической базой для создания ион-селективных электродов [22].

Исходя из этих соображений, был построен набор параметров потенциала 6–12, представленный в табл. 3. С подобранными параметрами выполняли минимизацию свободной энергии Гиббса для элементарных ячеек выбранных кристаллов. Расчеты показали, что построенный потенциальный набор вполне воспроизводит параметры элементарных ячеек (табл. 1), чьего воспроизведения в первую очередь добивались при построении набора параметров потенциала 6–12.

Таблица 3.

Разработанные параметры потенциалов Леннард–Джонса 6–12 для взаимодействий ионов в неорганических кристаллах. Используемые заряды ионов показаны верхними индексами

		σ, нм	ε, кДж/моль
O^{–1, 2}	O^{–1, 2}	0.286606	4.681993
P^{3, 0}	O^{–1, 2}	0.152660	2.744795 × 10¹
Si^{2, 4}	O^{–1, 2}	0.155210	7.015628 × 10¹
Ge^{2, 4}	O^{–1, 2}	0.199319	1.762161
Ti^{2, 4}	O^{–1, 2}	0.295934	7.548651 × 10^–2
Zr^{2, 4}	O^{–1, 2}	0.223406	1.269286 × 10¹
B^{1, 8}	O^{–1, 2}	0.124560	2.615126 × 10²
Al^{1, 8}	O^{–1, 2}	0.173320	1.596995 × 10¹
Ga^{1, 8}	O^{–1, 2}	0.329149	1.314018 × 10^–2
Cr^{1, 8}	O^{–1, 2}	0.264629	2.067842 × 10^–1
Fe^{1, 8}	O^{–1, 2}	0.349606	9.729074 × 10^–3
Sc^{1, 8}	O^{–1, 2}	0.314746	7.227098 × 10^–2
Y^{1, 8}	O^{–1, 2}	0.303543	3.150214 × 10^–1
Mn^{1, 2}	O^{–1, 2}	0.293912	2.004903 × 10^–1
Fe^{1, 2}	O^{–1, 2}	0.226718	6.394907
Co^{1, 2}	O^{–1, 2}	0.363822	8.450698 × 10^–3
Ni^{1, 2}	O^{–1, 2}	0.276104	1.827676 × 10^–1
Be^{1, 2}	O^{–1, 2}	0.205086	6.370660 × 10^–2
Mg^{1, 2}	O^{–1, 2}	0.315128	3.730477 × 10^–2
Ca^{1, 2}	O^{–1, 2}	0.326660	1.425480 × 10^–1
Sr^{1, 2}	O^{–1, 2}	0.313530	5.449972 × 10^–1
Ba^{1, 2}	O^{–1, 2}	0.309053	1.544622
Na^{0, 6}	O^{–1, 2}	0.247361	2.657484
K^{0, 6}	O^{–1, 2}	0.285044	2.892574

При построении потенциального набора также старались воспроизвести объемные модули кристаллов, используя эту величину как характеристику интенсивности межатомных взаимодействий. Выбор объемного модуля основывался также на том обстоятельстве, что именно эту величину чаще всего измеряли для любых неорганических кристаллов, тогда как модуль Юнга или упругие константы для многих рассмотренных кристаллов найти не представлялось возможным. Для некоторых кристаллов – оксидов скандия и иттрия, периклаза, рутила и гематита – не удалось одновременно воспроизвести параметры ячейки и объемный модуль. В таком случае воспроизведением объемного модуля жертвовали в пользу воспроизведения структуры (табл. 1).

Случай потенциала B–O требует частного обсуждения. В структуре тригонального оксида бора (α-B₂O₃), стабильного при нормальных условиях, взаимодействия бор–кислород неравны между собой: каждый атом бора окружен пятью атомами кислорода, расположенными в вершинах полярно растянутой тригональной бипирамиды, так что три экваториальных атома кислорода ковалентно связаны с бором и поэтому сближены с ним сильнее, чем два полярных атома кислорода (рис. 1). Поэтому оксид бора бессмысленно использовать для подгонки параметров потенциала Леннард–Джонса 6–12. Вместо этого использовался минерал малинкоит NaBSiO₄, в структуре которого атомы бора имеют тетраэдрическую координационную сферу с равноценными контактами бор–кислород. Исходя из тех же соображений минерал берлинит AlPO₄ использовался для подбора параметров фосфор–кислород вместо оксида фосфора P₄O₁₀.

Рис. 1.

Координационная сфера атома бора в его оксиде B₂O₃ (COD ID 1510796): бор показан синим, кислород – красным, грани воображаемой искаженной тригональной бипирамиды показаны черными линиями. Черным пунктиром показаны контакты между атомом бора и полярными атомами кислорода.

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

Выбранные кристаллические фазы были вовлечены в молекулярно-динамическое моделирование в пакете GROMACS при комнатной температуре и постоянном анизотропном давлении. В ходе расчетов длины ребер и объем моделируемой суперячейки стабилизировались между 1 и 4 нс траекторий, так что участок с 4 по 10 нс был выбран для расчета средних величин. Расхождение средних расчетных длин ребер суперячеек и соответствующих экспериментальных значений для большинства смоделированных фаз не превышает 3%, причем наилучший результат наблюдается для магнезиохромита MgCr₂O₄ (табл. 2). Чтобы охарактеризовать отклонение структуры кристаллической фазы в ходе расчета от экспериментально определенного ее состояния, рассчитывался средний RMSD на том же участке траектории. Среднее отклонение атомов в большинстве фаз не превосходит 0.06 нм, причем наименее искаженными оказались структуры магнезиохромита MgCr₂O₄, смешанного оксида железа(III) и скандия FeScO₃ и уваровита Ca₃Cr₂Si₃O₁₂. Однако структуры некоторых фаз оказывается ощутимо нарушаются в процессе расчета: для стронциевого граната Sr₃Al₂Si₃O₁₂, фосфата никеля-калия KNi₄P₃O₁₂, бората бария-марганца BaMnB₂O₅ и смешанного оксида бария, железа и иттрия Ba₂YFe₃O₈ расхождение средних длин ребер моделируемой ячейки с экспериментальными значениями весьма велико. Тем не менее, разработанный нами потенциальный набор успешно воспроизводит структурные параметры для большинства смоделированных весьма разнообразных кристаллических фаз сложного состава.

ЗАКЛЮЧЕНИЕ

В настоящей работе был построен набор параметров потенциала Леннард–Джонса 6–12, удовлетворительно воспроизводящий геометрические характеристики и объемные модули стабильных кристаллических соединений ряда химических элементов, составляющих употребительные неорганические материалы. Предлагаемый набор параметров потенциала 6–12 позволяет с высокой производительностью моделировать методами классической молекулярной динамики кристаллические и стекловидные материалы оксидного и алюмоборосиликатного характера.

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

Ma M., Li H., Xiong Y., Dong F. Rational design, synthesis, and application of silica/graphene-based nanocomposite: A review // Materials & Design, 2021. V. 198. P. 109367.
van Beest B.W.H., Kramer G.J., van Santen R.A. Force fields for silicas and aluminophosphates based on ab initio calculations // Physical Review Letters. 1990. V. 64. P. 1955–1958.
Hu Y.-J., Zhao G., Zhang M., Bin B., Del Rose T., Zhao Q., Zu Q., Chen Y., Sun X., de Jong M., Qi L. Predicting densities and elastic moduli of SiO₂-based glasses by machine learning // npj Computational Materials, 2020. V. 6. P. 25.
Goodman B.J. A Study of Vitrified Nuclear Wasteforms by Molecular Dynamics, Electron Microscopy and Raman Spectroscopy. University of Kent, 2015. Master thesis. 143 p.
Pedone A., Malavasi G., Menziani M.C., Cormack A.N., Segre U. A new self-consistent empirical interatomic potential model for oxides, silicates, and silica-based glasses // J. Physical Chemistry B. 2006. V. 110. P. 11780–11795.
Mishnev M., Korolev A., Bartashevich E., Ulrikh D. Effect of long-term thermal relaxation of epoxy binder on thermoelasticity of fiberglass plastics: multiscale modeling and experiments // Polymers. V. 14. P. 1712.
Wang J., Wolf R.M., Caldwell J.W., Kollman P.A., Case D.A. Development and testing of a general amber force field // J. Comput. Chem. 2004. V. 25. P. 1157–1174.
Soares T.A., Hünenberger P.H., Kastenholz M.A., Kräutler V., Lenz T., Lins R.D., Oostenbrink C., van Gunsteren W.F. An improved nucleic acid parameter set for the GROMOS force field // J. Comput. Chem., 2005. V. 26. P. 725–737.
Vanommeslaeghe K., Raman E.P., MacKerell A.D. Jr. Automation of the CHARMM General Force Field (CGenFF) II: Assignment of bonded parameters and partial atomic charges // J. Chemical Information and Modeling. 2012. V. 52. P. 3155–3168.
Wennberg C.L., Murtola T., Hess B., Lindahl E. Lennard–Jones lattice summation in bilayer simulations has critical effects on surface tension and lipid properties // J. Chem. Theory Comput. 2013. V. 9. P. 3527–3537.
Heinz H., Lin T.-J., Mishra R.K., Emami F.S. Thermodynamically consistent force fields for the assembly of inorganic, organic, and biological nanostructures: The INTERFACE force field // Langmuir. 2013. V. 29. P. 1754–1765.
Tsuneyuki S., Tsukada M., Aoki H., Matsui Y. First-principles interatomic potential of silica applied to molecular dynamics // Physical Review Letters, 1988. V. 61. P. 869–872.
Vaitkus A., Merkys A., Gražulis S. Validation of the Crystallography Open Database using the Crystallographic Information Framework // J. Applied Crystallography, 2021. V. 54. P. 661–672.
Smyth J.R., Jacobsen S.D., Hazen R.M. Comparative Crystal Chemistry of Dense Oxide Minerals // Reviews in Mineralogy and Geochemistry, 2000. V. 41. P. 157–186.
Gale J.D., Rohl A.L. The General Utility Lattice Program (GULP) // Molecular Simulation, 2003. V. 29. P. 291–341.
Abraham M., Murtola T., Schulz R., Páll S., Smith J., Hess B., Lindahl E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers // SoftwareX, 2015. V. 1–2. P. 19–25.
Bussi G., Donadio D., Parrinello M. Canonical sampling through velocity rescaling // J. Chem. Phys. 2007. V. 126. P. 014101.
Berendsen H.J.C., Postma J.P.M., van Gunsteren W.F., DiNola A., Haak J.R. Molecular dynamics with coupling to an external bath // J. Chem. Phys. 1984. V. 81. P. 3684–3690.
Diego Gatta G., Angel R.J., Rotiroti N., Carpenter M.A. High-pressure and low-temperature behaviour of trigonal kalsilite // Geophysical Research Abstracts, 2010. V. 12, EGU2010-12321.
Darden T., York D., Pedersen L. Particle mesh Ewald: An N⋅log(N) method for Ewald sums in large systems // J. Chem. Phys. 1993. V. 98. P. 10 089–10 092.
Wennberg C.L., Murtola T., Hess B., Lindahl E. Lennard–Jones Lattice Summation in Bilayer Simulations Has Critical Effects on Surface Tension and Lipid Properties // J. Chem. Theory Comput. 2013. V. 9. P. 3527–3537.
Shelby J.E. Introduction to Glass Science and Technology // Royal Society of Chemistry, 2 ed. 2007.
Lipinska-Kalita K.E., Kalita P., Hemmers O., Hartmann T. Equation of state of gallium oxide to 70 GPa: Comparison of quasihydrostatic and nonhydrostatic compression // Physical Review B. 2008. V. 77. P. 094123.
Barzilai S., Halevy I., Yeheskel O. Bulk modulus of Sc₂O₃: Ab initio calculations and experimental results // J. Applied Physics, 2011. V. 110. P. 043 532.
Palko J.W., Waltraud W.M., Sinogeikin S.V., Bass J.D., Sayir A. Elastic constants of yttria (Y₂O₃) monocrystals to high temperatures // J. Applied Physics, 2001. V. 89. P. 7791–7796.
Materials Data on Na₂SiO₃ by Materials Project // LBNL Materials Project; Lawrence Berkeley National Laboratory (LBNL), Berkeley, CA (United States), 2020.
Bass J.D. Elasticity of Minerals, Glasses, and Melts // In: Mineral Physics and Crystallography: A Handbook of Physical Constants. 1995. Eds. Ahrens T. J. Washington: American Geophysical Union. P. 45–63.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Физика и химия стекла

Архивы выпусков Информация о журнале Отправить рукопись в журнал