Топология и динамика сеток водородных связей в воде. Результаты компьютерного моделирования

УДК 540.1:532.7

ТОПОЛОГИЯ И ДИНАМИКА СЕТОК ВОДОРОДНЫХ СВЯЗЕЙ В ВОДЕ. РЕЗУЛЬТАТЫ КОМПЬЮТЕРНОГО МОДЕЛИРОВАНИЯ

В.Е. Петренко1, Н.А. Абакумова2, М.Л. Антипова1,

А.В. Боровков1, О.В. Ведь1, Д.Л. Гурина1

Институт химии растворов РАН, г. Иваново

Представлена членом редколлегии профессором А.Б. Килимником

Ключевые слова и фразы: водородные связи; межмолекулярные взаимодействия; метод молекулярной динамики; структурные и термодинамические свойства жидкой воды.

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

Свойства жидкой воды во многом определяются разветвленной сеткой водородных связей (ВС), поэтому в последнее время исследователи уделяют значительное внимание изучению сетки ВС как с топологической [1, 2], так и с динамической [3-5] точки зрения. Топология характеризует взаимное расположение и ориентацию молекул - иными словами, саму форму сетки ВС. Динамический подход рассматривает длительность происходящих в жидкости процессов: движения и переориентации молекул, разрушения и восстановления связей между ними. В компьютерном эксперименте динамического типа, каким является классическая молекулярная динамика, грань между этими точками зрения практически стирается, поэтому в настоящей работе применены оба подхода.

Молекулярно-динамический (МД) эксперимент проведен в каноническом КУТ-ансамбле, в кубической ячейке, содержащей 216 молекул, с использованием периодических граничных условий. Радиус обрезания межмолекулярных взаимодействий равен половине ребра ячейки. Температура и плотность моделируемой системы соответствуют нормальным условиям (Т = 300 К, Р = 0,1 МПа). Для интегрирования уравнений Ньютона применялась разностная схема Верле, а для поддержания постоянной температуры - термостат Берендсена [6]. Длительность моделирования после выхода на равновесие составляла от 200 пс до 1 нс с шагом по времени 2 фс.

На результаты компьютерного моделирования сетки ВС в жидкой воде оказывают весьма существенное влияние два момента: выбор потенциальной функции парного взаимодействия и выбор критерия, по которому оценивается существование ВС между молекулами.

В МД эксперименте пара молекул считается связанной, если она удовлетворяет заранее сформулированным требованиям (критериям образования ВС). Широко используемые определения ВС основаны либо на геометрическом, либо на энергетическом критерии [7]. Однако их категоричность позволяет получить не реальную сетку ВС, а ее упрощенную схему. Поэтому исследователями были разработаны альтернативные типы критериев образования ВС: комбинированный (представляющий собой сочетание энергетического и геометрического) [8, 9] и пространственный - «эффективный» [1]. Есть также мнение, что для окончательного признания существования ВС выполнения выбранного критерия в текущий момент недостаточно. Необходимо знать состояние рассматриваемой пары молекул в течение всего компьютерного эксперимента. На этом принципе основан динамический критерий образования ВС [10, 11].

Для иллюстрации зависимости результатов моделирования от выбранного критерия ВС мы выполнили расчет топологических характеристик жидкой воды: среднего числа ВС на молекулу <«Нв>, а также долей/„ «-связанных частиц в системе с использованием классического потенциала Т1Р4Р [12] для разных критериев и их пороговых значений (табл. 1). Как можно видеть, полученные результаты существенно зависят как от типа критерия, так и от выбора порогового значения энергии парного взаимодействия в случае энергетического и комбинированного критериев. При увеличении (по модулю) порогового значения энергии резко уменьшается число четырехкоординированных молекул за счет роста количества димеров и тримеров в системе. Между тем, широко используемый геометрический критерий (Коо £3,3 А, КоН £ 2,6А), несмотря на свою упрощенность и схематичность, показывает результаты, вполне согласующиеся с устоявшимися представлениями о тетраэдрической природе сетки ВС в воде, и, кроме того, позволяет интерпретировать поведение частиц в ближайшем окружении центральной молекулы, проследить их перемещение между координационными сферами.

Таблица 1

Доля связанных п водородными связями молекул в системе (/„), среднее число водородных связей на молекулу <пНВ> в зависимости от критерия водородной связи на примере модели ТІР4Р

Критерий f№ % f1, % f2, % f3, % f4, % f5, % f& % <пнв>

Г еометрический «ОО < 3,3 А; «ОИ < 2,6 А 0 4,3 12,5 29,6 44,2 9,4 0 3,4

Энергетический ЕИВ < -9,4 кДж/моль 2,0 3,1 13,2 32,5 40,7 8,5 0 3,3

Энергетический ЕИВ < -10,5 кДж/моль 0,5 6,2 18,3 29,4 37,6 8,0 0 3,2

Энергетический ЕИВ < -14,0 кДж/моль 0 12,1 28,3 33,7 18,6 7,3 0 2,8

Комбинированный ЕИВ < -9,4 кДж/моль; ЯОО < 3,3 А, ЯОИ < 2,6 А 1,0 4,6 15,9 30,9 40,1 7,5 0 3,3

Комбинированный ЕИВ < -14,0 кДж/моль; ЯОО < 3,3 А, ЯОИ < 2,6 А 0 11,8 30,1 34,2 17,2 6,7 0 2,8

Мы провели сравнение, как описывают топологию сетки ВС в воде известные парные потенциалы с различной геометрией модельной молекулы. В табл. 2 приведено среднее число ВС на молекулу, а также распределение молекул в системе по числу ВС для различных моделей [7]. Для корректности анализа во всех случаях использован одинаковый геометрический критерий: Лоо £ 3,3 А,

Лон £ 2,6 А. Как можно видеть, для всех моделей максимум на распределениях молекул по числу ВС соответствует четырехкоординированным молекулам, но самую высокую долю четырехсвязанных частиц в системе дает ТІР4Р, что свидетельствует о наиболее упорядоченной тетраэдрической структуре сетки ВС, несмотря на то, что значение <инв> для ТІР4Р мало отличается от величин, полученных для других потенциалов.

С нашей точки зрения, потенциал, претендующий на корректность описания структуры жидкой воды, должен быть способен воспроизвести такую важную структурную характеристику, как функции радиального распределения (ФРР). Разработке такого потенциала мы посвятили ряд работ [13-15], в результате получили модель, ФРР которой подтверждены экспериментально и удовлетворяют следующим требованиям [16]:

lmax 1max

gOmax < 3,0; rOoax ~ 0,45nm ; gOiT > 1; Pi = < 1; P2 = ^ > 1,

2 max gOH

2 max gHH

где Гц тах - положение к-го максимума на ФРР атомов у относительно атома і

к тах

(у = О или Н), - его высота.

Этот четырехточечный эффективный потенциал, которому мы дали название Т1Р4Р-НВ, в явном виде учитывает вклад неэлектростатических взаимодействий кислород - водород в энергию образования ВС

U (r) = UOO (r) + X Uf (r) + X U OH (r) = 4soo i,j OH

aoo

v rOO

12

X M-+X 14,93

и -

OH

U

OH

( „ oh А

12 („ oh А10

'OH

'OH

(aooА roo

6

+

а

а

+

Т аблица 2

Доля связанных п водородными связями молекул в системе

Модель /0, % /1, % /2, % /3, % /4, % /5, % /6, % <nHB>

SPC 0 4,0 18,1 29,9 43,2 4,8 0 3,27

TIP4P 0 4,3 12,5 29,6 44,2 9,4 0 3,41

BF 0 2,0 15,8 33,6 39,2 8,9 0,5 3,39

SRWK2 0 2,0 14,6 32,7 40,8 9,6 0,3 3,42

TIPS2 0 4,1 14,9 33,1 38,9 8,4 0,6 3,34

TIP4P-HB 0 3,5 10,6 33,2 46,5 6,0 0,2 3,42

Параметры модели, а также рассчитанные с ее использованием значения высот и положений пиков ФРР, термодинамические и динамические свойства приведены в табл. 3.

Результаты моделирования топологии сетки ВС в воде с новым потенциалом сопоставлены в табл. 2 с данными, полученными для известных моделей. Как и следовало ожидать, явный учет в парном потенциале вклада специфических взаимодействий OH в энергию образования ВС приводит к дополнительной стабилизации тетраэдрической координации молекул воды, что подтверждает увеличение количества четырехсвязанных молекул в системе.

С динамической точки зрения сетку водородных связей начали изучать лишь около двадцати лет назад. Здесь главными характеристиками ВС считают вероятность ее образования и время ее жизни.

В подавляющем большинстве случаев под временем жизни водородной связи понимают период, в течение которого в данной паре молекул соблюдается выбранный критерий существования ВС. Различают два типа времени жизни связи [17]. Непрерывное (continious) или короткое время длится до первого нарушения в данном димере выбранного критерия существования связи. Если данная связь разорвалась, а потом восстановилась, то она рассматривается как новая вне зависимости от длительности периода разрыва. Мерцающее (intermittent) или длинное время жизни водородной связи длится до тех пор, пока один из участников данной ВС не найдет себе партнера для образования новой связи. В этом случае промежуточные нарушения критерия существования связи не учитываются вне зависимости от их длительности.

Для расчета времени жизни ВС разработано несколько подходов, один из которых основан на исследовании изменения вероятности образования ВС с течением времени [18], другой - на анализе самих МД траекторий движения частиц [19].

Т аблица 3

Параметры потенциала TIP4P-HB и результаты моделирования

r0H, а ZHOH, град r , А 0M’ 4h, e D S00, кДж/моль °00, А \u 0H| min > кДж/моль r0H, А

0,9572 104,52 0,15 0,4415 1,85 0,6487 3,1536 10,46 1,40

^вн? кДж/моль rlmax А r00 , А 1 max g00 r 2max а r00 , А ax <N О P1 P2 -5 Dx10 , см2/с сУ, Дж/

-39,6 2,72 2,81 4,23 1,03 0,95 1,05 3,05 78,8

Примечание . гон - длина связи О-Н в молекуле воды; ^НОН - валентный угол; гом - расстояние между центром атома кислорода и смещенным по биссектрисе угла НОН избыточным отрицательным зарядом; дн - избыточный заряд на атоме водорода; т - дипольный момент молекулы; еоо и аоо - константы потенциала Леннард-Джонса;

|^ш1н | - абсолютная величина глубины минимума функции ^о™"61 (г) ; Г°н - абсцисса

точки, в которой функция ^оН1"61 (г) обращается в ноль; ивн - внутренняя энергия; В - коэффициент самодиффузии; сг - изохорная теплоемкость.

Наибольшую популярность в настоящее время приобрел метод, основанный на построении специфической автокорреляционной функции (АКФ) Снв(/) параметра существования ВС между молекулами i и у Sij(t) [5, 17], которая имеет вид

CHB(t) -

Sj

где параметр Sij(t) = 1, если ВС между молекулами i и у существовала в начальный момент времени, продолжает существовать в данный момент времени t и за истекший период не разрывалась более чем на заранее определенный интервал времени t; Si,(t) = 0 в противном случае. Промежуточные разрывы, если их продол-^ * * жительность не превышала t, игнорируются. В зависимости от значения t получают различные времена жизни ВС. Если взять t = 0, то результатом будет непре-

С *

рывное время жизни Тнв , если значение t не ограничивать - мерцающее время жизни тНв . Иными словами, при построении АКФ для расчета короткого времени жизни допускается только один переход параметра Sij от единицы к нулю, то есть разорванная связь больше не учитывается, даже если в ходе моделирования она восстанавливается. В случае длинного времени жизни учитывается восстановление любой ВС, вне зависимости от длительности периода ее отсутствия. Вероятность восстановления существовавшей связи очень высока, поэтому длинное время жизни ВС на порядок превышает короткое. При расчете длинного времени жизни существует вероятность (к счастью, очень малая) того, что рассматриваемая связь разорвется в начале моделирования, а восстановится в конце, и тогда временем ее жизни будет считаться весь период наблюдения за системой. Во избежание подобных неопределенностей на интервал промежуточного отсутствия ВС стремились наложить ограничения [5, 11]. В одной из работ [11] в качестве порогового было взято значение t = 0,08 пс, соответствующее точке перегиба на ранговом распределении времен жизни ВС. Другой вариант [5] - использовать * с

значение f = Тнв . Полученные с таким ограничением значения времени жизни

о

Тнв оказались близки к соответствующим значениям длинного времени жизни

тнв . Это показывает, что вероятность восстановления ВС, период отсутствия которой превысил короткое время жизни, очень мала.

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

Что же характеризует полученное время жизни? Для ответа на этот вопрос мы проанализировали распределения по времени жизни ВС для всех молекул системы. Они отражают, с какими именно частицами эти молекулы образуют ВС в каждый момент времени, а также роль, которую они играют: являются донором или акцептором протона для одной, двух или трех молекул. Кроме того, МД эксперимент позволяет проследить, как меняется с течением времени энергия каждой ВС. На рис. 1 в качестве примера приведены такие распределения для произвольно выбранной молекулы на протяжении 0,5 пс (расчет проведен с потенциалом Т1Р4Р-нв), которые отражают изменение характера ВС со временем (аналогичные распределения были получены для изменения энергии ВС). Можно видеть, что характер ВС с течением времени меняется незначительно и переходы

Таблица 4

Среднее время жизни водородной связи, полученное для различных потенциалов типа 12-6-1 и нового модельного потенциала ТГР4Р-ИБ с использованием геометрического критерия водородной связи

Модель Короткое время жизни тНв , пс Длинное время жизни тНв , пс

SPC 0,89 4,84

TIP4P 0,41 4,17

BF 0,38 4,10

SRWK2 0,22 5,08

TIPS2 0,21 5,45

TIP4P-HB 0,46 4,08

C

E

D

G

H

G,G

G,1

G,2

G,3

G,4

—I

G,5 t, пс

Рис. 1. Время существования водородных связей, образуемых одной из 216 молекул системы в течение 0,50 пс наблюдения с молекулами А, В, С, Б, Е, ґ, Є и Н первой координационной сферы. Изменение со временем типа связи - центральная молекула является: О - акцептором первого протона;  - акцептором второго протона;

А - донором первого протона; V - донором второго протона

молекулы из одного статуса в другой (когда, например, донор протона становится акцептором) относительно редки. Но, в это же время, энергия ВС существенно меняется, связь периодически становится то более сильной, то более слабой, хотя образующие ее молекулы сохраняют свои роли по отношению друг к другу. Усреднение таких распределений по числу связанных частиц позволило получить

визуальное представление о пребывании среднестатистической молекулы воды в водородносвязанном состоянии. Рис. 2 наглядно иллюстрирует изменение с течением времени взаимной ориентации молекул в первой координационной сфере, с которыми выбранная частица образует ВС. Разрыв ВС и, следовательно, замещение молекул в первой координационной сфере молекулами из объема происходит постепенно. Выяснено, что каждая частица может образовывать ВС не более чем с шестью молекулами одновременно. Была проведена оценка времени, необходимого для последовательного разрыва каждой из ВС (рис. 3). Если связи с 6-й и 5-й

V' ,,

Рис. 2. Мгновенные конфигурации, образованные водородными связями одной из 216 молекул системы с молекулами ближайшего окружения в течение 0,5 пс наблюдения: а) t = 0 пс; б) t = 0,05 пс; в) t = 0,11 пс; г) t = 0,20 пс; д) t = 0,30 пс; е) t = 0,50 пс

t, пс

Рис. 3. Распределение частиц по времени

частицами существуют непродолжительный промежуток времени (меньше

0,05 пс) и быстро разрушаются, то конфигурация с 4-мя ВС значительно более устойчива. При этом обращает на себя внимание тот факт, что третья частица быстро уходит вслед за четвертой, а одна из связей значительно прочнее всех остальных (продолжает существовать еще 0,7 пс после разрыва предпоследней). Полученные результаты полностью согласуются с рассчитанным нами средним значением короткого времени жизни ВС 0,462 пс и являются яркой иллюстрацией взаимной согласованности движения молекул в жидкой воде.

Сопоставив распределения частиц по временам жизни ВС и мгновенные конфигурации (см. рис. 1-3), мы пришли к выводу, что за интервал, равный короткому времени жизни ВС, в воде происходит полная перестройка ближайшего окружения данной молекулы, а за период, соответствующий длинному времени жизни, наблюдаются как входы частиц в первую координационную сферу выбранной молекулы, так и выходы из нее. Данные события можно объяснить совместным вращением ближайших соседей данной молекулы со своими координационными сферами. При таком типе движения изменяется не только взаимная ориентация молекул, но и расстояния между ними. То есть вращение координационных сфер как единого целого объединяет поступательное и вращательное движения и позволяет объяснить поступательную и вращательную диффузии (а значит и все связанные с ней свойства, например, вязкость жидкости) с микроскопической точки зрения.

Список литературы

1. Bergman, D.L. Topological properties of the hydrogen - bond network in liquid water / D.L. Bergman // J. Chem. Phys. - 2000. - Vol. 253, part 2/3. - P. 267-282.

2. Water in Biology, Chemistry and Physics. / W.G. Robinson and ets. // Singapore : World Scientific, 1996. - Р. 509.

3. Matsumoto, M. A new approach to the dynamics of hydrogen bond network in liquid water / M. Matsumoto, I. 0hmine // J. Chem. Phys. -1996. - Vol. 104, part 7. -P. 2705-2712.

4. Luzar, A. Resolving the hydrogen bond dynamics conundrum / A. Luzar // J. Chem. Phys. - 2000. - Vol. 113, part 23. - P. 10663-10675.

5. Guardia, E. Dynamics in hydrogen bonded liquids: water and alcohols / E. Guardia, J. Marti, J. A. Padro and ets. // J. Mol. Liq. - 2002, part 96-97. - P. 3-17.

6. Allen, M.P. Computer simulation of liquids / M.P. Allen, D.J. Tildesley. - London : Claren don Press, 1987. - 387 p.

7. Вода: структура, состояние, сольватация. Достижения последних лет (Water: Structure, State, Solvation. Achievements of the Last Years) / под ред. А.М. Куте-пова (A.M. Kutepov). - М. : Наука, 2003. - 403 c.

8. Sciortino, F. Hydrogen bond cooperativity in simulated water: Time dependence analysis of pair interactions / F. Sciortino, S.L. Fornili // J. Chem. Phys. -1989. -Vol. 90, part 5. - P. 2786-2792.

9. Kalinichev, A.G. Hydrogen bonding in supercritical water. 2. Computer Simulations / A.G. Kalinichev, J.D. Bass // J. Phys. Chem. A. - 1997. - Vol. 101, part 50. -P. 9720-9727.

10. Malenkov, G.G. Hydrogen bonds in computer - simulated water / G.G. Malenkov, D.L. Tytik, E.A. Zheligovskaya // J. Mol. Liq. - 1999. - Vol. 82. - P. 27-38.

11. Маленков, Г.Г. Динамический критерий водородной связи для анализа структуры водных кластеров / Г.Г. Маленков, Д.М. Тытик // Известия АН. - Сер. Физика. - 2000. - Т. 64, вып. 8. - С. 1469-1472.

12. Jorgensen, W.L. Comparison of simple potential functions for simulating liquid water / W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, and ets. // J. Chem. Phys. -1983. - Vol. 79, part. 2. - P. 926-935.

13. Петренко, В.Е. Вода в компьютерном эксперименте: противоречие в параметризации потенциалов / В.Е.Петренко, М.Л. Дубова, Ю.М. Кесслер // Журн. физ. химии. - 2000. - Т. 74, вып. 11. - С. 1957-1961 .

14. Антипова, М.Л. Параметры моделей и условные коэффициенты чувствительности функций радиального распределения воды по данным молекулярной динамики / М.Л. Антипова, В.Е. Петренко, Ю.М. Кесслер // Известия АН. Сер. Химия. - 2003. - Вып. 2. - С. 316-322.

15. Боровков, А.В. Влияние на структуру модельной жидкости учета водородных связей в потенциале парного взаимодействия вода - вода / А. В. Боровков, М.Л. Антипова, В.Е. Петренко // Журн. структ. химии. - 2004. - Т. 45, вып. 4. -

C. 678-682.

16. Soper, A.K. The radial distribution functions of water and ice from 220 to 673 K and pressures up to 400 MPa / A.K. Soper // Chem. Phys. - 2000. - Vol. 258. -P. 121-137.

17. Rapaport, D.C. Hydrogen bonds in water. Network organization and lifetimes /

D.C. Rapaport // Mol. Phys. - 1983. - Vol. 50, part 5. - P. 1151-1162.

18. Bertollini, D. Diffusion effects of hydrogen bond fluctuations. The long time regime of the translation and rotational diffusion of water / D. Bertollini, M. Cassettari., M. Ferrario and ets. // J. Chem. Phys. - 1989. - Vol. 91, part. 2. - P. 1179-1190.

19. Метод молекулярной динамики в физической химии (Method of molecular dynamics in physical chemistry) / под ред. Ю.К. Товбина. - М. : Наука, 1996. - 308 c.

Topology and Dynamics of Hydrogen Bond Nets in the Water. Computer Modeling Results

V.E. Petrenko1, N.A. Abakumova2, M.L. Antipova1, A.V. Borovkov1,

O.V. Ved1, D.L. Gurina1

Institute of Chemistry of Solutions RAS, Ivanovo (1);

Department of Chemistry, TSTU

Key words and phrases: hydrogen bonds; molecular interactions; molecular dynamics method; structural and thermo-dynamic properties of liquid water.

Abstract: Topological and dynamic characteristics of hydrogen bond nets in the water are calculated through the method of classical molecular dynamics using different model potentials. It is found out that taking into account specific molecular interactions of oxygen and hydrogen atoms leads to the increase in average number of hydrogen bonds per water molecule owing to extra stabilization of tetrahedrally coordinated molecules. It is shown that life cycle of hydrogen bonds characterizes correlation of motion of closely located molecules in coordinated spheres of central molecule.

Topologie und Dynamik der Netze der Wasserstoffverbindungen im Wasser. Ergebnisse der Computermodellierung

Zusammenfassung: Durch die Methode der klassischen Molekulardynamik sind die topologischen und dynamischen Charakteristiken der Netze der Wasserstoffverbindungen im Wasser mit der Nutzung der verschiedenen Modellpotentiale berech-net. Es ist erhalten, dass die offenbare Berucksichtigung der spezifischen Molekular-wechselwirkungen der Atome des Sauerstoffes und des Wasserstoffes zur Erhohung der mittleren Zahl der Wasserstoffverbindungen auf das Wassermolekul infolge der zusatz-lichen Stabilisierung der tetraedrisch koordinierten Molekule fuhrt. Es ist gezeigt, dass die Lebensdauer der Wasserstoffverbindungen die Koordination der Bewegungen der Molekule der nachsten Umgebung in den Koordinationsspharen des zentralen Molekuls charakterisiert.

Topologie et dynamique des reseaux des liens d’hydrogene dans l’eau. Resultats du modelage informatique

Resume: Par la methode de la dynamique moleculaire sont calculees les caracteristiques topologiques et dynamiquesl des reseaux des liens d’hydrogene dans l’eau avec l’emploi de differents potentiels du modelage. Est regu que l’enregistrement evident des relations intermoleculaires des atomes de l’oxygene et de l’hydrogene aboutit a l’elevation du nombre moyen des liens d’hydrogene sur une molecule de l’eau pour le compte de la stabilisation supplementaire des molecules coordonnees de la fagon tetraedrique. Est montre que les temps de vie des liens d’hydrogene caracterisent la correlation du mouvement des molecules de l’environnement le plus proche dans les spheres de coordination de la molecule centrale.