Построение потенциала ЕАМ-типа по структурным данным жидкости и кристалла Текст научной статьи по специальности «Физика»

Физика

УДК 536.421

ПОСТРОЕНИЕ ПОТЕНЦИАЛА ЕАМ-ТИПА ПО СТРУКТУРНЫМ ДАННЫМ ЖИДКОСТИ И КРИСТАЛЛА

А.Г. Воронцов, А.А. Мирзоев

В статье предлагается методика построения межчастичного ЕАМ-потенциала по данным дифракционных экспериментов в жидком и кристаллическом состоянии для молекулярно-динамического моделирования процессов плавления и кристаллизации металлов. Предлагаемая методика использована для построения ЕАМ-потенциала для моделирования расплава меди в широком диапазоне температур.

Введение

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

Самая простая форма взаимодействия, которая лежит в основе большинства моделей - парное взаимодействие. Наиболее распространенными являются парные потенциалы самого простого вида, потенциал Леннарда-Джонса, Морзе, Юкавы и др., параметры которых могут быть подобраны для описания той или иной системы. Простая форма потенциалов позволяет получить только качественное согласие моделей с реальными системами. Для количественного описания реальных систем были предложены методы Шоммерса [1] и Реатто [2], восстановления парного потенциала взаимодействия из данных эксперимента (экспериментальной функции радиального распределения атомов). Это методы, в которых реализована итерационная процедура получения парного потенциала, приводящего в процессе моделирования МД к заданной функции радиального распределения. Данные методы оказались полезны для получения потенциалов взаимодействия одно и несколько компонентных жидких и аморфных систем [3].

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

Реалистичное взаимодействие между атомами можно получить только с учетом электронной подсистемы и ее изменения при изменении температуры и плотности. Распределение электронов можно определить методом функционала электронной плотности или квантового Монте-Карло, а затем найти силы, действующие на атомы системы, но нахождение равновесной конфигурации электронов занимает на порядок больше времени, чем моделирование движения ионов (атомов). Успехи, достигнутые в распространении метода ab-initio молекулярной динамики на большие системы крайне скромны. Она позволяет проследить движение порядка 200 ионов на интервале времени около 100 пс. Для моделирования больших систем она оказывается непригодной из-за огромных затрат вычислительных ресурсов.

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

идея получила свое развитие в потенциалах погруженного атома, разработанных впервые Доу и Баскесом [4].

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

U, - Fi(pi) + '^J фу (г) (1)

парного взаимодействия ионов фу(г) и взаимодействия иона с остальной системой i5)(/?,■) посредством электронной плотности pi = где ц/ц{г) - электронная плотность, создаваемая

j*i

j атомом в точке, где находится атом /. Параметры функций взаимодействия, по мнению авторов, необходимо определять из каких-либо дополнительных источников. Для потенциала ЕАМ типа скорость вычисления сил, действующих на атомы, сравнима по быстродействию с таковой для парных потенциалов, но ЕАМ-потенциал позволяет решить фундаментальные противоречия имеющиеся для парного взаимодействия.

В серии статей Баскес, Доу, Фоилс [5-7] предложили форму потенциалов, позволяющих правильно воспроизводить параметры кристаллических Си, Ag, Аи, и других металлов. Основными свойствами, под которые подгонялся потенциал, были: постоянная решетки кристалла, упругие константы, уравнение состояния Росе [8], энергия сублимации, энергия образования вакансии. Долгое время эти потенциалы служили «рабочими» потенциалами для получения интересных результатов, но бурное развитие методов расчета из первых принципов и возрастающие требования к точности модельных расчетов привели к тому, что появилось множество работ, в которых получены новые потенциалы в форме погруженного атома.

Последние работы [9, 10], направленные на получение потенциалов ЕАМ-типа рассматривают все функции, входящие в потенциал, как подгоночные, направленные на правильное воспроизведение экспериментальных свойств. Предложены схемы подгонки потенциала для точного описания температуры плавления [11], сил, полученных методом ab-initio [12], структурного фактора жидкости [13], структурных и термодинамических свойств жидкости [10]. Стоит отметить, что все потенциалы удовлетворительно работают только в небольшом диапазоне температур и испытывают значительные трудности при решении задач, связанных с фазовыми переходами в системе. Например, потенциалы, основанные только на свойствах кристалла, не воспроизводят с нужной точностью свойства жидкости и наоборот. Для получения потенциалов погруженного атома хорошо описывающих и жидкость и кристалл наиболее приемлемой является комбинация классического метода для кристалла [4] с методом, недавно предложенным для жидкости [10, 14]. При таком подходе можно точно воспроизвести структурные свойства жидкости, энергию сублимации, изотермический модуль сжатия, кроме того, оставшиеся параметры позволяют получить правильный параметр решетки кристалла и упругие константы.

Метод

Определения функций, входящих в ЕАМ-потенциал происходило в два этапа. Сначала выбиралась жидкая система недалеко от точки плавления с известной из эксперимента функцией радиального распределения. При этих термодинамических условиях подбирался парный потенциал, при проведении молекулярно-динамического моделирования с которым максимально точно воспроизводились структурные свойства системы (ФРР и давление близкое к нулевому). Затем подбиралась функция погружения, правильно воспроизводящая уравнение состояния в области твердой и жидкой фазы, при различных термодинамических условиях. Подробности метода приведены ниже.

Построение парной части потенциала

Для построения парной части ЕАМ-потенциала меди нами использовалась схема, предложенная Белащенко [10]. Эффективная сила, действующая на атом i со стороны атома j равна

fy=-

ЯL

dp

dF. + —L dlV,j ёф

р-р, dP p=pj dr r=ru dr

(2)

В работах [10, 14] было показано, что на этапе выбора потенциалов в уравнении (2) можно

приближенно заменить сумму производных 5 =--------

ар

средним значением этой сум-

р=р,

мы по всем частицам 5 = 2-------

др

. На самом деле, сумма 5 зависит от выбранной пары частиц и

р=<р>

от времени, однако если система достаточно плотная и флуктуации плотности не слишком велики, то сумма 5 может флуктуировать с относительно небольшой амплитудой и ее можно приближенно принять постоянной. Это существенно упрощает вычисления. Если это предположение принять, то для определенных термодинамических условий на основе уравнения (2) можно свести все взаимодействия в системе к парным взаимодействиям с суммарной эффективной силой:

/.=-2^

СІГ

Щ . (3)

аг

Р=<Р> <Г=ГУ 'Г=Гу

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

N °°

<р>=Ал—^{г)г2у/{г)(іг. (4)

о

Выражение (3) позволяет сразу найти парную часть ЕАМ-потенциала, если рассматривать

аг

систему в точке, где — = 0. Для этого, как и в работе [10], можно воспользоваться методом

ар

Шоммерса [1].

Выбор электронной плотности

Зависимость безразмерной электронной плотности от расстояния до центра атома можно задать в простом виде [10]:

у/(г) = р1ехр(-р2-г). (5)

Параметр р2 выбирается из условия подобия - эмпирического наблюдения [3, 10, 14] о том, что спадание электронной плотности у атома в металле пропорционально положению первого

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

средней плотности (выражение (4)) на единицу при определенной температуре.

Определение функции погружения

Известно, что функция погружения Р(р) по своему смыслу равна 0 при р = 0 и имеет минимум при равновесном значении р [4], кроме того, существуют масштабные преобразований для ЕАМ-потенциала, не изменяющих силы, действующие в системе:

р{г) -> 5* р{г),

р{р)^р{р!8), ф(r)~*ф{r)-2g* р{г).

Подбором параметров і и g можно привести потенциал к каноническому виду, когда энергия погружения имеет минимальное значение в точке термодинамического равновесия. Кроме того, можно нормировать электронную плотность (что мы и сделали) таким образом, чтобы положение минимума функции Р(р) приходилось на значение электронной плотности равное 1.

В разных работах выбиралась различная функция погружения, например в работе [10] функция вида:

ір(р) = ал[р + /3р р < 0,8р0,

[Р( р) = ао + а2 (р - Ро )2 + «3 (р - Ро )3 Р > 0,8р0, где коэффициенты а и (3 определяются из условия гладкой сшивки функции погружения в точке /7 = 0.8р(). Выбор функции погружения в виде (7) позволяет удовлетворить требованиям ра-

венства Р(р) = 0 и минимума Р в точке р- р0 = 1. Параметры щ, а2, аг определялись по полной энергии, коэффициенту изотермической сжимаемости и высокотемпературному уравнению состояния [10].

Нами предлагается использовать экспериментальное уравнение состояния, как в кристаллической, так и в области жидкости. Можно заметить, что физически допустимому для кристалла и расплава изменению средней электронной плотности < р > соответствует диапазон от 0,80 (нагретая жидкость) до 1,2 (кристалл при комнатной температуре). Запишем давления, которые создают парная и ЕАМ части потенциала:

где производная берется при среднем значении плотности в системе, полученной из выражения (4). Проводя моделирование и приравнивая к 0 давление в системе можно получить условие для производной функции погружения:

где точка взятия производной определяется выражением (4). Уравнения (4), (8), (10) записаны в предположении, что можно пренебречь флуктуациями атомной плотности, что, конечно, проверяется при моделировании.

Результаты и обсуждение

Мы проводили построение ЕАМ-потенциала меди. Температура 1423 К была выбрана в качестве опорной точки построения потенциала т.к. для этой температуры известны данные дифракционных экспериментов [15]. Парный потенциал взаимодействия вычислялся методом Шоммерса [1], который позволяет построить парный потенциал, воспроизводящий при моделировании экспериментальную ФРР. Точность воспроизведения данных эксперимента оценивалась при помощи невязки:

где ^г) и go(r) - модельная и экспериментальная функции радиального распределения, заданные в точках г,. При моделировании невязка составила 0,02, при этом две функции радиального распределения практически неразличимы. Экспериментальная и модельная функции радиального распределения приведены на рис. 1. Моделирование проводилось в МУТ-ансамбле, но коррекцией потенциала мы добивались обращения давления в 0. При моделировании было получено значение давления 0,04 ГПа, что является малой величиной в сравнении с флуктуациями давления при моделировании. Полученный парный потенциал приведен на рис. 2.

Для меди параметр р2 был найден из условия подобия [14] и равнялся 1,6 А~'. Коэффициент ръ найденный из условия < р>= 1, равен 4,85. Для определения функции погружения строились модели различного объема и температуры, в которых вычислялась g(r), и из условия 11 находилась производная функции погружения. График функции погружения в области равновесных значений электронной плотности приведен на рис. 3.

Оптимизация функции погружения проводилась путем подгонки модели под реальное уравнение состояния меди. Для этого строился ряд моделей кристаллической меди при температурах 300-1300 К и моделей расплава меди при температурах 1423-1873 К. Модели содержали 4000 атомов и строились в ЫУТ-ансамбле при экспериментальных значениях объема и температуры.

'ч

(9)

(8)

откуда, заменяя суммирование интегрированием, получаем

(10)

(И)

(12)

Моделирование с чистым парным потенциалом, подобранным для температуры 1423 К приводит к появлению значительных давлений в модельной системе, значения которых приведены во втором столбце таблицы.

Рис. 1. Модельная и экспериментальная функции радиального распределения для меди при температуре 1423 К

Рис. 3. Функция погружения. Вертикальными линиями показаны средние значения электронной плотности при различных температурах

Ф,эВ,

0,4 -0,3-0,2-0,1 -0,0---0,1 --0,2--0,3--0,4-

-0,5 ^--,----;---,---,---,----,---,---,---,----;--=~

0 2 4 6 8 10 д

Г, А

Рис. 2. Парный потенциал взаимодействия для расплава меди при температуре 1423 К

3- -о-эксперимент

1 -----------1-------------!-------------'-------------1-------------'------------Г**-

2 4 6 8 0

Г, А

Рис. 4. Экспериментальные [15] и модельные функции радиального распределения расплава меди для ряда температур

Подбор функции погружения производился с целью приведения давления в системе при всех значениях температуры к экспериментальному (близкому к 0) значению. Значения давлений в моделях кристаллической и жидкой меди для полного ЕАМ-потенциала приведено во втором столбце таблицы. Как видно из таблицы, не удается привести к 0 давления при всех температурах, но удается сделать его отклонение от 0 несущественным.

Давление в модельной системе при разной температуре и разном потенциале взаимодействия

Т, К Парный потенциал Р, Ьаг ЕАМ потенциал Р, Ьаг

кристалл 300 -5 497 6 788

500 -18 051 -2 879

700 -24 742 -6 292

900 -28 859 -6 862

I 100 -29 152 -4 433

1 300 -25 420 35

жидкость 1 473 -399 6 133

1 573 1 239 2 045

1 773 3 584 -5 434

1 873 6 244 -5 185

Проверка полученного потенциала для кристалла проводилась моделированием в NPT-ансамбле и проверкой равновесного параметра решетки. Найдено, что построенный ЕАМ-потенциал воспроизводит экспериментальный параметр решетки для всех температур 300-1300 К с точностью более 0,3 %. Поскольку одновременно достигается близкая к нулю (и, следовательно, к равновесному значению) величина давления в системе, то построенный межчастич-ный потенциал успешно воспроизводит уравнение состояния кристалла. Сравнение моделей жидкости для температур 1573-1873 К с данными эксперимента проводилось по функциям радиального распределения. На рис. 4 приведены графики модельных и экспериментальных функций радиального распределения Васеды [15]. Можно заметить, что полученные модели достаточно точно воспроизводят g(r), что свидетельствует о способности предлагаемого потенциала одновременно воспроизводить и широкую область жидкой фазы меди.

Заключение

Предложена схема построения ЕАМ-типа на основе структурных данных жидкости и кристалла. Моделирование меди с полученным потенциалом показало хорошее согласие с фактическим уравнением состояния для меди, как в жидкой, так и в кристаллической фазе.

Работа выполнена при поддержке РФФИ-DFG, грант 07-03-91558-ННИО а и гранта губернатора Челябинской области, урчел_07-03-96038.

Литература

1. Schommers, W. Pair potentials in disordered many particle systems: A study for liquid gallium / W. Schommers // Phys. Rev. A. - 1983. -V. 28. - P. 3599-3605.

2. Reatto, L. Iterative predictor-corrector method for extraction of the pair interaction from structural data for dense classical liquids / L. Reatto, D. Levesque, J.J. Weis // Phys. Rev. A. - 1986. -V. 33.-P. 3451-3465.

3. Белащенко, Д.К. Компьютерное моделирование жидких и аморфных веществ / Д.К. Бе-лащенко. - М.: МИСИС, 2005. - 408 с.

4. Daw Murray, S. Embedded-atom method: derivation and application to impurities, surfaces, and other defects in metals / S. Daw Murray, M.I. Baskes // Phys. Rev. B. - 1984. - V. 29. - P. 6443-6453.

5. Foiles, S.M. Embedded-atom-method functions for the fee metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys / S.M. Foiles, M.I. Baskes, M.S. Daw // Phys. Rev. B. - 1986. - V. 33. - P. 7983-7991.

6. Daw Murray, S. Model of metallic cohesion: The embedded-atom method / S. Daw Murray // Phys. Rev. B. - 1989. - V. 39. - P. 7441-7452.

7. Baskes, M.I. Modified embedded-atom potentials for cubic materials and impurities / M.I. Baskes // Phys. Rev. B. - 1992. - V. 46. - P. 2727-2742.

8. Universal features of the equation of states of metals / H. Rose James, R. Smith John, F. Guinea, J. Ferrante // Phys. Rev. B. - 1984. - V. 29. - P. 2963-2969.

9. Nam, H.S. Solid-liquid phase diagrams for binary metallic alloys: Adjustable interatomic potentials / H.S. Nam, M.I. Mendelev, D.J. Srolovitz // Phys. Rev. B. - 2007. - V. 75. - P. 014204.

10. Белащенко, Д.К. Применение модели погруженного атома к жидким металлам. Жидкое железо / Д.К. Белащенко // Журнал физической химии. - 2006. - Т. 80. - С. 1-11.

11. Sturgeon Jess, В. Adjusting the melting point of a model system via Gibbs-Duhem integration: Application to a model of aluminum / B. Sturgeon Jess, B. Laird Brian // Phys. Rev. B. - 2000. -V. 62.-P. 14720-14727.

12. Ercolessi, F. Interatomic Potentials from First-Principles Calculations: the Force-Matching Method / F. Ercolessi, J.B. Adams // Europhys. Lett. - 1994. - V. 26. - P. 583-588.

13. Mendelev, M.I. Determination of alloy interatomic potentials from liquid-state diffraction data / M.I. Mendelev, D.J. Srolovitz // Phys. Rev. B. - 2002. -V. 66. - P. 1-9.

14. Белащенко, Д.К. Применение модели погруженного атома к жидким металлам. Жидкая ртуть / Д.К. Белащенко // Теплофизика высоких температур. - 2006. - Т. 44. - С. 1-44.

15. Waseda, Y. The Structure of Non-crystalline Materials, Liquid and Amorphous Solids / Y. Waseda. - McGraw-Hill, New York, 1980. - P. 324.

Поступила в редакцию 18 сентября 2007 г.