CC BY
8НАНОСИСТЕМЫ
DOI: 10.17725/rensit2021.13.149
Автокорреляционная функция скорости и коэффициент самодиффузии в больших молекулярно-динамических моделях жидкого аргона и воды Наберухин Ю.И., Аникеенко А. В., Волошин В.П.
Институт химической кинетики и горения им. В.В. Воеводского СО РАН, http://www.kinetics.nsc.ru
Новосибирск 630090, Российская Федерация
E-mail: [email protected], [email protected], [email protected]
Поступила 28.04.2021,рецензирована 07.05.2021, принята 31.05.2021
Аннотация: В моделях жидкого аргона и воды рассчитаны автокорреляционные функции скорости частиц Z(t) методом молекулярной динамики. Большой размер моделей (больше сотни тысяч частиц) позволил проследить эти функции до 50 пикосекунд в аргоне и до 10 пикосекунд в воде и добиться точности расчёта, достаточной для аналитического анализа их формы. Проанализировано различие в определении коэффициента самодиффузии, даваемое законом Эйнштейна и интегралом от Z(t) (интеграл Грина-Кубо), которое в лучшем случае составляет 3%. Асимптота функции Z(t) в аргоне близка по форме к предсказуемому гидродинамикой виду at—3/2, но с амплитудой, зависящей от рассматриваемого интервала времени. В воде асимптота Z(t) не имеет ничего общего с таковой в аргоне: у неё а < 0 и показатель близок к -5/2, а не к -3/2.
Ключевые слова: законы диффузии Эйнштейна и Грина-Кубо, асимптота автокорреляционной функции скорости
УДК 538.91
Для цитирования: Наберухин Ю.И., Аникеенко А. В., Волошин В.П. Автокорреляционная функция скорости и коэффициент самодиффузии в больших молекулярно-динамических моделях жидкого аргона и воды. РЭНСИТ, 2021, 13(2):149-156. DOI: 10.17725/rensit.2021.13.149._
Velocity autocorrelation function and self-diffusion coefficient in
large molecular dynamics models of liquid argon and water
Yuri I. Naberukhin, Alexey V. Anikeenko, Vladimir P. Voloshin
Voevodsky Institute of Chemical Kinetics and Combustion SB RAS, http://www.kinetics.nsc.ru/ Novosibirsk 630090, Russian Federation
E-mail: [email protected], [email protected], [email protected] Received 27 April, 2021, peer-reviewed 10 May, 2021, accepted 24 May, 2021
Abstract: Autocorrelation function of the particle velocity Z(t) is calculated using the molecular dynamics method in the models of liquid argon and water. The large size of the models (more than a hundred thousand particles) allowed us to trace these functions up to 50 picoseconds in argon and up to 10 picoseconds in water, and to achieve a calculation accuracy sufficient for analytical analysis of their shape. The difference in the determination of the self-diffusion coefficient using Einstein's law and the integral of Z(t) (Green-Kubo integral) is analyzed and it is shown to be 3% at best when t is of the order of several picoseconds. The asymptote of the function Z(t) in argon is close to the power law at-3/2 predicted by hydrodynamics, but with an amplitude that depends on the time interval under consideration. In water, the asymptote of Z(t) has nothing in common with that in argon: it has а < 0 and the exponent is close to -5/2, and not to -3/2. Keywords: diffusion laws of Einstein and Green-Kubo, asymptote of the velocity autocorrelation function UDC 538.91
НАБЕРУХИН Ю.И., АНИКЕЕНКО А.В., ВОЛОШИН В.П.
НАНОСИСТЕМЫ
For citation'. Yuri I. Naberukhin, Alexey V. Anikeenko, Vladimir P. Voloshin.Velocity autocorrelation function and self-diffusion coefficient in large molecular dynamics models of liquid argon and water. KENSIT, 2021, 13(2):149-156. DPI: 10.17725/rensit.2021.13.149._
Содержание
1. Введение (150)
2. Компьютерные модели (151)
3. Расчёт коэффициента самодиффузии
по автокорреляционной функции
скорости (152)
4. Асимптота автокорреляционной функции скорости (152)
5. Момент пересечения
автокорреляционной функции скорости оси абсцисс (154)
6. Заключение (154) Литература (155)
1. ВВЕДЕНИЕ
Коэффициент самодиффузии в жидкостях, D, в реальных экспериментах определяется из среднего квадратичного смещения частицы ^ДR2(tf) на основе закона диффузии Эйнштейна
(.AR\t)) = 6Dt.
частиц (АКФС). Основное соотношение между этими характеристиками имеет вид [13]
г
AR2(t)} = 6$(t-s)Z(s)ds,
(2)
или
AR¿(t))/6t =
t t 6 jz(s)ds— jsZ(s)ds = A(t) + B(t).
(3)
\-----/ --■ С1)
Этот закон не является строгим и выполняется только прп макроскопическом масштабе времени диффузии /, когда соотношение (1) является следствием закона диффузии Фика. Более фундаментальная связь среднего квадратичного смещения с микроскопическими характеристиками теплового движения выражается через автокорреляционную функцию скорости
Здесь Z(f) — автокорреляционная функция скорости частиц v(7), которая задаётся формулой
Z(0 = (v(0)v(0)/3, (4)
где скобки <... > означают усреднение по ансамблю. Закон диффузии Эйнштейна (1) получается из (3) как асимптота при очень большом /, когда A(f) » Wí) \. Коэффициент самодиффузии есть теперь
оо
D = \Z{t)dt. (5)
о
Эта формула выражает связь макроскопических и молекулярных величин и представляет собой один пз вариантов так называемых уравнений Грпна—Кубо [2, 3].
Универсальное соотношение (2) известно уже давно, но практически не использовалось,поскольку автокоррелятор
Water
Рис. 1. Но;
для жидкого аргона (слева) и воды (справа).
НАНОСИСТЕМЫ
АВТОКОРРЕЛЯЦИОННАЯ ФУНКЦИЯ СКОРОСТИ И КОЭФФИЦИЕНТ 151 САМОДИФФУЗИИ В БОЛЬШИХ МОЛЕКУЛЯРНО-ДИНАМИЧЕСКИХ МОДЕЛЯХ...
скорости Z(t) нельзя определить из эксперимента. Ситуация изменилась, когда метод молекулярной динамики (МД) сделал возможным рассчитывать Z(t) в жидкостях. Функции Z(t), типичный вид которых показан на рис. 1, очень быстро спадают со временем. При рассмотрении рисунка 1 создаётся впечатление, что хорошую оценку интеграла Грина—Кубо (5) можно получить, если в качестве верхнего предела в нём взять время 2-3 пикосекунды. Таким образом, в молекулярно-динамическом моделировании можно использовать два рецепта для расчёта коэффициента самодиффузии: на основе закона Эйнштейна (1) и при помощи уравнения Грина—Кубо (5). Точнее, уравнение (3) можно записать в виде:
Einstein ^ ' Green-Kubo
Первый рецепт задаётся функцией слева:
DEinstein (О = (М2(0)Ш. функция DGreen-Kubo(t) справа задаёт второй рецепт, в котором подразумевается, что интеграл Грина-Кубо
(5) берётся в конечных пределах, то есть это
(•t
теперь функция A(t) = I Z(s)ds. Возникает вопрос, насколько согласуется применение этих рецептов в конкретных МД расчётах. В данной работе мы предприняли исследование этой проблемы на основе универсальных соотношений (2) и (3). При этом, мы воспользовались большими МД моделями жидкого аргона, построенными ранее нами в [4,5], и специально рассчитали МД модель воды из 157 464 частиц.
2. КОМПЬЮТЕРНЫЕ МОДЕЛИ
Были построены две молекулярно-динамические (МД) модели жидкого аргона из 500 000 и 3 200 000 атомов в кубической ячейке с периодическими граничными условиями. Использовался потенциал Леннард-Джонса с параметрами о = 0.3405 нм, s/kB = 119.8 K, он обрезался на расстоянии 2.5о и сдвигался вверх для устранения разрыва на радиусе
обрезания. Молекулярно-динамическое моделирование проводилось в пакете GROMACS при шаге интегрирования 2 фс с применением стандартных коррекций для энергии и давления. Основное моделирование проводилось при постоянном объёме и температуре T = 101.83 К. В приведённых единицах у обоих моделей температура, плотность и давление составляли T* = 0.85, р* = 0.84, P* = 0.85, а ребро ячейки было a = 84о = 28.6 нм и a = 156о = 53.2 нм для малой и большой модели, соответственно.
Модели воды содержали 8 000 и 157 464 молекул, взаимодействующих с потенциалом TIP4P/2005 [6]. Электростатические взаимодействия рассчитывались методом Particle Mesh Ewald, и радиус обрезания короткодействующих взаимодействий
составлял 1 нм. В моделировании поддерживались постоянная температура T = 300 K и плотность р = 996 кг/м3. Ребро ячейки в большой модели составляло a = 16.8 нм. АКФС воды рассчитывалась для атомов кислорода.
Из-за большого числа частиц в некоторых моделях и ограниченного размера оперативной памяти компьютера рассчитать АКФС для полной траектории не всегда представляется возможным. Тогда расчёт проводился для коротких последовательных участков траектории, а полученные для них
Параметры МД вычисления
Таблица 1
траекторий, использованных для АКФС в моделях аргона и воды.
Тип Обозна- Чмсло Длина Длина Кол. Шаг
моде- чение молекул основ- отрезка учас- со-
ли модели ной усредн. тков хран,
траек- АКФС, ус- фс
тории, пс ред-
нс не-
ния
аргон Argon500-I 500 000 20 1000 20 500
аргон Argon500-II 500 000 1 20 50 20
аргон Argon3200 3200000 0.2 20 10 500
вода Water157-I 157 464 11 500 22 200
вода Water157-II 157 464 1.1 50 22 20
вода Water8 8 000 0.12 1 120 2
НАБЕРУХИН Ю.И., АНИКЕЕНКО А.В., ВОЛОШИН В.П.
НАНОСИСТЕМЫ
Е
о
со тз с го
А№
ад
о
ю О
00 тз с го
Рис. 2. Функции А^) = ^7и = из уравнения (3) для моделей Аг^опЗОО-П (слева)
\VaterS (справа). 0
функции АКФС усреднялись между собой. Длины траекторий, количество участков усреднения АКФС и параметры сохранения конфигураций перечислены в Таблице 1.
3. РАСЧЁТ КОЭФФИЦИЕНТА САМОДИФФУЗИИ ПО АВТОКОРРЕЛЯЦИОННОЙ ФУНКЦИИ СКОРОСТИ На рис. 2 показано поведение функций -Л(Т) и Б(/) для аргона и воды. Мы видим, что если ограничиться при интегрировании в формулах (3) пределом t — 2.5 пс (что зачастую делается при оценках коэффициента самодпффузпп в МД), то интегралом Б(/) вряд ли можно пренебречь по сравнению с интегралом А(Ч). Для аргона отношение В(2.5 пс)/ А(2.5 пс) составляет 4.5%, а для воды даже 17.2%. Более того, отношение 5(/)/^4(/) убывает со временем столь медленно, что для аргона оно при / = 20 пс составляет всё ещё 3.0%. Для воды функция спадает значительно быстрее, чем у аргона (см. раздел 4), но и здесь отношение В (20 пс) /А(20 пс) составляет заметную величину 2.6%.
Итак, можно утверждать, что два рецепта расчёта коэффициента самодиффузии согласуются между собой с точностью не более чем несколько процентов, если ограничить верхний предел интегрирования в (3) двумя-тремя ппкосекундамп.
4. АСИМПТОТА АВТОКОРРЕЛЯЦИОННОЙ ФУНКЦИИ СКОРОСТИ
Отмеченное выше медленное убывание автокорреляционной функции скорости обусловлено долговременным «хвостом» в её асимптоте. Особый интерес к поведению асимптоты АКФС возник после важной работы Олдера и Вэйнрайта [7]. Они обнаружили, что асимптота функции убывает в плотных жидкостях не экспоненциально (как в разреженных газах), а по степенному закону аГ3/2. Это приводит к тому, что в АКФС появляется длинный хвост, который хотя п имеет очень малую интенсивность, но зато длится очень долго. Именно этот долговременной хвост и приводит к результатам, изложенным в предыдущем разделе. Выявлению закона Г3/2 в МД моделях разных жидкостей посвящено много работ (упомянем для примера только [8, 9]), п во многих случаях такой закон как будто обнаруживается. Однако для такого рода исследований нужен расчёт очень длинных траекторий частиц и тщательное усреднение шумов, что связано с большими трудностями в МД моделировании. Кроме того теоретики выяснили, что закон Г3/2 является лишь приближением, и асимптота должна описываться бесконечным рядом степенных членов [10]. Истинное поведение асимптоты АКФС в жидкостях
НАНОСИСТЕМЫ
АВТОКОРРЕЛЯЦ1ЮННАЯ ФУНКЦ11Я СКОРОСП111КОЭФФ1Щ1IEHT 15 3 САМОД11ФФУ31IIIВ БОЛЫШIX МОАЕКУЛЯРНО-Д1IHAMII4ECKIIX МОДЕЛЯХ...
i
о £
Argon500-I
V
У = o.oos^2**^''"'
10 15 20 25 30 35
40 45
t, ps
0,0014
Li_ 0,0012-
O
^ 0,0010 0,00080,0006 0,00040,0002
0,0000-
у = 0,0084*Г у = 0,00727
6 7
t, ps
Рис. 3. Автокорреляционная функция скорости в моделях аргона, а) Три модели с двумя вариантами степенной асимптоты. Ь) Асимптота модели Аг£оп500-1. с) Модель Аг^опЗ200 с двумя вариантами степенной асимптоты.
наблюдение уже было сделано в работе
нельзя сейчас считать выясненным. Поэтому исследование асимптоты АКФС в наших больших моделях представляется полезным.
На рис. 3 показаны разные фрагменты АКФС в исследованных намп моделях аргона. Рпс. Ъа демонстрирует, что параметры расчёта траекторий в трёх моделях удаётся подобрать так, чтобы АКФС в них былп вполне совместимыми друг с другом. Рисунки Ъа и ЪЬ показывают, что асимптоту в модели \rgon500-I можно удовлетворительно натянуть на функцию [7) = 0.0084-Г3/2 в интервале времени 13-50 ис. Хотя при / > 35 пс заметно возрастает нерегулярность рассчитанных данных, эта функция как будто работает и здесь. Однако успех этого фитинга не означает, что определён вид асимптоты. Рис. Ъс показывает, что в очень большой модели .\rgon3200 асимптота на интервале 5.510 пс неплохо натягивается на функцию с другой амплитудой ^ (/) = 0.00727-Г3/2; но она плохо работает при / > 10 пс (см. рпс. 2а). Всё это означает, что асимптота вряд лп описывается простой степенной функцией типа
Заметим, что момент времени ~35
[8]. Время возвращения — это время, за которое флуктуация плотности проходит длину бокса п может исказить МД расчёт. Время / пропорционально кубическому корню от числа частиц в боксе. В работе [8] было 11 000 атомов и / ~ 10 пс, у нас 500 000 частиц п, значит, должно быть / ~ 35 пс. (При этой оценке мы использовали значение скорости звука пз [8] с параметрами состояния системы, близкими к нашим).
Автокорреляционная функция
скорости в моделях воды представлена на рис. 4. Рпс. 4а показывает, что модель ^7а1ег8 не даёт достоверных данных уже начиная с 2 пс, когда появляется нерегулярное поведение коррелятора с примерно одинаковой амплитудой хаотических колебаний вокруг осп абсцисс, которое длится до 20 пс. Модель ^7а1ег157-П специально сконструирована так, чтобы обеспечить достоверный расчёт асимптоты АКФС вплоть до 5 пс. Модель ^7а1ег157-1 позволяет ещё лучше усреднить шумы, так что удаётся проследить асимптоту до —10 пс (рпс. АН). Мы видим, что асимптота АКФС пс, при котором шумы в МД расчётах воды ведёт себя совершенно иначе, чем
автокоррелятора скорости резко возрастают (см. рисунки 2а и 2Ь), совпадает с так называемым «временем возвращения»
(recurrence time)
'r-
Аналогичное
в аргоне. Во-первых, автокоррелятор во всей исследованной области оказывается отрицательным. Только при / > 10 пс он, возможно, становится положительным
НАБЕРУХИН Ю.И., АНИКЕЕНКО А.В., ВОЛОШИН В.П.
НАНОСИСТЕМЫ
о.
N
у = - 0.002*Г
V\éter157-11 V\éter157-I
0,5 1,0 1,5 2,0 2,5 3,0 3,5 4,0 t, ps
2 3 4 5 6 7 8 910 20 30 40
t, ps
2 3 4 5
f, ps
b
Рис. 4. Нормированная автокорреляционная функция
\Vater157-1 и \Vater157-
(рис. АН), но из-за шумов в этой области достоверных выводов сделать нельзя. Во-вторых, аналитический вид асимптоты отнюдь не описывается функцией аГ3/2, а скорее соответствует функции аГ5/2 (рисунки 4Ь и Ас). Известно, что хвост асимптоты вида аГ5/2 встречается в автокорреляторах вращательной скорости в суспензиях Броуновских частиц [11,12]; но какое это имеет отношение к воде? Эти удивительные результаты были совершенно неожиданными для нас, и пока мы можем только констатировать их как факт.
5. МОМЕНТ ПЕРЕСЕЧЕНИЯ АКФС ОСИ АБСЦИСС
Проделанные расчёты позволяют вычислить значения при / < 3 пс с большой точностью. Это даёт возможность
Рис. 5..
аргона Argon500-II вблизи пресечения оси абсцисс.
скорости в моделях воды, а) Модель WaterS, b) Модели II, с) Модель Waterí 57-II.
выяснить некоторые отдельные моменты поведения автокоррелятора, на которые обычно не обращают внимания. Примечательной его характеристикой является, например, момент пересечения осп абсцисс, т.е. момент перехода пз области отрицательных значений в положительную область, где коррелятор и приобретает асимптоту вида аГ3/2. В моделях аргона это пересечение происходит при / = 2.14 пс (рис. 5).
В автокорреляторе скорости для воды такое пересечение нам установить не удалось. Возможно, что п в воде коррелятор в конце концов выйдет в положительную область. Но это может произойти только при t > 10 пс (см. рис. 4с), значительно позже, чем у аргона. Будет лп асимптота тогда подчиняться закону at—3/2 — это очередная загадка в списке аномальных свойств воды.
6. ЗАКЛЮЧЕНИЕ
Если под коэффициентом самодпффузпп понимать величину D, входящую в закон диффузии Эйнштейна (^ÁR2(tfj = 6Dt, то он имеет точный смысл только в макроскопическом масштабе времени. На небольших временах / это соотношение определяет коэффициент диффузии лишь приближённо. На микроскопическом, молекулярном уровне закономерности теплового перемещения частиц п, в частности, среднее квадратичное смещение
ÁR~(t)), строго описываются на языке
АВТОКОРРЕЛЯЦИОННАЯ ФУНКЦИЯ СКОРОСТИ И КОЭФФИЦИЕНТ 155 САМОДИФФУЗИИ В БОЛЬШИХ МОЛЕКУЛЯРНО-ДИНАМИЧЕСКИХ МОДЕЛЯХ...
автокорреляционной функции скорости 7(1) уравнениями (2) и (3), и в эти законы величина О не входит. Однако сейчас стало часто применяться так называемое уравнение Грина-Кубо (5), которое определяет О как интеграл от 7(1) в бесконечных пределах. На практике зачастую этот интеграл вычисляют до некоторого конечного времени. Ясно, что это также даёт приближённое значение О. Возникает вопрос о соотношении двух способов определения О. При больших 1 оба метода должны давать макроскопический коэффициент диффузии, и можно было бы ожидать, что разница между ними быстро стремится к нулю. Мы показали, что эта разница есть медленно спадающая функция, и различие составляет ~4% в аргоне и ~17% в воде от полной величины О, если интегрирование 7(1) обрывать при 2.5 пикосекунд.
Убывание коррелятора 7.(1) при больших имеет неэкспоненциальный вид. Предсказываемая гидродинамикой асимптота вида а Г3/2, действительно, может с какой-то точностью описать далёкий хвост. Однако значение амплитуды а при фитинге разных участков асимптоты получается разным. Поэтому нельзя сказать, что мы имеем ясное представление о форме асимптоты автокорреляционной функции скорости. Единственное, что можно надёжно утверждать, это что асимптота похожа на гиперболу а 1в с положительной амплитудой (а > 0). Это так в аргоне. А вода принципиально отличается от аргона: асимптота 7(1) в ней не выходит в положительную область даже при 1 ~ 10 пс, а в отрицательной области коррелятор, хотя и похож на гиперболу а 1в (с а < 0), но с показателем в = 5/2, а не 3/2. Итак, можно утверждать, что мы обнаружили ещё одно аномальное свойство воды, отличающее её от «обычных» простых жидкостей: асимптота автокорреляционной функции скорости в ней принципиально отличается
от асимптоты в аргоне.
Мы полагаем, что аномальные свойства воды обусловлены существованием в ней пространственной квази-тетраэдрической сетки водородных связей. Эта сетка ярко проявляется и в автокорреляторе скорости. В нем наличествуют модуляционные явления в области t < 0.5 пс (см. рис. 1), которые отражают заторможенные либрации и трансляции молекул в сетке водородных связей.
ЛИТЕРАТУРА
1. Egelstaff PA. An Introduction to the Liquid State. London, Academic Press, 1967.
2. Boon JP, Yip S. Molecular Hydrodynamics. Dover Publ. New York, 1991.
3. Hansen JP, McDonald IR. Theory of Simple Liquids. London, Academic Press, 2008.
4. Волошин ВП, Маленков ГГ, Наберухин ЮИ. Исследование коллективных движений в компьютерных моделях воды. Крупномасштабные и долговременные корреляции. Журн. структ. химии, 2013, 54, Приложение №2:S239-S257. [Voloshin VP, Malenkov GG, Naberukhin YuI. Collective motions in computer models of water. Large-scale and long-time correlations. J. Structur. Chem, 2013, 54, Supplement 2: S233-S251.]
5. Anikeenko AV, Malenkov GG, Naberukhin YuI. Visualization of the collective vortexlike motions in liquid argon and water: Molecular dynamics simulation. J. Chem. Phys., 2018, 148:094508. DOI: 10.1063/1.5018140.
6. Abascal JLF, Vega C. A general purpose model for the condensed phases of water: TIP4P/2005. J. Chem. Phys., 2005, 123:234505.
7. Alder BJ, Wainwright TE. Decay of the velocity autocorrelation function. Phys. Rev. A, 1970, 1:18-21.
8. Bellissima S, Neumann M, Guarini E, Bafile U, Barocchi F. Density of states and dynamical crossover in a dense fluid revealed by exponential mode analysis of the velocity autocorrelation function. Phys. Rev. E, 2017,
НАБЕРУХИН Ю.И., АНИКЕЕНКО А.В., ВОЛОШИН В.П.
НАНОСИСТЕМЫ
95:012108.
9. Kondratyuk ND, Norman GE, Stegailov VV Self-consistent molecular dynamics calculation of diffusion in higher n-alkanes. J. Chem. Phys., 2016, 145:204504.
10. Pomeau Y. Asymptotic properties of autocorrelation functions and the Enskog expansion in three-dimensional simple classical fluids. Phys. Rev. A, 1973, 7:11341147.
11. Hauge EH, Martin-Löf A. Fluctuating hydrodynamics and Brownian motion. J. Statist. Phys, 1973, 7:259-281.
12. Hermanns H-G, Felderhof BU. Long-time tails of translational and rotational Brownian motion in a suspension of hard spheres. J. Chem. Phys, 2007, 126:044902.
Наберухин Юрий Исаевич
д.х.н., проф.
Институт химической кинетики и горения им. В.В. Воеводского СО РАН Новосибирск 630090, Россия [email protected]
Аникеенко Алексей Владимирович к.ф.м.н.
Институт химической кинетики и горения им. В.В. Воеводского СО РАН
Новосибирск 630090, Россия [email protected] Волошин Владимир Павлович к.ф.-м.н, с.н.с.
Институт химической кинетики и горения им. В.В. Воеводского СО РАН
Новосибирск 630090, Россия [email protected].
