CC BY
10УДК 536.421
С. О. Ю р ч е н к о, Н. П. Крючков
НЕУПОРЯДОЧЕННЫЕ СОСТОЯНИЯ И ФУНКЦИИ ПАРНОГО РАСПРЕДЕЛЕНИЯ РАССТОЯНИЙ В ЛЕННАРД-ДЖОНСОВСКОЙ СИСТЕМЕ
Рассмотрены процессы плавления в гранецентрированной кубической и гексагональной плотноупакованной решетках. Предполагается, что взаимодействия между частицами описываются потенциалом Леннарда-Джонса. Вычислительные эксперименты проведены методом молекулярной динамики. Показано, что функции парного распределения могут быть построены при использовании приближенного метода s-функций. Установлено, что в точке плавления параметр неупорядоченности изменяется скачкообразно, в то время как температурная зависимость эффективного параметра решетки определяется видом решетки.
E-mail: [email protected]
Ключевые слова: конденсированное состояние, неупорядоченные структуры, фазовый переход, молекулярная динамика.
В современных инженерно-технических задачах большую роль играют процессы, протекающие в конденсированных средах, в частности в жидкостях. Для эффективного применения этих процессов требуется понимание свойств конденсированных состояний, однако в настоящее время пока не существует законченной теории описания неупорядоченных структур, жидкостей и фазовых переходов. Важная роль в описании свойств различных структур отводится функции радиального распределения частиц.
Построение функции парного распределения (ФПР) расстояний между частицами часто выполняют либо для анализа структур [1], либо для проверки пригодности новых моделей молекулярных систем [2]. Существует ряд работ, посвященных теоретическому построению функций парного распределения. В частности, в [3] была предпринята попытка построения ФПР для модели твердых сфер. Однако общего подхода к данной проблеме сегодня не существует.
В работе [4] предложен новый способ описания неупорядоченных структур, в котором ФПР играет одну из ключевых ролей. Суть нового метода состоит в использовании функции плотности вероятности для одного из ближайших узлов кристалла — s(r)-функции.
В данной работе ФПР системы построены для гранецентрированной кубической (ГЦК) и гексагональной плотноупакованной (ГПУ) решеток. Проведены сравнения результатов построения с ФПР, которые получены методами молекулярной динамики (МД) для систем
Леннарда - Джонса (Ы), находящихся в тех же кристаллических решетках.
Детали вычислительного эксперимента. Моделирование методами МД выполнено с помощью пакета ЬАММРБ [5]. В качестве модели была выбрана 3Б-система частиц (~ 25 • 103), парное взаимодействие между которыми описывается потенциалом Леннарда - Джонса
12 / , \ 6 4
и (г) = 4е
Вычислительные эксперименты выполнены в системе единиц е = 1, а = 1. Моделирование движения атомов проведено в ГЦК-и ГПУ-решетках, которые в ходе вычислительных экспериментов имели плотность р =1 и температуру Т = 0,1... 2,0 с шагом в 0,05 единиц Леннарда - Джонса.
Для идентификации момента плавления системы на каждом шаге моделирования рассчитывали параметр
а2 =
Е (Гк (*) - Гк (0))2
к
представляющий собой квадрат смещения частиц от их исходного положения, усредненный за время моделирования. Также выполнено построение ФПР расстояний между частицами системы р = р (г,Т).
Для построения ФПР расстояний методом з-функций, описанным в работе [4], необходимо выбрать некоторую аппроксимацию реальной з-функции. В данной работе з-функция взята в виде
/ |г — а|2 \ §§ ^ п, а) = ехр I--1 — 0,1
, §(г,п,а) §(г,п,а) > 0,
з (г, п, а) = < _
10 в (г, п, а) < 0,
где значение 0,1 в выражении для § введено для локализации области пребывания ближайшего узла. Построения ФПР расстояний между частицами методом з-функций (строили функцию р3 = р3 (г, а, п)) выполняли для таких наборов параметров (а, а), что а = 1,0... 1,2 с шагом 10-3 и а = 0,01... 0,21 с шагом 2 • 10-4. Функции р8 и р нормировали на единицу
J р3 (г,А,п)>Лг = J р(г,Т )д,г = 1,
где интегрирование проводили на отрезке (0, 5). Верхняя граница интегрирования обусловлена тем, что в результате вычислительных экс-
периментов как p, так и ps-функции строили в интервале 0 < r < 5. На больших расстояниях потенциал Леннарда - Джонса быстро убывает, поэтому взаимодействия со столь удаленными узлами не приводят к значимым физическим эффектам.
Сравнение функции ps (r, T, п) (фактически это теоретическая кривая) c p(r, T) (функция, наблюдаемая в вычислительном эксперименте) выполняли с помощью минимизации функционала ошибки
Ф = У (ps(r,a,n) — p(r,T))2 dr ^ min.
Здесь минимизация достигается посредством варьирования параметров а, п, откуда следуют температурные зависимости эффективного параметра решетки а = а (T) и параметра неупорядоченности
П = П (T)•
Результаты вычислительных экспериментов. На рис. 1 представлены зависимости o2(T) для ГЦК- и ГПУ-конфигураций. В момент плавления значение о2 начинает резко расти, что свидетельствует о резком увеличении подвижности узлов нагретой решетки. Последний эффект является характерным признаком жидкого состояния. Согласно проведенным расчетам, ГЦК-решетка плавится при T ~ 1,6, а ГПУ-решетка — при T ~ 1,55.
На рис.2 представлены температурные зависимости a(T) и п2(T) для ГЦК- и ГПУ-решеток. Видно, что у ГЦК-решетки значение параметра а в твердой фазе постоянно, а при плавлении скачкообразно изменяется. Зависимость n2(T) близка к линейной и может быть представлена как
n2(T) = 9,9 • 10-3T
практически до момента плавления, а при плавлении отмечается скачок.
Скачкообразное изменение параметра а может быть обусловлено тем, что при плавлении становятся неразличимыми второй и третий пики на функции p(r).
В отличие от случая ГЦК-решетки параметр а ГПУ-решетки изменяется линейно и не претерпевает скачков в момент плавления. Поведение зависимости n2(T) аналогично случаю ГЦК-решетки и в твердой фазе может быть представлено как
n2(T) = 9,5 • 10-3T.
В данном случае ввиду непрерывности параметра а говорить о принципиальном несовпадении правил суммирования для ГПУ-решетки и жидкости уже нельзя. Тогда существует два способа объяснения скачка зависимости п2 (T):
0^5 1,0 1,5 2,0
б
Рис. 1. Зависимости а2 от Т для ГЦК- (а) и ГПУ-решеток (б)
1) структура жидкости (по которой ведется суммирование) совпадает с ГПУ-решеткой и скачок зависимости принципиально неизбежен;
2) структура жидкости незначительно отличается от ГПУ-решетки, что приводит к расхождениям р8 и р, которые могут быть уменьшены за счет резкого увеличения значения параметра п. Однако незначительное изменение правил суммирования при плавлении может позволить избежать возникновения скачка п2 при плавлении.
Отметим, что в аналогичных расчетах для двумерной гексагональной решетки каких-либо заметных скачков параметра п при плавлении не наблюдается.
0,5 * ' 1,0 ' ' 1,5 ' ' 2,0
б
Рис. 2. Зависимости параметров а и п2 от Т для ГЦК- (*) и ГПУ-решеток (•)
На рис. 3 приведены функции р3 (г, а (Т), п (Т)) и р (г, Т) для ГЦК-решетки. Теоретические и экспериментальные кривые практически совпадают вплоть до момента плавления. Однако видно несколько незначительных различий:
• при низких температурах дальние пики теоретической кривой "размываются" сильнее экспериментальных, но это не очень существенно, так как наибольшее значение имеет совпадение ближних пиков (основной вклад в свободную энергию структуры дают взаимодействия с ближайшими узлами);
• при высоких температурах наблюдаются некоторые расхождения в области между первым и вторым пиками. Скорее всего они обусловлены неточностью выбранного приближения з-функции (усеченная гауссоида) при высоких температурах.
На рис. 4 для сравнения представлены функции р3 (г, а (Т), п (Т)) и р (г, Т) для ГПУ-решетки. Как и в случае ГЦК-решетки, отчетли-
Рис. 3. Функция парного распределения для ГЦК-решетки, построенные методом й-функции (сплошные линии — рв) и методом МД (пунктирные линии — р):
а — Т = 0,1, А = 1,12, а = 0,0326; б — Т = 0,5, А = 1,12, а = 0,0706; в — Т = 0,9, А = 1,12, а = 0,0942; г — Т = 1,3, А = 1,12, а = 0,1114; д — Т = 1,6, А = 1,11, а = 0,1338; е — Т = 1,65, А = 1,09, а = 0,1564
во видно хорошее совпадение ФПР расстояний, полученных разными методами. Характер различий аналогичен тем, которые описаны выше для ГЦК-решетки.
В отличие от ГЦК-решетки выше точки плавления различия между кривыми ФПР расстояний для ГПУ-решетки существенно меньше. Последнее говорит о том, что структурно жидкое состояние ближе к ГПУ-решетке с делокализованными узлами.
Возвращаясь к вопросу о необходимости смены правил суммирования при плавлении системы, можно отметить некоторые характерные особенности изменения ФПР расстояний узлов ГПУ-решетки при плавлении. На рис. 5 видно, что в момент плавления ФПР расстояний системы претерпевает скачок. На первом пике, однако, существенных изменений не наблюдается.
Рис. 4. Функция парного распределения для ГЦК-решетки, построенные методом й-функции (сплошные линии — рв) и методом МД (пунктирные линии — р):
а — Т = 0,2, А = 1,093, а = 0,04800; б — Т = 0,70, А = 1,092, а = 0,0814; в — Т = 1,1, А = 1,092, а = 0,102; г — Т = 1,55, А = 1,092, а = 0,122; д — Т = 1,6, А = 1,09, а = 0,158; е — Т = 1,9, А = 1,09, а = 0,162
При построении ФПР расстояний методом з-функций обнаружено, что форма первого пика практически полностью определяется видом з-функции и слабо чувствительна к малым изменениям правил суммирования (т.е. слабо зависит от конфигурации узлов решетки). В то же время вид остальных пиков крайне чувствителен к изменениям правил суммирования. Таким образом, можно предположить, что скачок параметра п при плавлении во многом обусловлен необходимостью смены правил суммирования как в ГЦК-, так и в ГПУ-решетках.
Проведенное построение ФПР расстояний между частицами системы методом з-функций и с помощью метода МД для ГЦК- и ГПУ-решеток позволяет сделать несколько важных выводов.
В области твердой фазы наблюдается очень хорошее совпадение ФПР расстояний, построенных двумя разными методами, что подтверждает возможность использования метода з-функций. Незначительные расхождения могут быть обусловлены выбором приближения
р
0,6 0,5 0,4 0,3 0,2 ОД
1
2
3
4
г
Рис. 5. Экспериментально найденное изменение ФПР системы ниже (отлошная линия) и выше (пунктирная линия) точки плавления при изменении температуры на AT = 0,05. Метод МД, ГПУ-решетка
s-функции. После плавления ГЦК-решетки наблюдаются существенные расхождения с результатами метода МД при сравнении с плавлением ГПУ-решетки. Последнее говорит о том, что характерная структура жидкости ближе к ГПУ-решетке, чем к ГЦК-решетке. Зависимости a(T) и n2(T) в области твердой фазы для обеих решеток линейны. В случае ГЦК-решетки обе зависимости претерпевают скачок в момент плавления. В ГПУ-решетках отмечаются скачкообразные изменения только параметра n2 (T).
Работа выполнена при поддержке РФФИ (проекты № 12-08-31104 мол_а, 12-08-33112 мол_а_вед).
СПИСОК ЛИТЕРАТУРЫ
1. Yiannourakou M., Economou I. G., Bitsanis I. A. Structural and dynamical analysis of monodisperse and polydisperse colloidal systems // Journal of Chemical Physics. - 2010. - Vol. 133, no. 22. - P. 224901.
2. Modeling noncrystalline materials: Use of vibrational spectra as a protocol for validation / F. Gaspari, I. M. Kupchak, A. I. Shkrebtii, and J. M. Perz // Phys. Rev. B. - 2009. - Jun. - Vol. 79. - P. 224203.
3.Cochran T. W., Chiew Y. C. Radial distribution function of freely jointed hard-sphere chains in the solid phase // Journal of Chemical Physics. - 2006. -Vol. 124, no. 7. - P. 074901.
4. S t r u c t u r e of the nanobubble clusters of dissolved air in liquid media / N.V. Su-yazov, V. Kozlov, N.F. Bunkin, S.O. Yurchenko // Journal of Biological Physics. -2012. - Vol. 38, no. 1 - P. 121-152.
5. http://lammps.sandia.gov/.
Статья поступила в редакцию 05.07.2012
