CC BY
13
Труды БГТУ, 2016, № 6, с. 71-75
71
УДК 531.19
Г. С. Бокун
Белорусский государственный технологический университет
АЛГОРИТМИЗАЦИЯ ДВУХУРОВНЕВОГО МОЛЕКУЛЯРНО-СТАТИСТИЧЕСКОГО ПОДХОДА ДЛЯ РАСЧЕТА ПАРАМЕТРОВ КРИСТАЛЛИЧЕСКИХ НАНОЧАСТИЦ СО СТРУКТУРОЙ АЛМАЗНОЙ РЕШЕТКИ
Для изучения структуры и термодинамических характеристик наночастиц используются общие статистико-механические формулы, полученные в рамках двухуровневого молекулярно-статистического подхода, который является модификацией метода условных распределений и используется при описании свойств неоднородных молекулярных систем с помощью потенциалов средних сил. Эти потенциалы являются функционалами от искомых полей плотности. С помощью одночастичных потенциалов рассчитываются младшие условные функции распределения. По последним определяются среднеквадратичные отклонения, описывающие распределение частиц возле узлов решетки. Найден способ приближенного решения полученной ранее сложной системы интегральных уравнений, требующей выполнения большого объема численных расчетов. В результате получены аналитические выражения для одночастичных потенциалов, которые использованы для исследования изменения степени локализации частиц в узлах решетки вдоль радиальных направлений в сферических наночастицах. Численные расчеты проведены для кристаллических наночастиц со структурой алмаза. Для этого осуществлено компьютерное построение алмазной решетки. Применение общей статистической теории неоднородных систем позволило разработать алгоритмы анализа взаимного расположения частиц в локальном и дальнем окружениях. Исследована структура углеродных нанообразований с алмазоподобной решеткой.
Ключевые слова: коррелятивные функции, потенциалы средних сил, неоднородная система, химический потенциал, функционал свободной энергии, наночастица.
G. S. Bokun
Belarusian State Technological University
THE ALGORITHMIZATION OF THE TWO-LEVEL MOLECULAR-STATISTICAL APPROACH FOR CALCULATION OF PARAMETERS OF CRYSTAL NANOPARTICLES WITH DIAMOND LATTICE
To describe the structure and thermodynamic characteristics of nanoparticles the general statistical-mechanical expressions are used that are obtained in the framework of the two-level molecular-statistical approach, which is a modification of the method of conditional distributions and is used to describe properties of non-uniform molecular systems using the potentials of average forces. The single-particle potentials are functionals of the density fields. The potentials are used for calculation of younger conditional distribution functions. The latter are used for calculating the parameters describing the particle distribution near the lattice sites. The way to approximately solve the system of integral equations is found for that, which requires the implementation of a large amount of numerical calculations. As a result, the analytical expressions for the single-particle average potentials used to determine changes in the degree of localization of particles around lattice sites along the radial directions of spherical nanoparticles are found. The numerical calculations for crystalline nanoparticles with the diamond structure are performed. To do this, a computer construction of the lattice with a diamond-like structure is carried out. The application of the general statistical theory of inhomogeneous systems allowed the development of algorithms for the analysis of the mutual arrangement of the particles in the local and distant environments. The structure of the carbon nanostructures with a diamond lattice is investigated.
Key words: correlation functions, approximation of the potential medium-strength, homogeneous system, chemical potential, free energy functional, limited system, nanoparticle.
Введение. Наряду со все более широким применением численных методов Монте-Карло и молекулярной динамики для моделирования и прогнозирования различных свойств вещества сохраняют значение подходы, основанные на сокращенном статистическом описании, к которым, в частности, относятся методы коррелятивных функций [1, 2]. Если в первом случае
коррелятивные функции находятся при усреднении результатов вычислений, то здесь они играют роль исходных объектов, для которых записывается система соответствующих уравнений. Как правило, это бесконечная цепочка интегро-дифференциальных уравнений, решение которой возможно при внесении приближений, позволяющих замкнуть цепочку на
первых уравнениях. Первоначально для замыкания применялись соотношения между младшими и старшими функциями распределения. В последующем цепочки записывались и аппроксимировались для некоторых характеристик этих функций. Конструктивной оказалась аппроксимация ячеечных потенциалов средних сил. На использовании этих одночастичных потенциалов основана настоящая работа, так как эта аппроксимация передает основные особенности макроскопической системы в конденсированном состоянии вещества. Идея дробления потенциалов, предопределившая успехи в рассмотрении однородных систем, была перенесена на описание неоднородных систем путем дополнительного дробления пространства конфигураций на подпространства. На каждом из подпространств вводятся наборы парциальных функций распределения, так что статистико-механи-ческое описание приобретает двухуровневый характер [2]. Каждое отдельное подпространство характеризуется соответствующим полем плотности числа частиц во всем объеме системы.
Было показано, что суммарные равновесные характеристики системы в целом (при пренебрежении флуктуациями) соответствуют такому полю плотности, при котором функционал свободной энергии минимален [2, 3]. Очевидно, что рассматриваемый двухуровневый метод наиболее удобен для описания систем, содержащих малое число частиц (кластеров). При наличии центра симметрии структурные и термодинамические характеристики таких систем (например, наночастиц) зависят от расстояния до их центра.
Основная часть. В развиваемой теории для описания структуры нанообразований применяются три параметра: средняя заселенность узлов решетки, смещение узлов в связи с ее деформацией (пространственной релаксацией [4]) и среднеквадратичное отклонение молекул относительно узлов. В данной работе исследуется изменение величины среднеквадратичного отклонения молекул по мере удаления от центра нано-частицы. В области кристаллического состояния функция распределения имеет резко выраженный пик. Для описания распределения частиц в этом случае достаточно ограничиться одночастичны-ми ячеечными потенциалами. Рассмотрим процедуру их вычисления, моделируя неоднородность в системе с помощью микроячеек разных объемов юг-, удовлетворяющих условию
N
IИ = V,
(1)
'=1
где V - полный объем системы; ' = 1, 2, ..., N N - общее количество частиц в системе.
Рассмотрим одночастичное распределение частиц (атомов или молекул), отвечающее пер-
вому приближению метода условных распределений, которое описывается унарной ненормированной функцией распределения (в = 1 / кТ):
( N \
рп{я') = ехр -вIФ; (я')
^ ](')
(2)
Для вычисления одночастичных потенциалов ф;(я') используем вариационный принцип Н. Н. Боголюбова, согласно которому
QN > О0ехр(-в ^ (Н-Н0) у0).
(3)
Угловые скобки в (3) означают усреднение с помощью базисной функции распределения (2); QN, Q0 - конфигурационные интегралы исследуемой и базисной систем соответственно.
В соотношении (3) энергия взаимодействия молекул между собой определяется формулой
N N
Н = 2i1ф (я'
'), ' * 1.
(4)
' =1 1 =1
Энергия базисной системы представляется средними потенциалами в виде
Н 0 = 11ф 1 (Я),
'=1 1=1
' * 1.
(5)
Конфигурационный интеграл базисной системы записывается в факторизованной форме:
Q0N =п Qi,
'=1
( N \
Q' = | ехр -вI ф; (Я) йЯ.
(6)
1 (')
Варьирование (3) позволяет найти одночас-тичные потенциалы, которые определяют наилучшее приближение к конфигурационному интегралу QN. В результате получаем интегральное уравнение для расчета этих потенциалов:
ф 1 (я)=ОI ф( Я1 - Я' Ж1 (Я1 )йЯ. (7)
Преобразуем это уравнение в соотношение, связывающее между собой среднюю энергию системы с ее свободной энергией. Для этого умножим уравнение (7) на нормированную од-ночастичную функцию распределения (2) и, выполнив суммирование по 1, получим:
1 N
оiф1 ( я' ш я' ) =
^ 1 (-)
1 (- ) 1 N 1
=ОIОт I ф(Я' - Я1 ЖЛЯ')Ри(Я1)йя1. (8)
О 1 (') и,
Из соотношения (8) находим, что внутренняя энергия системы определится выражением
N
<H >0=-1 £ 1"П а2-
(9)
Из соотношения (9) следует выражение для конфигурационного интеграла в форме
QN (а г
(10)
г =1
Нормированные функции в ^ ^приближении метода условных распределений имеют резкие пики в окрестности узлов, поэтому заменим их на функции ^ с равномерным распределением молекул внутри сфер с радиусами Ьг, центры которых совпадают с узлами решетки-Параметры Ьг подбираем так, чтобы среднеквадратичные отклонения сг от узлов решетки были равны для функций и Е11:
с, = | г2(х,у,г)с1щ = Л Ь. (11)
В этом случае для среднего потенциала (7) получается аналитическое выражение вида [5]:
ф*(р)= 4
1
+-
7,2Ь2
14,4Ь4
(р2 -Ь2)6 (р2 -Ь2)7 (р2 -Ь2)8
+
,(12)
где р - расстояние от молекулы в ячейке юг до узла ячейки ф и р обезразмерены с помощью параметров используемого здесь потенциала Леннард-Джонса.
Задача, решаемая в работе, состояла в том, чтобы найти зависимость параметра Ь от расстояния до центра наночастицы. Для этого необходимо выполнить классификацию узлов используемой решетки. Рассмотрим случай решетки алмаза, которая состоит из двух гране-центрированных подрешеток, одна из которых смещена относительно основной подрешетки на четверть пространственной диагонали ее элементарной ячейки. Базис образовывается восьмью атомами с координатами (0, 0, 0; 0, 0,5, 0,5; 0,5, 0,5, 0; 0,5, 0, 0,5) для четырех атомов основной подрешетки и с координатами (0,25, 0,25, 0,25; 0,75, 0,25, 0,75; 0,25, 0,75, 0,75; 0,75, 0,75, 0,75) для атомов второй смещенной подрешетки. Молекулярный объем V для алмазной решетки определяется формулой V = а3 / 8, где а - параметр основной решетки, равный длине ребра куба элементарной ячейки. Базисными векторами в каждой подрешетке являются векторы а1, а2, а3 и Ь1, Ь2, Ь3 со следующими одинаковыми координатами (0,5, 0,5, 0,5), (0,5, 0, 0,5), (0, 0,5, 0,5). По координатам базисных векторов построены массивы координат узлов
первой и второй подрешеток в основной (глобальной) системе координат:
хп = 0,5(т + «2); хт = хп + 0,25, уп = 0,5( тз + т2); Ут = Уп + 0,25, (13)
гп = 0,5(т1 + т3); гт = гп + 0,25.
В выписанных здесь формулах т1, т2, т3 - целочисленные переменные, изменяющиеся в интервале [-М, М], определяющем размеры наночастицы. Тогда общее число частиц в первой и второй подрешетках составляет величину N = = 2(2М + 1)3. Индекс п изменяется в интервале 1 < п < (2М +1)3, следовательно, индекс т в (13) принимает значения в интервале (2М +1)3 < т < < 2(2М +1)3. Расстояния от узлов первой и второй подрешеток до начала координат основной системы находились по формуле
=4
х2 + у 2 + г2.
(14)
Для рассчитанных по формулам (13) и (14) координат и расстояний сформирована матрица А. Номера столбцов этой матрицы совпадают с номерами узлов. В первой, второй, третьей и четвертой строках матрицы записаны расстояния и координаты узлов решетки соответственно.
При выполнении суммирования в (2) удобно осуществлять сначала суммирование по узлам, принадлежащим координационной сфере с номером I, относительно центра ячейки юг, а затем по I (0 < I < Ь, здесь Ь - число учитываемых в расчетах ближайших координационных сфер). Для этого массив А был переформатирован в новый массив В, где узлы решетки перенумерованы так, чтобы с возрастанием их номера расстояние до центра наночастицы увеличивалось при переходе от одной координационной сферы к другой либо сохранялось в пределах координационной сферы.
Для дальнейшей классификации узлов решетки из массива В был сформирован новый массив С из трех строк и с числом столбцов, равных числу координационных сфер. В первой строке этого массива указывается номер последнего узла из матрицы В, принадлежащего выбранной координационной сфере. Во второй строке записывается количество узлов, принадлежащих этой сфере, а в третьей - расстояния до этих узлов.
В связи с тем, что искомый параметр Ь зависит от радиуса координационной сферы относительно центра наночастицы, необходимо установить, к какой координационной сфере принадлежит выбранный узел решетки. Указанную связь находили с помощью специального одномерного массива В, в котором номера узлов записаны в первом столбце, а номера соответствующих координационных сфер - во втором.
Для расчета коррелятивной функции я' ) в окрестности узла с номером ' необходимо выполнить суммирование по всем узлам 1 в локальном окружении узла '. Параметры Ь этих узлов зависят от их положения в основной (глобальной) системе координат. Для их определения по координатам узлов 1, окружающих избранный узел ', устанавливались их номера в матрице В и соответствующие им номера координационных сфер относительно центра нано-частицы. По номерам сфер определялись параметры Ь и вычислялись потенциалы средних сил, входящие в определение функции распределения (2). Найденные таким образом функции распределения использовались для вычисления среднеквадратичных отклонений и новых значений параметров Ь (радиусов сфер), моделирующих сферически симметричное распределение частицы в окрестности своего узла. В связи со сказанным представляется возможным интегрирование по молекулярному объему
V = а3 / 8 заменить интегрированием по объему шара, радиус которого вычисляется из условия равенства этих объемов, т. е.
4 пт? = а8Т ^ Т = 0,308. (15)
Все расчеты выполнялись в единицах параметра решетки а, соответственно, безразмерное расстояние между ближайшими соседними узлами составляло 0,432. Радиус тч, определенный по (15), оказался больше радиуса сферы, вписанной в ячейку Вигнера - Зейтца, который равен 0,216.
Интегрирование выполнялось по 10 лучам, направленным на четыре ближайшие узлы, и вдоль трех координатных осей в положительном и отрицательном направлениях. Поскольку расчеты проводились в единицах параметра решетки а, то входящие в потенциал безразмерные расстояния пересчитывались путем умножения на параметр
^ = 2,59 ^ ) , (16)
где vЬ - молекулярный объем, при котором расстояние между ближайшими узлами решетки соответствует минимуму потенциала Леннард-Джонса с параметрами, определенными на основе характеристик монокристалла алмаза
V = 5,1 А3).
Парный межатомный потенциал Леннард-Джонса учитывает зависимость энергии взаимодействия атомов углерода только от расстояния между ними. Поскольку атомы углерода, расположенные в ближайших ячейках, образуют прочную ковалентную связь, то необходимо учесть зависимость их энергии от всех валент-
ных углов в тетраэдрах, образующих сферическую наночастицу. Это взаимодействие описывалось трехчастичными добавками к парному потенциалу с помощью формулы А и = с(Аф) . Коэффициент угловой жесткости связей с в тетраэдрах принимался в соответствии с данными работы [6]. При этом энергия угловых деформаций в тетраэдрах ограничивалась максимальным значением, которое достигалось при смещении атома из узла на величину, равную ширине ямы потенциала Леннард-Джонса.
В результате расчетов (с учетом вкладов только от взаимодействия с ближайшими соседями) установлено, как изменяется микроструктура анализируемых наночастиц с ростом их размеров, и проведено сопоставление с аналогичными данными для случая неограниченного кристалла.
Полученные результаты свидетельствуют о высокой степени локализации колебаний атомов возле узлов в макрокристалле до безразмерной температуры 9, равной 0,1. При росте температуры вплоть до этого значения параметр Ь медленно увеличивается до 0,058. Как только температура превышает указанное значение, параметр Ь испытывает резкий скачок, связанный с делокализацией функции распределения.
Это можно интерпретировать как фазовый переход, вызванный плавлением кристалла. Рассчитанная таким образом температура близка к экспериментальному значению температуры плавления алмаза (9 ~ 0,05).
Полученные характеристики микроструктуры для сферической наночастицы изменяются при удалении от ее центра и зависят от общего количества атомов в наночастице. При температурах, близких к температуре плавления, на-ночастицы сохраняют устойчивость при высокой степени делокализации поверхностных атомов за счет удержания их атомами, находящимися во внутренней области наночастицы.
Дальнейшие расчеты, с учетом вкладов во взаимодействие от вторых и последующих соседей, показали, что алмазоподобная кристаллическая структура не может существовать, если взаимодействие со вторыми и последующими соседями описывается тем же парным потенциалом, что и для первых соседей. Это следует из того, что при использовании одного и того же парного потенциала Леннард-Джонса локализованных решений не существует, по крайней мере при температурах выше 0,004, для которых удалось провести вычисления по применяемой в этой работе программе. При описанном выше взаимодействии между атомами реализуется кристалл с гранецентрированной решеткой.
Алмазоподобная структура наночастицы может существовать при образовании прочной
(ковалентной) связи с первыми соседями при условии, что энергия связи со вторыми и последующими соседями будет, по крайней мере, на четыре порядка меньше, чем с первыми. Это означает, что взаимодействие со вторыми и последующими соседями носит лондоновский (ван-дер-ваальсовый) характер. Оценки, проведенные с использованием справочных данных по энергиям связи атомов углерода в молекулах и кристаллах, показывают, что для первых соседей глубина ямы должна соответствовать температуре порядка 60 000 К, а для остальных -60 К, т. е. даже меньше, чем для случая взаимодействия между атомами аргона (110 К).
Заключение. Выполненная в работе алгоритмизация двухуровневого статистического описания наноструктур позволяет исследовать их свойства с помощью коррелятивных функций при различных термодинамических условиях.
Исследования проведены в соответствии с «Координационным планом работ, выполняемых в Объединенном институте ядерных исследований (г. Дубна, РФ) с участием организаций и учреждений Республики Беларусь в 2015 г.». Автор благодарен главному научному сотруднику лаборатории теоретической физики (ОИЯИ) В. Б. Приезжеву за рекомендации, высказанные при обсуждении проведенных исследований.
Литература
1. Rott L. A., Vikhrenko V. S. Statistical method of conditional distributions // Fortschr. Phys. 1975. Vol. 23, no. 3. P. 133-164.
2. Narkevich I. I. Statistical theory of nonuniform systems and reduced description in the density fluctuation theory // Physica. 1982. Vol. 112A. P. 167-192.
3. Бокун Г. С., Вихренко В. С., Наркевич И. И. Применение вариационных методов для описания структурных и термодинамических характеристик наночастиц // Автоматический контроль и автоматизация производственных процессов: материалы Междунар. науч. конф., 22-24 окт. 2015 г. Минск, 2015. С. 239-243.
4. Narkevich I. I. A statistical study of the defect crystal lattice relaxation // Physica A. 1988. Vol. 150. P.659-671.
5. Наркевич И. И., Немцов В. Б., Ротт Л. А. Вычисление на ЭВМ коррелятивных функций молекулярного кристалла // Известия высших учебных заведений. Физика. 1973. № 4. С. 95-100.
6. Schröder C., Schwarzer D., Vikhrenko V. Molecular Dynamics Simulation of Heat Conduction through a Molecular Chain // Journ. Phys. Chem. A. 2009. Vol. 113, no. 51. P. 14039-14051.
References
1. Rott L. A., Vikhrenko V. S. Statistical method of conditional distributions. Fortschr. Phys., 1975, vol. 23, no. 3, pp. 133-164.
2. Narkevich I. I. Statistical theory of nonuniform systems and reduced description in the density fluctuation theory. Physica, 1982, vol. 112A, pp. 167-192.
3. Bokun G. S., Vikhrenko V. S., Narkevich I. I. [Application of variational methods to describe the structural and thermodynamic properties of nanoparticles] Materialy mezhdunarodnoy nauchnoy konferentsii (Avtomaticheskiy kotrol' i avtomatizatsiya proizvodstvennykh protsessov) [Materials of the Internat. Sci. Conf. (Automatic control and automation of production processes)]. Minsk, 2015, pp. 239-243 (In Russian).
4. Narkevich I. I. A statistical study of the defect crystal lattice relaxation. Physica A, 1988, vol. 150, pp.659-671.
5. Narkevich I. I. Nemtsov V. B., Rott L. A. Computer calculations of correlation function of molecular crystals. Izvestiya vysshikh uchebnykh zavedeniy. Fizika [Izvsetiya of institution of higher education. Physics], 1973, no. 4, pp. 95-100 (In Russian).
6. Schröder C., Schwarzer D., Vikhrenko V. Molecular Dynamics Simulation of Heat Conduction through a Molecular Chain. Journ. Phys. Chem. A., 2009, vol. 113, no. 51, pp. 14039-14051.
Информация об авторе
Бокун Георгий Станиславович — кандидат физико-математических наук, доцент, доцент кафедры теоретической механики. Белорусский государственный технологический университет (220006, г. Минск, ул. Свердлова, 13а, Республика Беларусь). E-mail: gBokun@mail.ru
Information about the author
Bokun George Stanislavovich — PhD (Physics and Mathematics), Assistant Professor, Assistant Professor, the Department of Theoretical Mechanics. Belarusian State Technological University (13a, Sverdlova str., 220006, Minsk, Republic of Belarus). E-mail: gBokun@mail.ru
Поступила 07.03.2016
