Рис. 1.1. Двумерная модель перио- дических граничных условий. Мо- лекула 5
’
является «образом» моле- кулы 5.
21 ния энергию взаимодействий полагают равной нулю (метод минимального об- раза, метод сферического ограничения) [6].
Согласно методу минимального образа расчет энергии взаимодействия мо- лекул основной ячейки проводят следующим образом. Рассматривают только те молекулы и образы молекул, которые попадают в куб с центром на выделенной молекуле. Ребра куба параллельны ребрам основной ячейки. Затем рассчитывают энергию взаимодействия центральной молекулы с молекулами, попавшими в куб.
В методе сферического ограничения области действия потенциала при опре- делении энергии взаимодействия молекулы учитывают только частицы, попав- шие в сферу определенного радиуса. Как правило, радиус сферы полагают рав- ным половине длины ребра элементарной ячейки.
Данные методы расчета энергии вносят определенные погрешности, по- скольку пренебрегают учетом дальнодействующей части потенциала взаимодей- ствия. В случае электрически нейтральных молекул на больших расстояниях по- тенциал взаимодействия быстро стремится к нулю, и погрешности расчета не ока- зывают существенного влияния на свойства жидкости. Здесь достаточно часто используют метод реакционного поля [13, 14]. Однако при изучении сольватации ионов электростатические взаимодействия описываются медленно изменяющим- ся потенциалом, и ограничение области его действия может оказать существенное влияние на конечный результат. В этом случае, как правило, используют более сложный и медленный способ расчета энергии квазибесконечной системы – ме- тод суммирования по Эвальду. Для уменьшения времени расчета применяют и другие методы, например, метод GENB (general expansion neighbor-box) имеет существенные преимущества перед методом суммирования по Эвальду [15].
Для выяснения степени влияния периодических граничных условий и спо- соба учета дальних взаимодействий проводят дополнительные исследования изу- чаемой системы. Наиболее интенсивно методические аспекты компьютерного моделирования исследовались в 80-е годы [16, 17, 18, 19, 20].
1.4. Оценка эргодичности и сходимости результатов моделирования
Суть проблемы эргодичности состоит в том, что не для всякой системы ус- реднение по цепи Маркова совпадает с реальным средним по ансамблю Гиббса.
22
Условие эргодичности может быть сформулировано следующим образом. Веро- ятность перехода за n шагов из любого возможного состояния i в любое возмож- ное состояние j не равна нулю, т.е. из любого состояния принципиально возмож- но попасть во все остальные [3, 5]. При проведении расчетов методом Монте-
Карло это означает, что любая цепь Маркова, независимо от выбора стартовой конфигурации, приведет в итоге к одному и тому же среднему по ансамблю. Это условие, однако, выполняется полностью только для систем, у которых потенциа- лы взаимодействия нигде не обращаются в бесконечность. В случае же межмоле- кулярных взаимодействий потенциальная энергия обращается в бесконечность при сильном сближении атомов, т.е. существуют области конфигурационного пространства, в принципе недостижимые в процессе расчета. При больших плот- ностях, характерных для жидкостей, конфигурационное пространство разбивает- ся на отдельные области, внутри которых условие эргодичности выполняется, од- нако переходы между ними маловероятны или невозможны. Существует пробле- ма квазиэргодичности: даже если система в целом является эргодичной, возмож- но появление состояний, вероятность перехода между которыми очень мала, хотя и отлична от нуля.
Применяют различные способы решения проблем эргодичности и квазиэр- годичности. Так в работах [21, 22] предлагается метод определения эргодичности путем проверки зависимости значений рассчитанных энергий от стартовых кон- фигураций. Поскольку этот путь ведет к резкому удорожанию расчетов, наряду с другими условиями рекомендуется использовать так называемый "горячий старт", т.е. проведение расчетов в начальной стадии эксперимента при очень вы- сокой температуре.
Проблема квазиэргодичности в принципе может быть решена путем удлине- ния цепи Маркова. В частности проблему эргодичности рекомендуют решать удачным выбором граничных условий. В одной из наших работ [23] исследована зависимость структурных свойств модели воды от стартовой конфигурации и длины цепи Маркова.
Проблема сходимости тесно связана с эргодичностью метода. Если восполь- зоваться обычной формулой для математической дисперсии, то можно считать
23 сходимость среднего значения
f
произвольной функции
f R
( )
r
, достигнутой при малости стандартного отклонения [24]:
( )
[
]
2 1
)
1
(
1
∑
=
2
−
−
=
M
m
m
f
R
f
M
M
r
ο
, (1.21) где М — число испытаний.
Для корректного использования (1.21) требуется, чтобы соседние
f R
m
(
)
r были независимы (некоррелированы). Очевидно, что соседние конфигурации, ге- нерируемые в методе Монте-Карло, не удовлетворяют этому требованию. Поэто- му вместо отдельных значений
f R
m
(
)
r используют «контрольные функции» — средние значения функций в заданном интервале шагов. Тогда М в выражении
(1.21) — число контрольных функций, или, другими словами, число интервалов, на каждом из которых ведется независимое усреднение.
Выбор необходимой длины интервала усреднения и числа таких интервалов в первую очередь определяется спецификой изучаемой жидкости, а также рассчи- тываемой характеристикой. Согласно [25], число контрольных функций должно быть не менее пяти.
Размеры изучаемой системы, используемые потенциалы взаимодействия
[26], схема генерации новых состояний, а также граничные условия [27] оказы- вают существенное влияние на сходимость результатов. Таким образом, пути улучшения сходимости могут быть различными.
Чаще всего уменьшения среднеквадратичного отклонения
σ добиваются двумя способами: 1) увеличением числа сгенерированных конфигураций; 2) оп- тимизацией алгоритма получения новых конфигураций [28]. Первый путь требует больших затрат компьютерного времени, второй может оказаться более эконом- ным.
1.5. Современные представления о межмолекулярных взаимодействиях
Для определения энергии термодинамической системы необходимо уметь рассчитывать энергию межмолекулярных взаимодействий. Теории межмолеку- лярных взаимодействий посвящено большое количество монографий и статей [29,
30, 31, 32].
24
С формальной точки зрения все достаточно просто - для определения энергии взаимодействия молекул необходимо решить уравнение Шредингера. Однако аналитическим способом решение найдено только для ограничен- ного круга простейших атомно-молекулярных сис- тем. Поэтому приходится прибегать к значитель- ным упрощениям. Уравнение Шредингера заме- няют уравнением Хартри – Фока, полагая, что об- щую волновую функцию можно представить в ви- де произведения одноэлектронных спин- орбиталей. Молекулярные орбитали записывают в виде линейной комбинации атомных орбиталей (МО-ЛКАО). В качестве базиса разложения используют базис атомных орбиталей, который заведомо неполон. В этом приближении уравнение
Хартри – Фока переходит в систему уравнений Рутаана решаемую методом ите- раций.
В методе Хартри-Фока полагают, что электроны с разными спинами движутся независимо друг от друга. Величина разности между экспериментальным и тео- ретическим значением энергии существенно зависит от метода расчета, базиса, способа учета электронной корреляции. Как следует из представленной на рис.
1.2 схемы энергетических состояний, в рамках метода Хартри-Фока принципи- ально невозможно точно определить энергию молекулы.
Энергию взаимодействия двух молекул рассчитывают как разность между энергией комплекса и энергиями изолированных молекул. Однако ее величина крайне мала по сравнению с электронной энергией молекул, поэтому погреш- ность расчета существенно зависит от точности определения трех относительно больших величин. Дисперсионное взаимодействие возникает вследствие мгно- венной корреляции флуктуаций электронной плотности молекул. Оно не учиты- вается в методе самосогласованного поля. Для учета энергии электронной корре- ляции проводят вычисления в широком базисе и применяют метод конфигураци- онного взаимодействия или теорию возмущений Меллера-Плессета [33]. Если энергию димера и молекул определять в различных базисах, то в результате су-
E Х-Ф расчет в минимальном базисе в расширенном базисе
Х-Ф предел с учетом конфиг. взаимодейств ия
Точная нереляти- вистская энергия
Эксперимент
Рис. 1.2
Зависимость энергии молекулы от метода расчета .
25 перпозиции базисных рядов возникнет погрешность (BSSE), для минимизации которой применяют метод уравнивания.
Действуя таким способом, в принципе, можно достаточно точно определить энергию взаимодействия, распределение зарядов, геометрические характеристики и другие свойства димера. Однако технические затруднения, возникающие при расчете сложных молекул, существенно ограничивают круг «точно» решаемых задач.
Альтернативный способ решения опирается на соотношения теории возму- щений Релея-Шредингера. Она основана на разложении волновых функций и га- мильтониана взаимодействий в ряд по малому параметру. Область ее применения ограничена относительно большими межмолекулярными расстояниями, слабыми взаимодействиями. Исторически сложилось так, что длительное время оценка энергии взаимодействия проводилась с привлечением теории возмущений. В пер- вом порядке теории рассчитывается вклад от прямых электростатических взаимо- действий. Здесь чаще всего используют разложение энергии в ряд по мульти- польным моментам, что естественно вносит свои погрешности в величину энер- гии. Поправки к энергии взаимодействия во втором порядке теории носят назва- ние поляризационной и дисперсионной энергии. Поляризационный вклад обу- словлен переходом одной из молекул в возбужденное состояние под воздействи- ем поля другой молекулы, дисперсионный – переходом в возбужденное состоя- ние обеих молекул.
Энергия взаимодействия двух молекул является функцией пространственных и угловых переменных, U(R,
Ω
1
,
Ω
2
), где R(r,
θ,ϕ) – вектор, соединяющий молеку- лярные системы координат,
Ω(α, β, γ) - углы Эйлера, характеризующие ориента- цию молекулярных систем координат. Энергия зависит от шести переменных, но расчет удобно проводить в лабораторной системе координат. Пространственные и угловые переменные можно разделить, записав U(R,
Ω
1
,
Ω
2
) в виде ряда [32]:
∑
Ω
Ω
=
Ω
Ω
)
,
и
,
,
(
)
(
)
,
и
,
,
,
(
2 1
2 1
2 1
2 1
2 1
2 1
ϕ
ϕ
k
k
j
l
l
k
k
j
l
l
S
r
f
r
U
, (1.22) где суммирование ведется по повторяющимся индексам. Функция S предстает в виде ряда:
26
)
)(
,
(
)]
,
,
(
[
)]
,
,
(
[
)
000
(
2 1
2 1
*
2 2
2
*
1 1
1 1
2 1
2 2
2 1
1 1
2 1
2 1
2 1
m
m
m
j
l
l
С
D
D
j
l
l
i
S
jm
l
k
m
l
k
m
j
l
l
k
k
j
l
l
ϕ
θ
γ
β
α
γ
β
α
×
×
=
∑
−
−
−
(1.23)
Здесь D – матрица Вигнера, С – сферические гармоники, коэффициенты Клебша-
Гордана записаны в матричном виде. Подобные разложения очень часто исполь- зуют в квантовой механике, теории взаимодействий. Детальному описанию мето- да посвящено множество монографий и статей [4, 34, 35, 36, 37].
Метод, основанный на разложении энергии межмолекулярных взаимодейст- вий в ряд по бимолекулярным базисным функциям, не получил широкого распро- странения при компьютерном моделировании жидкостей вследствие сложности расчета. Значительно проще определить энергию взаимодействия, используя атом-атомную схему разложения потенциала.
)
(
)
,
(
,
ik
B
k
A
i
ik
r
U
R
U
∑
∈
∈
=
Ω
, (1.24) где r
ik
– расстояние между атомами i и k, принадлежащими молекулам A, B.
В настоящее время для описания взаимодействий используют модельные по- тенциалы, состоящие, как правило, из двух слагаемых. Первое имеет универсаль- ный характер и определяется отталкиванием атомов на малых расстояниях за счет обменного взаимодействия и дисперсионным притяжением, доминирующим в интервале средних и дальних расстояний. Для описания этих универсальных ван- дер-ваальсовых взаимодействий применяют различные виды функциональной за- висимости.
В качестве тестового примера можно рассмотреть способы описания взаимо- действий атомов Ar. Это наиболее простая задача. Тем не менее, только в конце
70-х годов был разработан модельный парный потенциал, воспроизводящий экс- периментальные данные о свойствах аргона в газовой, жидкой и твердой фазах.
Он имеет следующую функциональную форму [32]:
∑
∑
=
+
+
=
+
+
−
−
=
2 0
6 2
*
6 2
5 0
)
)
/
(
(
1 1
exp
)
(
j
j
m
j
i
m
i
i
m
R
r
C
R
r
A
R
r
r
U
δ
α
ε
, (1.25)
27 где кроме большого набора подгоночных параметров используют дисперсионные константы C
6
, C
8
и C
10
. Данное выражение слишком громозко и для расчета энер- гии взаимодействия более сложных молекул необходимо знать множество пара- метров, которые трудно определить для многих атомов или групп атомов. Поэто- му используют упрощенные формулы.
Физически наиболее обосновано описывать отталкивание молекул экспонен- циальной зависимостью, как в потенциале Букингема. Однако при низких темпе- ратурах и давлениях, когда расстояния между атомами молекул не сильно отли- чаются от равновесных значений, степенная функция также неплохо аппроксими- рует теоретическую кривую. В диполь-дипольном приближении дисперсионное взаимодействие имеет ассимптотику (1/r
6
). Поэтому наиболее широкое распро- странение для описания ван-дер-ваальсовых взаимодействий получил потенциал
Леннард-Джонса (12-6):
U(r) = 4
ε [(σ/r)
12
- (
σ/r)
6
] =
ε[(R
m
/r)
12
– 2 (R
m
/r)
6
]
, (1.26) где
ε - глубина потенциальной ямы, σ - значение r, при котором энергия равна нулю, U(R
m
) = -
ε. Выбор показателя степени равным 12 обусловлен математиче- ским удобством.
Второе слагаемое в потенциале описывает электростатические взаимодейст- вия, которые можно рассчитать по формуле :
U(r
ij
) =
Σ q
i
q
j
/ r
ij
, (1.27) где q
i
, q
j
– заряды атомов или специально выделенных центров, а r
ij
.- расстояния между атомами (центрами взаимодействия) молекул. Электростатический вклад отражает специфические особенности взаимодействия молекул.
В атом-атомной схеме расчета энергии необходимо использовать параметры взаимодействия для каждой пары атомов. Возникает проблема их определения.
Количество параметров можно существенно уменьшить, если использовать ком- бинационные правила:
ε
ij
= (
ε
ii
ε
jj
)
0.5
,
σ
ij
= (
σ
ij
+
σ
ij
)/2 .
(1.28)
Следующая проблема связана с переносимостью параметров. В простейшем варианте параметризации можно подбирать параметры для каждого атома, вхо- дящего в состав молекулы. Однако более предпочтительной выглядит процедура
28 аппроксимации, когда найденные параметры можно использовать для целого класса веществ. С этой проблемой разработчики потенциалов столкнулись при построении силового поля в методе молекулярной механики [38].
1.6. Водородная связь
Водородную связь (Н-связь) можно рассматривать как частный случай специ- фического взаимодействия. Она образуется между атомом водорода, ковалентно связанным с атомом А в одной молекуле, и атомом В, принадлежащим той же или соседней молекуле. В качестве атомов А, В обычно выступают электроотрица- тельные атомы O, N, F, Cl, S. Выделяют внутри- и межмолекулярную Н-связь.
Более подробно рассмотрим физико-химические проявления межмолекулярных
Н-связей, которые подразделяются на связи средней прочности, образованные нейтральными молекулами, и сильные ионные Н-связи.
Образование межмолекулярной Н-связи сопровождается сдвигом в сторону длинных волн и уширением полос поглощения в ИК- и комбинационных спек- трах. Согласно данным спектроскопии ЯМР на протоне, участвующим в Н-связи всегда наблюдается уменьшение электронной плотности. Физико-химические свойства веществ, молекулы которых образуют межмолекулярные Н-связи, суще- ственно отличаются от свойств остальных веществ. Температура плавления и ки- пения, теплота испарения, вязкость, диэлектрическая проницаемость, как прави- ло, имеют более высокие значения в случае образования Н-связи.
В работе [31] предложена классификация по силе межмолекулярной Н-связи: различают слабую (
∆Е =0.4÷4 кДж моль
-1
), среднюю (нейтральную) (
∆Е =20÷60 кДж моль
-1
) и сильную (ионную) (
∆Е =80÷240 кДж моль
-1
) связь (
∆Е – энергия разрыва Н-связи). Сильная Н-связь образуется при взаимодействии иона с моле- кулой, содержащей функциональные группы ОН, ОN, FН. Приведенная класси- фикация достаточно условна, т.к. невозможно каким-либо физическим обосно- ванным способом найти границы деления. Со стороны слабых Н-связей отличия от ван-дер-ваальсовых комплексов становятся малосущественными. Со стороны сильных Н-связей характеристики комплексов близки величинам, наблюдаемым при образовании ковалентных химических связей.
29
Квантово-химические расчеты позволили выявить природу Н-связи. Показа- но, что с точки зрения энергетических характеристик комплексы с Н-связью не обладают никакими преимуществами по сравнению с другими молекулярными комплексами [31]. Относительные вклады различных типов взаимодействий в энергию комплексов с Н-связью примерно такие же, как и в энергию донорно- акцепторных комплексов. Спецификой нейтральной Н-связи считают образова- ние водородного мостика, содержащего умеренно полярную и сильную химиче- скую связь. В отличие от нейтральной, при образовании сильной ионной Н-связи энергетические и структурные характеристики мономеров существенно изменя- ются. Результаты квантово-химического анализа комплексов с Н-связями показа- ли, что особой уникальной природой Н-связи не обладают.
Авторы статьи [39] замечают, что невозможно дать точное определение Н- связи и даже указать, какие взаимодействия, ковалентные или электростатиче- ские, играют основную роль при ее образовании. В некоторых случаях (бесконеч- ные цепочки из молекул карбамида) основными являются электростатические и поляризационные взаимодействия, в других
(енольная форма 1,3- циклогександиона) – ковалентные взаимодействия. Обычно считают [31], что СН
–группы не образуют Н-связей. Однако показано [40], что их взаимодействие с атомами О соседних молекул играет важную роль в формировании структуры жидкой муравьиной кислоты.
Взаимодействие молекул в жидкой фазе обладает той особенностью, что в рамках молекулярного ансамбля реализуются как относительно сильные, так и слабые Н-связи. При этом геометрические и энергетические характеристики Н- связей могут существенно отличаться от величин, наблюдаемых в кристалле или в наиболее стабильных конфигурациях димеров. Как правило, они описываются плавными монотонными кривыми, и в рамках I-ансамбля невозможно однозначно установить наличие связи между двумя выделенными молекулами. Степень на- дежности определения существенно возрастает, если ввести динамический крите- рий Н-связи, рассматривая систему в рамках колебательно - усредненного V- ансамбля. Однако это требует проведения дополнительных расчетов в рамках ме- тода молекулярной динамики и не применимо к конфигурациям, полученным ме- тодом Монте-Карло. Наряду с линейными Н-связями в жидкости могут сущест-
30 вовать бифуркатные [41, 42], трифуркатные и циклические связи. Они имеют раз- ный статистический вес, но оказывают большое влияние на динамические и структурные характеристики жидкости.
Авторы работы [43], посвященной исследованию сверхкритической жидкой воды методами компьютерного моделирования, наряду с линейными Н-связями рассматривают циклические, бифуркатные и дважды согнутые (twofold). (см. рис.
1.3). Линейные связи доминируют при нормальных и сверхкритических условиях, но с ростом температуры доля остальных типов связей существенно повышается.
1 2 3 4 5 6 7 8 9 ... 33