1.9 Экспериментальные сведения о функциях
радиального распределения атомов в жидкой воде
Практически во всех работах по компьютерному моделированию воды прово- дят сравнение расчетных и экспериментальных атом-атомных функций радиаль- ного распределения (ФРР). Решение задачи экспериментального определения
ФРР имеет крайне важное значение для дальнейшего развития теории жидкого состояния. Однако даже в простых случаях не удается получить надежные сведе- ния о поведении ФРР. Для лучшего понимания сути проблемы необходимо обра- титься к основам теории рассеяния элементарных частиц.
Дифференциальное поперечное сечение нейтронов, рассеивающихся на ядрах атомов m, l , описывается выражением [108]
∑
∑
+
=
Ω
n
m
l
m
j
l
m
iQr
l
m
ml
e
b
b
d
d
,
,
,
4
δ
π
σ
σ
, (1.32) где b
i
, b
j
– длины когерентного рассеяния, σ
j
– поперечное сечение некогерентно- го рассеяния. Дифференциальное сечение можно представить в виде суммы от рассеяния на одинаковых (l=m) и разных ядрах (l≠m). В случае D
2
O первый вклад
(self) записывают следующим образом
+
+
=
Ω
π
σ
σ
4 2
)
2
(
2 2
D
D
O
self
b
b
d
d
(1.33) где b
O
, b
D
– длины для когерентного рассеяния атомов O и D, σ
D
– сечение некоге- рентного рассеяния атомов D (σ
O
= 0). Второй вклад (distinct) можно представить в виде суммы от рассеяния на атомах одной молекулы (intra) и атомах разных мо- лекул (inter). Структурный фактор S
M
(Q) записывают в виде
2 2
2
inter intra
)
2
/(
)
2
(
)
(
)
(
)
(
D
O
D
O
M
b
b
b
b
d
d
d
d
Q
S
+
+
+
Ω
+
Ω
=
σ
σ
, (1.34) где Q – волновой вектор.
Можно выделить две составляющие структурного фактора: f(Q) – молекуляр- ный форм-фактор и D
M
(Q) – функцию, содержащую все межмолекулярные вкла- ды.
f(Q)=[b
O
2
+4b
O
b
D
j
0
(Q r
OD
) exp(-γ
OD
Q
2
)+
45
+2 b
D
2
j
0
(Qr
DD
) exp(-γ
DD
Q
2
)]/(b
O
+2b
D
)
2
, (1.35) где j
0
– функция Бесселя, exp(-γ
OD
Q
2
), exp(-γ
DD
Q
2
) – дебаевские факторы, учиты- вающие изменения внутримолекулярных расстояний вследствие колебания ато- мов.
При больших значениях Q D
M
(Q) стремится к нулю. Здесь главный вклад в структурный фактор дает функция f(Q). Однако для жидкостей, в которых суще- ствует система межмолекулярных Н-связей, межмолекулярный вклад убывает очень медленно, поэтому его сложно определить и оценить.
Полную парную корреляционную функцию G(r) определяют с помощью фу- рье-преобразования структурного фактора
4πrρ(G(r)-1) = 2/π∫Q(S
M
(Q)- S
M
(
∝)) sin Qr dQ
, (1.36) где интеграл берется от 0 до
∝, а асимптотическое значение для D
2
O S
M
(
∝) =
f(Q)= 0.3346,
ρ - плотность жидкости. Функция G(r) является суммой атом- атомных ФРР и включает вклады от внутримолекулярных и межмолекулярных расстояний.
Если из функции S
M
(Q) вычесть молекулярный форм-фактор, то получим функцию D
M
(Q), фурье-преобразование которой дает возможность определить межмолекулярную парную корреляционную функцию
4πrρ(G
inter
(r)-1) = 2/π ∫ Q D
M
(Q) sin Qr dQ ,
(1.37)
Можно показать, что для D
2
O[109]
G
inter
(r) = 0.485 g
DD
(r) +0.423 g
OD
(r)+0.092 g
OO
(r)
(1.38)
Следовательно, основной вклад в экспериментально определяемую функцию да- ют корреляции в пространственном положении атомов D относительно атомов D и O других молекул.
Рентгеновское излучение рассеивается на электронах, а не на ядрах атомов.
Определив угловую зависимость когерентного рассеяния можно определить пар- ную корреляционную функцию G
X
(r), основной вклад в которую дает ФРР g
OO
(r).
Ее рассчитывают по соотношению
G
X
(r) = 1+ 1/(2
π
2
ρ r) ∫ Q [S
M
(Q) –1] M(Q) sin Qr dQ
, (1.39) где M(Q) –функция, демпфирующая осцилляции структурного фактора при больших Q. В этой формуле необходимо интегрировать по всем значениям вол-
46 нового вектора Q. Однако в реальном эксперименте доступный интервал Q суще- ственно ограничен. В связи с этим, на кривой G
X
(r) появляются ложные макси- мумы.
Таким образом, процедура определения атом- атомных ФРР на основании экспериментальных сведений об угловой зависимости интенсивности рассеяния сопряжена с необходимостью выделе- ния отдельных составляющих, проведения боль- шого количества расчетов и привлечения сведе- ний о дополнительных константах и функциях. В результате, даже в случае относительно простых жидкостей, молекулы которых состоят из не- большого числа атомов, невозможно надежно ус- тановить поведение атом-атомных ФРР. Сущест- вует большое количество методов уменьшения погрешностей определения корреляционных функций [110, 111, 112, 113, 114, 115, 116, 117,
118], однако устранить экспериментальные и расчетные погрешности не удается. Например, существенно отличаются ФРР g
OO
(r), определен- ные для воды разными авторами, или одной груп- группой авторов, но в разное время (см. рис. 1.5). Функции отличаются по высоте, положению и количеству пиков, что крайне затрудняет процедуру разработки модельных потенциалов.
Первоначально тестирование потенциалов проводили по данным NDIS-86
(Neutron Diffraction with Isotope Substitution) [112]. Затем были опубликованы ре- зультаты исследований жидкой воды, выполненных в широком интервале пара- метров состояний (NDIS-93 [119]). В литературе разгорелась дисскуссия [120,
121, 122, 123], поскольку новые экспериментальные данные существенно отлича- лись от предыдущих. Естественно, что результаты многочисленных компьютер- ных моделирований стали не соответствовать экспериментальным данным [124].
К счастью, последующие исследования показали, что данные NDIS-93 являются
200 400 600
r/пм
0 1
2 3
g
2
3
1
Рис. 1.5.
Функции радиаль- ного распределения атомов
О молекул воды, определен- ные по экспериментальным данным. 1 – [112], 2- [113],
3- [114].
47 артефактом. При обработке экспериментальных результатов для H
2
O была сдела- на неоправданно большая поправка на неупругое рассеяние. Новый набор атом- атомных ФРР (NDIS-97 [125]), принятый в качестве стандартного, существенно отличается от предыдущих. Амплитуды первых пиков g
OO
, g
OH
(NDIS-86) оказа- лись завышенными. Новые результаты (NDIS-97) по функциям g
OO
в большей степени соответствуют данным Горбатого и Демьянца [115, 116]
С помощью методов компьютерного моделирования получают набор молеку- лярных конфигураций, при этом известны координаты каждого атома. Достаточ- но просто можно рассчитать ФРР g(r), которая по определению равна отношению локальной численной плотности
ρ(r) точек в сферическом слое толщиной dr, на- ходящихся на расстоянии r от фиксированной точки, к среднему значению плотности
ρ.
g(r) =
ρ(r)/ρ = [dN(r)/dV]/ρ=(1/(4πr
2
ρ)) [dN(r)/dr]
, (1.40) где dN(r) - число точек в сферическом слое объемом dV.
Функции, полученные экспериментальными и теоретическими методами, можно непосредственно сравнивать. Однако неустранимые погрешности опреде- ления экспериментальной функции не позволяют делать однозначные выводы о степени адекватности модельного потенциала. Иногда сравнивают структурные факторы. Для этого с помощью фурье-преобразования ФРР переходят в про- странство волновых векторов. В этом случае необходимо интегрировать по всем расстояниям 0 < r <
∝. Однако в компьютерном эксперименте ФРР можно рас- считать только при r < L/2, где L –длина ребра кубической ячейки. В результате ограничения предела интегрирования структурный фактор будет определен с большими погрешностями.
В настоящее время все большее число исследователей при тестировании по- тенциалов для воды обращается к экспериментально определяемым изохорным температурным дифференциалам (ИТД). Если обратиться к формуле (1.40), то можно отметить, что в числителе содержится структурная информация
ρ(r), а в знаменателе стоит среднее значение плотности
ρ, характеризующее свойства жидкости на макроуровне. С изменением внешних условий плотность воды меня- ется, поэтому для получения информации о структурных свойствах воды некор- ректно непосредственно сравнивать ФРР, определенные при разных значениях
48 плотности. Как известно, для температурной зависимости плотности жидкой во- ды характерно аномальное поведение. В окрестности максимума функции ρ(T) можно выбрать две точки, отвечающие одному значению плотности, но разным температурам. В этом случае разность ФРР,
∆g
OO
(r,
∆T), будет характеризовать изменение закономерностей расположения молекул в сферических слоях, вы- званное изменением температуры.
∆g
OO
(r,
∆T) = g
OO
(r,T
2
) - g
OO
(r,T
1
)
(1.41)
Легко заметить, что в этом случае уменьшаются систематические погрешно- сти определения ФРР на основе экспериментальных данных. При расчете ИТД не обязательно вычитать ФРР. Можно первоначально рассчитать функцию
∆D
M
(Q,∆T), а затем сделать фурье-преобразование по формуле (1.34). Такая про- цедура расчета позволяет обходиться без демпфирующих функций и существенно уменьшает ошибки, связанные с ограничением верхнего предела интегрирования
[109].
Изохорные температурные дифференциалы были определены по данным о рассеянии рентгеновского излучения [113] и нейтронов [109]. На рис. 1.6 пред- ставлены ФРР и ИТД рассчитанные нами для модели воды с потенциалом Мален- кова-Полтева (МП) и определенные по экспериментальным данным [113]. В пар- но аддитивном неполяризуемом эффективном потенциале МП задается строение жесткой молекулы [126, 127]. Первый пик g
OO
вдвое выше пика «эксперимен- тальной» функции. Можно отметить, что за пределами первой координационной сферы расчетные и экспериментальные данные близки. Внутри сферы наблюда- ется качественное соответствие поведения кривых. Подобные результаты были получены и для модели воды (BSV), учитывающей поляризацию молекул [128].
Следует отметить, что в этой работе соответствие экспериментальных и расчет- ных функций достигнуто в меньшей степени. По результатам воспроизведения экспериментальной температурной зависимости локальной плотности атомов в сферических слоях модель МП вполне адекватна, а эффекты поляризации без по- тери точности можно учесть в параметрах эффективного потенциала.
Изохорные температурные дифференциалы, полученные по результатам ис- следования рассеяния нейтронов [109], отражают изменения функции G(r). Авто- ры работы отмечают несколько интересных закономерностей поведения ИТД.
49
При увеличении температуры на 45К наблюдаются достоверные изменения про- странственной структуры воды на межмолекулярных расстояниях, достигающих
1200 пм. Ранее было принято считать, что осцилляции ФРР g
DD
(r) и g
OD
(r), даю- щих основной вклад в G(r), а следовательно и корреляции в положении этих ато- мов, наблюдаются на расстояниях до 600 пм. Сопоставление двух функций, полу- ченных в результате деления ИТД на разность температур (
∆G(r, ∆T)/∆T) при ∆T
=45K и
∆T =25.5К показало, что они практически идентичны за пределами облас- ти малых расстояний (r<250 пм).
Сейчас не приходится сомневаться в том, что все широко используемые и вновь разрабатываемые потенциалы должны тестироваться на степень воспроиз- ведения экспериментальных ИТД. Систематические погрешности, связанные с постановкой эксперимента и обработкой экспериментальных данных, в наимень- шей степени оказывают влияние на эти функции. Они характеризуют темпера- турную зависимость пространственной структуры воды.
200 400 600 0
2 4
200 400 600
-0.8
-0.4 0.0
1
2
g g
∆
r/ r/
пм пм
a b
Рис. 1.6.
Функции радиального распределения атомов O (a) 1 – модель воды с потенциалом МП, 2- эксперимент [113], и изохорные температурные дифферен- циалы (b). 1 – расчет при
∆T=55K, модель МП, 2- эксперимент при ∆T=51K [113].
50
1.10 Тестирование компьютерных моделей
по экспериментальным данным
Очень сложно перечислить все экспериментальные данные, по которым про- водят проверку адекватности компьютерных моделей жидкостей. Парные корре- ляционные функции чаще всего дают ограниченные сведения о закономерностях пространственного расположения атомов. Остальные экспериментально опреде- ляемые характеристики связаны более сложными функциональными зависимо- стями с потенциалами взаимодействий. При определении любого эксперимен- тально измеряемого свойства жидкости всегда проводят усреднение по фазовому пространству изучаемой системы. Следовательно, они являются характеристика- ми D-структуры жидкости. Поскольку обратная задача восстановления потенциа-
0 . 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 0 . 8 0 . 9 1 . 0 1 . 1 1 . 2 1 . 3 1 . 4 2 5 0 3 0 0 3 5 0 4 0 0 4 5 0 5 0 0 5 5 0 6 0 0 6 5 0 7 0 0 7 5 0 8 0 0 8 5 0 9 0 0 9 5 0
1 3
П а р
Ж и д к о с т ь + П а р
Ж и д к о с т ь
2
4
5
1 2
1 1
1 0
9
8
7
6
3
1
2 7 3 . 1 6 K
1 M P a
0 . 1 M P a
4 0 0
Д а в л е н и е , M П a
3 0 5 0 1 0 0 2 0 0 6 0 0 1 0 0 0
К р и т и ч е с к а я т о ч к а
( T = 6 4 7 . 1 0 K , P = 2 2 . 0 6 4 M П a )
Т
е м
п ер ат у
р а
, K
П л о т н о с т ь , г р / м л
Рис. 1.7
Фазовая диаграмма воды. Сполошными линиями обозначены изо- бары, пунктирной – кривая равновесия жидкость-пар. Точками обозначены обсуждаемые состояния воды.
51 ла по экспериментальным данным является математически некорректно постав- ленной [56], близость расчетных и экспериментальных значений лишь косвенно свидетельствует о корректности параметризации.
Первоначальные исследования свойств жидкой воды проводили при условиях близких к нормальным. Поэтому параметры потенциалов подбирали по степени воспроизведения экспериментальных данных в максимально исследованной об- ласти параметров состояния. Было опубликовано большое количество работ [6,
26, 53, 61, 77, 129, 130, 131, 132, 133, 134], посвященных сопоставлению расчет- ных и экспериментальных: ФРР, термодинамических функциий, данных о плот- ности, втором вириальном коэффициенте пара, диэлектричекой постоянной, дан- ным о свойствах льдов.
В настоящее время проводят тестирования вновь разрабатываемых поляри- зуемых и гибких потенциалов в широкой области параметров состояния. На рис.
1.7 приведена фазовая диаграмма воды, на которой обозначены обсуждаемые да- лее состояния. Остановимся более подробно на результатах двух недавно опуб- ликованных работ [107, 135]. Chialvo с соавторами [107] провел рассчеты струк- турных и термодинамических свойств для двух гибких (BHJ [136], SPC-mTR[88]) и двух поляризуемых (TIP4P-FQ, PPC [78]) моделей воды. В табл. 1.4 приведены результаты расчета термодинамических характеристик моделей воды.
Таблица 1.4
Термодинамические характеристики поляризуемых и гибких моделей во- ды по данным работы [107].
TIP4P-FQ PPC SPC-mTR BJH
N
T
К
P
эксп кбар
P -U P -U P -U P -U
1 2
3 4
5 6
7 298 423 423 573 573 573 673 0.001 1.9 0.1 2.8 0.5 0.1 0.8
-0.52 3.02 0.95 4.74 1.62 1.02 1.79 41.5 33.9 32.4 27.6 24.4 23.1 19.6 0.7 3.2 1.18 4.3 1.6 1.04 1.7 43.1 36.6 35.6 30.2 27.6 26.4 22.6 0.23 0.423 0.8 3.37 0.96 0.69 1.2 41.5 33.8 32.9 26.1 23.6 22.4 17.5 13.4 15.6 11.6 14.5 8.4 6.7 6.3 40.8 35.4 34.5 29.5 26.9 25.6 21.8
Примечание: N – номер точки на рис. 1.7, P –давление (кбар), U – конфи- гурационная потенциальная энергия (кДж моль
-1
)
52
Как следует из приведенных данных, ни одна из исследованных современных моделей не воспроизводит экспериментальные значения давления. Отличия зна- чений внутренней энергии в ряду моделей достигают 20 %. Авторы работы [107] также провели сопоставление атом-атомных ФРР и пришли к выводу, что поляри- зационные модели более предпочтительны, по сравнению с гибкими.
Jedlovszky и Richardi [135] сопоставили экспериментальные данные с данны- ми, полученными в результате моделирования трех поляризуемых моделей: Brod- holt-Sampoli-Vallauri (BSV), Chialvo – Cummings (CC), Dang- Chang (DC) [94] и двух парно-аддитивных (TIP4P, SPC/E). В табл. 1.5 приведены результаты расчета термодинамических характеристик моделей воды, рассчитанные в NPT-ансамбле.
Авторы пришли к выводу, что поляризуемые модели, в противоположность
Таблица 1.5
Термодинамические характеристики поляризуемых и парно- аддитивных моделей воды по данным работы [135]
Свойство BSV
CC DC SPC/E TIP4P Эксп.
N=1, T=298K, P=1 бар
-U
ρ
C
P
42.72 3.622 114.04 39.38 2.956 88.93 41.27 3.542 80.38 47.48 3.542 80.38 42.52 3.52 67.13 41.53 3.34 75.31
N=13, T=573К, Р=1970 бар
-U
ρ
C
P
26.57 2.776 84.85 22.54 2.304 76.10 24.97 2.608 63.52 33.98 2.96 78.83 29.96 2.887 43.77 30.05 2.96 73.41
N=6, T=573К, Р=100 бар
-U
ρ
C
P
2.96 0.152 42.65 2.23 0.140 34.12 2.76 0.164 30.85 29.82 2.220 74.18 8.21 0.242 60.53 26.72 2.41 102.33
N=7, T=673К, Р=800 бар
-U
ρ
C
P
16.22 1.647 88.37 12.17 1.207 86.93
-15.07 1.539 76.52
-26.53 2.138 64.84 20.55 1.857 80.33 23.09 2.17 94.94
Примечание:
U – конфигурационная потенциальная энергия (кДж моль
-1
),
ρ – численная плотность (10
-2
Ǻ
-3
),
C
P
– теплоемкость (Дж моль
-1
К
-1
)
53 неполяризуемым, не описывают температурные изменения внутренней энергии корректно, даже при фиксированном экспериментальном значении плотности
(NVT-ансамбль). С увеличением температуры расчетные значения плотности уменьшаются быстрее экспериментальных при моделировании в NPT-ансамбле.
Быстрое уменьшение плотности влияет и на другие термодинамические свойства.
Авторы были удивлены неудачей в описании термодинамических свойств во- ды в широкой области параметров состояния с помощью поляризуемых моделей.
Термодинамические свойства оказались очень чувствительными к параметрам моделей. Напротив, традиционные модели SPC/E и TIP4P дали очень неплохие результаты. Авторы рассчитывали теплоемкость по флуктуационным формулам, поэтому достоверность ее определения крайне невелика во всех исследованных случаях.
В рамках поляризуемых моделей экспериментальные атом-атомные ФРР вос- производятся лучше. Удлинение Н-связи, происходящее при увеличении темпе- ратуры и давления, наблюдалось только в поляризуемых моделях. Однако ни в одной из моделей форма первого пика ФРР g
OO
не воспроизводится.
Большое количество работ по исследованию свойств воды в широкой области параметров состояния выполнено А. Калиничевым с соавторами [122, 137, 138,
47, 139, 140, 141, 142]. Они исследовали более 40 термодинамических состояний модели TIP4P (273K < T <1273K и 1 бар < P < 10000 бар), при этом плотность из- менялась от 0.02 до 1.67 г см
-3
. В последнее время они часто используют гибкую модель BJH. Авторам не только удалось воспроизвести известные термодинами- ческие и структурные данные о свойствах воды, но и получить важную информа- цию о системе Н-связей, особенностях межмолекулярных и внутримолекулярных движений молекул воды в широком диапазоне параметров состояния.
Интересные результаты об атом-атомных ФРР сверхкритической воды полу- чены международной группой авторов [43]. Они использовали вариант метода молекулярной динамики, разработанный Car и Parrinello (CPMD). Метод основан на первых принципах. Для расчета энергетических и структурных характеристик молекул воды и молекулярных конфигураций авторы использовали теорию функционала плотности. Обладая самой современной компьютерной техникой, они смогли промоделировать систему из 32 молекул воды с периодическими гра-
54 ничными условиями на протяжении 9 пс. При параметрах критической точки воды получено хорошее согласие с экспериментальными ФРР g
OO
и g
HH
. Однако функции g
ОH
заметно отличаются. На расчетной функции наблюдается плечо в области 250 пм, отсутствующее на экспериментальной кривой. Вопрос о причине рассогласования данных остался открытым. Достоверность метода подтверждена и близостью расчетных (46.2 10
-5
, 103.5 10
-5
см
2
с
-1
) и экспериментальных значе- ний (47.4 10
-5
, 112 10
-5
см
2
с
-1
) коэффициента диффузии воды, определенных, со- ответственно, при T=653K,
ρ=0.73 гр см
-3
и T=647K,
ρ=0.32 гр см
-3
Тестирование потенциалов проводят не только в области параметров со- стояния, соответствующих жидкой и газообразной фазам воды, но и в области существования переохлажденной воды и льдов. В настоящее время опубликовано множество результатов экспериментальных и теоретических исследований воды в этой области параметров состояния и выявлено большое количество аномалий.
Тем не менее, множество вопросов о структурных свойствах, особенностях дви-
1 2 3 4 5 6 7 8 9 ... 33
скачати
радиального распределения атомов в жидкой воде
Практически во всех работах по компьютерному моделированию воды прово- дят сравнение расчетных и экспериментальных атом-атомных функций радиаль- ного распределения (ФРР). Решение задачи экспериментального определения
ФРР имеет крайне важное значение для дальнейшего развития теории жидкого состояния. Однако даже в простых случаях не удается получить надежные сведе- ния о поведении ФРР. Для лучшего понимания сути проблемы необходимо обра- титься к основам теории рассеяния элементарных частиц.
Дифференциальное поперечное сечение нейтронов, рассеивающихся на ядрах атомов m, l , описывается выражением [108]
∑
∑
+
=
Ω
n
m
l
m
j
l
m
iQr
l
m
ml
e
b
b
d
d
,
,
,
4
δ
π
σ
σ
, (1.32) где b
i
, b
j
– длины когерентного рассеяния, σ
j
– поперечное сечение некогерентно- го рассеяния. Дифференциальное сечение можно представить в виде суммы от рассеяния на одинаковых (l=m) и разных ядрах (l≠m). В случае D
2
O первый вклад
(self) записывают следующим образом
+
+
=
Ω
π
σ
σ
4 2
)
2
(
2 2
D
D
O
self
b
b
d
d
(1.33) где b
O
, b
D
– длины для когерентного рассеяния атомов O и D, σ
D
– сечение некоге- рентного рассеяния атомов D (σ
O
= 0). Второй вклад (distinct) можно представить в виде суммы от рассеяния на атомах одной молекулы (intra) и атомах разных мо- лекул (inter). Структурный фактор S
M
(Q) записывают в виде
2 2
2
inter intra
)
2
/(
)
2
(
)
(
)
(
)
(
D
O
D
O
M
b
b
b
b
d
d
d
d
Q
S
+
+
+
Ω
+
Ω
=
σ
σ
, (1.34) где Q – волновой вектор.
Можно выделить две составляющие структурного фактора: f(Q) – молекуляр- ный форм-фактор и D
M
(Q) – функцию, содержащую все межмолекулярные вкла- ды.
f(Q)=[b
O
2
+4b
O
b
D
j
0
(Q r
OD
) exp(-γ
OD
Q
2
)+
45
+2 b
D
2
j
0
(Qr
DD
) exp(-γ
DD
Q
2
)]/(b
O
+2b
D
)
2
, (1.35) где j
0
– функция Бесселя, exp(-γ
OD
Q
2
), exp(-γ
DD
Q
2
) – дебаевские факторы, учиты- вающие изменения внутримолекулярных расстояний вследствие колебания ато- мов.
При больших значениях Q D
M
(Q) стремится к нулю. Здесь главный вклад в структурный фактор дает функция f(Q). Однако для жидкостей, в которых суще- ствует система межмолекулярных Н-связей, межмолекулярный вклад убывает очень медленно, поэтому его сложно определить и оценить.
Полную парную корреляционную функцию G(r) определяют с помощью фу- рье-преобразования структурного фактора
4πrρ(G(r)-1) = 2/π∫Q(S
M
(Q)- S
M
(
∝)) sin Qr dQ
, (1.36) где интеграл берется от 0 до
∝, а асимптотическое значение для D
2
O S
M
(
∝) =
f(Q)= 0.3346,
ρ - плотность жидкости. Функция G(r) является суммой атом- атомных ФРР и включает вклады от внутримолекулярных и межмолекулярных расстояний.
Если из функции S
M
(Q) вычесть молекулярный форм-фактор, то получим функцию D
M
(Q), фурье-преобразование которой дает возможность определить межмолекулярную парную корреляционную функцию
4πrρ(G
inter
(r)-1) = 2/π ∫ Q D
M
(Q) sin Qr dQ ,
(1.37)
Можно показать, что для D
2
O[109]
G
inter
(r) = 0.485 g
DD
(r) +0.423 g
OD
(r)+0.092 g
OO
(r)
(1.38)
Следовательно, основной вклад в экспериментально определяемую функцию да- ют корреляции в пространственном положении атомов D относительно атомов D и O других молекул.
Рентгеновское излучение рассеивается на электронах, а не на ядрах атомов.
Определив угловую зависимость когерентного рассеяния можно определить пар- ную корреляционную функцию G
X
(r), основной вклад в которую дает ФРР g
OO
(r).
Ее рассчитывают по соотношению
G
X
(r) = 1+ 1/(2
π
2
ρ r) ∫ Q [S
M
(Q) –1] M(Q) sin Qr dQ
, (1.39) где M(Q) –функция, демпфирующая осцилляции структурного фактора при больших Q. В этой формуле необходимо интегрировать по всем значениям вол-
46 нового вектора Q. Однако в реальном эксперименте доступный интервал Q суще- ственно ограничен. В связи с этим, на кривой G
X
(r) появляются ложные макси- мумы.
Таким образом, процедура определения атом- атомных ФРР на основании экспериментальных сведений об угловой зависимости интенсивности рассеяния сопряжена с необходимостью выделе- ния отдельных составляющих, проведения боль- шого количества расчетов и привлечения сведе- ний о дополнительных константах и функциях. В результате, даже в случае относительно простых жидкостей, молекулы которых состоят из не- большого числа атомов, невозможно надежно ус- тановить поведение атом-атомных ФРР. Сущест- вует большое количество методов уменьшения погрешностей определения корреляционных функций [110, 111, 112, 113, 114, 115, 116, 117,
118], однако устранить экспериментальные и расчетные погрешности не удается. Например, существенно отличаются ФРР g
OO
(r), определен- ные для воды разными авторами, или одной груп- группой авторов, но в разное время (см. рис. 1.5). Функции отличаются по высоте, положению и количеству пиков, что крайне затрудняет процедуру разработки модельных потенциалов.
Первоначально тестирование потенциалов проводили по данным NDIS-86
(Neutron Diffraction with Isotope Substitution) [112]. Затем были опубликованы ре- зультаты исследований жидкой воды, выполненных в широком интервале пара- метров состояний (NDIS-93 [119]). В литературе разгорелась дисскуссия [120,
121, 122, 123], поскольку новые экспериментальные данные существенно отлича- лись от предыдущих. Естественно, что результаты многочисленных компьютер- ных моделирований стали не соответствовать экспериментальным данным [124].
К счастью, последующие исследования показали, что данные NDIS-93 являются
200 400 600
r/пм
0 1
2 3
g
2
3
1
Рис. 1.5.
Функции радиаль- ного распределения атомов
О молекул воды, определен- ные по экспериментальным данным. 1 – [112], 2- [113],
3- [114].
47 артефактом. При обработке экспериментальных результатов для H
2
O была сдела- на неоправданно большая поправка на неупругое рассеяние. Новый набор атом- атомных ФРР (NDIS-97 [125]), принятый в качестве стандартного, существенно отличается от предыдущих. Амплитуды первых пиков g
OO
, g
OH
(NDIS-86) оказа- лись завышенными. Новые результаты (NDIS-97) по функциям g
OO
в большей степени соответствуют данным Горбатого и Демьянца [115, 116]
С помощью методов компьютерного моделирования получают набор молеку- лярных конфигураций, при этом известны координаты каждого атома. Достаточ- но просто можно рассчитать ФРР g(r), которая по определению равна отношению локальной численной плотности
ρ(r) точек в сферическом слое толщиной dr, на- ходящихся на расстоянии r от фиксированной точки, к среднему значению плотности
ρ.
g(r) =
ρ(r)/ρ = [dN(r)/dV]/ρ=(1/(4πr
2
ρ)) [dN(r)/dr]
, (1.40) где dN(r) - число точек в сферическом слое объемом dV.
Функции, полученные экспериментальными и теоретическими методами, можно непосредственно сравнивать. Однако неустранимые погрешности опреде- ления экспериментальной функции не позволяют делать однозначные выводы о степени адекватности модельного потенциала. Иногда сравнивают структурные факторы. Для этого с помощью фурье-преобразования ФРР переходят в про- странство волновых векторов. В этом случае необходимо интегрировать по всем расстояниям 0 < r <
∝. Однако в компьютерном эксперименте ФРР можно рас- считать только при r < L/2, где L –длина ребра кубической ячейки. В результате ограничения предела интегрирования структурный фактор будет определен с большими погрешностями.
В настоящее время все большее число исследователей при тестировании по- тенциалов для воды обращается к экспериментально определяемым изохорным температурным дифференциалам (ИТД). Если обратиться к формуле (1.40), то можно отметить, что в числителе содержится структурная информация
ρ(r), а в знаменателе стоит среднее значение плотности
ρ, характеризующее свойства жидкости на макроуровне. С изменением внешних условий плотность воды меня- ется, поэтому для получения информации о структурных свойствах воды некор- ректно непосредственно сравнивать ФРР, определенные при разных значениях
48 плотности. Как известно, для температурной зависимости плотности жидкой во- ды характерно аномальное поведение. В окрестности максимума функции ρ(T) можно выбрать две точки, отвечающие одному значению плотности, но разным температурам. В этом случае разность ФРР,
∆g
OO
(r,
∆T), будет характеризовать изменение закономерностей расположения молекул в сферических слоях, вы- званное изменением температуры.
∆g
OO
(r,
∆T) = g
OO
(r,T
2
) - g
OO
(r,T
1
)
(1.41)
Легко заметить, что в этом случае уменьшаются систематические погрешно- сти определения ФРР на основе экспериментальных данных. При расчете ИТД не обязательно вычитать ФРР. Можно первоначально рассчитать функцию
∆D
M
(Q,∆T), а затем сделать фурье-преобразование по формуле (1.34). Такая про- цедура расчета позволяет обходиться без демпфирующих функций и существенно уменьшает ошибки, связанные с ограничением верхнего предела интегрирования
[109].
Изохорные температурные дифференциалы были определены по данным о рассеянии рентгеновского излучения [113] и нейтронов [109]. На рис. 1.6 пред- ставлены ФРР и ИТД рассчитанные нами для модели воды с потенциалом Мален- кова-Полтева (МП) и определенные по экспериментальным данным [113]. В пар- но аддитивном неполяризуемом эффективном потенциале МП задается строение жесткой молекулы [126, 127]. Первый пик g
OO
вдвое выше пика «эксперимен- тальной» функции. Можно отметить, что за пределами первой координационной сферы расчетные и экспериментальные данные близки. Внутри сферы наблюда- ется качественное соответствие поведения кривых. Подобные результаты были получены и для модели воды (BSV), учитывающей поляризацию молекул [128].
Следует отметить, что в этой работе соответствие экспериментальных и расчет- ных функций достигнуто в меньшей степени. По результатам воспроизведения экспериментальной температурной зависимости локальной плотности атомов в сферических слоях модель МП вполне адекватна, а эффекты поляризации без по- тери точности можно учесть в параметрах эффективного потенциала.
Изохорные температурные дифференциалы, полученные по результатам ис- следования рассеяния нейтронов [109], отражают изменения функции G(r). Авто- ры работы отмечают несколько интересных закономерностей поведения ИТД.
49
При увеличении температуры на 45К наблюдаются достоверные изменения про- странственной структуры воды на межмолекулярных расстояниях, достигающих
1200 пм. Ранее было принято считать, что осцилляции ФРР g
DD
(r) и g
OD
(r), даю- щих основной вклад в G(r), а следовательно и корреляции в положении этих ато- мов, наблюдаются на расстояниях до 600 пм. Сопоставление двух функций, полу- ченных в результате деления ИТД на разность температур (
∆G(r, ∆T)/∆T) при ∆T
=45K и
∆T =25.5К показало, что они практически идентичны за пределами облас- ти малых расстояний (r<250 пм).
Сейчас не приходится сомневаться в том, что все широко используемые и вновь разрабатываемые потенциалы должны тестироваться на степень воспроиз- ведения экспериментальных ИТД. Систематические погрешности, связанные с постановкой эксперимента и обработкой экспериментальных данных, в наимень- шей степени оказывают влияние на эти функции. Они характеризуют темпера- турную зависимость пространственной структуры воды.
200 400 600 0
2 4
200 400 600
-0.8
-0.4 0.0
1
2
g g
∆
r/ r/
пм пм
a b
Рис. 1.6.
Функции радиального распределения атомов O (a) 1 – модель воды с потенциалом МП, 2- эксперимент [113], и изохорные температурные дифферен- циалы (b). 1 – расчет при
∆T=55K, модель МП, 2- эксперимент при ∆T=51K [113].
50
1.10 Тестирование компьютерных моделей
по экспериментальным данным
Очень сложно перечислить все экспериментальные данные, по которым про- водят проверку адекватности компьютерных моделей жидкостей. Парные корре- ляционные функции чаще всего дают ограниченные сведения о закономерностях пространственного расположения атомов. Остальные экспериментально опреде- ляемые характеристики связаны более сложными функциональными зависимо- стями с потенциалами взаимодействий. При определении любого эксперимен- тально измеряемого свойства жидкости всегда проводят усреднение по фазовому пространству изучаемой системы. Следовательно, они являются характеристика- ми D-структуры жидкости. Поскольку обратная задача восстановления потенциа-
0 . 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 0 . 8 0 . 9 1 . 0 1 . 1 1 . 2 1 . 3 1 . 4 2 5 0 3 0 0 3 5 0 4 0 0 4 5 0 5 0 0 5 5 0 6 0 0 6 5 0 7 0 0 7 5 0 8 0 0 8 5 0 9 0 0 9 5 0
1 3
П а р
Ж и д к о с т ь + П а р
Ж и д к о с т ь
2
4
5
1 2
1 1
1 0
9
8
7
6
3
1
2 7 3 . 1 6 K
1 M P a
0 . 1 M P a
4 0 0
Д а в л е н и е , M П a
3 0 5 0 1 0 0 2 0 0 6 0 0 1 0 0 0
К р и т и ч е с к а я т о ч к а
( T = 6 4 7 . 1 0 K , P = 2 2 . 0 6 4 M П a )
Т
е м
п ер ат у
р а
, K
П л о т н о с т ь , г р / м л
Рис. 1.7
Фазовая диаграмма воды. Сполошными линиями обозначены изо- бары, пунктирной – кривая равновесия жидкость-пар. Точками обозначены обсуждаемые состояния воды.
51 ла по экспериментальным данным является математически некорректно постав- ленной [56], близость расчетных и экспериментальных значений лишь косвенно свидетельствует о корректности параметризации.
Первоначальные исследования свойств жидкой воды проводили при условиях близких к нормальным. Поэтому параметры потенциалов подбирали по степени воспроизведения экспериментальных данных в максимально исследованной об- ласти параметров состояния. Было опубликовано большое количество работ [6,
26, 53, 61, 77, 129, 130, 131, 132, 133, 134], посвященных сопоставлению расчет- ных и экспериментальных: ФРР, термодинамических функциий, данных о плот- ности, втором вириальном коэффициенте пара, диэлектричекой постоянной, дан- ным о свойствах льдов.
В настоящее время проводят тестирования вновь разрабатываемых поляри- зуемых и гибких потенциалов в широкой области параметров состояния. На рис.
1.7 приведена фазовая диаграмма воды, на которой обозначены обсуждаемые да- лее состояния. Остановимся более подробно на результатах двух недавно опуб- ликованных работ [107, 135]. Chialvo с соавторами [107] провел рассчеты струк- турных и термодинамических свойств для двух гибких (BHJ [136], SPC-mTR[88]) и двух поляризуемых (TIP4P-FQ, PPC [78]) моделей воды. В табл. 1.4 приведены результаты расчета термодинамических характеристик моделей воды.
Таблица 1.4
Термодинамические характеристики поляризуемых и гибких моделей во- ды по данным работы [107].
TIP4P-FQ PPC SPC-mTR BJH
N
T
К
P
эксп кбар
P -U P -U P -U P -U
1 2
3 4
5 6
7 298 423 423 573 573 573 673 0.001 1.9 0.1 2.8 0.5 0.1 0.8
-0.52 3.02 0.95 4.74 1.62 1.02 1.79 41.5 33.9 32.4 27.6 24.4 23.1 19.6 0.7 3.2 1.18 4.3 1.6 1.04 1.7 43.1 36.6 35.6 30.2 27.6 26.4 22.6 0.23 0.423 0.8 3.37 0.96 0.69 1.2 41.5 33.8 32.9 26.1 23.6 22.4 17.5 13.4 15.6 11.6 14.5 8.4 6.7 6.3 40.8 35.4 34.5 29.5 26.9 25.6 21.8
Примечание: N – номер точки на рис. 1.7, P –давление (кбар), U – конфи- гурационная потенциальная энергия (кДж моль
-1
)
52
Как следует из приведенных данных, ни одна из исследованных современных моделей не воспроизводит экспериментальные значения давления. Отличия зна- чений внутренней энергии в ряду моделей достигают 20 %. Авторы работы [107] также провели сопоставление атом-атомных ФРР и пришли к выводу, что поляри- зационные модели более предпочтительны, по сравнению с гибкими.
Jedlovszky и Richardi [135] сопоставили экспериментальные данные с данны- ми, полученными в результате моделирования трех поляризуемых моделей: Brod- holt-Sampoli-Vallauri (BSV), Chialvo – Cummings (CC), Dang- Chang (DC) [94] и двух парно-аддитивных (TIP4P, SPC/E). В табл. 1.5 приведены результаты расчета термодинамических характеристик моделей воды, рассчитанные в NPT-ансамбле.
Авторы пришли к выводу, что поляризуемые модели, в противоположность
Таблица 1.5
Термодинамические характеристики поляризуемых и парно- аддитивных моделей воды по данным работы [135]
Свойство BSV
CC DC SPC/E TIP4P Эксп.
N=1, T=298K, P=1 бар
-U
ρ
C
P
42.72 3.622 114.04 39.38 2.956 88.93 41.27 3.542 80.38 47.48 3.542 80.38 42.52 3.52 67.13 41.53 3.34 75.31
N=13, T=573К, Р=1970 бар
-U
ρ
C
P
26.57 2.776 84.85 22.54 2.304 76.10 24.97 2.608 63.52 33.98 2.96 78.83 29.96 2.887 43.77 30.05 2.96 73.41
N=6, T=573К, Р=100 бар
-U
ρ
C
P
2.96 0.152 42.65 2.23 0.140 34.12 2.76 0.164 30.85 29.82 2.220 74.18 8.21 0.242 60.53 26.72 2.41 102.33
N=7, T=673К, Р=800 бар
-U
ρ
C
P
16.22 1.647 88.37 12.17 1.207 86.93
-15.07 1.539 76.52
-26.53 2.138 64.84 20.55 1.857 80.33 23.09 2.17 94.94
Примечание:
U – конфигурационная потенциальная энергия (кДж моль
-1
),
ρ – численная плотность (10
-2
Ǻ
-3
),
C
P
– теплоемкость (Дж моль
-1
К
-1
)
53 неполяризуемым, не описывают температурные изменения внутренней энергии корректно, даже при фиксированном экспериментальном значении плотности
(NVT-ансамбль). С увеличением температуры расчетные значения плотности уменьшаются быстрее экспериментальных при моделировании в NPT-ансамбле.
Быстрое уменьшение плотности влияет и на другие термодинамические свойства.
Авторы были удивлены неудачей в описании термодинамических свойств во- ды в широкой области параметров состояния с помощью поляризуемых моделей.
Термодинамические свойства оказались очень чувствительными к параметрам моделей. Напротив, традиционные модели SPC/E и TIP4P дали очень неплохие результаты. Авторы рассчитывали теплоемкость по флуктуационным формулам, поэтому достоверность ее определения крайне невелика во всех исследованных случаях.
В рамках поляризуемых моделей экспериментальные атом-атомные ФРР вос- производятся лучше. Удлинение Н-связи, происходящее при увеличении темпе- ратуры и давления, наблюдалось только в поляризуемых моделях. Однако ни в одной из моделей форма первого пика ФРР g
OO
не воспроизводится.
Большое количество работ по исследованию свойств воды в широкой области параметров состояния выполнено А. Калиничевым с соавторами [122, 137, 138,
47, 139, 140, 141, 142]. Они исследовали более 40 термодинамических состояний модели TIP4P (273K < T <1273K и 1 бар < P < 10000 бар), при этом плотность из- менялась от 0.02 до 1.67 г см
-3
. В последнее время они часто используют гибкую модель BJH. Авторам не только удалось воспроизвести известные термодинами- ческие и структурные данные о свойствах воды, но и получить важную информа- цию о системе Н-связей, особенностях межмолекулярных и внутримолекулярных движений молекул воды в широком диапазоне параметров состояния.
Интересные результаты об атом-атомных ФРР сверхкритической воды полу- чены международной группой авторов [43]. Они использовали вариант метода молекулярной динамики, разработанный Car и Parrinello (CPMD). Метод основан на первых принципах. Для расчета энергетических и структурных характеристик молекул воды и молекулярных конфигураций авторы использовали теорию функционала плотности. Обладая самой современной компьютерной техникой, они смогли промоделировать систему из 32 молекул воды с периодическими гра-
54 ничными условиями на протяжении 9 пс. При параметрах критической точки воды получено хорошее согласие с экспериментальными ФРР g
OO
и g
HH
. Однако функции g
ОH
заметно отличаются. На расчетной функции наблюдается плечо в области 250 пм, отсутствующее на экспериментальной кривой. Вопрос о причине рассогласования данных остался открытым. Достоверность метода подтверждена и близостью расчетных (46.2 10
-5
, 103.5 10
-5
см
2
с
-1
) и экспериментальных значе- ний (47.4 10
-5
, 112 10
-5
см
2
с
-1
) коэффициента диффузии воды, определенных, со- ответственно, при T=653K,
ρ=0.73 гр см
-3
и T=647K,
ρ=0.32 гр см
-3
Тестирование потенциалов проводят не только в области параметров со- стояния, соответствующих жидкой и газообразной фазам воды, но и в области существования переохлажденной воды и льдов. В настоящее время опубликовано множество результатов экспериментальных и теоретических исследований воды в этой области параметров состояния и выявлено большое количество аномалий.
Тем не менее, множество вопросов о структурных свойствах, особенностях дви-
1 2 3 4 5 6 7 8 9 ... 33