CC BY
30
р
УДК 541.12.011.2/3:541.135.4 Вестник СПбГУ. Сер. 4, 2006, вып. 4
О. С. Терёхина, Е. Н. Бродская, Е. М. Пиотровская
КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ СВОЙСТВ И СТРУКТУРЫ КЛАСТЕРОВ ФТОРИДА БАРИЯ
Введение. Целью настоящей работы является молекулярное моделирование кластеров фторида бария, которые, с одной стороны, представляют самостоятельный интерес, а с другой - могут быть составными частями наноструктурных композиционных материалов, широко применяемых в настоящее время. Методы компьютерного моделирования дают подробную информацию о структуре исследуемой системы, сопоставление которой с данными спектральных экспериментов позволяет уточнить характер взаимодействия составляющих систему компонентов и механизмы происходящих в ней процессов. Однако нам известно весьма ограниченное число работ, посвященных изучению ионных кластеров методами компьютерного моделирования.
Наибольший интерес среди них представляют работы [1-3], в которых были исследованы фторидные кластеры ^Р, СаР2, ЬаРз. В работе [1] моделировалось плавление нанокристаллитов NaF и СаР2. Было показано, что в расплавах кластеров КаР наблюдается новая структура - стехиометрические петли типа ... .. Р~.... Ка+... Р-, тогда как в системах с СаР2 в расплаве такие структуры не образуются, и обсуждаются причины такого поведения систем. Разработан алгоритм расчета температуры плавления малых ионных кластеров, который и был использован в настоящей работе. Изучению кластеров ЬаР3 разных размеров (160, 552 и 3120 ионов) методом молекулярной динамики (МД) посвящены работы [2, 3]. В них моделирование проводилось при температурах от нескольких градусов до температуры испарения. Особое внимание уделялось исследованию механизма диффузии ионов, при этом рассматривались отличия в поведении ионов, лежащих на поверхности и внутри кластера. Было получено хорошее согласие рассчитанных температур плавления и коэффициентов диффузии с экспериментальными данными, что позволило сделать вывод о возможности применения простых модельных потенциалов для изучения малых кластеров со структурой типа тисонита.
В настоящей работе ставились задачи моделирования кластеров фторида бария разных размеров с использованием простейших потенциалов межионного взаимодействия, применяемых при описании макроскопических систем [4], и оценивалась пригодность этих потенциалов для расчета свойств малых систем. При использовании таких потенциальных моделей была оценена температура плавления, получены термодинамические свойства и структурные характеристики систем, определены коэффициенты диффузии ионов. Исследовалась зависимость всех перечисленных выше свойств от температуры и размера кластера, полученные результаты сравнивались с данными работы [4], посвященной моделированию макроскопической системы фторида бария.
Описание модели и методики численного эксперимента. Как и во многих более ранних исследованиях по моделированию флюоритоподобных твердых фаз, применялась простая модель «жестких ионов» с потенциалом Борна-Майера-Хаггинса. С учетом результатов [5-7] есть основания считать, что она вполне пригодна для моделирования данных систем. В ней кристалл рассматривается как система, образованная сферически симметричными ионами; потенциальная энергия взаимодействия ионов г и э описывается в виде
© О. С. Терёхина, Е. Н. Бродская, Е. М. Пиотровская, 2006 60
иф) = +Ац ехр(—) - % (1)
где г — расстояние между взаимодействующими ионами; — их эффективные заряды; е
- элементарный заряд; Ац, р^ и Сц - постоянные параметры потенциальной модели; первое слагаемое в (1) описывает кулоновское взаимодействие ионов, второе - энергию отталкивания за счет перекрывания электронных оболочек, третье - ван-дер-ваальсовское притяжение.
Кристалл рассматривался как «чисто ионный», т. е. эффективный заряд катионов бария считался равным (+2), анионов фтора - (-1). При этом предполагалось, что двухзарядные катионы Ва2+ из-за достаточно сильного отталкивания не могут находиться вблизи друг от друга, поэтому для них можно не учитывать ван-дер-ваальсовское взаимодействие и эффект перекрывания электронных оболочек, т. е. >1++= С++= 0. Поскольку разноименные ионы сильно притягиваются и находятся на малых расстояниях друг от друга, то для них ван-дер-ваальсовское взаимодействие оказывается «замороженным» и проявляется слабо, т. е. С+- = 0.
Приведем использовавшиеся значения параметров и С^ для взаимодействий
катион-анион и анион-анион для ВаГг [6, 7] .
А---1018, Дж 700,1
р— • Ю10, м 0,2753
С— ■ 1078, Дж-м6 23,07
А+_ • 1018, Дж 841,4
р+_ ■ Ю10, м 0,2792
Число частиц N в кластерах - 768 (512 анионов и 256 катионов), 1500 (1000 анионов и 500 катионов) и 2592 (1728 анионов и 864 катионов). Это означает, что кластер состоит из 64 (4x4x4), 125 (5 х 5 х 5) и 216 (6x6x6) кубических элементарных ячеек флюоритовой решетки, каждая из которых содержит 4 формульных единицы ВаГ2. Параметр элементарной ячейки кристаллической решетки оо взят из [5] и был равен 6,2 А.
Расчеты проводились методом МД, который в применении к исследованным системам подробно изложен в [8]. При решении ньютоновских уравнений движения в численном эксперименте использовался алгоритм Бимана [9]. На расстоянии порядка 30 А от поверхности кластера на систему накладывалась непроницаемая оболочка с «отталкивательным» полем. При построении начальной конфигурации ионы располагались в положениях, соответствующих идеальной решетке.
МД-моделирование осуществлялось в три этапа. На первом происходило уравновешивание системы при температуре 100 К в течение 5000 пс. При этом временной шаг устанавливался такой, чтобы достигнуть полного сохранения энергии с любой выбранной точностью. Далее моделировали процесс равномерного нагревания системы с целью определения температуры плавления. На третьем этапе проводился расчет равновесной локальной структуры и термодинамических свойств кластеров при температурах 1700, 1900 и 2500 К (что соответствует температуре до точки плавления, в ней и после нее) с использованием пакета программ молекулярно-динамических вычислений БЬ_РОЬУ [10] и разработанной нами программы, описанной выше. Временной шаг при движении по фазовой траектории системы был равен 1 фс, общая длина траектории, на протяжении которой изучались свойства системы, составляла 100 пс.
Фазовый переход. Как сказано выше, расчеты начинались с уравновешивания системы при низкой температуре. На рис. 1 показаны флуктуации температуры для кластера из 768 ионов в процессе уравновешивания. Видно, что они не превышают 3°.
После приведения системы к равновесию при 100 К начинали процесс нагревания системы. Каждую пикосекунду скорости ионов системы увеличивались таким образом, чтобы температура системы возрастала на 100 К. На рис. 2 приведена зависимость кинетической энергии ионов Ва+2 от времени для кластера из 768 ионов. При проведении
61
т, к
103 г
102 101 100 99
98 97
1000 2000 3000 4000 5000
I пс 768 ионов
Рис. 1. Зависимость температуры Т от времени t для кластера из
Е- 10Ш, эрг
Рис. 2. Зависимость кинетической энергии ионов Ва+2 {Е) от времени t в кластере из 768 ионов.
первых двух этапов расчета был применен пакет компьютерных программ, разработанный ранее Е. Н. Бродской, согласно алгоритму, использованному в [1].
Процессу плавления микрокристалла должен соответствовать определенный участок на зависимости Т = /(<), где температура колеблется около некоторого постоянного значения. Благодаря малой молекулярной массе фтор обладает высокой подвижностью. В связи с этим участок, соответствующий плавлению ядра кластера, следует искать на зависимости температуры системы от времени для ионов бария (см. рис. 2).
Анализ изменения температуры системы во времени для кластеров различных размеров позволил установить, что температура их плавления составляет примерно 1900 К. Заметной зависимости температуры плавления от размера кластера не наблюдается. В работе [4] по моделированию макроскопической фазы ВаР2 температура 62
Рис. 3. Зависимость локальной парциальной плотности ионов от радиуса кластера (1500 ионов).
1 - для ионов Ва+2 при Т = 1700 К] 2 -для ионов Ва+2 при Т = 2500 К; 3 - для ионов Р" при Т = 1700 К.
плавления была определена в интервале 2080-2140 К. В пределах точности измерения (±100 К) полученные нами результаты находятся в удовлетворительном согласии с данными [4]. Экспериментальная температура плавления кристаллического ВаРг заметно ниже, чем в численном эксперименте: так, в [11] она равняется 1641 К, а в [12] -1628 К. Расхождение между реальным и численным экспериментами объясняется недостатками потенциальной модели.
Локальная структура. Для расчета локальных характеристик в равновесном состоянии (третий этап) микрокристалл разбивали на сферические слои толщиной 0,1 А. Далее в течение 1 пс рассчитывали среднее количество ионов в каждом слое и значение энергии, приходящейся на один ион. В качестве примера на рис. 3 приведены профили парциальных локальных плотностей для ионов Ва+2 и Р- для кластера из 1500 ионов.
Вид профиля локальной плотности Ва+2 в центральной области кластера при 1700 К характерен для высокоструктурированного вещества, что говорит о наличии кристаллического состояния. При 2500 К эта область значительно сокращается, и уже нельзя говорить о наличии дальнего порядка в системе. При движении к поверхности профиль сглаживается и становится характерным для расплава.
В отличие от структурообразующего иона Ва+2 профиль плотности для ионов фтора уже при 1700 К является характерным для разупорядоченной фазы. Это свидетельствует о том, что при этой температуре катионная подрешетка еще кристаллическая, а анионная подрешетка практически полностью разупорядочена, что может косвенно подтверждать наличие суперионной проводимости для исследуемой системы. Такое поведение системы иллюстрируют и приведенные ниже корреляционные функции распределения.
В работе рассчитывались парные радиальные корреляционные функции катион-катион, катион-анион, анион-анион. На рис. 4 приведены функции Р](г)д^{г), где Рз(г) - локальная плотность ионов, - корреляционная функция распределения,
для кластера из 1500 ионов (предплавление - при 1700 К, начало плавления - при 1900 К, расплав - при 2500 К).
При исследованных температурах вид всех корреляционных функций, включая 9Ва-Ва) типичен для расплава, хотя, согласно профилям локальной парциальной плотности В а, в центральной области кластера при 1700 К сохраняется еще кристаллическая структура. Это связано с тем, что основной вклад в корреляционные функции
Р. А
■
6 ' 8 10 12 14
г-10 , м
Рис. 4• Парные радиальные корреляционные функции для кластера ВаРг из 1500 ионов.
а - катион-катион, б - анион-анион, в - анион-катион. 1 - Т = 1700 К; 2 - Т = 1900 К; 3 - Т = 2500 К.
0,020г
происходит от разупорядоченной области и поверхностного слоя кластеров. Признаком присутствия кристаллического ядра в кластере при 1700 К можно считать наличие четко выраженного третьего максимума для рва-Ва, который исчезает при нагревании до 2500 К (рис. 4, а, кривая 1).
Описывая данные функции в целом, можно сказать, что по виду и положению максимумов корреляционные функции аналогичны полученным в численном эксперименте для объемного фторида бария [4]. При 1700 К (рис. 4, кривые 1) ближние корреляции катион-катион и катион-анион характеризуются острыми узкими пиками. Положения максимумов отвечают расстояниям между соответствующими позициями ионов в идеальной флюоритовой решетке. Число ближайших соседей в первой координационной сфере приведено в табл. 1. Как и следовало ожидать, оно уменьшается с увеличением температуры, поскольку система становится более разупорядоченной. Для объемной системы в численном эксперименте [4] число ближайших соседей практически не зависит от температуры в температурном интервале от 1500 до 1900 К и соответствует значениям для идеального кристалла: 12 - для Ва-Ва, 4 - для Ва-Р и б - для Р-Р. Полученные значения первых координационных чисел для кластера оказались существенно ниже, чем в идеальной кристаллической решетке. Это связано с тем, что большая доля ионов в кластере находится в поверхностном слое, для которого кристаллическая решетка в значительной степени разрушена.
Таблица 1. Первые корреляционные коэффициенты для кластера из 1500 ионов
Центральный ион Соседний ион Т, К
1700 1900 2500
Р~ Р" 4,2 4,0 3,5
Р~ Ва+2 3,1 3,0 2,8
Ва+2 Ва+2 4,8 4,5 4,0
При температуре плавления, 1900 К (рис. 4, кривые 2), для катион-катионной и катион-анионной корреляционных функций происходит некоторое снижение первых и уширение вторых пиков. В то же время положение максимумов этих функций меняется незначительно, потому можно предположить, что и в вблизи этой точки значительная часть катионов еще находится около своих «стандартных» кристаллографических положений. В фазе расплава (рис. 4, кривая 3) для катион-анионных корреляций особенно хорошо видно, что второй и третий максимумы сливаются в один, такое сглаживание подтверждает исчезновение дальнего порядка. Корреляции же анион-анион существенно сглаживаются уже при переходе к плавлению. Однако и при температуре 1700 К видно, что пики достаточно уширены, а анион-анионная функция даже в первом минимуме значительно выше нуля, что характеризует подрешетку Р~ в кластере как достаточно разупорядоченную.
Термодинамические свойства и коэффициент диффузии. В качестве термодинамической характеристики кластеров для анионов и катионов рассчитывали парциальную энергию на один ион. Следует указать, что для разнородных пар потенциальная энергия взаимодействия делилась пополам при учете энергии каждого иона. Как и следовало ожидать, полная энергия ионов с ростом температуры снижается по абсолютной величине в связи с большей подвижностью ионов и разупорядоченностью решетки. При увеличении размера кластера изменение полной энергии решетки также носит вполне предсказуемый характер: при заданной температуре энергия становится
65
N 768 1500 2592
Т, К 1700 2500 1700 2500 1700 2500
Ва+2 - 1,5 1Д 1,4 0,4 1Д
Р" 3,5 4,9 2,4 4,7 2,7 3,4
более отрицательной из-за уменьшения относительного вклада поверхности. Полная энергия двухзарядных ионов бария почти в три раза ниже энергии фтора, что объясняет большую устойчивость катионной подрешетки в процессе нагревания.
Подвижность ионов характеризовали коэффициентами диффузии, которые рассчитывали на основе соотношения Эйнштейна для временной зависимости среднеквадратичного смещения
б£>в* + А, = (Дг*(*)>,
где (Дг^(¿)) - средний квадрат смещения частицы сорта ос за время Ба - коэффициент диффузии частиц этого сорта; слагаемое Аа описывает термодинамические флуктуации частиц, не приводящие к систематическому смещению; отсюда
П^шЩР^. (2)
Юо ос
Соотношение (2) применимо в численном эксперименте только в случае «флюидо-подобного» движения частиц, в рассматриваемой системе - для анионов и катионов в расплаве. Результаты расчета коэффициента диффузии приведены в табл. 2.
Для флюидного состояния кластеров следует ожидать увеличения общего коэффициента диффузии при повышении температуры из-за роста подвижности ионов и, наоборот, уменьшение величины с ростом размера кластера из-за сокращения относительного вклада более подвижных ионов, находящихся в поверхностном слое. Эти тенденции для коэффициента диффузии можно наблюдать для ионов фтора, так как данные температуры оказываются выше температуры суперионного перехода для ВаРг [4] и анионная подсистема находится уже в разупорядоченном состоянии. Некоторое нарушение с изменением N от 1500 до 2592 при 1700 К в поведении Бр, которое в пределах погрешности (10%) можно считать постоянным, обусловлено, скорее всего, недостаточной продолжительностью счета для кластера с N = 2592.
Такое же поведение демонстрирует и средний коэффициент диффузии ионов Ва. Это связано с подавляющим вкладом тех катионов, которые расположены в поверхностном слое и в разупорядоченной, находящейся во флюидном состоянии области.
Заключение. По результатам работы можно сделать вывод, что моделирование методом МД кластеров фторида бария с использованием простого модельного потенциала межионного взаимодействия позволяет качественно воспроизвести основные свойства этих систем. Так, получено значение температуры плавления кластера, в пределах погрешности расчета согласующееся с температурой плавления объемного фторида бария, определенной в компьютерном эксперименте. Однако установить температуру плавления кластеров удалось только по динамической зависимости температуры от времени, на которой плавлению отвечает постоянство температуры на некотором временном интервале. Как правило, для определения температуры фазового перехода используются три основных критерия - изменение наклона калорической кривой, скачок величины коэффициента диффузии и изменение вида парных корреляционных 66
функций. Однако для изучаемых малых кластеров они оказались непригодны. Вероятно, этого удастся избежать при применении потенциалов межионного взаимодействия, разработанных специально для моделирования малых фторидных систем.
Согласно полученным результатам, температура плавления кластеров фторида бария практически не зависит от размера кластера, в то время как парциальная энергия ионов и их коэффициенты диффузии с увеличением числа ионов в кластере уменьшаются, а с повышением температуры растут, но такие изменения невелики.
Summary
Terekhina О. S., Brodskaya Е. N., Piotrovskaya Е. М. Computer simulation of properties and structure of baxium fluorides clusters.
Molecular dynamics simulation of small barium fluoride clusters were carried out using a simple Born-Mayer-Huggins model potential. It was shown that such simple model potentials may be used only for qualitative description of small ionic clusters. The obtained data correlates well with the properties calculated earlier for macroscopic BaF2. It was shown that the melting temperature does not almost depend on the cluster size, while partial energy of ions and diffusion coefficients decrease with the increase of the number of ions in the cluster.
Литература
1. Dornford-Smith A., Grimes R. W. 11 Philosophical Magazine. B. 1995. Vol. 72, N 6. P. 563-576. 2. Bulatov V. L., Grimes R. W., Harker A. H. // Phil. Mag. Letters. 1998. Vol. 77. P. 267- 273. 3. Bulatov V. L., Grimes R. W., Harker A. H. 11 J. Materials-e. 1997. Vol. 49, N 2 (http://www.tms.org/pubs/journals/JOM/9704/Bulatov). 4. Готлиб И. Ю., Мурин И. В., Пиотровская Е. М., Бродская Е. Н. // Вестн. С.-Петерб. ун-та. Сер. 4: Физика, химия. 2000. Вып. 2 (№ 12). С. 62-80. 5. Brass А. М. // Phil. Mag. А. 1989. Vol. 59, N 4. P. 843-859. 6. Айтьян С. X., Иванов-Шиц А. К., Бухштаб А. С., Лошманов А. А. Суперионное состояние флюорита BaF2: молекулярно-динамические расчеты и нейтронодифракционный эксперимент. - М., 1990. 38 с. - Деп. в ВИНИТИ 23 мая 1990 г., № 2804-В90. 7. Ivanov-Shitz А. К., Buchtab A. S., Aityan S. К., Kohler Н.-Н. // Appl. Phys. A: Solids and Surfaces. 1992. Vol. 54, N 3. P. 251-257. 8. Allen M. P., Tildesley D. J. Computer simulations in liquids. Oxford, 1987. 9. Schofield P. // Сотр. Phys. Commun. 1973. Vol. 5. P. 17-23. 10. Пакет программ DL.POLY (http://www.cse.clrc.ac.uk/msi/software/DL-POLY/). 11. Термодинамические свойства индивидуальных веществ: В 4 т. / Отв. ред. В. П. Глушко. М., 1978-1982. 12. Воронин Б. М., Присяжный В. Д. // Электрохимия. 1980. Т. 16, № 2. С. 131-137.
Статья поступила в редакцию 23 мал 2006 г.
