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

Ивановский государственный химико-технологический университет

На правах рукописи

БУШУЕВ Юрий Гениевич

СТРУКТУРНЫЕ СВОЙСТВА ЖИДКОСТЕЙ С РАЗЛ ИЧНЫМИ

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

ПО ДАННЫМ КОМПЬЮТЕРНОГО МОДЕЛ ИРОВ АНИЯ

02.00.04 – физическая химия Диссертация на соискание ученой степени доктора химических наук

Иваново 2001

ОГЛАВЛЕНИЕ

стр.

ВВЕДЕНИЕ 7

1. ПРИМ ЕНЕНИЕ МЕТОДА МОНТЕ-КАРЛО ДЛЯ МОДЕЛИРОВАНИЯ СТРУКТУРЫ ЖИДКОСТЕЙ

1.1. Общие теоретические положения 1.2. Алгоритм Метрополиса 1.3. Периодические граничные условия 1.4. Оценка эргодичности и сходимости результатов моделирования 1.5. Современные представления о межмолекулярных взаимодействиях 1.6. Во дородная связь 1.7. Методы аппроксимации межмолекулярных взаимодействий 1.8. Описание взаимодействий молекул воды 1.9. Экспериментальные сведения о функциях радиального распределения атомов в жидкой воде 1.10. Тестирование компьютерных моделей по экспериментальным данным 1.11. Полуфеноменологические модели строения воды 1.12. Изучение структуры воды компьютерными методами 1.13. Современные представления о структуре жидкого ацетона 1.14. Современные представления о структуре и свойствах жидких амидов 1.14.1. Формамид 1.14.2. N-метилформамид 1.14.3. N,N-диметилформамид 1.15. Структурные свойства жидкого метанола 1.16. Современные представления о структурных свойствах водных растворов 1.17. Некоторые результаты исследования водных растворов методами компьютерного моделирования 2. ЭКСПЕРИМЕНТАЛЬНАЯ ЧАСТЬ 2.1. Общие положения 2.2. Методические особенности моделирования индивидуальных жидкостей 2.3. Методика моделирования растворов Cs +, Хе I- в диметилфор- мамиде и метаноле.

2.4. Методика моделирования растворов атомов благородных газов в воде 2.5. Методика моделирования бинарных смесей 2.6. Оценка корректности расчетов 2.7. Методика определения структурных свойств жидкостей 2.8. Топологические свойства сеток связей 2.9. Перколяционные свойства сеток связей 3. ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ 3.1. Структурные свойства индивидуальных жидкостей. 3.1.1. Ацетон 3.1.1.1. Свойства ацетона на субмолекулярном и молекулярном структурных уровнях 3.1.1.2. Свойства надмолекулярной структуры жидкого ацетона. 3.1.1.3. Выводы 3.1.2. Диметилформамид 3.1.2.1. Свойства ДМФА на субмолекулярном и молекулярном структурных уровнях 3.1.2.2. Свойства надмолекулярной структуры жидкого ДМФА 3.1.2.3. Выводы 3.1.3. N-метилформамид 3.1.3.1. Свойства МФА на субмолекулярном и молекулярном структурных уровнях 3.1.3.2. Свойства надмолекулярной структуры жидкого МФА 3.1.3.3. Выводы 3.1.4. Метанол 3.1.4.1. Общие черты и специфические особенности структур ме- 3.1.4.2. Общие черты и специфические особенности структур ме- 3.1.5.1. Свойства ФА на субмолекулярном и молекулярном 3.1.6.1. Структурные свойства воды на субмолекулярном и 3.1.6.2. Свойства сетки Н-связей воды вблизи порога 3.1.6.3. Свойства сетки Н-связей воды в области средних энергий межмолекулярного взаимодействия 3.1.6.4. Дальние корреляции во взаимном расположении молекул 3.2.1. Структурные особенности гидратации гидрофобных частиц 3.2.2. Структурные свойства бинарных смесей вода – органический 3.2.2.2. Закономерности структурных преобразований на субмо- 3.2.2.3. Закономерности структурных преобразований на надмо- 3.2.3. Структурные свойства жидкой смеси формамида с диметилформамидом 3.2.3.1. Закономерности структурных преобразований на субмо- лекулярном и молекулярном уровнях 3.2.3.2. Закономерности структурных преобразований на надмо- 3.2.4 Структурные особенности сольватации частиц в N,N-диметилформамиде и метаноле 3.2.4.2. Закономерности структурных преобразований на субмо- лекулярном и молекулярном уровнях 3.2.4.3. Закономерности структурных преобразований на надмо-

ОСНОВНЫЕ РЕЗУЛЬТАТЫ И ВЫВОДЫ

1. Алгоритм и программа расчета молекулярных конфигураций жидкой воды методом Монте-Карло в NVT-ансамбле 2. Алгоритм и программа расчета числа замкнутых циклов и функций радиального распределения их геометрических центров в сетках Н-связей жидкой воды Метод нейтронной дифракции с изотопным Функции радиального распределения ФРР Изо хорный температурный дифференциал ИТД

ВВЕДЕНИЕ

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

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

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

Однако он имеет множество очевидных недостатков.

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

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

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

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

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

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

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

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

В качестве объектов исследования выбраны жидкости, состоящие из молекул разной формы, описываемой отталкивательной частью потенциала. В ряду жидкостей существенно варьируются и притягивательные взаимодействия. Жидкая вода обладает уникальными физико-химическими свойствами. Невозможно переоценить ее роль и значение в биологических и технологических процессах. Ацетон, метанол, амиды алифатических карбоновых кислот находят широкое практическое применение. Исследование взаимодействий и определение структурных свойств молекул формамида и N-метилформамида способствуют пониманию закономерностей структурной организации многих биологически значимых веществ.

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

Диссертационная работа выполнена в соответствии с координационным планом РАН по направлению «Химические науки и науки о материалах» (разделы 3.1, 3.12) и представляет собой часть исследований по одному из основных научных направлений Ивановского государственного химико-технологического университета "Термодинамика, строение растворов и кинетика жидкофазных реакций". На различных этапах работа получала финансовую поддержку Министерства образования (грант 1994 г.), шестого конкурса-экспертизы научных проектов молодых ученых РАН (грант 1999 г.), Российского фонда фундаментальных исследований (проекты №95-03-08426а, №96-03-33642а, №99-03-32413а, №99-03а) Цель работы состояла в установлении закономерностей структурной организации индивидуальных жидкостей, бинарных растворов и разбавленных растворов электролитов и неэлектролитов; в разработке методов анализа молекулярных конфигураций; в определении влияния отдельных составляющих межмолекулярных взаимодействий на свойства жидкостей на субмолекулярном, молекулярном и надмолекулярном структурных уровнях; в проверке многочисленных гипотез, объясняющих экспериментальные свойства жидкостей.

В связи с этим были определены основные задачи исследования:

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

- установить закономерности формирования надмолекулярной структуры, определить топологические свойства систем Н-связей жидкостей;

- определить наборы наиболее характерных молекулярных ассоциатов, выявить особенности их пространственного расположения;

- установить закономерности концентрационных изменений структурных, топологических и энергетических характеристик в бинарных растворах;

- исследовать влияние незаряженных и заряженных частиц растворенного вещества на структурные и энергетические характеристики индивидуальных растворителей;

- установить взаимосвязи между структурными свойствами жидкости и ее экспериментальными характеристиками.

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

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

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

VI Всесоюзная менделеевская дискуссия "Результаты экспериментов и их обсуждение на молекулярном уровне". Харьков. (1983). III Всесоюзное совещание "Проблемы сольватации и комплексообразования в растворах" Иваново. (1984), II Всесоюзная конференция "Химия и применение неводных растворов". Харьков.

(1989). IV Всесоюзное совещание «Проблемы сольватации и комплексообразования в растворах». Иваново (1989), VI Всесоюзная конференция "Термодинамика органических соединений". Минск. (1990), Х Менделеевская дискуссия "Периодический закон и свойства растворов" С.-Петербург. (1993), III Российская конференция "Химия и применение неводных растворов". Иваново. (1993), EuchemConference on statistical mechanical simulations of complex liquids. Lund, Sweden (1994), YI Conference on Calorimetry and Thermal Analysis. Zakopane, Poland (1994). VI Международная конференция "Проблемы сольватации и комплексообразования в растворах". Иваново. (1995), Международная конференция. "Теория и практика процессов сольватации и комплексообразования в смешанных растворителях", Красноярск (1996), I международная научно-техническая конференция "Актуальные проблемы химии и химической технологии", Иваново (1997). VII Международная конференция «Проблемы сольватации и комплексообразования в растворах». Иваново (1998), II Международная научно-техническая конференция "Актуальные проблемы химии и химической технологии "Химия-99" (весенняя сессия)" Иваново (1999), XIX Всероссийское Чугаевское совещание по химии комплексных соединений. Иваново (1999), Международная научная конференция "Жидкофазные системы и нелинейные процессы в химии и химической технологии". Иваново (1999), Международная научная конференция «Кинетика и механизм кристаллизации». Иваново. (2000), VIII Международная конференция «Проблемы сольватации и комплексообразования в растворах» Иваново (2001).

1. ПРИМЕНЕНИЕ МЕТОДА МОНТЕ-КАРЛО ДЛЯ

МОДЕЛИРОВАНИЯ СТРУКТУРЫ ЖИДКОСТЕЙ.

Метод Монте-Карло, применяемый в статистической физике, является частным случаем общего метода статистического моделирования, который используют для решения широкого круга задач в различных областях науки. В рамках метода Гиббса термодинамические характеристики вещества получают в результате усреднения по ансамблю, т.е. по совокупности очень большого числа идентичных по природе систем, находящихся в одинаковых внешних условиях и различающихся только по микросостоянию [1, 2, 3, 4, 5, 6]. Принято выделять три типа ансамблей, состояние которых задается тремя типами функций распределения.

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

Большинство физических систем не являются полностью изолированными.

Они могут обмениваться энергией и частицами с окружающей средой. При этом полагают, что рассматриваемая система мала по сравнению с окружающей ее системой, и любое изменение характеристик малой системы не сказывается на состоянии большой. Большая система действует как тепловой резервуар или тепловая баня с заданной абсолютной температурой Т. В большом каноническом ансамбле (T,V,µ) системы способны обмениваться и энергией, и частицами. Его состояние задается температурой Т, объемом V и химическим потенциалом µ. Расчет термодинамических характеристик, как правило, проводится в рамках канонического ансамбля (NVT, NPT). Каноническое распределение Гиббса – статистическое распределение для систем, содержащих заданное число частиц N, объем V (или давление P) и способных обмениваться энергией с окружением.

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

где k – постоянная Больцмана, Z – статистическая сумма, сумма по состояниям системы:

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

В классической статистической термодинамике микросостояние определяется заданием обобщенных координат q и обобщенных импульсов p. Для канонического ансамбля N частиц вероятность иметь значения импульсов в интервале (p, p+p), значения координат в интервале (q, q+q) определяется как:

Для классической системы статистическая сумма (1.2) заменяется статистическим интегралом fi - число степеней свободы молекул сорта i.

Здесь где m - масса молекулы. Функцию Гамильтона, отсчитываемую от нулевой энергии молекул, можно представить как сумму энергии внутренних молекулярных движений (электронные состояния, колебания, вращения и т.д.) Hвнутр, энергии поступательного движения центров масс и потенциальной энергии межмолекулярных взаимодействий.

С учетом этого статистическая сумма Q может быть представлена в виде:

где E0 – энергия молекулы в самом низком энергетическом состоянии, Qпост - статистическая сумма, связанная с поступательным движением молекулы, Qвнутр статистическая сумма, связанная с внутренними молекулярными движениями.

Подставляя выражение для Гамильтониана в (1.3) и, учитывая условие нормировки, получим выражение для статистического интеграла в виде где Zконф =... e kT d q1... d q N – конфигурационный интеграл. Если статистический интеграл известен, то для рассматриваемой системы можно найти все термодинамические величины; так свободная энергия системы равна С помощью известных термодинамических соотношений могут быть найдены давление, энтропия и химический потенциал системы:

Конфигурационный интеграл Zконф как функция температуры T, объема V и числа частиц N дает полную статистико-механическую информацию о системе и позволяет по формулам (1.7)-(1.9) рассчитать термодинамические свойства системы.

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

Особенностям применения метода Монте-Карло для моделирования жидкостей посвящено большое количество работ. Укажем лишь основные монографии, в которых излагаются теоретические положения и даются рекомендации к практическому применению метода [6, 7, 8, 9, 10].

Как отмечалось выше, в случае сложных жидкостей практически невозможно точно рассчитать значение конфигурационной статистической суммы системы.

Однако его можно оценить с помощью метода статистических испытаний, генерируя конечный набор молекулярных конфигураций и определяя вероятность их появления wi Если в некоторый фиксированный объем помещать случайным образом молекулы, энергия взаимодействия между которыми задается набором потенциальных функций, то в зависимости от конфигурации системы больцмановский множитель exp (-Ui /kT) может принимать различные значения. Некоторые конфигурации дают значительный вклад в канонические средние, а некоторые – практически нулевой (например, когда две частицы сближены настолько, что между ними имеется сильное отталкивание).

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

Среднее по ансамблю от любой физической величины M рассчитывают по формуле где i – номер конфигурации (среднее берется по всем рассмотренным конфигурациям системы). Поскольку усреднение (1.11.) проводят по конечному числу конфигураций m со смещенной выборкой, для исключения влияния смещения каждую конфигурацию необходимо брать с весом 1/i :

Метрополис с соавторами [7] предложил в качестве i взять распределение Больцмана В результате среднее значение любой физической величины M можно записать в виде Ансамбль, состоящий из m конфигураций, получают путем задания вероятностей перехода от одной конфигурации к другой. Вероятность перехода от i-й конфигурации к j-й pij считают зависящей от энергий этих конфигураций, а точнее от величины (Uj –Ui)/kT:

Таким образом, строят простые цепи Маркова, т.е. последовательности случайных событий, в которых вероятность определенного события зависит от исхода предыдущего испытания. В соответствии с условием микроскопической обратимости вероятности pij должны удовлетворять условиям:

Обычно полагают, что при где wi – вероятность появления некоторой конфигурации при случайном выборе (с использованием последовательности равномерно распределенных случайных чисел). Используя центральную предельную теорему теории вероятностей может быть доказано [8], что рассматриваемая цепь Маркова задает распределение, асимптотически стремящееся к каноническому.

На практике алгоритм Метрополиса реализуют следующим образом. Пусть заданы потенциал взаимодействия, конфигурация системы (начальное расположение частиц в элементарной ячейке моделирования) и температура Т. Рассчитывают потенциальную энергию системы Ui и вносят случайное изменение в конфигурацию (случайным образом выбирают k-ю частицу в ячейке и смещают ее).

При этом энергия системы становится равной Uj. Если Uj < Ui, то считают, что система перешла в новое состояние. Если Uj > Ui, то сравнивают величину exp[Uj –Ui )/kT ] со случайным числом (0,1). Если exp[ -(Uj –Ui )/kT ], то считают, что система перешла в j-е состояние. Если же > exp[-(Uj –Ui )/kT ], то переход в новое состояние не происходит, k-я частица сохраняет свои прежние координаты. При этом j-ую конфигурацию в цепи не учитывают, а рассматривают прежнее расположение частиц, соответствующее энергии Ui. Таким образом, чем больше значение энергии имеет система при случайном изменении конфигурации, тем с меньшей вероятностью она переходит в это состояние.

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

В результате генерирования цепей Маркова длиной в несколько миллионов конфигураций отбрасывают начальный неравновесный участок цепи, а на равновесном участке отбирают m статистически независимых молекулярных конфигураций, по которым рассчитывают средние значения физических величин. Очевидно, что расчет средних значений сопряжен с большим количеством вычислений, выполнение которых стало возможным только при появлении быстродействующих ЭВМ.

Одно из преимуществ метода Монте-Карло состоит в том, что алгоритм легко адаптировать к любому статистическому ансамблю. Например, при моделировании системы в NPT-ансамбле необходимо периодически изменять объем ячейки, а в уравнении (1.17) вместо разности потенциальных энергий использовать разность энтальпий:

где P – давление, V – объем ячейки, V – изменение объема ячейки.

Усредняя термодинамические функции по конфигурациям на равновесном участке цепи Маркова по уравнению (1.14), легко рассчитать конфигурационную энтальпию и мольный объем :

Для расчета других термодинамических функций – изобарической теплоемкости СP, изотермической сжимаемости k T и коэффициента температурного расширения, необходимо проводить моделирования системы при различающихся параметрах состояния, а затем находить конечные разности. Оценку указанных термодинамических функций можно сделать и по результатам одного моделирования, используя флуктуационные формулы:

Для решения большинства задач вполне достаточно проведения вычислений в каноническом ансамбле. Однако если требуется, то можно использовать и большой канонический ансамбль [11, 12].

Очевидно, что какими бы мощными не были компьютеры, невозможно решать уравнения для макро объемов жидкости, содержащих порядка 10 23 молекул.

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

Для минимизации влияния указанных эффектов используют периодические граничные условия. В элементарную ячейку моделирования, которую выбирают чаще всего в форме куба, помещают N частиц. Длину ребра ячейки L рассчитывают по экспериментальному значению плотности жидкости. Все • • • • • • бесконечное пространство заполняют аналогичными ячейками – образами основной ячейки. Так, на • • • • • • а в трехмерном пространстве – двадцать шесть. Процедура моделироРис. 1.1. Двумерная модель периодических граничных условий. Мо- вания по методу Монте-Карло училекула 5 является «образом» моле- тывает независимые смещения часкулы 5.

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

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

Согласно методу минимального образа расчет энергии взаимодействия молекул основной ячейки проводят следующим образом. Рассматривают только те молекулы и образы молекул, которые попадают в куб с центром на выделенной молекуле. Ребра куба параллельны ребрам основной ячейки. Затем рассчитывают энергию взаимодействия центральной молекулы с молекулами, попавшими в куб.

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

Данные методы расчета энергии вносят определенные погрешности, поскольку пренебрегают учетом дальнодействующей части потенциала взаимодействия. В случае электрически нейтральных молекул на больших расстояниях потенциал взаимодействия быстро стремится к нулю, и погрешности расчета не оказывают существенного влияния на свойства жидкости. Здесь достаточно часто используют метод реакционного поля [13, 14]. Однако при изучении сольватации ионов электростатические взаимодействия описываются медленно изменяющимся потенциалом, и ограничение области его действия может оказать существенное влияние на конечный результат. В этом случае, как правило, используют более сложный и медленный способ расчета энергии квазибесконечной системы – метод суммирования по Эвальду. Для уменьшения времени расчета применяют и другие методы, например, метод GENB (general expansion neighbor-box) имеет существенные преимущества перед методом суммирования по Эвальду [15].

Для выяснения степени влияния периодических граничных условий и способа учета дальних взаимодействий проводят дополнительные исследования изучаемой системы. Наиболее интенсивно методические аспекты компьютерного моделирования исследовались в 80-е годы [16, 17, 18, 19, 20].

1.4. Оценка эргодичности и сходимости результатов моделирования Суть проблемы эргодичности состоит в том, что не для всякой системы усреднение по цепи Маркова совпадает с реальным средним по ансамблю Гиббса.

Условие эргодичности может быть сформулировано следующим образом. Вероятность перехода за n шагов из любого возможного состояния i в любое возможное состояние j не равна нулю, т.е. из любого состояния принципиально возможно попасть во все остальные [3, 5]. При проведении расчетов методом МонтеКарло это означает, что любая цепь Маркова, независимо от выбора стартовой конфигурации, приведет в итоге к одному и тому же среднему по ансамблю. Это условие, однако, выполняется полностью только для систем, у ко торых по тенциалы взаимодействия нигде не обращаются в бесконечность. В случае же межмолекулярных взаимодействий потенциальная энергия обращается в бесконечность при сильном сближении атомов, т.е. существуют области конфигурационного пространства, в принципе недостижимые в процессе расчета. При больших плотностях, характерных для жидкостей, конфигурационное пространство разбивается на отдельные области, внутри которых условие эргодичности выполняется, однако переходы между ними маловероятны или невозможны. Существует проблема квазиэргодичности: даже если система в целом является эргодичной, возможно появление состояний, вероятность перехода между которыми очень мала, хотя и отлична от нуля.

Применяют различные способы решения проблем эргодичности и квазиэргодичности. Так в работах [21, 22] предлагается метод определения эргодичности путем проверки зависимости значений рассчитанных энергий от стартовых конфигураций. Поскольку этот путь ведет к резкому удорожанию расчетов, наряду с другими условиями рекомендуется использовать так называемый "горячий старт", т.е. проведение расчетов в начальной стадии эксперимента при очень высокой температуре.

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

Проблема сходимости тесно связана с эргодичностью метода. Если воспользоваться обычной формулой для математической дисперсии, то можно считать сходимость среднего значения f произвольной функции f (R ), достигнутой при малости стандартного отклонения [24]:

где М — число испытаний.

Для корректного использования (1.21) требуется, чтобы соседние f (R m ) были независимы (некоррелированы). Очевидно, что соседние конфигурации, генерируемые в методе Монте-Карло, не удовлетворяют этому требованию. Поэтоr му вместо отдельных значений f (R m ) используют «контрольные функции» — средние значения функций в заданном интервале шагов. Тогда М в выражении (1.21) — число контрольных функций, или, другими словами, число интервалов, на каждом из которых ведется независимое усреднение.

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

Размеры изучаемой системы, используемые потенциалы взаимодействия [26], схема генерации новых состояний, а также граничные условия [27] оказывают существенное влияние на сходимость результатов. Таким образом, пути улучшения сходимости могут быть различными.

Чаще всего уменьшения среднеквадратичного отклонения добиваются двумя способами: 1) увеличением числа сгенерированных конфигураций; 2) оптимизацией алгоритма получения новых конфигураций [28]. Первый путь требует больших затрат компьютерного времени, второй может оказаться более экономным.

1.5. Современные представления о межмолекулярных взаимодействиях Для определения энергии термодинамической системы необходимо уметь рассчитывать энергию межмолекулярных взаимодействий. Теории межмолекулярных взаимодействий посвящено большое количество монографий и статей [29, 30, 31, 32].

С формальной точки зрения все достаточно просто - для определения энергии минимальном уравнение Шредингера. Однако аналитическим Х-Ф предел ного круга простейших атомно-молекулярных сисс учетом конфиг.

Точная нерелятиным упрощениям. Уравнение Шредингера замевистская энергия Эксперимент няют уравнением Хартри – Фока, полагая, что обРис. 1.2 Зависимость щую волновую функцию можно представить в виэнергии молекулы от метода расчета. де произведения одноэлектронных спинорбиталей. Молекулярные орбитали записывают в виде линейной комбинации атомных орбиталей (МО-ЛКАО). В качестве базиса разложения используют базис атомных орбиталей, который заведомо неполон. В этом приближении уравнение Хартри – Фока переходит в систему уравнений Рутаана решаемую методом итераций.

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

1.2 схемы энергетических состояний, в рамках метода Хартри-Фока принципиально невозможно точно определить энергию молекулы.

Энергию взаимодействия двух молекул рассчитывают как разность между энергией комплекса и энергиями изолированных молекул. Однако ее величина крайне мала по сравнению с электронной энергией молекул, поэтому погрешность расчета существенно зависит от точности определения трех относительно больших величин. Дисперсионное взаимодействие возникает вследствие мгновенной корреляции флуктуаций электронной плотности молекул. Оно не учитывается в методе самосогласованного поля. Для учета энергии электронной корреляции проводят вычисления в широком базисе и применяют метод конфигурационного взаимодействия или теорию возмущений Меллера-Плессета [33]. Если энергию димера и молекул определять в различных базисах, то в результате суперпозиции базисных рядов возникнет погрешность (BSSE), для минимизации которой применяют метод уравнивания.

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

Альтернативный способ решения опирается на соотношения теории возмущений Релея-Шредингера. Она основана на разложении волновых функций и гамильтониана взаимодействий в ряд по малому параметру. Область ее применения ограничена относительно большими межмолекулярными расстояниями, слабыми взаимодействиями. Исторически сложилось так, что длительное время оценка энергии взаимодействия проводилась с привлечением теории возмущений. В первом порядке теории рассчитывается вклад от прямых электростатических взаимодействий. Здесь чаще всего используют разложение энергии в ряд по мультипольным моментам, что естественно вносит свои погрешности в величину энергии. Поправки к энергии взаимодействия во втором порядке теории носят название поляризационной и дисперсионной энергии. Поляризационный вклад обусловлен переходом одной из молекул в возбужденное состояние под воздействием поля другой молекулы, дисперсионный – переходом в возбужденное состояние обеих молекул.

Энергия взаимодействия дву х молекул является функцией пространственных и угловых переменных, U(R, 1, 2 ), где R(r,,) – вектор, соединяющий молекулярные системы координат, (,, ) - углы Эйлера, характеризующие ориентацию молекулярных систем координат. Энергия зависит от шести переменных, но расчет удобно проводить в лабораторной системе координат. Пространственные и угловые переменные можно разделить, записав U(R, 1, 2 ) в виде ряда [32]:

где суммирование ведется по повторяющимся индексам. Функция S предстает в виде ряда:

Здесь D – матрица Вигнера, С – сферические гармоники, коэффициенты КлебшаГордана записаны в матричном виде. Подобные разложения очень часто используют в квантовой механике, теории взаимодействий. Детальному описанию метода посвящено множество монографий и статей [4, 34, 35, 36, 37].

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

где rik – расстояние между атомами i и k, принадлежащими молекулам A, B.

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

В качестве тестового примера можно рассмотреть способы описания взаимодействий атомов Ar. Это наиболее простая задача. Тем не менее, только в конце 70-х годов был разработан модельный парный потенциал, воспроизводящий экспериментальные данные о свойствах аргона в газовой, жидкой и твердой фазах.

Он имеет следующую функциональную форму [32]:

где кроме большого набора подгоночных параметров используют дисперсионные константы C6, C8 и C10. Данное выражение слишком громозко и для расчета энергии взаимодействия более сложных молекул необходимо знать множество параметров, которые трудно определить для многих атомов или групп атомов. Поэтому используют упрощенные формулы.

Физически наиболее обосновано описывать отталкивание молекул экспоненциальной зависимостью, как в потенциале Букингема. Однако при низких температурах и давлениях, когда расстояния между атомами молекул не сильно отличаются от равновесных значений, степенная функция также неплохо аппроксимирует теоретическую кривую. В диполь-дипольном приближении дисперсионное взаимодействие имеет ассимптотику (1/r6 ). Поэтому наиболее широкое распространение для описания ван-дер-ваальсовых взаимодействий получил потенциал Леннард-Джонса (12-6):

где - глубина потенциальной ямы, - значение r, при котором энергия равна нулю, U(Rm ) = -. Выбор показателя степени равным 12 обусловлен математическим удобством.

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

где qi, qj – заряды атомов или специально выделенных центров, а rij.- расстояния между атомами (центрами взаимодействия) молекул. Электростатический вклад отражает специфические особенности взаимодействия молекул.

В атом-атомной схеме расчета энергии необходимо использовать параметры взаимодействия для каждой пары атомов. Возникает проблема их определения.

Количество параметров можно существенно уменьшить, если использовать комбинационные правила:

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

Водородную связь (Н-связь) можно рассматривать как частный случай специфического взаимодействия. Она образуется между атомом водорода, ковалентно связанным с атомом А в одной молекуле, и атомом В, принадлежащим той же или соседней молекуле. В качестве атомов А, В обычно выступают электроотрицательные атомы O, N, F, Cl, S. Выделяют внутри- и межмолекулярную Н-связь.

Более подробно рассмотрим физико-химические проявления межмолекулярных Н-связей, которые подразделяются на связи средней прочности, образованные нейтральными молекулами, и сильные ионные Н-связи.

Образование межмолекулярной Н-связи сопровождается сдвигом в сторону длинных волн и уширением полос поглощения в ИК- и комбинационных спектрах. Согласно данным спектроскопии ЯМР на протоне, участвующим в Н-связи всегда наблюдается уменьшение электронной плотности. Физико-химические свойства веществ, молекулы которых образуют межмолекулярные Н-связи, существенно отличаются от свойств остальных веществ. Температура плавления и кипения, теплота испарения, вязкость, диэлектрическая проницаемость, как правило, имеют более высокие значения в случае образования Н-связи.

В работе [31] предложена классификация по силе межмолекулярной Н-связи:

различают слабую (Е =0.44 кДж моль-1 ), среднюю (нейтральную) (Е = кДж моль-1 ) и сильную (ионную) (Е =80240 кДж моль-1 ) связь (Е – энергия разрыва Н-связи). Сильная Н-связь образуется при взаимодействии иона с молекулой, содержащей функциональные группы ОН, ОN, FН. Приведенная классификация достаточно условна, т.к. невозможно каким-либо физическим обоснованным способом найти границы деления. Со стороны слабых Н-связей отличия от ван-дер-ваальсовых комплексов становятся малосущественными. Со стороны сильных Н-связей характеристики комплексов близки величинам, наблюдаемым при образовании ковалентных химических связей.

Квантово-химические расчеты позволили выявить природу Н-связи. Показано, что с точки зрения энергетических характеристик комплексы с Н-связью не обладают никакими преимуществами по сравнению с другими молекулярными комплексами [31]. Относительные вклады различных типов взаимодействий в энергию комплексов с Н-связью примерно такие же, как и в энергию донорноакцепторных комплексов. Спецификой нейтральной Н-связи считают образование водородного мостика, содержащего умеренно полярную и сильную химическую связь. В отличие от нейтральной, при образовании сильной ионной Н-связи энергетические и структурные характеристики мономеров существенно изменяются. Результаты квантово-химического анализа комплексов с Н-связями показали, что особой уникальной природой Н-связи не обладают.

Авторы статьи [39] замечают, что невозможно дать точное определение Нсвязи и даже указать, какие взаимодействия, ковалентные или электростатические, играют основную роль при ее образовании. В некоторых случаях (бесконечные цепочки из молекул карбамида) основными являются электростатические и поляризационные взаимодействия, в других (енольная форма 1,3циклогександиона) – ковалентные взаимодействия. Обычно считают [31], что СН –группы не образуют Н-связей. Однако показано [40], что их взаимодействие с атомами О соседних молекул играет важную роль в формировании структуры жидкой муравьиной кислоты.

Взаимодействие молекул в жидкой фазе обладает той особенностью, что в рамках молекулярного ансамбля реализуются как относительно сильные, так и слабые Н-связи. При этом геометрические и энергетические характеристики Нсвязей могут существенно отличаться от величин, наблюдаемых в кристалле или в наиболее стабильных конфигурациях димеров. Как правило, они описываются плавными монотонными кривыми, и в рамках I-ансамбля невозможно однозначно установить наличие связи между двумя выделенными молекулами. Степень надежности определения существенно возрастает, если ввести динамический критерий Н-связи, рассматривая систему в рамках колебательно - усредненного Vансамбля. Однако это требует проведения дополнительных расчетов в рамках метода молекулярной динамики и не применимо к конфигурациям, полученным методом Монте-Карло. Наряду с линейными Н-связями в жидкости могут существовать бифуркатные [41, 42], трифуркатные и циклические связи. Они имеют разный статистический вес, но оказывают большое влияние на динамические и структурные характеристики жидкости.

Авторы работы [43], посвященной исследованию сверхкритической жидкой воды методами компьютерного моделирования, наряду с линейными Н-связями рассматривают циклические, бифуркатные и дважды согнутые (twofold). (см. рис.

1.3). Линейные связи доминируют при нормальных и сверхкритических условиях, но с ростом температуры доля остальных типов связей существенно повышается.

Рис. 1.3. Типы Н-связей, рассматриваемых в молекулярных конфигурациях воды. a – линейная, b – циклическая, c – бифуркатная, d – дважды согнутая.

Для преодоления затруднений, связанных с отсутствием физически обоснованного критерия Н-связи, принято исследовать зависимость свойств системы Нсвязей жидкости от параметров критерия [44]. Наиболее часто используют энергетический пороговый критерий Н-связи. Полагают, что между молекулами образовалась связь, если их энергия взаимодействия меньше заранее заданной фиксированной величины EHB. В некоторых случаях, например в амидах и сверхкритической воде, данного условия недостаточно, поскольку энергия электростатических взаимодействий может оказаться достаточно низкой и без образования Нсвязи. Здесь в определение критерия приходится вводить дополнительные геометрические ограничения, позволяющие разделить различные состояния молекул [45, 46, 47].

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

При исследовании жидкостей и растворов методами компьютерного моделирования необходимо проводить расчет потенциальной энергии взаимодействия 102 -104 частиц в десятках миллионов конфигураций. В данном случае абсолютно невозможно использовать квантово-химические методы определения энергии для группы молекул. Необходим упрощенный расчет энергий межмолекулярных взаимодействий. Его можно осуществить с использованием простых аналитических функций, аппроксимирующих истинную поверхность потенциальной энергии небольшой группы молекул (в простейшем парно-аддитивном варианте – димера). Для расчета термодинамических свойств растворов точность определения межмолекулярных потенциалов не менее важна, чем методические особенности реализации методов компьютерного моделирования.

Условно методы построения потенциалов межмолекулярных взаимодействий можно подразделить на три группы: эмпирические, неэмпирические (ab initio) и полуэмпирические. В пределах каждой группы потенциалы различают: по способу задания геометрии молекулы, по величине и характеру размещения эффективных зарядов, по способу учета внутренних степеней свободы (гибкие или жесткие молекулы) и поляризации молекул. Различаются потенциалы и выбором вида аналитической зависимости [21].

Эмпирические методы основаны на подборе параметров межмолекулярных потенциалов на основании известных свойств вещества. В наиболее строгом варианте - проводится уточнение параметров потенциала путем проведения компьютерного эксперимента и согласования полученных данных с экспериментальными свойствами. Таким способом получен потенциал Дьяконовой – Маленкова для моделирования воды [48] и потенциалы, используемые в расчетах методом молекулярной механики [38].

Неэмпирические методы основаны на построении потенциальной поверхности взаимодействия двух или более молекул с помощью квантово-химических расчетов и дальнейшей аппроксимации этой поверхности аналитическими функциями, используемыми в качестве потенциалов взаимодействия. Здесь следует отметить работы, выполненные под руководством Е. Clementi [49, 50] Полуэмпирические методы калибровки параметров потенциалов основаны на учете результатов теоретических расчетов и экспериментальных данных. Один из наиболее широко используемых наборов потенциальных функций – OPLS (Optimized Potentials for Liquid Simulations), - создан таким способом [46, 51, 52].

Недостатки неэмпирических методов связаны с необходимостью проведения высокоточных расчетов энергии взаимодействия для сотен и тысяч возможных способов взаимного расположения двух молекул. Это требует больших затрат компьютерного времени. Поверхность потенциальной энергии взаимодействующих молекул может иметь очень сложную форму, которую не всегда удается хорошо аппроксимировать аналитическими функциями. Приближение ищется численными методами в интегральной норме L2. В результате в отдельных областях пространства значения энергии, рассчитанные "точным" и приближенным методами, могут сильно различаться [53].

Если молекулы не являются сферически симметричными, энергия взаимодействия двух молекул является функцией шести переменных. Допустим, что пространство сканируют по каждой координате, разбивая числовую ось на частей. В результате получают миллион точек. Однако сетка все равно является слишком редкой. Действительно, шаг по пространственной координате составит ~100 пм, а по углу - 360 или 180. Недавно методом функционала плотности аппроксимировали формулой (1.16). Были предприняты попытки по возможности формулой (1.16). Были предприняты попытки по возможности наиболее точного определения потенциала для димера воды [55]. В результате он оказался настолько сложным, что его нельзя использовать в компьютерном моделировании. Набор http://fandango.ch.cam.ac.uk).

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

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

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

Среди множества существующих потенциалов следует отметить систему функций OPLS, аппроксимирующую взаимодействия углеводородов, амидов, пептидов, воды [46] и спиртов [51]. Она разработана на основе воспроизведения как результатов квантово-химических расчетов взаимодействия димеров молекул, так и физико-химических и термодинамических свойств жидкостей. Авторы системы добились выполнения условия переносимости параметров взаимодействия.

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

Однако необходимо отметить, что при любом способе построения потенциала функциональная зависимость лишь аппроксимирует реальную поверхность потенциальной энергии простыми функциями. В разных областях пространства точность аппроксимации различна. Следовательно, всегда существует большая вероятность, что при удовлетворительном описании ряда экспериментальных свойств жидкости, точность расчета других свойств окажется крайне низкой. Хорошее согласие расчетных и экспериментальных величин может быть обусловлено балансом между несколькими конкурирующими эффектами. Так, например, наиболее широко используемые потенциалы, воспроизводящие многие термодинамические и физико-химические свойства воды, не позволяют правильно рассчитать вибрационные спектры льда [57, 58]. Поэтому определение влияния различных составляющих потенциала взаимодействия на рассчитываемые характеристики жидкости является крайне актуальным.

Рассмотрим более подробно проблему описания взаимодействий на примере молекул воды. Молекула воды имеет относительно простое строение. Имеется большой объем экспериментальных данных о свойствах воды и водных систем.

При разработке потенциалов авторы выбирают воду в качестве тестовой модели.

1.8 Описание взаимодействий молекул воды Исследованию воды и водных растворов методами компьютерного моделирования посвящено огромное количество публикаций. Первые работы по моделированию воды [59, 60] стимулировали интенсивные исследования водных систем.

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

Существует несколько подходов к описанию межмолекулярных взаимодействий молекул Н2 О. Практически во всех случаях отделяют универсальную ван-дерваальсову составляющую потенциала от электростатической. Универсальные взаимодействия аппроксимируют простыми аналитическими функциями типа потенциалов Леннард-Джонса, Букингема или Морзе. Энергию электростатических взаимодействий рассчитывают разными методами.

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

Как правило, в первом случае требуется значительно меньше затрат компьютерного времени для расчета энергии молекулярной конфигурации. Но затраты времени будут еще меньше, если использовать один центр взаимодействия, а угловую зависимость энергии аппроксимировать специально подобранными функциями. Такие модели редко используют в расчетах и, как правило, в тех случаях, когда решение поставленной задачи сопряжено с необходимостью проведения большого объема вычислений. Среди моделей этого типа можно отметить модель 3D [61] и модель MB [62]. Несмотря на существенные упрощения, данные модели во многих аспектах правильно характеризуют термодинамические и структурные свойства воды и водных растворов.

На первом этапе работы по созданию потенциала необходимо задать геометрические характеристики молекулы воды. Согласно экспериментальным данным в газовой фазе длина связи О-Н составляет 95.72 пм, а валентный угол Н-О-Н – 104.520 [63], дипольный момент – 6.17 10-30 К м (1.85 Д) [64]. В жидкой фазе средняя длина связи О-Н увеличивается до 99.6 пм, а валентный угол немного уменьшается и составляет 101±50 [65].

Важные данные для проверки адекватности модельного потенциала и точности квантово-химических расчетов дают экспериментальные сведения о строении Рис. 1.4. Геометрические характеристики димера воды димера (Н2 О)2. В оптимальной конфигурации, схема которой приведена на рис.

1.4 расстояние между атомами О (ROO ) составляет 298±1 пм, угол = 57±100, угол = 6±100 [66]. Энергия взаимодействия молекул в димере равна –22.6±2. кДж моль-1 [67].

Ранее полагали [68], что кластеры, содержащие небольшое число молекул воды, образуют циклические структуры. В этом случае увеличивается число Нсвязей. После исследования тримеров [69], тетрамеров [70], пентамеров [71] и гексамеров [72] данное мнение в значительной степени было пересмотрено. Найдено, что с увеличением размера кластера расстояние между атомами О молекул, участвующих в Н-связи, сокращается. В пентамерах и гексамерах оно составляет 270 пм. Структура кластеров существенно изменяется при переходе от пентамеров к гексамерам. Пентамер имеет циклическую структуру, а гексамер является трехмерным объектом.

Недавно Е. Крячко провел интенсивные квантово-химические исследования структуры гексамеров [73]. Он выделил 15 структур с близкими значениями энтальпии (Н < 11.17 кДж моль-1 ). Минимальной энергией обладает структура в форме призмы, где девять Н-связей образуют треугольные основания и четырехугольные боковые грани. В большинстве структур можно выделить трех-, четырех- и пятичленные циклы. Шестичленные циклы имеют сложную пространственную форму, и большинство из них можно разбить на более простые циклы.

Структуры, характерные для гексагонального льда Ih, имеют высокую энтальпию относительно минимума: Н=7.11 («кресло»), 11.17 кДж моль-1 («ванна»).

Полную потенциальную энергию конфигурации из N молекул можно представить в виде ряда где V2, V3 …означают дву-, трех- … частичные взаимодействия, r – координаты молекулы. Основной вклад в энергию взаимодействия молекул воды дают первые две суммы ряда. Попытка прямой параметризации потенциальной энергии, основанная на использовании методов квантовой химии [74, 75], не привела к успеху.

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

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

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

Полную потенциальную энергию можно представить в виде где V/ eff Veff. Поляризационные взаимодействия рассчитывают с помощью иттерационной процедуры. Естественно, что такой способ аппроксимации энергии взаимодействия требует значительно больше затрат компьютерного времени по сравнению с парно аддитивной схемой расчета.

Энергетические, геометрические и электростатические характеристики димера воды Эмпирический (PPC) Veff + Vpol -20.9 280 53 19 2.6 0. В табл. 1.1 приведены экспериментальные энергетические [67] и геометрические [76] характеристики димера воды, и рассчитанные с помощью модельных потенциалов (SPC [77], PPC [78]) и с помощью квантовохимических методов [79].

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

Эффективные неполяризуемые потенциалы воспроизводят энергетические характеристики малых кластеров с большими погрешностями. В работе [80] на примере квантово-химического расчета малых циклических кластеров показано, что трехчастичные взаимодействия играют значительную роль при описании кооперативных эффектов, и они должны быть включены в схему расчета при моделировании структуры кластеров или жидкой воды. Эффекты от других многочастичных взаимодействий оказались незначительными. Однако, автор подчеркивает, что рассмотрены только простейшие циклические кластеры. Для окончательного вывода о роли многочастичных взаимодействий необходимо исследовать и другие структуры кластеров.

К аналогичным выводам приходят и авторы другого квантово-химического исследования [81] водных кластеров. Электростатические взаимодействия дают основной вклад в энергиию взаимодействия димера воды. Однако поляризация и перенос заряда являются доминирующими компонентами трехчастичной неаддитивной составляющей энергии тримера. Четырехчастичные взаимодействия оказались настолько малы, что ими можно пренебречь.

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

Молекула воды гибкая. Атомы находятся в постоянном движении, которое можно представить в виде суммы трех типов колебаний, обусловленных симметричным и антисимметричным растяжением связей О-Н, и изменением величины валентного угла Н-О-Н. Если в модельном потенциале, основанном на точечном распределении зарядов, допустить колебания заряженных центров, то дипольный момент молекулы будет линейно зависеть от степени растяжения связей. Однако Параметры эффективных потенциалов для моделирования воды Примечание: и – параметры потенциала Леннард-Джонса, расстояние от атома О до точки М, находящейся на биссектрисе угла Н-О-Н в направлении атомов Н.

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

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

Возникает вопрос - по какому закону должен изменяться дипольный момент молекулы при колебаниях атомов? Для корректной параметризации необходимо проводить квантово-химические расчеты.

Существует несколько причин для использования гибких моделей: исследование спектроскопических свойств водных систем, исследование переноса энергии между нормальными модами колебаний в растворах, исследование свойств водных систем в экстремальных условиях. При решении других задач жесткие модели более предпочтительны, поскольку они являются более простыми [82].

В литературе существуют противоречивые мнения по поводу влияния и роли гибкости молекул. Так, модель [83], основанная на потенциале SPC-ТR, по поведению атом-атомных ФРР оказалась очень близкой к исходной жесткой модели (SPC). В другом исследовании [84] показано, что допущение подвижности атомов в SPC потенциале приводит к существенному изменению свойств воды. Авторы работы [85] показали, что гибкость молекулы оказывает малое воздействие на свойства воды. Результаты МД моделирования [86] «гибкой» (CF) и «жесткой»

воды (RCF) указывают на то, что при нормальной температуре внутримолекулярная подвижность атомов оказывает крайне малое воздействие на физические свойства воды и растворов 1-1 электролитов. Напротив, диэлектрические свойства очень чувствительны к типу используемого потенциала. Время дебаевской релаксации при переходе от RCF к CF модели снизилось почти вдвое с 32 до 18 пс..

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

Гибкие модели постоянно совершенствуются. В работе [87] предложена простая гибкая модель F3C и исследованы ее структурные и динамические свойства.

В серии работ японских ученых [88, 89, 90] рассмотрены относительно старые гибкие модели SPC-TR, SPC-ZW, разработаны новые модели SPC-mTR, cm4P и определены их свойства. Основное внимание авторы уделили исследованию свойств воды в области суб- и надкритических состояний, где преимущества гибких моделей проявляются наиболее ярко. С введением в схему расчета внутримолекулярных степеней свободы возникает необходимость подбора дополнительных параметров. Для решения этой задачи используют данные ИКспектроскопии. О степени соответствия расчетных и экспериментальных значений частот внутримолекулярных колебаний можно судить по данным, приведенным в табл 1.3.

Частоты колебаний (см ) атомов молекул гибких моделей воды [88] SPC-TR SPC-ZW SPC-mTR Эксперим. Эксперим.

(2 ) Bending (1 ) Sym.str.

(3 ) Asym.str.

В настоящее время широко применяют поляризуемые потенциалы [91, 78, 92, 93, 94, 95, 96, 97], поскольку прогресс технических средств позволил за разумное время рассчитывать свойства водных систем, используя иттерационные процедуры. В тех случаях, когда задачи исследования связаны с определением состояния молекул в сильных электростатических полях (гидратация ионов, вода вблизи поверхности), с определением диэлектрических свойств растворов, поляризуемые потенциалы имеют несомненные преимущества перед неполяризуемыми.

Были созданы модели учитывающие одновременно и гибкость и поляризуемость молекулы воды. Некоторые из них основаны на известных жестких парно аддитивных моделях SPC-FP [98], MST-FP [99], NCC-vib [100]. В отдельных моделях молекула воды может диссоциировать [101].

Другой тип поляризуемых моделей основан на идее о внутримолекулярной флуктуации заряженных центров под воздействием внешнего электрического поля. Так в модели Sprik и Klein [102] четыре флуктуирующих заряда расположены в центрах взаимодействия модели молекулы TIP4P. Авторы ввели в схему расчета дополнительный лагранжиан взаимодействия, который описывает движение зарядов. Другой подход развит в работе [103]. В этой модели заряды связаны гармоническими силами с атомами. В модели PPC [78, 104] величина точечного заряда изменяется в зависимости от локального значения напряженности электрического поля.

В серии работ, выполненных под руководством A. Chialvo [105, 106, 107], проведено тестирование нескольких поляризуемых и гибких моделей. Авторы предложили новую поляризуемую модель в которой распределение электронной плотности задается гауссианами. Вместо простого точечного распределения, принятого в других моделях, здесь заряд «размазан» (smear) в окрестности центров взаимодействия. С помощью новой модели авторам удалось описать структурные свойства жидкой воды в широкой области параметров состояния.

Таким образом, в настоящее время существует множество потенциалов, аппроксимирующих взаимодействия молекул воды. При их разработке было использовано большое количество идей о роли и характере отдельных типов взаимодействий и внутримолекулярных движений. Каждая модель воды имеет свои преимущества и недостатки и предназначена для решения определенного круга задач. В рамках отдельных работ или серий работ определяют ограниченный набор параметров, характеристик жидкости. Более детально исследован относительно небольшой набор широко известных моделей (например, ST2, MCY, SPC/E, TIP4P). Задача, связанная с выявлением особенностей и общих структурных, динамических, диэлектрических, спектроскопических закономерностей различных моделей, несомненно является актуальной, а ее решение способствует пониманию механизмов обуславливающих специфических особенностей водных систем.

радиального распределения атомов в жидкой воде Практически во всех работах по компьютерному моделированию воды проводят сравнение расчетных и экспериментальных атом-атомных функций радиального распределения (ФРР). Решение задачи экспериментального определения ФРР имеет крайне важное значение для дальнейшего развития теории жидкого состояния. Однако даже в простых случаях не удается получить надежные сведения о поведении ФРР. Для лучшего понимания сути проблемы необходимо обратиться к основам теории рассеяния элементарных частиц.

Дифференциальное поперечное сечение нейтронов, рассеивающихся на ядрах атомов m, l, описывается выражением [108] где bi, bj – длины когерентного рассеяния, j – поперечное сечение некогерентного рассеяния. Дифференциальное сечение можно представить в виде суммы от рассеяния на одинаковых (l=m) и разных ядрах (lm). В случае D2 O первый вклад (self) записывают следующим образом где bO, bD – длины для когерентного рассеяния атомов O и D, D – сечение некогерентного рассеяния атомов D (O = 0). Второй вклад (distinct) можно представить в виде суммы от рассеяния на атомах одной молекулы (intra) и атомах разных молекул (inter). Структурный фактор SM(Q) записывают в виде где Q – волновой вектор.

Можно выделить две составляющие структурного фактора: f(Q) – молекулярный форм-фактор и DM(Q) – функцию, содержащую все межмолекулярные вклады.

где j0 – функция Бесселя, exp(-OD Q2 ), exp(-DD Q2 ) – дебаевские факторы, учитывающие изменения внутримолекулярных расстояний вследствие колебания атомов.

При больших значениях Q DM(Q) стремится к нулю. Здесь главный вклад в структурный фактор дает функция f(Q). Однако для жидкостей, в которых существует система межмолекулярных Н-связей, межмолекулярный вклад убывает очень медленно, поэтому его сложно определить и оценить.

Полную парную корреляционную функцию G(r) определяют с помощью фурье-преобразования структурного фактора где интеграл берется от 0 до, а асимптотическое значение для D2 O SM () = f(Q)= 0.3346, - плотность жидкости. Функция G(r) является суммой атоматомных ФРР и включает вклады от внутримолекулярных и межмолекулярных расстояний.

Если из функции SM (Q) вычесть молекулярный форм-фактор, то получим функцию DM(Q), фурье-преобразование которой дает возможность определить межмолекулярную парную корреляционную функцию Можно показать, что для D2 O [109] Следовательно, основной вклад в экспериментально определяемую функцию дают корреляции в пространственном положении атомов D относительно атомов D и O других молекул.

Рентгеновское излучение рассеивается на электронах, а не на ядрах атомов.

Определив угловую зависимость когерентного рассеяния можно определить парную корреляционную функцию GX(r), основной вклад в которую дает ФРР gOO (r).

Ее рассчитывают по соотношению где M(Q) –функция, демпфирующая осцилляции структурного фактора при больших Q. В этой формуле необходимо интегрировать по всем значениям волнового вектора Q. Однако в реальном эксперименте доступный интервал Q существенно ограничен. В связи с этим, на кривой GX(r) появляются ложные максимумы.

рассеяния сопряжена с необходимостью выделения отдельных составляющих, проведения большого количества расчетов и привлечения сведений о дополнительных константах и функциях. В жидкостей, молекулы которых состоят из небольшого числа атомов, невозможно надежно установить поведение атом-атомных ФРР. Существует большое количество методов уменьшения Рис. 1.5. Функции радиального распределения атомов функций [110, 111, 112, 113, 114, 115, 116, 117, О молекул воды, определеноднако устранить экспериментальные и ные по экспериментальным данным. 1 – [112], 2- [113], расчетные погрешности не удается. Например, 3- [114].

существенно отличаются ФРР gOO (r), определенные для воды разными авторами, или одной групгруппой авторов, но в разное время (см. рис. 1.5). Функции отличаются по высоте, положению и количеству пиков, что крайне затрудняет процедуру разработки модельных потенциалов.

Первоначально тестирование потенциалов проводили по данным NDIS- (Neutron Diffraction with Isotope Substitution) [112]. Затем были опубликованы результаты исследований жидкой воды, выполненных в широком интервале параметров состояний (NDIS-93 [119]). В литературе разгорелась дисскуссия [120, 121, 122, 123], поскольку новые экспериментальные данные существенно отличались от предыдущих. Естественно, что результаты многочисленных компьютерных моделирований стали не соответствовать экспериментальным данным [124].

К счастью, последующие исследования показали, что данные NDIS-93 являются артефактом. При обработке экспериментальных результатов для H2 O была сделана неоправданно большая поправка на неупругое рассеяние. Новый набор атоматомных ФРР (NDIS-97 [125]), принятый в качестве стандартного, существенно отличается от предыдущих. Амплитуды первых пиков gOO, gOH (NDIS-86) оказались завышенными. Новые результаты (NDIS-97) по функциям gOO в большей степени соответствуют данным Горбатого и Демьянца [115, 116] С помощью методов компьютерного моделирования получают набор молекулярных конфигураций, при этом известны координаты каждого атома. Достаточно просто можно рассчитать ФРР g(r), которая по определению равна отношению локальной численной плотности (r) точек в сферическом слое толщиной dr, находящихся на расстоянии r от фиксированной точки, к среднему значению плотности.

g(r) = (r)/ = [dN(r)/dV]/=(1/(4r2)) [dN(r)/dr], (1.40) где dN(r) - число точек в сферическом слое объемом dV.

Функции, полученные экспериментальными и теоретическими методами, можно непосредственно сравнивать. Однако неустранимые погрешности определения экспериментальной функции не позволяют делать однозначные выводы о степени адекватности модельного потенциала. Иногда сравнивают структурные факторы. Для этого с помощью фурье-преобразования ФРР переходят в пространство волновых векторов. В этом случае необходимо интегрировать по всем расстояниям 0 < r <. Однако в компьютерном эксперименте ФРР можно рассчитать только при r < L/2, где L –длина ребра кубической ячейки. В результате ограничения предела интегрирования структурный фактор будет определен с большими погрешностями.

В настоящее время все большее число исследователей при тестировании потенциалов для воды обращается к экспериментально определяемым изохорным температурным дифференциалам (ИТД). Если обратиться к формуле (1.40), то можно отметить, что в числителе содержится структурная информация (r), а в знаменателе стоит среднее значение плотности, характеризующее свойства жидкости на макроуровне. С изменением внешних условий плотность воды меняется, поэтому для получения информации о структурных свойствах воды некорректно непосредственно сравнивать ФРР, определенные при разных значениях плотности. Как известно, для температурной зависимости плотности жидкой воды характерно аномальное поведение. В окрестности максимума функции (T) можно выбрать две точки, отвечающие одному значению плотности, но разным температурам. В этом случае разность ФРР, gOO (r,T), будет характеризовать изменение закономерностей расположения молекул в сферических слоях, вызванное изменением температуры.

Легко заметить, что в этом случае уменьшаются систематические погрешности определения ФРР на основе экспериментальных данных. При расчете ИТД не обязательно вычитать ФРР. Можно первоначально рассчитать функцию DM (Q,T), а затем сделать фурье-преобразование по формуле (1.34). Такая процедура расчета позволяет обходиться без демпфирующих функций и существенно уменьшает ошибки, связанные с ограничением верхнего предела интегрирования [109].

Изохорные температурные дифференциалы были определены по данным о рассеянии рентгеновского излучения [113] и нейтронов [109]. На рис. 1.6 представлены ФРР и ИТД рассчитанные нами для модели воды с потенциалом Маленкова-Полтева (МП) и определенные по экспериментальным данным [113]. В парно аддитивном неполяризуемом эффективном потенциале МП задается строение жесткой молекулы [126, 127]. Первый пик gOO вдвое выше пика «экспериментальной» функции. Можно отметить, что за пределами первой координационной сферы расчетные и экспериментальные данные близки. Внутри сферы наблюдается качественное соответствие поведения кривых. Подобные результаты были получены и для модели воды (BSV), учитывающей поляризацию молекул [128].

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

Изохорные температурные дифференциалы, полученные по результатам исследования рассеяния нейтронов [109], отражают изменения функции G(r). Авторы работы отмечают несколько интересных закономерностей поведения ИТД.

При увеличении температуры на 45К наблюдаются достоверные изменения пространственной структуры воды на межмолекулярных расстояниях, достигающих Рис. 1.6. Функции радиального распределения атомов O (a) 1 – модель воды с потенциалом МП, 2- эксперимент [113], и изохорные температурные дифференциалы (b). 1 – расчет при T=55K, модель МП, 2- эксперимент при T=51K [113].

щих основной вклад в G(r), а следовательно и корреляции в положении этих атомов, наблюдаются на расстояниях до 600 пм. Сопоставление двух функций, полученных в результате деления ИТД на разность температур (G(r, T)/T) при T =45K и T =25.5К показало, что они практически идентичны за пределами области малых расстояний (r Рис. 3.53. Экспериментальная изохорная дифференциальная функция [113] g(r,51) (1) и определенная по уравнению (3.5) функция g(r), (2).

ру молекулы воды - 0.31, 0.63, 0.95, 1.25 нм. Напротив, для молекул с низкой энергией функция gL(r) осциллирует значительно сильнее и имеет максимумы при r 0.275, 0.45, 0.68, 0.91 и 1.1 нм, что соответствует межмолекулярным расстояниям в линейных цепочках тетраэдрически координированных молекул. Кривая g(r) для частиц с промежуточными значениями потенциальной энергии практически не отличается от кривой, рассчитанной с учетом всех молекул.

На рис. 3.53 представлена экспериментальная изохорная дифференциальная функция [113] g(r,51) и функция g(r), определенная по уравнению Несмотря на разную физическую природу явлений, описываемых этими функциями, наблюдается хорошая корреляция положений их экстремумов (существенное расхождение в величине дифференциальных эффектов при r < 0.6 нм связано с условностью процедуры разделения молекул на два подмножества). Это указывает на возможность интерпретации экспериментальных результатов с помощью привлечения представлений о существовании двух типов пространственной координации молекул в жидкой воде. Повышение температуры жидкости сдвигает динамическое равновесие в сторону увеличения доли молекул, обладающих более высокой энергией и менее выраженной тетраэдрической упорядоченностью [114, 115, 116].

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

Этот вывод соответствует совокупности представлений [157, 158, 159, 160], на которых основывается интерпретация многих экспериментально наблюдаемых свойств воды, и позволяет объяснить зависимость ФРР от температуры и давления.

Для определения состава и геометрических параметров подструктуры, образованной молекулами с низкой потенциальной энергией, необходимо получить сведения о кластерах - группах молекул, объединенных сильными Н-связями. Для этого мы использовали методы, разработанные в физике полимеров [8, 10, 218, 219, 361], классическую теорию самосогласованного поля и теорию перколяций.

Результаты компьютерного моделирования [44, 216, 217] и экспериментально определенное [379] аномальное скейлинговое поведение зависимости нейтронной вибрационной плотности состояний от энергии указывают на правомерность использования математического аппарата теории перколяции для описания свойств сетки связей воды. Однако результаты экспериментов [379] по малоугловому рассеянию рентгеновского (SAXS) и нейтронного (SANS) излучения не выявили структурной неоднородности и наличия дальних корреляций во взаимном расположении молекул воды. По мнению авторов [379], причинами отсутствия проявлений микрогетерогенности могут быть малые различия свойств областей с разными типами координации молекул, или, что более вероятно, динамический характер пространственной корреляции молекул в кластерах.

Ранее [44] расчет перколяционных параметров системы Н-связей жидкой воды проводили методом молекулярной динамики при 284 К для 216 молекул воды, взаимодействие которых описывали потенциалом ST2. Для выяснения влияния особенностей методики постановки компьютерного эксперимента на значения параметров необходимо сопоставить результаты расчетов, полученных с использованием различных моделей воды. Логично предположить, что увеличение количества частиц и конфигураций, по которым проводится усреднение, должно повысить надежность определения перколяционных параметров.

В молекулярных конфигурациях, рассчитанных методом Монте-Карло, мы выделяли кластеры и определяли зависимости их характеристик от вероятности образования связи p. Варьирование значения энергетического порогового критерия связи приводило к изменению величины p. Полученные зависимости свойств кластеров аппроксимировали функциями (2.10). Необходимо отметить сложности, возникающие при подборе оптимальных параметров методами нелинейной регрессии [380].

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

В табл. 3.11 приведены результаты аппроксимации зависимостей S() и P() Результаты аппроксимации по соотношениям (2.10) зависимостей среднего числа молекул в конечных кластерах (S ) и доли молекул в инфинитных кластерах (P) от параметра = (p-pc)/pc Для связей О-Н молекул воды характерно тангенциальное расположение по отношению к вектору, соединяющему молекулу с гидратируемой частицей. Только небольшая часть атомов Н может находиться ближе атомов О к неполярной Рис. 3.65. Функции радиального расобласти расстояний, соответствующих пределения атомов молекул воды в первой гидратной оболочке (1) и в положению первых пиков gАO. Это объеме (2) для растворов инертных газов а – He, b – Ar, c – Xe. Графики разнесены по оси ординат. объема». Замещая молекулу воды, атом численной плотности атомов О.

Отличия в поведении функций gOН проявляются за пределами первых пиков, отражающих закономерности расположения атомов соседних молекул воды образующих Н-связи (Н1..О2 на рис. 3.60). Второму пику функции отвечают расстояния между атомами О1…Н3, О1…Н4 и О2…Н2 представленными на рис.

3.60. Вследствие стерических ограничений, создаваемых атомами газов не реализуются некоторые ориентации молекул воды находящихся в ПГО.

Можно заключить, что в ПГО атомов инертных газов число Н-связей, образуемых молекулами воды, меняется незначительно, а изменение поведения ФРР объясняется понижением локальной численной плотности атомов Н и изменением распределения молекул воды по параметрам, характеризующим их взаимные ориентации.

На рис. 3.66 представлены графики функций W(E), характеризующих распределение молекул воды в зависимости от их полной энергии взаимодействия с Рис. 3.66. Функции W(E), характемами газов от расстояния до атомов, ризующие распределения молекул воды от полной энергии взаимодей- представленные на рис. 3.67.a, имеют вид ствия с другими молекулами воды.

1 – для ПГО частиц: a – He, b – Ar, c- Xe, 2 – для молекул воды в объе- Леннард-Джонса. На рис. 3.67.b предме раствора. ставлены зависимости полных потенциальных энергий молекул воды от расстояния до атома газа, являющиеся суммой энергий взаимодействия молекул воды друг с другом и с атомами газа. Погрешность их определения относительно высока, но можно заметить, что вблизи атомов инертных газов наблюдается небольшая энергетическая стабилизация части молекул, находящихся в ПГО. Высота энергетического барьера, отделяющего молекулы на ближних и дальних расстояниях составляет ~ 1 кДж моль-1, что не может заметно влиять на подвижность молекул воды в гидратных оболочках. Этот Рис. 3.67. Зависимость средней тицами (а) и зависимость полпм, угол О1..Н1..О2 > 1300 (см. рис. 3.60);

ной энергии молекул (b) от расстояния до атома: 1 – He, 2 и энергетический - EHB < -10.5 кДж моль-1.

– Ar, 3 – Xe.

Это позволило повысить достоверность выводов. Распределения молекул воды по числу образуемых ими Н-связей представлены на рис. 3.68. Они дают качественно одинаковую картину изменения свойств сетки связей вблизи гидратируемой частицы. Во всех случаях гидратации в ПГО частиц уменьшается доля молекул воды с бифуркатными связями (i=5).

Ограничение вращательной подвижности молекул ведет к уменьшению количества искаженных, «несовершенных» Н-связей. Доля молекул с четырьмя Нсвязями практически одинакова в ПГО и объеме растворителя. Наиболее заметно по числу образуемых ими Н-связей для молекул в ПГО неполярных частиц (1) средняя энергия образования Н-связи и молекул в объеме растворителя (2). a, b – He, c, d – Ar, e, f – Xe. a, c, e – энергетический, b, d, f - геометрический говорить о стабилизации сетки связей.

критерий связи Следовательно, классификация растворенных веществ, основанная на характере воздействия вещества на структуру растворителя (“structure maker”, “structure breaker”) абсолютно неверна, если не указано, что понимается под термином «структура», какие элементы и отношения между ними рассматриваются.

Полученные нами данные подтверждаются результатами независимых исследований, выполненных методами МД и МК для растворов метана и этана в воде [382-386]. Несмотря на то, что авторы использовали для моделирования другой ансамбль, метод расчета, потенциал взаимодействия, количество частиц в ячейке и способ учета дальних взаимодействий, они получили такую же картину процесТаблица 3. Среднее число Н-связей на молекулу воды (n) и средняя энергия Н-связи (E) для молекул воды находящихся в ПГО частиц и в объеме растворителя (bulk) сов, происходящих при растворении относительно небольших неполярных частиц в воде. Этот факт свидетельствует о достоверности выводов, сделанных на основании анализа расчетных характеристик растворов.

Результаты компьютерного моделирования гидратации трех неполярных частиц подтверждают основные положения феноменологических моделей гидрофобной гидратации (HSHB, LCW, решеточной модели). Неполярные частицы оказывают небольшое пространственно локализованное воздействие на структурные свойства воды. Характер изменения ФРР атомов воды в ПГО частиц объясняется «эффектом исключенного объема» - локальным понижением численной плотности атомов. Вид ФРР атомов воды относительно гидратируемых частиц определяется размером частицы. Кривые имеют сходную форму. На графиках ФРР наблюдаются широкие первые пики. Вторые пики практически отсутствуют.

Для молекул воды в ПГО свойственна преимущественно тангенциальная ориентация векторов дипольных моментов по отношению к поверхности неполярной частицы. Стерические препятствия, создаваемые частицей в растворе ограничивают вращательную подвижность молекул воды, что приводит к уменьшению числа искаженных, бифуркатных Н-связей, уменьшению среднего числа Н-связей и понижению их энергии. В зависимости от расстояния до частицы полная потенциальная энергия молекул воды осциллирует с небольшой амплитудой колебаний ~ 1 кДж моль-1. Вблизи частицы наблюдается энергетическая стабилизация молекул воды, приводящая к небольшому понижению их коэффициента диффузии, что характерно для положительно гидратируемых частиц.