CC BY
71
УДК 533.1
ОПРЕДЕЛЕНИЕ КОЭФФИЦИЕНТА САМОДИФФУЗИИ ВОДЫ В ПАКЕТЕ GROMACS
И.В. Булдашев, А.А. Мирзоев
Проведены молекулярно-динамическое моделирование температурных зависимостей коэффициента самодиффузии воды вблизи критической области при помощи пакета Gromacs с использованием трехточечного меж-частичного SPC-потенциала и сравнение результатов моделирования с опытными данными. Показано, что выбранная методика позволяет обеспечить хорошее согласие результатов моделирования с данными других авторов и эксперимента. Показано, что температурная зависимость коэффициента самодиффузии в критической области одинаково хорошо аппроксимируется как кривыми типа Аррениуса, так и уравнением Коэна-Турнбулла.
Ключевые слова: диффузия, вода, Gromacs, потенциал SPC, молекулярнодинамическое моделирование.
Вода имеет ключевое значение во многих физических процессах, поэтому изучение ее свойств является важной научной задачей на протяжении многих лет. И хотя к настоящему времени по данной тематике проведено множество исследований [1,2], некоторые вопросы все еще остаются без ответа [3]. В частности, до сих пор не закончено изучение структурных свойств воды в критической области, поэтому эта область не затрагивается во многих обзорах [4]. Целью настоящей работы является определение границ применимости потенциала SPC при моделировании поведения воды в широком температурном диапазоне критических состояний в пакете Gromacs (рис. 1). Исследуемая область
выбрана с учетом возможного применения в технике, например, паросиловых установках и ВВЭР (Водоводяной энергетический реактор).
Gromacs - многофункциональный пакет, использующий метод молекулярной динамики для систем, состоящих из большого числа частиц, в состав которого входит большая база полу-эмпирических парных межчастичных потенциалов для наиболее часто используемых веществ.
В работе проведено моделирование системы, состоящей из 216 молекул воды с использованием трехточечного потенциала, называемого в текущей литературе SPC [5], с периодическими граничными условиями. Модельный SPC-потенциал представляет молекулу воды в виде трех точечных зарядов: двух с зарядом q\ = 0,4le (атомы водорода) и одного с ¿72 = “0,8le (атом кислорода), для которого кроме электростатического потенциала добавлен потенциал типа Леннарда-Джонса (сг= 3,016 Á, энергия е- 15,319 кДж/моль). Приближенный характер SPC-потенциала несколько искажает равновесные параметры одиночной молекулы. В молекуле воды угол между О—Н связями и их длина составляют соответственно в = 104,52° и /] = 0,9584 Â против в = 109,47° и /= 1 Â в модели SPC.
Стартовая конфигурация для моделирования задавалась следующим образом: координаты атомов были равномерно распре-
300
ТОО 800 SG0 1000
400 500 600
Температура. К Рис. 1. Фазовая диаграмма воды, жирной линией отмечена исследуемая область
Рис. 2, Схема модели воды потенциала SPC
Булдашев И.В., Мирзоев А.А.
делены по кубической ячейке, а их начальные скорости инициализировались с помощью встроенного в пакет генератора случайных чисел, и удовлетворяли распределению Максвелла. Временной шаг моделирования был выбран равным 0,002 пс
[6], а так как число шагов 20 ООО, то полное время моделирования составило 40 пс
[7]. Тестовые прогоны показали, что для установления равновесия требуется не более 2 пс. Для моделирования системы в ЫРТ-ансамбле использовались термостат и баростат Берендсена [8], обеспечивающие сходимость температуры и давления системы к установленным значениям Т0 и Р0 в соответствии с динамическими уравнениями:
с1Т _Т0-Т с1Р _Р0-Р
Л т ’ Л тр
где тр, т - константы времени, являющиеся входными параметрами программы моделирования.
Моделирование проводилось при постоянном давлении и температуре (в т. н. ИРТ-ансамбле), в диапазоне изменений давления Р от 500 бар до 3000 бар и температуры Т от 300 К до 700 К. Число молекул N сохранялось постоянным и равным 216, что как показано в работе [9] является достаточным для получения результатов с погрешностью не хуже 5 %. Дополнительные моделирования с 820 молекулами подтвердили, что погрешность действительно лежит в этих пределах. По усредненным данным определялись величины объема, давления и температуры, а также парная корреляционная функция g(r). Отклонение средних величин давления и температуры от заданных баростатом и термостатом лежит в пределах 6 %. Коэффициент самодиффузии вычислялся по формуле
П /Л 2\
и =------, где (Аг ) - среднии квадрат
6т
полного смещения молекулы воды за время г . Использовался метод наименьших квадратов при его определении. Полученные при моделировании результаты сопоставлялись с экспериментальными данными работ [7, 10].
Сравнение результатов моделирования температурной зависимости коэффициента самодиффузии О воды с данными эксперимента (рис. 3) показало, что в области
Температура, К
Рис. 3. Сравнение результатов моделирования температурной зависимости коэффициента самодиффузии воды с данными эксперимента [10] при Р = 1000 бар
1/Т, К'1
Рис. 4. Температурные зависимость коэффициента диффузии воды в координатах
Температура, К
Рис. 5. Аппроксимация температурной зависимости коэффициента диффузии
температур 300-450 К расхождение находится в пределах расчетной погрешности, однако в области более высоких температур моделирование несколько занижает величину IX Существует альтернативный безакти-вационный подход к рассмотрению процессов диффузии в жидкости. Для многих из них зависимость коэффициента само-диффузии О в широком диапазоне температур описывается уравнением Аррениуса, которое при учете приложенного давления, имеет вид
0,45
0,40 -
0,35-
О - Д) ехр
Ев + РУ,
и
ЯТ
[11], где О0
«
а
0,25
постоянная слабо зависящая от темпера-
туры,
молярная энергия активации,
Я - универсальная газовая постоянная, V,
~ш— Т=700 К - Т=600 к Т=500 К
£>
так называемый активационныи
18
-----I----
20
22
24
—I—
26
—і—
28
32
34
Молярный объем, см /моль Рис. 6. Зависимость ОТ~°'1вз от молярного объема, построенная по данным МД-моделирования
объем. Зависимость 1п(£>) при высоком давлении, как видно из рис. 4, ведет себя прямо пропорционально 1/Т в диапазоне от 300 К до 700 К. Это позволяет нам выделить актива-
ционные параметры: ЕОР = -
д\п{Р) д 1п(1 / Т)
К. и V,
о
-КГ
д\п{Р)
дР
. Результаты расчета указан-
ных величин приведены в табл. 1. Для сравнения там же приведены определенные экспериментально значения энергии активации при постоянном объеме ЕВУ, которые определяются аналогично и согласно [10] слабо отличаются от величин . Результаты сравнения показывают, что найденные при моделировании энергии активации находятся в неплохом согласии с данными эксперимента. При увеличении температуры, при постоянной плотности, энергия активации убывает, а при увеличении плотности, при постоянной температуре, - увеличивается. Расхождение значений энергии активации при моделировании и эксперименте, не выходит за границы допустимой погрешности, и возникает в первую очередь из-за отличий в определении ЕОР и ЕГ)У . Активационный объем У0 меняется следующим образом: при 7’= 700 К - от 9,6 см3/моль при 2,5 кбар до 33,7 см3/моль при 0,95 кбар, при Г=600К - от 6,2 см3/моль при 2,5 кбар до 14,8 см3/моль при 0,95 кбар, что согласуется с экспериментальными данными[10].
Таблица 1
Моделирование Эксперимент
р, кг/м3 Г,К ЕОР, кДж/моль р, кг/м3 Г, К Еоу , кДж/моль
928 425 11,5 940 423 12,2
964 387 10,7 9,6 383 14,3
955 425 11,9 940 423 12,2
987 387 11,7 980 383 13,6
1016 350 11,4 1020 343 14,1
977 425 10,6 980 423 7,9
1007 387 10,1 1000 383 11,8
Существует еще одна точка зрения на природу диффузии, согласно которой основную роль играют не активационные процессы, а случайное перераспределение свободного пространства внутри жидкости [13]. При такой трактовке зависимость коэффициента диффузии от температуры и давления имеет вид [10, 13]
В(Р)
£> = А(Р)Т1/2 ехр
Т-Т0(Р)
Булдашев И.В., Мирзоев А.А.
Как видно из рис. 5, указанная зависимость также прекрасно описывает температурное поведение коэффициента диффузии воды в исследуемом диапазоне состояний, как и уравнение Аррениуса. Получаемые при такой подгонке значения коэффициентов А, В, Т0 приведены в табл. 2 для ряда значений давления.
Отметим также, что в работе [13] была теоретически выведена и подтверждена экспериментально зависимость И ос ТпУт, где Ут - молярный объем вдоль изотерм в сверхкритической области. Полученные нами результаты моделирования подтверждают указанную зависимость при значении п = 0,763 в диапазоне температур 500-700 К (рис. 6).
Таблица 2
Зависимость коэффициентов А,В,Т0 от давления
Давление, кбар 0,8 0,9 1,0 1,5 3
А 2273,3 147,9 94,0 27,0 9,8
В 11196,7 4833,0 4125,3 2487,8 1618,7
То -947,0 -460,7 -398,5 -227,6 -130,4
Таким образом, проведенное нами моделирование самодиффузии в широком температурном диапазоне вдоль нескольких изобар в сверхкритической области состояний, позволяет сделать следующие выводы:
а) потенциал SPC позволяет получать количественно точные значения коэффициента диффузии воды в диапазоне 300-500 К, однако в области более высоких температур приводит к заниженным оценкам указанной величины;
б) температурная зависимость коэффициента самодиффузии воды при различных значениях давления одинаково хорошо аппроксимируется как кривыми типа Аррениуса, так и уравнением Коэна-Турнбулла [13];
в) показано, что предложенная в работе [10] эмпирическая зависимость коэффициента самодиффузии D<x.TnVm, где Vm - молярный объем, при значении « = 0,763 вдоль изотерм в сверхкритической области, выполняется и в исследованном нами диапазоне температур 500-700 К.
Литература
1. Malenkov, G.G. Liquid water and ices: understanding the structure and physical properties /
G.G. Malenkov // J. Phys.: Condens. Matter. - 2009. - V. 21. - P. 283101(35 pp.).
2. Toukan, K. Molecular-dynamics study of atomic motions in water / K. Toukan, A. Rahman // Phys. Rev. B. - 1985. -V. 31. - P. 2643-2648.
3. Harrington, S. Liquid-Liquid Phase Transition: Evidence from Simulations / S. Harrington, R. Zhang, P.H. Poole et al. // Phys. Rev. Lett. - 1997. - V. 78. - P. 2409-2412.
4. Саркисов, Г.Н. Структурные модели воды / Г.Н. Саркисов // УФН. - 2006. - Т. 176. -Вып. 8.-С. 833-845.
5. Berendsen, H.J.C. Interaction models for water in relation to protein hydration /
H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren et al. // Intermolecular forces. - Dordrecht: D. Reidel Publishing Company, 1981. - P. 331-342.
6. Киселев, М.Г. Влияние отталкивательного взаимодействия на структурные и динамические особенности жидкой воды. Роль молекулярной поляризуемости / М.Г. Киселев, Ю.П. Пуховский, Д.В. Ивлев и др. // Журнал структурной химии. - 1999. - Т. 40. - Вып. 2. -С. 296-303.
7. Star, F.W. Dynamics of simulated water under pressure / F.W. Star, F. Sciortino, H.E. Stanley // Phys. Rev. E - 1999. - V. 60. - P. 6757-6768.
8. Berendsen, H.J.C. Molecular dynamics with coupling to an external bath / H.J.C. Berendsen, J.P.M. Postma, A. DiNola et al // J. of Chem. Phys. - 1984. - V. 81. - P. 3886-3892.
9. Spoel, D. A systematic study of water models for molecular simulation: Derivation of water models optimized for use with a reaction field / D. Spoel, P.J. van Maaren, H.J.C. Berendsen // J. of Chem. Phys.- 1998.-V. 108.-P. 10220-10229.
10. Krynicki, K. Pressure and temperature dependence of self-diffusion in water / K. Krynicki, C.D. Green, D.W. Sawyer // Faraday Discuss. Chem. Soc. - 1978. - V. 66. - P. 199-208.
11. Жирифалько, JI.А. Статистическая физика твердого тела / Jl.А. Жирифалько; пер. с англ. А.В. Ведяева, Ю.Г. Рудого. - М.: Мир, 1975. - 382 с.
12. Cohen, М.Н. Molecular Transport in Liquids and Glasses / M.H. Cohen, D. Turnbull I I J. of Chem. Phys. - 1959. - V. 31. - P. 1164-1169.
13. Lamb, W.J. Self-diffusion in compressed supercritical water / W.J. Lamb, G.A. Hoffman, J. Jonas // J. of Chem. Phys. - 1981. -V. 74. - P. 6875-6880.
Поступила в редакцию 28 декабря 2010 г.
WATER SELF-DIFFUSION COEFFICIENT CALCULATED BY GROMACS SOFTWARE
We performed the molecular dynamics simulation of the temperature dependence of self-diffusion coefficient of water in the supercritical area by means of Gromacs software using three-point interparticle SPC-potential. Simulation results were compared with experimental data. It is shown that the chosen method ensures good agreement of simulation results with those of other authors and experimental data. It is shown that the temperature dependence of the self-diffusion coefficient in the supercritical area is equally well approximated with Arrhenius type curves and the Cohen-Tumbull equation.
Keywords: diffusion, water, Gromacs, SPC potential, molecular dynamic simulation.
Buldashev Ivan Vladimirovich is student, Optics and Spectral Analysis Department, South Ural State University.
Булдашев Иван Владимирович - студент, кафедра оптики и спектроскопии, Южно-Уральский государственный университет.
e-mail: buldashev.ivan@mail.ru
Mirzoev Aleksandr Aminulaevich is Dr. Sc. (Physics and Mathematics), Professor, General and Theoretical Physics Department, South Ural State University.
Мирзоев Александр Аминулаевич - доктор физико-математических наук, профессор, кафедра общей и теоретической физики, Южно-Уральский государственный университет.
e-mail: mirzoev@physics.susu.ac.ru
