Физика и химия стекла. T. 49, Номер 3, 2023

Физика и химия стекла, 2023, T. 49, № 3, стр. 294-312

Моделирование расплавов SiO2 и процессов образования стекол методом молекулярной динамики

А. С. Мысовский 1*, А. С. Паклин 1

1 Институт геохимии им. А.П. Виноградова СО РАН
664033 Иркутск, ул. Фаворского, 1а, Россия

* E-mail: mysovsky@gmail.com

Поступила в редакцию 13.10.2022
После доработки 10.11.2022
Принята к публикации 08.02.2023

Полный текст (PDF)

Аннотация

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

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

ВВЕДЕНИЕ

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

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

Общепринятая гипотеза структуры кварцевого стекла в виде трехмерной непрерывной случайной сетки, предложенная Захариасеном [1], не предполагает каких-либо отличий в структуре стекла, полученных из различного фазового состава SiO2, поскольку считается, что базовой структурной единицей стекла является кремний-кислородный тетраэдр, а дальний порядок нарушен случайным распределением угла Si–O–Si в широком диапазоне (от 120° до 180°). Что было впоследствии подтверждено измерениями структуры кварцевого стекла методом рентгеновской дифракции в работе Уоррена [2].

На практике стекла, синтезированные из кварца и кристобалита, обладают рядом отличий. Такие стекла были, например, получены в нашем коллективе, и в работе [3] наблюдались отличия в парамагнитных свойствах радиационных дефектов. Впоследствии также были обнаружены ощутимые различия и в полосах оптического поглощения и люминесценции точечных дефектов [4]. В этой связи вопрос о структурных различиях и особенностях стекол, полученных из разных кристаллических модификаций SiO2, приобретает большую практическую значимость. В данной работе мы исследовали методом молекулярной динамики особенности плавления кварца и кристобалита, а также различия в полученных из них аморфных структурах, которые мы воспринимаем как модели стекол.

Есть, однако, и более фундаментальный и общий вопрос физики конденсированного состояния, на который данная работа может попытаться ответить: “помнит” ли стекло, как оно было кристаллом? А именно, сохраняются ли в стекле какие-либо структурные особенно той кристаллической фазы вещества, из которой оно было выплавлено? Разумеется, ответ на данный вопрос может оказаться совершенно различным для стекол разной химической композиции и с разной термической историей. В рамках данной краткой работы мы вряд ли дадим исчерпывающий ответ на этот весьма интересный вопрос для всех стекол из диоксида кремния, но, возможно, какие-то предположительные контуры ответов у нас и получится наметить.

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

Начиная с самой ранней работы (Woodcock, 1976 [5]), а затем в работах 1988–1990 гг. Feuston & Garofalini [6, 7], Vessal [8], Vashishta [9] молекулярная динамика успешно справлялась с генерацией аморфной структуры SiO2. В частности, полученная в этих работах радиальная функция распределения (Radial Distribution Function, RDF) хорошо согласовывалась с экспериментальной, однако в том, что касается механических свойств материала, согласие с экспериментом оказалось существенно хуже. При этом для короткодействующей части взаимодействия ионов использовались потенциалы Борна–Майера–Хаггинса (Born–Mayer–Huggins).

Ситуация улучшилась с появлением в 1990 г. потенциала BKS [10]. Данный потенциал не только успешно описывает различные кристаллические фазы SiO2 и реалистично предсказывают структуру аморфной фазы, но справляются также и с описанием механических свойств этих материалов. Короткодействующая часть BKS основана на потенциалах Бакингема (Buckingham). Забегая вперед, можно сказать, что BKS оказался на удивление удачным потенциалом, несмотря на свою простоту и древность. Не так давно в работе [11] проводилось сравнение ряда потенциалов, предназначенных для описания SiO2. Эти потенциалы включали BKS, PMMCS [12], а также потенциалы из работ [1316]. В описании статических свойств кристаллического и аморфного SiO2 потенциал PMMCS оказался точнее, чем BKS, хотя и BKS везде дает качественно верные и количественно удовлетворительные результаты. Однако PMMCS на удивление плохо воспроизводит фононную плотность состояний материала, с чем BKS, опять же, на уровне качественного согласия, вполне справляется. В работе [11] сделан вывод, что из рассмотренных потенциалов только BKS адекватно описывает как статические, так и динамические свойства SiO2.

Еще до работы [11] потенциалы для SiO2 были подогнаны на основе неэмпирических расчетов молекулярной динамики Кара-Паринелло [17]. Использовалось два подхода: FM (force matching) – подгонка по силам и CHIK [18] – подгонка по структурным данным. Многообещающий, на первый взгляд, подход FM дал удивительно плохие, качественно неверные результаты. Потенциал CHIK оказался, возможно, чуть лучше, чем BKS. Еще раньше аналогичные подходы и потенциалы сравнивались в работе [19]. Сравнивались потенциалы BKS, CHIK, Soules 1990 [20] и Takada 2004 [21]. Сделан вывод, что BKS и CHIK воспроизводят основные свойства SiO2, однако делают это при завышенных, по сравнению с экспериментом, температурах. Известны также попытки использовать потенциалы Терсофа [22], параметризованные на основе неэмпирического расчета [23], и подгонять параметры потенциалов по неэмпирически рассчитанным силам и их градиентам [24].

По истории молекулярной динамики кристаллического и аморфного SiO2 с подобными потенциалами имеются обзорные работы, например [25, 26]. Следует отметить также использование различных параметризаций потенциалов Китинга [27] для моделирования аморфной структуры SiO2 по методологии WWW [28]. Но данный подход не имеет ничего общего с молекулярной динамикой, а является разновидностью методов Монте-Карло.

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

Ситуация улучшилась с появлением потенциалов ReaxFF, предложенных впервые для углеводородов в работе [29]. Успех потенциалов ReaxFF не случаен. Идея применять в молекулярной динамике эквилибрацию зарядов, видимо, “носилась в воздухе”. Например, в 2003 г. в работе [30] предлагалось использовать для этого метод QE [31].

Для SiO2 потенциалы типа ReaxFF появились впервые также в 2003 г. в работе [32]. Затем они дополнялись и обобщались, например, для взаимодействия SiO2–вода [33]. С помощью таких потенциалов был исследован процесс гидроксилирования и растворения в воде поверхности аморфного SiO2. Затем потенциалы ReaxFF для этих систем улучшались и дорабатывались [34, 35], что позволило исследовать процесс растворения диоксида кремния и боросиликатных стекол через формирования пористых областей [36].

Более подробную информацию о развитии потенциалов ReaxFF и их применении для исследования стекол, и в частности аморфного SiO2, можно найти в недавнем обзоре [37]. Приведенные данные убеждают, что эти потенциалы позволяют правильно описывать многие свойства стекол и происходящие в них процессы. При этом все же отметим, что полноценное сравнение структур аморфного SiO2, полученных с помощью ReaxFF и BKS, с экспериментом и друг с другом, насколько нам известно, не проводилось.

Плавление кристаллического SiO2 в молекулярной динамике приходится моделировать при температурах, значительно превышающих экспериментальную температуру плавления, вплоть до значения 7500 К [25]. При этом начальная стадия плавления в молекулярной динамике едва ли сопоставима с экспериментом: экспериментально плавление начинается с поверхности, межзеренных границ и прочих дефектов [38] и является гетерогенным процессом [39]. Однако полученная в ходе симуляции аморфная структура уже, как правило, соответствует экспериментальным данным, о чем говорят функции парного и углового распределения [25], структурный фактор [40]. В структуре доминируют 5- и 6-членные кольца, есть в небольшом количестве и 3-членные, которые могут отвечать за пик комбинационного рассеяния 606 см1 [41]. Не все потенциалы одинаково удачны. Так, полученная в работе [42] с помощью потенциалов CMAS94 [43] плотность аморфного SiO2 величиной 2.39 г/см3 представляется существенно завышенной.

Наряду с классической молекулярной динамикой, для SiO2 велись и расчеты из первых принципов. В 1995 г. Sarnthein, Pasquarello и Car [44, 45] использовали для этого приближение LDA и ячейку из 72 атомов. После нагрева до 3500 К и охлаждения им удалось получить аморфную структуру без дефектов координации, при этом структурный фактор согласовывался с экспериментом. Молекулярную динамику из первых принципов использовали для расчета инфракрасных и рамановских спектров [46, 47], зависимости температуры плавления от давления [48], сравнения структуры расплавов силикатов и SiO2 [49] и других исследований.

В ряде работ исследовалась диффузия атомов кремния и кислорода [5058]. Эти данные мы приведем чуть позже, когда будем обсуждать диффузию в моделируемых в данной работе расплавах.

В серии основательных работ [5560] процесс перехода расплав–стекло исследован с фундаментальной точки зрения. Выдвинут ряд аргументов в пользу MCT (mode coupling theory, теория связанных мод) [61], в частности показано [60], что зависимость температуры стеклования от скорости охлаждения плохо описывается соотношением Вогеля–Фалчера–Тамманна (Vogel–Fulcher–Tammann) [62], и гораздо лучше описывается с помощью MCT.

В последнее время в молекулярной динамике SiO2 интерес к моделированию структуры объемного материала снизился в связи с большим количеством уже посвященных этому работ, зато исследуется поверхность SiO2 и ее взаимодействие с различными средами [6365], пористый оксид кремния [66, 64 ], формирование треков заряженных частиц при облучении материала [67]. Ряд работ посвящен SiO2 при высоких давлениях [48, 68, 69]. Несколько ранее была предсказана экзотическая фаза SiO2, существующая только в условиях неаксиального стресса [70]. Молекулярная динамика может использоваться также в сочетании с другими расчетно-теоретическими методами, например, для исследования электронных и оптических свойств разнообразных дефектов в SiO2 [7173].

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

МЕТОД РАСЧЕТА

Расчеты выполнялись методом молекулярной динамики (МД, MD) в программе LAMMPS [74]. Взаимодействие атомов кислорода и кремния моделировалось с помощью реактивных потенциалов ReaxFF [29, 32]. Как уже указывалось выше, данная методика успешно применялась в ряде работ [34, 36, 7577, и др]. Достоинством потенциалов ReaxFF является способность описывать реакции, например, разрыв и образование ковалентных связей. Был выбран шаг интегрирования по времени в размере 0.5 фс. Расчетная ячейка кристобалита содержала 3000 атомов и имела форму куба с ребром 37.3 Å. Ячейка кварца содержала 3024 атома и имела близкую к кубической форму и схожие размеры. На расчетную ячейку были наложены периодические граничные условия.

Для моделирования синтеза стекол с кислородно-дефецитными центрами (ODC, oxygen deficiency centers) из ячеек бездефектного кварца и кристобалита удалялось по 30 случайно выбранных атомов кислорода. Таким образом, содержание вакансий кислорода в этих ячейках составило около 1.5% от общего числа атомов кислорода, и мы будем называть их “кварц с 1.5% ODC” и “кристобалит с 1.5% ODC”.

Начальная эквилибрация системы производилась при температурах 300 или 500 K в течении 40 пс в ансамбле NVT, а затем еще 40 пс в ансамбле NPT. После этого в течении 200 пс производился нагрев системы до температур от 3500 до 5000 K. При таких температурах расчетная ячейка кварца или кристобалита выдерживалась в течении 500 пс. При необходимости это время могло быть увеличено до 750 пс, или же, наоборот, уменьшено, после чего во всех случаях фиксировалось разрушение кристаллической структуры и образование расплава. Охлаждение до низких температур (1 K) производилось в течении 50 пс.

Таким образом, были смоделированы 4 основные плавки: бездефектного кварца и кварца с 1.5% ODC при 4050 K, а также бездефектного кристобалита и кристобалита с 1.5% ODC при 4000 K. Кроме того, был смоделирован ряд вспомогательных плавок – для бездефектных ячеек при температурах 4200 и 4500 K, а для ячеек с ODC при более низких температурах: 3500 и 3900 K. В результате всех плавок были получены аморфные структуры SiO2, которые затем можно было анализировать. Однако, при данном подходе невозможно контролировать степень аморфизации расплава, и структурные особенности стекол из кварца и кристобалита оказываются утрачены.

Более полезной оказалась техника взятия мгновенных снимков системы, или “снапшотов”, в ходе перечисленных выше основных плавок. На протяжении всего плавления каждые 1000 шагов записывался файл рестарта, позволяющий возобновить расчет с данного момента. Для каждого из рассмотренных материалов (бездефектные кварц и кристобалит, они же с 1.5% ODC) было отобрано от 6 до 8 снимков. Затем для этих снимков была сделана очень быстрая заморозка (за 10 пс), а также отжиг при температурах 3000, 3250 и 3500 K длительностью в 250 пс с последующим охлаждением в течении 100 пс. Анализируя аморфные структуры, полученные из последовательно взятых снимков, удается шаг за шагом отслеживать, как именно изменяется структура SiO2 в ходе плавления и отжига.

ДИФФУЗИЯ В РАСПЛАВЕ SiO2

Для того, чтобы наблюдать за диффузией атомов на всем протяжении расчетов молекулярной динамики, был организован регулярный (каждые 100 шагов МД) вывод такой величины, как средний квадрат отклонения (СКО) атомов от своих начальных позиций. Усреднение производилось отдельно для атомов кислорода и кремния, что дало две величины: $\left\langle {\Delta r_{{\text{O}}}^{2}} \right\rangle $ и $\left\langle {\Delta r_{{{\text{Si}}}}^{2}} \right\rangle .$ Если бы коэффициент диффузии (КД) во время моделирования был постоянным, то эти величины подчинялись бы хорошо известному соотношению ∆r2 = 6Dt. Однако мы моделируем процесс плавления кристалла и застывания расплава, поэтому КД очевидным образом не должен быть постоянным. В этом случае для его вычисления нужно использовать соотношение

(1)
$D = \frac{1}{6}\frac{{d\left\langle {\Delta {{r}^{2}}} \right\rangle }}{{dt}},$

На врезке рис. 1а приведена зависимость $\left\langle {\Delta r_{{\text{O}}}^{2}} \right\rangle $ от времени для случая плавления кварца с ODC. Эта величина флуктуирует, как и следует ожидать, в силу хаотичного теплового движения атомов, что затрудняет ее дифференцирование по времени. Различные методики численного дифференцирования шумных данных не дают здесь хороших результатов, поэтому эти данные были подогнаны аналитической функцией, которая и была затем продифференцирована для нахождения коэффициента диффузии.

Рис. 1.

Поведение коэффициента диффузии атомов кислорода и плотности расплава при плавлении кварца с 1.5% ODC при различных температурах (a); бездефектного кристобалита при 4000 K (б). На врезках приведены соответствующие кривые СКО атомов кислорода.

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

$f\left( t \right) = \left\{ \begin{gathered} {{a}_{1}}t + {{b}_{1}},\,\,\,\,t < {{t}_{1}} \hfill \\ {{P}_{4}}\left( t \right),\,\,\,\,{{t}_{1}} \leqslant t < {{t}_{2}} \hfill \\ {{a}_{3}}t + {{b}_{3}},\,\,\,{{t}_{2}} \leqslant t \hfill \\ \end{gathered} \right..$
Здесь P4(t) – полином четвертой степени. Коэффициенты этого полинома, а также моменты времени t1 и t2 были подгоночными параметрами при нахождении данной функции. Параметры a1, b1, a3, b3 подобраны так, чтобы обеспечить непрерывность f(t) и ее производной. На рис. 1a показано поведение коэффициента диффузии при температурах 3500, 3900 и 4050 K (кривые 1–3), полученного из такой функции. На том же графике показано поведение плотности в ходе плавления (кривые 4–6), и корреляция между изменением плотности и ростом коэффициента диффузии видна очень хорошо. При плавлении кристобалита с 1.5% ODC поведение СКО, плотности и коэффициента диффузии выглядит аналогично.

Это позволяет разделить наблюдаемый процесс на три хорошо определенных стадии, которые мы будем называть в дальнейшем стадиями 1, 2 и 3:

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

Стадия активного плавления, которая характеризуется резким ростом КД и одновременным изменением плотности. При этом происходит разрушение кристаллической решетки.

Стадия равновесного расплава, когда все, что могло, уже расплавилось. КД и плотность больше не изменяются. Расплав достигает состояния термодинамического равновесия.

Следует отметить, что стадию 1 едва ли возможно напрямую сравнивать с экспериментальными наблюдениями. Как уже отмечалось во введении, экспериментально плавление начинается с поверхности, межзеренных границ и других дефектов. Однако, начиная со стадии 2 сопоставление симуляции с экспериментом вполне возможно, во всяком случае, никаких очевидных препятствий к этому нет.

Плавление кристобалита с 1.5% ODC качественно подобно плавлению кварца с ODC. Плавление же бездефектных кварца и кристобалита несколько отличается от них на начальной стадии, а именно, подогнать зависимость СКО от времени прямой линией не представляется возможным. Это означает, что на начальной стадии плавления бездефектного кварца и кристобалита КД не имеет постоянного значения, а увеличивается, начиная с близкого к нулю значения. Для подгонки этих кривых были использованы полиномы 5ой степени для стадий 1 и 2 и линейная функция для стадии 3. В моменты t1 и t2 эти функции сшивались по значению, а также по первой и второй производной. Результаты подгонки показаны на рис. 1б. Стадии 2 и 3 во всех случаях (плавление кристобалита и кварца; бездефектного и с 1.5% ODC) выглядят очень похоже. Также на рис. 1б показан вертикальными пунктирными линиями выбор моментов времени для взятия снимков – мы старались отобрать их таким образом, чтобы покрыть ими стадию активного плавления.

В результате проделанного анализа данных оказалось возможным оценить коэффициенты диффузии для: 1) бездефектного кварца на стадии 3 при температурах 4050, 4200, 4500 и 5000 K; 2) бездефектного кристобалита на стадии 3 при температурах 4000, 4200, 4500 и 5000 K; 3) кварца с 1.5% ODC на стадиях 1 и 3 при температурах 4050, 3900 и 3500 K; 4) кристобалита с 1.5% ODC на стадиях 1 и 3 при температурах 4000, 3900 и 3500 K. Это позволяет построить для всех этих систем диаграммы Аррениуса (рис. 2) и определить энергии активации, а также префактор для коэффициента диффузии.

Рис. 2.

Диаграммы Аррениуса для коэффициента диффузии при плавлении бездефектных кварца и кристобалита (a); кварца с 1.5% ODC (б); кристобалита с 1.5% ODC (в).

Все эти данные представлены в табл. 1. Для каждого из расплавов приведены значения энергии активации ($E_{{\text{O}}}^{A}$ и $E_{{{\text{Si}}}}^{A}$) и префактора ($D_{{\text{O}}}^{0}$ и $D_{{{\text{Si}}}}^{0}$) для атомов кислорода и кремния соответственно. Следует отметить несколько моментов касательно этих данных. Во-первых, определение КД на начальной стадии плавления, судя по всему, сопровождается значительной численной ошибкой. Об этом говорит большой разброс между $E_{{\text{O}}}^{A}$ и $E_{{{\text{Si}}}}^{A}$ в случае кварца и отличающееся на три порядка значение префактора для атомов кремния. Во-вторых, в конечной стадии плавления все энергии активации для расплава кристобалита оказываются выше соответствующих величин в расплаве кварца на 0.8–1 эВ. Это не совсем понятно и свидетельствует либо об, опять же, численной ошибке, либо о том, что данные расплавы в ходе симуляции все же не достигли полного термодинамического равновесия, иначе систематическое различие между ними трудно объяснить. В третьих, наличие ODC понижает энергию активации диффузии для атомов кремния (но не кислорода) в расплавах как кварца, так и кристобалита, примерно на 0.3–0.5 эВ. Данная закономерность вполне может реальной и имеет разумное объяснение – наличие вакансий кислорода может облегчать диффузию атомов кремния.

Таблица 1.  

Энергии активации и префакторы для коэффициента диффузии в моделируемых расплавах

Стадия плавления Кристалл Наличие ODC $E_{{\text{O}}}^{A}$, эВ $D_{{\text{O}}}^{0},$ Å2/пс $E_{{{\text{Si}}}}^{A}$, эВ $D_{{{\text{Si}}}}^{0}$, Å2/пс
Начальная Кварц 1.5% 3.48 1.7 × 102 6.15 1.8 × 105
  Кристобалит 1.5% 4.23 3.3 × 103 4.41 9.9 × 102
Конечная Кварц Нет 4.64 4.9 × 104 5.26 1.5 × 105
    1.5% 4.98 1.6 × 105 4.80 6.3 × 104
  Кристобалит Нет 5.48 5.3 × 105 6.29 3.1 × 106
    1.5% 5.70 1.4 × 106 5.95 1.9 × 106

Учитывая значительную численную ошибку представленных данных, делать окончательные выводы на их основании пока рано. Требуются дополнительные расчеты, с большими временами выдержки расплавов и с большим количеством значений температуры, для уточнения этих результатов. Однако, в пределах ошибки определения полученные значения энергии активации для диффузии в расплаве SiO2 удовлетворительно согласуются с экспериментальным значением 4.7 эВ [82], которая пользуется наибольшим доверием исследователей. Также они хорошо согласуются и с заключением работы [78] о том, что наиболее надежное значение для энергии активации вязкости расплава SiO2 составляет 5.34 эВ. Учитывая, что вязкость и коэффициент диффузии связаны друг с другом соотношением Стокса–Эйнштейна, их энергии активации должны быть, по крайней мере, близки друг к другу.

Затем были определены коэффициенты диффузии для всех снимков при различных температурах отжига. На рис. 3а показаны данные по СКО атомов кислорода для снимков кварца с 1.5% ODC при температурах 3000, 3250 и 3500 K. Эти зависимости хорошо и достоверно подгоняются прямыми линиями, что дает значения КД при данных температурах, а значение при 4050 K (или 4000 K для кристобалита) вычислено при помощи подогнанной функции (1). По всем этим данным были построены диаграммы Аррениуса (пример для кристобалита с ODC показан на рис. 3б, остальные три выглядят похоже) и были определены энергии активации и префакторы КД.

Рис. 3.

Определение коэффициентов диффузии для снимков расплава кристобалита с 1.5% ODC (a) и построенная из них диаграмма Аррениуса (б).

Как оказалось, все энергии активации диффузии, определенные таким образом, имеют величину около 3 эВ, которая существенно отличается от определенной выше величины около 5 эВ. При этом, если для каждого снимка индивидуально подгонять энергию активации, то они варьируются от 2.5 до 3.6 эВ. Однако последовательные снимки, взятые с одной плавки, хорошо подгоняются и одной общей энергией активации. Для кварца и кристобалита без дефектов данные энергии активации составляют 3.31 и 3.21 эВ, а для кварца и кристобалита с ODC – 3.12 и 2.85 эВ соответственно.

Самое подходящее время обсудить полученные результаты и сопоставить их с литературными данными. Различие между энергиями активации диффузии для расплавов (около 5 эВ) и для снимков (около 3 эВ), на первый взгляд, обескураживающее, легко объяснимо. Это давно известное различие между равновесными и изоструктурными характеристиками стекол, подробно разъясненное, например, в книге [79]. Дело в том, что расплавы SiO2, выдержанные при разных температурах, структурно отличаются друг от друга, поскольку их структура успевает подстроиться под заданную температуру. По ним получено равновесное значение $E_{{eq}}^{A} \approx 5$ эВ, которое характеризует изменение КД в разных структурах при разных температурах. В то же время снимки, почти мгновенно охлажденные до температуры 3000, 3250 или 3500 K, не успели претерпеть структурных изменений, характерных для их новой температуры. Поэтому можно сказать, что значение $E_{{iso}}^{A} \approx 3$ эВ получено для одной и той же структуры при разных температурах и является, таким образом, изоструктурным.

Данное обстоятельство было известно еще в 1951 г. [80], когда для наблюдаемого отличия энергий активации диффузии в расплаве и в стекле, а также для неаррениусового поведения КД в области перехода было предложено похожее объяснение. Однако, обращает на себя внимание то обстоятельство, что многие зарубежные авторы не акцентируют внимания на различии равновесных и изоструктурных характеристик стекол, а в отдельных случаях, возможно, и не подозревают о нем. Так, например, авторы расчетно-теоретической работы [81] рассчитывают энергию активации диффузии кислорода методом функционала плотности в локальном приближении (DFT-LDA), однозначно понимая ее как энергетический барьер для этого процесса. Получив значение 3.8 эВ (в неплохом согласии с изоструктурной величиной, причем, судя по графикам, эта величина должна быть еще меньше), они, тем не менее, находят согласие с экспериментом неудовлетворительным и начинают манипулировать положением уровня Ферми. В результате они все-таки добиваются значения 4.6 эВ, в согласии с результатами работы [82], причем эта величина равновесная и уже не имеет смысла энергетического барьера.

Экспериментальные данные по диффузии кислорода в кристаллическом и аморфном SiO2 имеют весьма широкий разброс энергий активации: от 1.25 [83] до 4.7 эВ [82]. Однако в действительности это не должно обескураживать. Значения в районе 1–2 эВ [83, 84] относятся к диффузии молекулярного кислорода. Затем следует группа значений в районе 3–3.5 эВ (например, [85, 86]). К этим экспериментальным данным примыкают и некоторые теоретические результаты, например, в работе [54] в ситуации, аналогичной нашим “снимкам”, получена энергия активации 3.09 ± 0.18 эВ. Мы полагаем, что все эти результаты правильные и согласуются между собой в пределах точности использованных методов.

Точно так же правильной мы считаем и величину 4.7 эВ, полученную в результате очень тщательных измерений в работе [82]. Важная для интерпретации результатов деталь: в этой работе Mikkelsen длительно выдерживал образец при каждой из промеренных температур. Несколько выше, но в том же интервале, лежит величина 5.34 эВ, полученная для энергии активации коэффициента вязкости расплава SiO2 [87]. В обзоре [88] именно эта величина выделена как наиболее надежная. Экспериментальное значение энергии активации диффузии атомов кремния составляет 5.2 эВ, по данным работы [89]. Большинство теоретических работ также дают близкие к 5 эВ величины: 4.9 эВ [50, 51], 5.2 эВ [52], 4.77 эВ [53]. В работах [5558] получены энергии активации 4.66 и 5.18 эВ для диффузии атомов кислорода и кремния соответственно.

Итак, мы полагаем, что все эти величины правильные. Просто величины в районе 3–3.5 эВ относятся к изоструктурной, а в районе 5 эВ – к равновесной энергии активации диффузии. Длительная выдержка образца при заданной температуре (работа [82], величина 4.7 эВ) как раз и означает, что структура успела перестроиться в состояние термодинамического равновесия при этой температуре. Наши результаты, таким образом, находятся в согласии с экспериментальными и теоретическими данными других авторов, и тоже ярко демонстрируют отличие равновесных и изоструктурных характеристик расплавов и стекол.

Разность равновесной $E_{{eq}}^{A} \approx 5$ эВ и изоструктурной $E_{{iso}}^{A} \approx 3$ эВ энергий активации диффузии можно назвать структурным или конфигурационным вкладом в энергию активации диффузии атомов кислорода $\Delta E_{{conf}}^{A} = E_{{eq}}^{A} - E_{{iso}}^{A}.$ Интересно здесь то, что эта величина позволяет оценить температурную зависимость конфигурационной энтропии, связанной с актом диффузии. Изоструктурная энергия активации EA представляет собой энергетический барьер для акта диффузии. Выразим равновесную и изоструктурную зависимости КД:

${{D}_{{eq}}} = D_{{eq}}^{0}{{e}^{{ - E_{{eq}}^{A}/kT}}},\,\,\,\,{{D}_{{iso}}} = D_{{iso}}^{0}{{e}^{{ - E_{{iso}}^{A}/kT}}}.$

С другой стороны, те же величины можно выразить через барьер свободной энергии Гиббса $\Delta G = \Delta H - T\Delta S.$ Энтальпия диффузии ∆H одинакова для равновесного и изоструктурного КД, и, по-видимому, должна совпадать с изоструктурной энергией активации:

${{D}_{{eq}}}\sim {{e}^{{ - \Delta H/kT}}}{{e}^{{\Delta S\left( T \right)/k}}},\,\,\,\,{{D}_{{iso}}}\sim {{e}^{{ - \Delta H/kT}}}{{e}^{{\Delta S\left( {{{T}_{0}}} \right)/k}}}.$

Энтропия диффузии разделяется на колебательный и конфигурационный вклады $\Delta S = \Delta {{S}_{{vib}}} + \Delta {{S}_{{conf}}},$ причем колебательный вклад в гармоническом приближении определяется отношением частот колебаний в переходном и в начальном состоянии и от температуры не зависит. Также не зависит от температуры величина $\Delta S_{{conf}}^{{iso}} = \Delta {{S}_{{conf}}}({{T}_{0}})$ – конфигурационная энтропия в изоструктурном случае, поскольку здесь сама структура не успевает подстроиться к изменяющейся температуре и соответствует начальной температуре эксперимента T0. Логарифмируя отношение двух КД, получим

$\Delta {{S}_{{conf}}}\left( T \right) = - \frac{{\Delta E_{{conf}}^{A}}}{T} + const.$

Данная функция является монотонно возрастающей и означает, что с ростом температуры структура расплава перестраивается таким образом, что диффузия становится более интенсивной. При этом барьер диффузии не изменяется, изменяется только число доступных для диффузии конфигураций. Эта же зависимость обеспечивает Аррениусовский характер поведения величины Deq(T). Отметим также, что в большинстве работ конфигурационная энтропия считается независящей от температуры, поскольку она определяется только числом доступных для диффузии конфигураций. Эта величина очевидным образом не зависит от температуры для кристаллов, и это соображение, по всей видимости, автоматически переносится некоторыми исследователями и на аморфные структуры, что уже не вполне верно.

ХАРАКТЕРИСТИКИ АМОРФНЫХ СТРУКТУР

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

Ближний порядок можно охарактеризовать с помощью RDF. Для мгновенно замороженных снимков кварца и кристобалита с 1.5% ODC радиальные функции распределения приведены на рис. 4. По ним очень хорошо видно, что если в начале стадии активного плавления они имеют характерные для соответствующей кристаллической фазы максимумы и минимумы, то к ее окончанию эта структурная информация полностью утрачивается. Для наглядности на врезке рис. 4a приведены кривые для кварца и кристобалита в начале плавления, и они существенно отличаются друг от друга. На врезке же рис. 4б аналогичные же кривые приведены в конце плавления, и они практически неразличимы. Частичные функции радиального распределения (т.е. построенные отдельно для пар Si–Si, Si–O, O–O), а также функции углового распределения в данной статье не приводятся, но они демонстрируют точно такое же поведение.

Рис. 4.

Радиальные функции распределения для снимков кварца с 1.5% ODC (a); кристобалита с 1.5% ODC (б). На врезках: сравнение RDF в начальных снимках кварца и кристобалита с ODC (a); то же в конечных снимках (б).

Данное поведение означает, что ближний порядок расплавов кварца и кристобалита к концу стадии активного плавления становится одинаковым. Это полностью согласуется с нашим представлением о том, что этому моменту расплавы достигают термодинамического равновесия. В таком состоянии между ними и не должно быть структурных различий, иначе говоря, к данному моменту структурная “память” расплава оказывается полностью утерянной.

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

Это становится видно, если проанализировать статистику колец в полученных аморфных структурах. С этой целью в них все были найдены К-кольца, или кольца Кинга [90]. Статистика колец для всех снимков кварца и кристобалита с 1.5% ODC показана на рис. 5. На каждой из гистограмм по оси абсцисс отложен размер, или число звеньев кольца (звеном считается элемент Si–O, т.е. атом кремния вместе с мостиковым атомом кислорода, осуществляющим ковалентную связь с соседним звеном). По оси ординат отложено число таких колец в структуре, приходящихся на один атом кремния.

Рис. 5.

Статистика колец для снимков кварца с 1.5% ODC (a); кристобалита с 1.5% ODC (б)

Первая особенность, на которую стоит обратить внимание – это избыток 8-звенных колец в структурах, полученных из кварца. В кристаллическом кварце на один атом кремния приходится одно 6-звенное и 7 штук 8-звенных колец. На гистограммах рис. 5a видно, что избыток 8-звенных колец уменьшается от снимка к снимку, но все же сохраняется даже еще в предпоследнем снимке № 455. При этом радиальные функции распределения практически сливаются с таковыми для равновесного расплава уже к снимку № 415. Это говорит о том, что топологические особенности расплава сохраняются дольше, чем его ближний порядок, и могут затем обнаруживаться и в стеклах.

В кристобалите же единственный имеющийся вид колец – это 6-звенные, их приходится по 2 штуки на один атом кремния. Пик 6-звенных колец разрушается быстро, и уже к третьему снимку (№ 430) он отсутствует. Возможно, это говорит о том, что потеря структурной “памяти” у расплава кристобалита происходит быстрее, чем у расплава кварца, и кристобалит при плавлении практически сразу же превращается в равновесный расплав.

Ближе к концу активного плавления (последние 3 снимка для кварца и последние 6 – для кристобалита) обнаруживается еще одно немаловажное обстоятельство. А именно, соотношение 4-звенных и 5-звенных колец сильно зависит от температуры, при которой был выдержан расплав. “Мгновенно замороженные” снимки (instant) до этого плавились при 4000 или 4050 K. В них 4-звенные кольца доминируют над 5-звенными. По мере убывания температуры отжига (3500, 3250 и 3000 K) количество 4-звенных колец монотонно убывает, а количество 5-звенных – растет. В какой-то мере одновременно с 5-звенными возрастает и число 6-ти и 7-звенных колец. Соотношение 4-звенных и 5/6/7-звенных колец может являться, пожалуй, самым ярким температурным маркером расплава и полученных из него стекол. Следует отметить, что это соотношение одинаковым образом зависит от температуры и в равновесном расплаве, и в расплаве со структурной памятью. Именно это обстоятельство наиболее отчетливо иллюстрирует то, что структурные особенности расплава меняются с ростом температуры, что и приводит к различию изоструктурных и равновесных энергий активации диффузии.

Наконец, мы рассмотрели дефектные, то есть недокоординированные, атомы кремния. За редкими исключениями дефекты во всех структурах с 1.5% ODC (из которых были удалены по 30 атомов кислорода) сводятся к 60 атомам кремния, координированным тремя атомами кислорода. Для каждого такого атома можно определить расстояние rmin до ближайшего соседнего атома кремния. Если это проделать, а потом усреднить по всей структуре, то между расплавами кварца и кристобалита не обнаруживается никакой разницы (рис. 6, нижние кривые).

Рис. 6.

Среднее расстояние между 3-координированным атомом кремния и ближайшим несвязанным с ним через кислородный мостик атомом кремния (〈rNB〉) и усредненное наименьшее расстояние между 3-координированным атомом кремния и другим атомом кремния (〈rmin〉).

Но если поставить вопрос несколько иначе – на каком расстоянии от данного 3-координированного атома кремния находится ближайший к нему атом кремния, не связанный с ним через кислородный мостик, то между расплавами кварца и кристобалита обнаружится различие. А именно, данное расстояние rNB, усредненное по всей структуре, для снимков кристобалита оказывается примерно на 0.2 Å больше, чем для снимков кварца. Это видно по двум верхним кривым на рис. 6.

Надо сказать, что именно эти, не связанные кислородным мостиком атомы кремния, и должны образовывать кислородно-дефицитный центр ODC. На том же рисунке показаны и характерные структуры ODC, взятые из этих снимков. Наименьшие значения rNB отвечают паре 3-координированных атомов кремния. Такие структуры оказываются ближе к модели ODC(I) центра [91].

Большие значения rNB отвечают или паре 3-координированных атомов, связанных короткой цепочкой звеньев Si–O, либо паре 3-х и 4-координированного атомов кремния. Все это существенно влияет на тип и свойства ODC. Таким образом, указанное различие может оказаться значимым для объяснения разницы в оптических свойствах стекол, синтезированных из кварца и из кристобалита.

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

ЗАКЛЮЧЕНИЕ

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

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

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

Соотношение количества 4-звенных колец к 5-звенным в равновесном расплаве монотонно растет с температурой расплава.

Структуры ODC в стеклах из кварца и кристобалита несколько отличаются. Среднее расстояние между атомами кремния, образующими ODC, в стеклах из кварца меньше, а кратчайшая цепочка звеньев Si–O, соединяющих эти атомы, наоборот, длиннее.

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

Исследование проведено в рамках выполнения государственного задания по Проекту № 0284-2021-0004. Расчеты проведены на вычислительном кластере “Академик В.М. Матросов” (Иркутский суперкомпьютерный центр СО РАН).

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

  1. Zachariasen W.H. The atomic arrangement in glass // J. American Chemical Society. 1932. V. 54. № 10. P. 3841–3851.

  2. Warren B.E. X-ray diffraction of vitreous silica // Zeitschrift für Kristallographie-Crystalline Materials. 1933. V. 86. № 1–6. P. 349–358.

  3. Mashkovtsev R.I. Nepomnyashchikh A.I., Zhaboedov A.P., Paklin A.S. EPR study of the E’defects in optical glasses and cristobalite // Europhysics Letters. 2021. V. 133. № 1. P. 14003.

  4. Garmysheva T.Y., Nepomnyashchikh A.I., Shalaev A., Kaneva E., Paklin A., Chernenko K., Kozlova A.P., Pankratov V., Shendrik R. Luminescence of ODC (II) in quartz and cristobalite glasses // J. Non-Crystalline Solids. 2022. V. 575. P. 121199.

  5. Woodcock L.V., Angell C.A., Cheeseman P. Molecular dynamics studies of the vitreous state: Simple ionic systems and silica // The J. Chemical Physics. 1976. V. 65. № 4. P. 1565–1577.

  6. Feuston B.P., Garofalini S.H. Empirical three-body potential for vitreous silica // The J. Chemical Physics. 1988. V. 89. № 9. P. 5818–5824.

  7. Feuston B.P., Garofalini S.H. Oligomerization in silica sols // J. Physical Chemistry. 1990. V. 94. № 13. P. 5351–5356.

  8. Vessal B., Amini, M., Fincham D., Catlow C.R.A. Water-like melting behaviour of SiO2 investigated by the molecular dynamics simulation technique // Philosophical Magazine B. 1989. V. 60. № 6. P. 753–775.

  9. Vashishta P.P., Kalia R.K., Rino J.P., Ebbsjö I. Interaction potential for SiO2: A molecular-dynamics study of structural correlations // Physical Review B. 1990. V. 41. № 17. P. 12197.

  10. 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. № 16. P. 1955.

  11. Afify N.D., Mountjoy G., Haworth R. Selecting reliable interatomic potentials for classical molecular dynamics simulations of glasses: The case of amorphous SiO2 // Computational Materials Science. 2017. V. 128. P. 75–80.

  12. 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 // The J. Physical Chemistry B. 2006. V. 110. № 24. P. 11780–11795.

  13. Tsuneyuki S., Tsukada M., Aoki H., Matsui Y. First-principles interatomic potential of silica applied to molecular dynamics // Physical Review Letters. 1988. V. 61. № 7. P. 869.

  14. Cormack A.N., Du J., Zeitler T.R. Alkali ion migration mechanisms in silicate glasses probed by molecular dynamics simulations // Physical Chemistry Chemical Physics. 2002. V. 4. № 14. P. 3193–3197.

  15. Flikkema E., Bromley S.T. A new interatomic potential for nanoscale silica // Chemical Physics Letters. 2003. V. 378. № 5–6. P. 622–629.

  16. Du J., Cormack A.N. The structure of erbium doped sodium silicate glasses // J. Non-Crystalline Solids. 2005. V. 351. № 27–29. P. 2263–2276.

  17. Carré A., Ispas S., Horbach J., Kob W. Developing empirical potentials from ab initio simulations: The case of amorphous silica // Computational Materials Science. 2016. V. 124. P. 323–334.

  18. Carre A., Horbach J., Ispas S., Kob W. New fitting scheme to obtain effective potential from Car-Parrinello molecular-dynamics simulations: Application to silica // EPL (Europhysics Letters). 2008. V. 82. № 1. P. 17001.

  19. Soules T.F., Gilmer G.H., Matthews M.J., Stolken J.S., Feit M.D. Silica molecular dynamic force fields – A practical assessment // J. Non-Crystalline Solids. 2011. V. 357. № 6. P. 1564–1573.

  20. Soules T.F. Computer simulation of glass structures // J. Non-Crystalline Solids. 1990. V. 123. № 1–3. P. 48–70.

  21. Takada A., Richet P., Catlow C.R.A., Price G.D. Molecular dynamics simulations of vitreous silica structures // J. Non-Crystalline Solids. 2004. V. 345. P. 224–229.

  22. Tersoff J. Empirical interatomic potential for carbon, with applications to amorphous carbon // Physical Review Letters. 1988. V. 61. № 25. P. 2879.

  23. Munetoh S., Motooka T., Moriguchi K., Shintani A. Interatomic potential for Si–O systems using Tersoff parameterization // Computational Materials Science. 2007. V. 39. № 2. P. 334–339.

  24. Tangney P., Scandolo S. An ab initio parametrized interatomic force field for silica // The J. Chemical Physics. 2002. V. 117. № 19. P. 8898–8904.

  25. Garofalini S.H. Molecular dynamics simulations of silicate glasses and glass surfaces // Reviews in Mineralogy and Geochemistry. 2001. V. 42. № 1. P. 131–168.

  26. Pedone A. Properties calculations of silica-based glasses by atomistic simulations techniques: a review // The J. Physical Chemistry C. 2009. V. 113. № 49. P. 20773–20784.

  27. von Alfthan S., Kuronen A., Kaski K. Realistic models of amorphous silica: a comparative study of different potentials // Physical Review B. 2003. V. 68. № 7. P. 073203.

  28. Wooten F., Winer K., Weaire D. Computer generation of structural models of amorphous Si and Ge // Physical Review Letters. 1985. V. 54. № 13. P. 1392.

  29. Van Duin A.C., Dasgupta S., Lorant F., Goddard W.A. ReaxFF: a reactive force field for hydrocarbons // The J. Physical Chemistry A. 2001. V. 105. № 41. P. 9396–9409.

  30. Wang C., Kuzuu N., Tamai Y. Molecular dynamics study on surface structure of a-SiO2 by charge equilibration method // J. Non-Crystalline Solids. 2003. V. 318. № 1–2. P. 131–141.

  31. Rappe A.K., Goddard III W.A. Charge equilibration for molecular dynamics simulations // The J. Physical Chemistry. 1991. V. 95. № 8. P. 3358–3363.

  32. Van Duin A.C., Strachan A., Stewman S., Zhang Q., Xu X., Goddard W.A. ReaxFFSiO reactive force field for silicon and silicon oxide systems // The J. Physical Chemistry A. 2003. V. 107. № 19. P. 3803–3811.

  33. Fogarty J.C., Aktulga H.M., Grama A.Y., Van Duin A.C., Pandit S.A. A reactive molecular dynamics simulation of the silica-water interface // The J. Chemical Physics. 2010. V. 132. № 17. P. 174 704.

  34. Rimsza J.M., Yeon J., Van Duin A.C.T., Du J. Water interactions with nanoporous silica: comparison of ReaxFF and ab initio based molecular dynamics simulations // The J. Physical Chemistry C. 2016. V. 120. № 43. P. 24803–24816.

  35. Yeon J., Van Duin A.C.T. ReaxFF molecular dynamics simulations of hydroxylation kinetics for amorphous and nano-silica structure, and its relations with atomic strain energy // The J. Physical Chemistry C. 2016. V. 120. № 1. P. 305–317.

  36. Rimsza J.M., Du J. Interfacial structure and evolution of the water–silica gel system by reactive force-field-based molecular dynamics simulations // The J. Physical Chemistry C. 2017. V. 121. № 21. P. 11534–11543.

  37. Musgraves J.D., Hu J., Calvez L. Springer handbook of glass // Cham: Springer; 2019. P. 326.

  38. Cahn R.W. Materials science: melting and the surface // Nature. 1986. V. 323. № 6090. P. 668–669.

  39. Wolf D., Yip S. MRS Bulletin. 1995. V. 20. Issue 1. P. 63.

  40. Nakano A., Kalia R.K., Vashishta P. First sharp diffraction peak and intermediate-range order in amorphous silica: finite-size effects in molecular dynamics simulations // J. Non-Crystalline Solids. 1994. V. 171. № 2. P. 157–163.

  41. Galeener F.L., Mikkelsen Jr J.C. Vibrational dynamics in O18-substituted vitreous SiO2 // Physical Review B. 1981. V. 23. № 10. P. 5527.

  42. Bin L., Jing-Yang W., Yan-Chun Z., Fang-Zhi L. Temperature dependence of elastic properties for amorphous SiO2 by molecular dynamics simulation // Chinese Physics Letters. 2008. V. 25. № 8. P. 2747.

  43. Matsui M. A transferable interatomic potential model for crystals and melts in the system CaO–MgO–Al2O3–SiO2 // Mineral. Mag. 1994. V. 58. P. 571–572.

  44. Sarnthein J., Pasquarello A., Car R. Structural and electronic properties of liquid and amorphous SiO2: An ab initio molecular dynamics study // Physical Review Letters. 1995. V. 74. № 23. P. 4682.

  45. Sarnthein J., Pasquarello A., Car R. Model of vitreous SiO2 generated by an ab initio molecular-dynamics quench from the melt // Physical Review B. 1995. V. 52. № 17. P. 12690.

  46. Spiekermann G., Steele-MacInnis M., Schmidt C., Jahn S. Vibrational mode frequencies of silica species in SiO2–H2O liquids and glasses from ab initio molecular dynamics // The J. Chemical Physics. 2012. V. 136. № 15. P. 154501.

  47. Spiekermann G., Steele-MacInnis M., Kowalski P.M., Schmidt C., Jahn S. Vibrational properties of silica species in MgO–SiO2 glasses obtained from ab initio molecular dynamics // Chemical Geology. 2013. V. 346. P. 22–33.

  48. Usui Y., Tsuchiya T. Ab initio two-phase molecular dynamics on the melting curve of SiO2 // J. Earth Science. 2010. V. 21. № 5. P. 801–810.

  49. Benoit M., Ispas S., Tuckerman M.E. Structural properties of molten silicates from ab initio molecular-dynamics simulations: Comparison between CaO–Al2O3−SiO2 and SiO2 // Physical Review B. 2001. V. 64. № 22. P. 224205.

  50. Litton D.A., Garofalini S.H. Vitreous silica bulk and surface self-diffusion analysis by molecular dynamics // J. Non-Crystalline Solids. 1997. V. 217. № 2–3. P. 250–263.

  51. Litton D.A., Garofalini S.H. Modeling of hydrophilic wafer bonding by molecular dynamics simulations // J. Applied Physics. 2001. V. 89. № 11. P. 6013–6023.

  52. Soules T.F. Molecular dynamic calculations of glass structure and diffusion in glass // J. Non-Crystalline Solids. 1982. V. 49. № 1–3. P. 29–52.

  53. Kubicki J.D., Lasaga A.C. Molecular dynamics simulations of SiO2 melt and glass; ionic and covalent models // American Mineralogist. 1988. V. 73. № 9–10. P. 941–955.

  54. Della Valle R.G., Andersen H.C. Molecular dynamics simulation of silica liquid and glass // The J. Chemical Physics. 1992. V. 97. № 4. P. 2682–2689.

  55. Horbach J., Kob W., Binder K. Molecular dynamics simulation of the dynamics of supercooled silica // Philosophical Magazine B. 1998. V. 77. № 2. P. 297–303.

  56. Horbach J., Kob W., Binder K. The dynamics of supercooled silica: acoustic modes and boson peak // J. Non-Crystalline Solids. 1998. V. 235. P. 320–324.

  57. Horbach J., Kob W., Binder K. Specific heat of amorphous silica within the harmonic approximation // The J. Physical Chemistry B. 1999. V. 103. № 20. P. 4104–4108.

  58. Horbach J., Kob W. Static and dynamic properties of a viscous silica melt // Physical Review B. 1999. V. 60. № 5. P. 3169.

  59. Binder K., Horbach J., Knoth H., Pfleiderer P. Computer simulation of molten silica and related glass forming fluids: recent progress // J. Physics: Condensed matter. 2007. V. 19. № 20. P. 205102.

  60. Vollmayr K., Kob W., Binder K. Cooling-rate effects in amorphous silica: A computer-simulation study // Physical Review B. 1996. V. 54. № 22. P. 15808.

  61. Gotze W., Sjogren L. Relaxation processes in supercooled liquids // Reports on Progress in Physics. 1992. V. 55. № 3. P. 241.

  62. Garca-Coln L.S., Del Castillo L.F., Goldstein P. Theoretical basis for the Vogel-Fulcher-Tammann equation // Physical Review B. 1989. V. 40. № 10. P. 7040.

  63. Quenneville J., Taylor R.S., Van Duin A.C.T. Reactive molecular dynamics studies of DMMP adsorption and reactivity on amorphous silica surfaces // The J. Physical Chemistry C. 2010. V. 114. № 44. P. 18894–18902.

  64. Tranh D.T.N., Van Hoang V. Molecular dynamics simulation of amorphous SiO2 thin films // The European Physical J. Applied Physics. 2015. V. 70. № 1. P. 10302.

  65. Athanasopoulos D.C., Garofalini S.H. Molecular dynamics simulations of the effect of adsorption on SiO2 surfaces // The J. Chemical Physics. 1992. V. 97. № 5. P. 3775–3780.

  66. Vo T., He B., Blum M., Damone A., Newell P. Molecular scale insight of pore morphology relation with mechanical properties of amorphous silica using ReaxFF // Computational Materials Science. 2020. V. 183. P. 109881.

  67. Pakarinen O.H., Djurabekova F., Nordlund K., Kluth P., Ridgway M.C. Molecular dynamics simulations of the structure of latent tracks in quartz and amorphous SiO2 //Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms. 2009. V. 267. № 8–9. P. 1456–1459.

  68. Takada A., Bell R.G., Catlow C.R.A. Molecular dynamics study of liquid silica under high pressure // J. Non-Crystalline Solids. 2016. V. 451. P. 124–130.

  69. Le V.V., Nguyen G.T. Molecular dynamics simulation of structural transformation in SiO2 glass under densification // J. Non-Crystalline Solids. 2019. V. 505. P. 225–233.

  70. Badro J., Barrat J.L., Gillet P. Numerical simulation of α-quartz under nonhydrostatic compression: memory glass and five-coordinated crystalline phases // Physical Review Letters. 1996. V. 76. № 5. P. 772.

  71. Szymanski M.A., Shluger A.L., Stoneham A.M. Role of disorder in incorporation energies of oxygen atoms in amorphous silica // Physical Review B. 2001. V. 63. № 22. P. 224207.

  72. Mukhopadhyay S., Sushko P.V., Stoneham A.M., Shluger A.L. Modeling of the structure and properties of oxygen vacancies in amorphous silica // Physical Review B. 2004. V. 70. № 19. P. 195203.

  73. El-Sayed A.M., Watkins M.B., Afanas’ev V.V., Shluger A.L. Nature of intrinsic and extrinsic electron trapping in SiO2 // Physical Review B. 2014. V. 89. № 12. P. 125201.

  74. https://lammps.sandia.gov

  75. Newsome D.A., Sengupta D., Foroutan H., Russo M.F., Van Duin A.C. Oxidation of silicon carbide by O2 and H2O: a ReaxFF reactive molecular dynamics study, Part I // The J. Physical Chemistry C. 2012. V. 116. № 30. P. 16111–16121.

  76. Yu Y., Wang B., Wang M., Sant G., Bauchy M. Revisiting silica with ReaxFF: towards improved predictions of glass structure and properties via reactive molecular dynamics // J. Non-Crystalline Solids. 2016. V. 443. P. 148–154.

  77. Yeon J., Van Duin A.C.T. ReaxFF molecular dynamics simulations of hydroxylation kinetics for amorphous and nano-silica structure, and its relations with atomic strain energy // The J. Physical Chemistry C. 2016. V. 120. № 1. P. 305–317.

  78. Doremus R.H. Viscosity of silica // J. Applied Physics. 2002. V. 92. № 12. P. 7619–7629.

  79. Мазурин О.В. Стеклование. Наука, Ленинград; 1986, 158 с.

  80. Johnson J.R., Bristow R.H., Blau H.H. Diffusion of ions in some simple glasses // J. American Ceramic Society. 1951. V. 34. № 6. P. 165–172.

  81. Roma G., Limoge Y., Martin-Samos L. Oxygen and silicon self-diffusion in quartz and silica: the contribution of first principles calculations // Defect and Diffusion Forum. Trans Tech Publications Ltd, 2006. V. 258. P. 542–553.

  82. Mikkelsen Jr J.C. Self-diffusivity of network oxygen in vitreous SiO2 // Applied Physics Letters. 1984. V. 45. № 11. P. 1187–1189.

  83. Williams E.L. Diffusion of oxygen in fused silica // J. American Ceramic Society. 1965. V. 48. № 4. P. 190–194.

  84. Kalen J.D., Boyce R.S., Cawley J.D. Oxygen tracer diffusion in vitreous silica // J. American Ceramic Society. 1991. V. 74. № 1. P. 203–209.

  85. Rodríguez-Viejo J., Sibieude F., Clavaguera-Mora M.T., Monty C. 18O diffusion through amorphous SiO2 and cristobalite // Applied Physics Letters. 1993. V. 63. № 14. P. 1906–1908.

  86. Sucov E.W. Diffusion of oxygen in vitreous silica // J. American Ceramic Society. 1963. V. 46. № 1. P. 14–20.

  87. Richet P., Bottinga Y., Denielou L., Petitet J.P., Tequi C. Thermodynamic properties of quartz, cristobalite and amorphous SiO2: drop calorimetry measurements between 1000 and 1800 K and a review from 0 to 2000 K // Geochimica et Cosmochimica Acta. 1982. V. 46. № 12. P. 2639–2658.

  88. Doremus R.H. Viscosity of silica // J. Applied Physics. 2002. V. 92. № 12. P. 7619–7629.

  89. Takahashi T., Fukatsu S., Itoh K.M., Uematsu M., Fujiwara A., Kageshima H., Takahashi Y., Shiraishi K. Self-diffusion of Si in thermally grown SiO2 under equilibrium conditious // J. Applied Physics. 2003. V. 93. 1 6. P. 3674–3676.

  90. King S.V. Ring configurations in a random network model of vitreous silica // Nature. 1967. V. 213. № 5081. P. 1112–1113.

  91. Skuja L. Optically active oxygen-deficiency-related centers in amorphous silicon dioxide // J. Non-Crystalline Solids. 1998. V. 239. № 1–3. P. 16–48.

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