CC BY
44
УДК 538.911, 538.913
ПРИМЕНЕНИЕ МЕТОДА ШОММЕРСА ДЛЯ РАСЧЕТА ПАРНОГО ПОТЕНЦИАЛА В ЫРТ АНСАМБЛЕ
В.А. Старухин, А.А. Мирзоев, А.Г. Воронцов
Предложен алгоритм расчета эффективного парного потенциала на основе схемы Шоммерса в №Т ансамбле. С помощью данного алгоритма рассчитан эффективный парный потенциал для жидкого железа вблизи температуры плавления. Особенность построенного алгоритма позволила получить потенциал, корректно воспроизводящий структуру и плотность исследуемой системы.
Ключевые слова: эффективный парный потенциал, компьютерное моделирование, неупорядоченные системы, структура жидкости, молекулярная динамика.
Введение
На сегодняшний день компьютерное моделирование представляет собой очень мощное средство для изучения различных веществ, которое дает возможность предсказания их свойств при условиях, недоступных в эксперименте. Существует два основных подхода в этой области - пер-вопринципное и полуэмпирическое моделирование. В первом случае расчет ограничивается несколькими сотнями атомов и малым масштабом времени, но при этом представляет собой наиболее точный квантово-механический способ описания структуры и свойств исследуемой системы. Второй вариант менее точен, поскольку основан на задании эмпирического потенциала межчас-тичного взаимодействия, но при этом существенно менее ресурсоемкий, что позволяет моделировать системы из большого числа частиц.
Среди эмпирических методов наиболее точными являются подходы, использующие для моделирования экспериментальные дифракционные данные о структуре исследуемого объекта. Из них можно выделить два вида алгоритмов. К первой группе относятся методы, которые позволяют построить только структурную модель системы по данным дифракционных экспериментов. Это обратный метод Монте-Карло [1] и силовой алгоритм Менделева Белащеико [2]. Вторую группу составляют методы, которые помимо структурной модели позволяют определить эффективный парный потенциал, воспроизводящий при моделировании известную из эксперимента парную корреляционную функцию (ПКФ). Наиболее удачными примерами этих методов являются метод Шоммерса [3, 4], метод Реатто [5], а также распространенный в настоящее время алгоритм Сопера [6].
С точки зрения возможностей моделирования, методы второй группы представляют значительно больший интерес, поскольку с помощью полученного эффективного потенциала появляется возможность рассчитать не только структурные характеристики системы, но и широкий круг физических свойств вещества. Кроме того, парный потенциал можно использовать для разработки более реалистических потенциалов межчастичного взаимодействия имеющих зависимость от объема, называемых потенциалами погруженного атома (ЕАМ) [7]. Наиболее простым методом, позволяющим рассчитать эффективный парный потенциал, воспроизводящий экспериментальную структуру, является метод Шоммерса. Методы подобные методу Шоммерса, в том числе и основанный на нем метод Сопера, хорошо воспроизводят парную корреляционную функцию и структурный фактор при проведении молекулярно-динамических расчетов в КУ'Г ансамбле (моделирование с фиксированным количеством частиц, объемом и температурой системы), однако при этом не обеспечивают адекватного значения давления в системе. Поэтому использование такого потенциала для расчета физических характеристик может приводить к большой погрешности. В настоящей работе предлагается метод построения эффективного парного потенциала в NРТ ансамбле (в системе постоянным остается количество частиц, давление и температура), основанный на методе Шоммерса и приводящий не только к заданному виду парной корреляционной функции, но и правильному значению равновесной плотности.
Алгоритм
Сохранив математическую основу схемы Шоммерса [3, 4], мы внесли в нее ряд изменений. Наиболее существенное из них состоит в переходе от итерационной схемы для потенциала к итерациям по парным силовым функциям. Для этого уравнения метода Шоммерса были продифференцированы и взяты с обратным знаком. Таким образом, мы пришли к следующей итерационной схеме для силовых функций:
где (г) - эффективная силовая функция на г -й итерации, g(r) - целевая ПКФ, к — постоянная Больцмана, Т - температура, gj(r) - расчетная ПКФ на г'-й итерации, функции со штрихом -производные соответствующих функций.
Переход к силовой итерационной процедуре был сделан по двум причинам. Во-первых, схема Шоммерса в таком представлении позволяет получать более гладкие кривые эффективного потенциала и силовой функции. Дело в том, что в оригинальном методе Шоммерса используется итерационная процедура для потенциала, в которую входит ПКФ полученная из модели. Последняя, в силу ограниченного числа частиц модели, достаточно сильно флуктуирует, что приводит к флуктуациям потенциала парного взаимодействия и, как следствие, немонотонному, скачкообразному поведению силовой функции. Данную проблему можно решить, используя методы регуляризации, например сглаживание кубическими сплайнами [8], но лучшие результаты получаются при использовании силового варианта алгоритма Шоммерса. Во-вторых, использование силового алгоритма Шоммерса позволяет избавиться от сдвига уровня потенциала при проведении итераций, поскольку потенциал рассчитывается интегрированием силовой функции при условии, что потенциал на бесконечности равен нулю.
Критерием близости модели к данным эксперимента служила невязка по ПКФ, которая на / -й итерации рассчитывалась по формуле
где g(r) - целевая ПКФ, gl(r) - расчетная ПКФ на г -й итерации, N - количество точек в целевой и модельной ПКФ. Для уменьшения случайных выбросов и получения стабильной сходимости при расчетах использовалась линейная схема примешивания. В этом случае новое приближение выбиралось по формуле
где Р^г) - эффективная силовая функция на г-й итерации; р - коэффициент примешивания, принимающий значения от 0 до 1.
Еще одним отличием предлагаемого алгоритма от оригинальной схемы Шоммерса является использование NРТ ансамбля вместо МУТ ансамбля при моделировании. Можно ожидать, что если невязка по ПКФ в №Т ансамбле будет уменьшаться, то и плотность системы должна приближаться к экспериментальному значению, что не происходит для NУТ ансамбля. К сожалению, итерационная процедура Шоммерса в N РТ ансамбле значительно менее устойчива. Например, если объем системы не ограничен, то ее случайное расширение может оказаться необратимым. Этот эффект объясняется тем, что притягивающая ветвь потенциала достаточно слаба и ее действие при моделировании ограничено, т.е. система медленно собирается, а если атомы разлетелись на расстояния большие радиуса обрыва потенциала, то сжатия может не происходить вовсе. Поэтому данная схема очень чувствительна к выбору стартового потенциала. В данной работе в качестве стартового потенциала и стартовой конфигурации атомов мы использовали результат рас-
Р](г) = кТ, ё(г)
Р1(г) = Е](г) + кТ
+ kтg,(r)g^(r)-g(r)g[(r)
g{r)g^ (г)
1 J:Tg'(r)gi(r)-g(r)g'j(r)
g{r)gi{r)
(1)
Рм(> Ї^Р^ + кТ
N
'E(g(rj)-gi(rj))2
(2)
ЪГ(г) = рРм(г) + (1-рЩ(г),
(3)
чета методом Шоммерса для ЫУТ ансамбля. Т.е. потенциал, полученный методом Шоммерса в ЫУТ ансамбле, доводится до искомого значения в КРТ ансамбле.
Расчет потенциала
В результате тестовых расчетов было установлено, что схема Шоммерса в КР'Г ансамбле сходится, если стартовый потенциал соответствует атомной плотности, превышающей экспериментальное значение, и одновременно воспроизводит ПКФ. Поэтому процедура расчета потенциала, воспроизводящего ПКФ и плотность в ОТТ ансамбле, была разбита на две стадии. На первом этапе рассчитывался потенциал методом Шоммерса в ИУТ ансамбле, а на втором он уточнялся в МРТ ансамбле.
В качестве целевой ПКФ были выбраны результаты для жидкого железа при температуре 1820 К [9]. Экспериментальная ПКФ показана на рис. 1 сплошной линией. Все парные корреляционные функции задавались в 200 точках, равномерно распределенных по интервалу от 0 до 9 А. Система состояла из 1000 атомов, при этом модель представляла собой куб с длиной грани 23,6497 А, что соответствует фактической атомной плотности 0,0756 А'3. Моделирование производилось с использованием периодических граничных условий. Шаг по времени составлял 1 фс.
г, А
Рис. 1. Парная корреляционная функция жидкого железа при температуре 1820 К [9].
Сплошная линия - экспериментальная кривая. Проколотые кружки - ПКФ, рассчитанная в ЫРТ ансамбле со стартовым потенциалом Леннард-Джонса. Закрашенные кружки - ПКФ, рассчитанная с тем же потенциалом в N>/1 ансамбле
В качестве стартового потенциала для моделирования в ЫУТ ансамбле использовался потенциал Леннард-Джонса. Положение ямы выбиралось по положению первого пика ПКФ равное приблизительно 2,5 А. Для удовлетворения условия о завышенном значении атомной плотности глубина ямы Леннард-Джонса выбиралась достаточно большой: ЕЫп = -\0кТ , где к - постоянная Больцмана, Т - 1820 К - температура. Радиус обрыва потенциала был выбран равным 9 А. Это было обусловлено предыдущими расчетами с большими значениями радиуса обрыва, по итогам которых было замечено, что потенциал после 9 А приблизительно равен нулю. Начальный потенциал показан на рис. 2 сплошной линией.
На рис. 1 приведены кривые ПКФ, рассчитанные для начальных моделей с потенциалом Леннард-Джонса в №3Т и ТчГУТ ансамблях. Расчет в ЫУТ ансамбле выполнялся при атомной плотности жидкого железа равной экспериментальному значению. Расчет в №Т ансамбле выполнялся при нулевом давлении в баростате. Как видно из рисунка, ПКФ для двух моделей немного отличаются друг от друга и совершенно не похожи на экспериментальную ПКФ. Это связано с тем, что используемый потенциал Леннард-Джонса соответствует атомной плотности отличной от экспериментальной, а именно 0,0912 А“3. Невязка (2) для результата в ЫУТ ансамбле
составила 4,52 ■ 10“2, в №Т ансамбле - 1,04 • 10“2.
г, А
Рис. 2. Эффективный парный потенциал. Сплошная линия - стартовый потенциал Леннард-Джонса. Проколотые кружки - результат моделирования в МТ ансамбле.
Закрашенные кружки - результат моделирования в №Т ансамбле
Для восстановления потенциала методом Шоммерса в ЫУТ ансамбле на каждой итерации после изменения потенциала выполнялась релаксация модели в течение 5000 шагов моделирования, после этого в течение 1000 шагов с интервалом в 10 шагов накапливалась информация для вычисления модельной ПКФ. Для достижения сходимости при моделировании в ЫУТ ансамбле потребовалось 30 итераций с коэффициентом примешивания 0,5. Был получен эффективный потенциал (рис. 2 проколотые кружки), хорошо воспроизводящий целевую ПКФ с невязкой
1,48-10_3. При таком значении невязки целевая ПКФ визуально совпадает с ПКФ модели (рис. 3).
г, А
Рис. 3. Парная корреляционная функция жидкого железа при температуре 1820 К [9]. Сплошная линия - экспериментальная кривая. Закрашенные кружки - ПКФ, рассчитанная в N>/1 ансамбле
Полученный при моделировании с ТЧУТ ансамблем потенциал использовался в качестве стартового для ОТТ ансамбля. Процедура применения метода Шоммерса в этом случае была практически идентичной описанной выше. Отличия состояли в следующем. Во-первых, для релаксации модели использовалось 10 000 шагов вместо 5000 шагов, т.к. релаксация объема происходит медленнее, чем давления. Во-вторых, усреднение ПКФ проводилось по 1000 кривым, что уменьшало статистическую ошибку. Флуктуации по объему при этом составляли около 50-100 А3. И, в-третьих, на каждой итерации в качестве стартовой конфигурации использовалась одна и та же модель, которая была получена для потенциала в ЫУТ расчете. Это позволило избе-
жать ситуации, когда при запуске расчета методом молекулярной динамики силы в системе с новым потенциалом оказываются слишком большими, что приводит к быстрому расширению системы. Баростат в №Т ансамбле устанавливался в нулевое значение. Коэффициент примешивания был задан равным 0,3.
Г, А
Рис. 4. Парная корреляционная функция жидкого железа при температуре 1820 К [9]. Сплошная линия - экспериментальная кривая. Проколотые кружки - ПКФ стартовой конфигурации атомов для №Т ансамбля. Закрашенные кружки - конечная
ПКФ для ЫРТ ансамбля
—2 —3
В результате расчета в ЫРТ ансамбле невязка по ПКФ опустилась с 1,04 10 до 2,34-10 . Модельная ПКФ при этом существенно приблизилась к целевой кривой (рис. 4). Для достижения этого результата потребовалась 21 итерация методом Шоммерса в ЫРТ ансамбле. Итоговый эффективный парный потенциал приведен на рис. 2 (закрашенный кружки). Атомная плотность, как и ожидалось, тоже изменилась, опустившись на первой итерации в ЫРТ ансамбле со значения
0,0861 до 0,0756 А”3 на итерации 21, достигнув экспериментальной величины. Для сравнения, потенциал Леннард-Джонса, использованный на первой итерации в ЫУТ ансамбле, соответствует значению атомной плотности 0,0912 А"3.
Анализ полученного потенциала
Как видно из рис. 2, существенно изменилась отталкивающая часть эффективного потенциала. Слишком крутая отталкивающая часть Леннард-Джонса стала намного мягче. При этом притягивающая осталась приблизительно того же порядка. В отличие от исходного потенциала Лен-нард-Джонса корректировка потенциала по ПКФ методом Шоммерса привела к появлению осциллирующего хвоста, который объясняют эффектом экранирования ионного остова атомов электронами проводимости. Доводка потенциала в ЫРТ ансамбле привела к незначительным поправкам к потенциалу, полученному в 1ЧУТ ансамбле. Хотя ПКФ при этом претерпела существенные изменения (рис. 4).
Изменение расчетной ПКФ при доводке в ЫРТ ансамбле особенно заметно в области второй и последующих координационных сфер (рис. 4). Все пики сдвинулись вправо и уменьшились по амплитуде. Также заметно изменение формы кривой. Например, самым заметным изменениям подверглась впадина за первым пиком. При всем этом потенциал изменился в основном только в области ямы, немного увеличив ее глубину и сместив ее положение вправо. Хвост потенциала практически не изменился.
В таблице приведены значения атомной плотности п, потенциальной энергии 11роЬ коэффициента самодиффузии Д и коэффициента теплового линейного расширения ои Доверительный интервал приведен для вероятности 95 %. Потенциальная энергия рассчитана по значению теплоты испарения [13]:
А11 -АН-ЯТ, (4)
где Аи - изменение внутренней энергии, АН - теплота испарения, Я - универсальная газовая постоянная, Т - температура. Коэффициент теплового линейного расширения рассчитан по значениям массовой плотности жидкого железа [14]:
1 ь-> - А
а=-------------Ц (5)
и
■Т,
(6)
где Ь - линейные размеры образца, т - масса, р - массовая плотность, Т - температура. Величина Ь рассчитывалась по формуле (6) при массе для одного моля вещества. Параметр а рассчитывался по формуле (5) при нагревании с 1823 К до 1873 К.
Для сравнения в таблице приводятся значения величин, полученных с парным потенциалом, использованным для разработки потенциала погруженного атома [10]. Несмотря на хорошее воспроизведение плотности, процедура получения этого потенциала требует ручной корректировки параметров, поэтому ее сложно сделать автоматической.
Расчетные свойства жидкого железа при температуре 1820 К с различными парными потенциалами
и, А3 - Vро,, кДж/моль £>, Н) 5 см2/с а о 1 С'
Эксперимент 0,0756 [11] 339,2 [12] 13,0 [13] 55,0 [14]
Парный потенциал [10] 0,07529±0,00008 175,7±0,2 3,61±0,15 31,6±0,8
Стартовый потенциал Леннарда - Джонса 0,094896±0,000002 1220,704±0,007 0,0018±0,0008 2,72±0,02
Результат ]ч!\Т 0,0860±0,0001 468,7±0,2 2,32±0,19 16,0±0,6
Результат ОТТ 0,07553±0,00004 458,0±0,1 3,26±0,13 8,2±0,6
Два существенно различных парных потенциала (потенциал из работы [10] и потенциал, полученный в данной работе) хорошо воспроизводят экспериментальную ПКФ и плотность системы, но при этом дают существенно различные результаты для коэффициента теплового линейного расширения и потенциальной энергии. Это явилось следствием неоднозначности воспроизведения эффективного парного потенциала по структурным данным. В работе [15] было показано, что парный потенциал определен с точностью до функции, использование которой в качестве потенциала при моделировании методом непрерывной статической релаксации [16] приводит к воспроизведению этой же структуры при нулевой температуре. Поэтому значение других физических характеристик тоже должно быть использовано в качестве критерия при расчете эффективного потенциала.
Сравнивая результат метода Шоммерса в ЫУТ ансамбле и его доводки в ОТТ ансамбле, видно, что помимо плотности изменились все остальные параметры. Самое существенное изменение получил коэффициент теплового линейного расширения, уменьшившись почти в два раза. Хотя при этом потенциалы очень близки друг к другу (рис. 2). Отсюда следует высокая чувствительность некоторых физических характеристик к форме эффективного потенциала.
Заключение
В результате работы был получен алгоритм расчета эффективного парного потенциала в ОТТ ансамбле на основе схемы Шоммерса. Получаемый потенциал позволяет не только хорошо воспроизводить экспериментальную кривую ПКФ, но и атомную плотность системы при расчете в ОТТ ансамбле с нулевым давлением. Тестовый расчет эффективного парного потенциала жидкого железа при температуре 1820 К привел к совпадению расчетной и экспериментальной атомной
плотности и достаточно низкой невязке по ПКФ (2,34 • 10-3), что при сравнении на графике соответствует визуальному совпадению кривых экспериментальной и модельной ПКФ.
Однако алгоритм оказался неустойчив по отношению к стартовому потенциалу. Расчеты показали, что расчетная схема сходится, если стартовый потенциал хорошо воспроизводит ПКФ в NУТ ансамбле и дает завышенную в сравнении с экспериментальной плотность. Поэтому предложенная схема на данный момент представляет собой не самостоятельный метод, а скорее способ тонкой настройки потенциала для корректного воспроизведения ПКФ в более естественных условиях ОТТ ансамбля.
Неоднозначность воспроизведения эффективного парного потенциала по данным структурного эксперимента, продемонстрированная в работе [15], приводит к возможности получения потенциалов соответствующих одной и той же структуре, но с различными физическими свойствами. Поэтому при расчете эффективного парного потенциала, помимо воспроизведения структуры, необходимо также накладывать условие на воспроизведение при моделировании других физических свойств системы.
Работа выполнена при поддержке грантов РФФИ 10-03-96088 р_урал_а, 11-08-00891-а.
Литература
1. Kaplow, R. Atomic arrangement in vitreous selenium / R. Kaplow, T.A. Rowe, B.L. Averbach // Phys. Rev. - 1968. -V. 168, № 3. _p. Ю68-1079.
2. Менделев, М.И. Новые алгоритмы реконструкции атомных моделей жидких и аморфных тел / М.И. Менделев, Д.К. Белащенко // Расплавы. - 1992. - № 4. - С. 60-65.
3. Schommers, W. A pair potential for liquid rubidium from the pair correlation function / W. Schommers//Physics Letters A. - 1973. - V. 43. - P. 157-158.
4. 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.
5. Levesque, D. Pair interaction from Structural Data for Dense Classical Liquids / D. Levesque, J.J. Weis, L. Reatto // Phys. Rev. Lett. - 1985. - V. 54. - P. 451-454.
6. Soper, A.K. Empirical potential Monte Carlo simulation of fluid structure / A.K. Soper // Chem. Phys. - 1996. - V. 202. - P. 295-306.
7. Murray, S.D. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals / S.D. Murray, M.I. Baskes // Phys. Rev. B. - 1984. - V. 29. - P. 6443-6453.
8. Шикин, E.B, Кривые и поверхности на экране компьютера / Е.В. Шикин, А.И. Плис. - М.: ДИАЛОГ-МИФИ, 1996.-240 с.
9. Il’inskii A. Structure of Liquid Fe-Al Alloys / A. Il’inskii, S. Slyusarenko, O. Slukhovskii et al. // Mater. Sci. Eng. A. - 2002. - V. 325. - P. 98-102.
10. Белащенко, Д.К. Применение модели погруженного атома к жидким металлам. Жидкое железо / Д.К. Белащенко // Журнал физической химии. - 2006. - Т. 80, № 5. - С. 872-883.
11.Waseda. Y. The structure of Non-crystalline Materials. Liquids and Amorphous Solids / Y. Waseda. - New York: McGraw-Hill, 1980. - 325 p.
12. Химическая энциклопедия / под. ред. И.Л. Кнунянца. - М.: Советская энциклопедия, 1990.-Т. 2.-671 с.
13. Yang, L. Self diffusion in molten Fe-C alloys / L. Yang, M.T. Simnad, G. Derge // J. Metals. Soc. Am. - 1956,-V. 206.-P. 1577-1580.
14. Еланский, Т.Н. Строение и свойства металлических расплавов / Г.Н. Еланский, Д.Г. Еланский. - М.: МГВМИ, 2006 - 228 с.
15. Belashchenko, D.K. Families of Intermolecular Potentials Corresponding to Identical Structure of Noncrystalline Substances in the Method of Molecular Dynamics / D.K. Belashchenko // Russian Journal of Physical Chemistry. - 2004. - V. 78, № 9. - P. 1621-1628.
16. Belashchenko, D. K. Structural Modeling of Amorphous Iron / D.K. Belashchenko // Fiz. Met. Metalloved. - 1985. -V. 60, № 6. - P. 1076-1080.
Поступила в редакцию 10 марта 2011 г.
THE APPLICATION OF SCHOMMERS SCHEME FOR THE CALCULATION OF EFFECTIVE PAIR POTENTIAL IN NPT ENSEMBLE
An algorithm for effective pair potential obtaining was elaborated. It is based on Schommers scheme and molecular dynamics simulation in NPT ensemble. An effective pair potential for liquid iron near its melting point was calculated. Peculiarity of the new scheme allowed us to obtain an effective potential which could reproduce both structure and density of the system.
Keywords: effective pair potential, computer simulation, disordered systems, liquid structure, molecular dynamics.
Staruhin Victor Andreevich is Post-graduate Student, General and Theoretical Physics Department, South Ural State University.
Старухин Виктор Андреевич - аспирант, кафедра общей и теоретической физики, Южно-Уральский государственный университет.
e-mail: victor_staruhin@mail.ru
Mirzoev Alexander Aminulaevicb is Dr. Sc. (Physics and Mathematics), Professor, General and Theoretical Physics Department, South Ural State University.
Мирзоев Александр Аминулаевич - профессор, доктор физико-математических наук, кафедра общей и теоретической физики, Южно-Уральский государственный университет.
e-mail: mirzoev@physics.susu.ac.ru
Vorontsov Alexander Gennadevich is Cand. Sc. (Physics and Mathematics), Associate Professor, General and Theoretical Physics Department, South Ural State University.
Воронцов Александр Геннадьевич - кандидат физико-математических наук, доцент, кафедра общей и теоретической физики, Южно-Уральский государственный университет.
e-mail: sas@physics.susu.ac.ru
