1. ПРИМЕНЕНИЕ МЕТОДА МОНТЕ-КАРЛО ДЛЯ
МОДЕЛИРОВАНИЯ СТРУКТУРЫ ЖИДКОСТЕЙ
.
1.1 Общие теоретические положения
Метод Монте-Карло, применяемый в статистической физике, является част- ным случаем общего метода статистического моделирования, который использу- ют для решения широкого круга задач в различных областях науки. В рамках ме- тода Гиббса термодинамические характеристики вещества получают в результате усреднения по ансамблю, т.е. по совокупности очень большого числа идентичных по природе систем, находящихся в одинаковых внешних условиях и различаю- щихся только по микросостоянию [1, 2, 3, 4, 5, 6]. Принято выделять три типа ан- самблей, состояние которых задается тремя типами функций распределения.
В микроканоническом ансамбле (N,V,E) рассматривают замкнутые изолиро- ванные системы, в которых фиксированы число частиц N, объем V и полная энер- гия E. На микроскопическом уровне существует бесконечное число различных способов, или конфигураций, в которых может быть реализовано данное макро- состояние.
Большинство физических систем не являются полностью изолированными.
Они могут обмениваться энергией и частицами с окружающей средой. При этом полагают, что рассматриваемая система мала по сравнению с окружающей ее системой, и любое изменение характеристик малой системы не сказывается на состоянии большой. Большая система действует как тепловой резервуар или теп- ловая баня с заданной абсолютной температурой Т. В большом каноническом ан- самбле (T,V,µ) системы способны обмениваться и энергией, и частицами. Его со- стояние задается температурой Т, объемом V и химическим потенциалом µ. Рас- чет термодинамических характеристик, как правило, проводится в рамках кано-
нического ансамбля (NVT, NPT). Каноническое распределение Гиббса – статисти- ческое распределение для систем, содержащих заданное число частиц N, объем V
(или давление P) и способных обмениваться энергией с окружением.
Вероятность нахождения системы в микросостоянии i с энергией E
i рассчи- тывают по формуле:
w
i
= (1/Z) exp (-E
i
/kT),
(1.1)
14 где k – постоянная Больцмана, Z – статистическая сумма, сумма по состояниям системы:
Z =
∑
i
exp(-E
i
/kT)
(1.2)
Квантовые статистические распределения для ансамблей фермионов и бозонов различны. В обычных флюидных системах эти различия не проявляются, и при решении задач теории молекулярных растворов практически всегда можно поль- зоваться классической статистикой. Для реальных систем квантовые закономер- ности требуется учитывать лишь при описании внутримолекулярных состояний, прежде всего электронных и колебательных. Вклад межмолекулярных взаимо- действий в термодинамические функции, структурные характеристики можно найти, пользуясь формулами классической статистической термодинамики, рас- сматривая молекулы как объекты, подчиняющиеся законам классической механи- ки.
В классической статистической термодинамике микросостояние определяет- ся заданием обобщенных координат q и обобщенных импульсов p. Для канониче- ского ансамбля N частиц вероятность иметь значения импульсов в интервале (p,
p+
∆p), значения координат в интервале (q, q+∆q) определяется как: dw(p,q) = (1/Z) exp[-H(p,q)/kT] dpdq
(1.3)
Для классической системы статистическая сумма (1.2) заменяется статистическим интегралом
∫ ∫
∏
−
=
⋅
⋅
=
q
p
e
N
h
Z
kT
q
p
H
i
l
i
f
N
i
i
d d
!
1
)
,
(
1
(1.4)
fi - число степеней свободы молекул сорта i.
Здесь
)
(
2
)
(
1 1
2 1
внутр
N
N
N
i
i
N
i
.....q
q
U
m
p
H
p,q
H
+
+
=
∑
∑
=
=
(1.5) где m - масса молекулы. Функцию Гамильтона, отсчитываемую от нулевой энер- гии молекул, можно представить как сумму энергии внутренних молекулярных
15 движений (электронные состояния, колебания, вращения и т.д.) H
внутр
, энергии поступательного движения центров масс и потенциальной энергии межмолеку- лярных взаимодействий.
С учетом этого статистическая сумма Q может быть представлена в виде: внутр пост
0
Q
Q
e
Q
kT
E
⋅
⋅
=
−
(1.6) где E
0
– энергия молекулы в самом низком энергетическом состоянии, Q
пост
- ста- тистическая сумма, связанная с поступательным движением молекулы, Q
внутр
- статистическая сумма, связанная с внутренними молекулярными движениями.
Подставляя выражение для Гамильтониана в (1.3) и, учитывая условие нор- мировки, получим выражение для статистического интеграла в виде конф внутр
2 3
2
!
2
Z
N
Q
h
mkT
Z
N
N
⋅
=
π
(1.7) где
∫ ∫
−
=
N
kT
U
q
q
e
Z
d d
1
конф
– конфигурационный интеграл. Если стати- стический интеграл известен, то для рассматриваемой системы можно найти все термодинамические величины; так свободная энергия системы равна
Z
kT
F
ln
−
=
(1.8)
С помощью известных термодинамических соотношений могут быть найдены давление, энтропия и химический потенциал системы:
T,V
V,N
T,N
N
F
T
F
S
V
F
p
=
=
=
d d
м d
d d
d
(1.9)
Конфигурационный интеграл Z
конф как функция температуры T, объема V и числа частиц N дает полную статистико-механическую информацию о системе и позволяет по формулам (1.7)-(1.9) рассчитать термодинамические свойства сис- темы.
Невозможность точного вычисления конфигурационного интеграла для ре- альных систем приводит к необходимости применения новых методов, в которых избегают непосредственного вычисления Z
конф
. Одним из таких методов расчета является метод Монте-Карло.
16
1.2. Алгоритм Метрополиса
Особенностям применения метода Монте-Карло для моделирования жидко- стей посвящено большое количество работ. Укажем лишь основные монографии, в которых излагаются теоретические положения и даются рекомендации к прак- тическому применению метода [6, 7, 8, 9, 10].
Как отмечалось выше, в случае сложных жидкостей практически невозможно точно рассчитать значение конфигурационной статистической суммы системы.
Однако его можно оценить с помощью метода статистических испытаний, гене- рируя конечный набор молекулярных конфигураций и определяя вероятность их появления
w
i
w
i
= exp (-
U
i
/
kT) /
∑
i
exp (-
U
i
/
kT)
(1.10)
Если в некоторый фиксированный объем помещать случайным образом моле- кулы, энергия взаимодействия между которыми задается набором потенциальных функций, то в зависимости от конфигурации системы больцмановский множи- тель exp(-U
i
/kT) может принимать различные значения. Некоторые конфигура- ции дают значительный вклад в канонические средние, а некоторые – практиче- ски нулевой (например, когда две частицы сближены настолько, что между ними имеется сильное отталкивание).
При случайной генерации конфигураций подавляющее их большинство будет давать вклад близкий нулю. Поэтому необходимо пользоваться методом сущест- венной выборки, в соответствии с которым конфигурации генерируют с заданной функцией распределения вероятностей
π
i
Среднее по ансамблю от любой физической величины M рассчитывают по формуле
∑
∑
∑
=
−
−
>=
<
i
i
i
i
i
i
i
i
w
M
kT
U
kT
U
M
M
)
/
(
exp
/
)
/
(
exp
(1.11) где i – номер конфигурации (среднее берется по всем рассмотренным конфигура- циям системы). Поскольку усреднение (1.11.) проводят по конечному числу кон- фигураций m со смещенной выборкой, для исключения влияния смещения каж- дую конфигурацию необходимо брать с весом 1/
π
i
:
17
∑
∑
=
=
−
−
>≈
<
m
i
m
i
i
i
i
i
i
kT
U
kT
U
M
M
1 1
)
/
(
exp
)
р
/
1
(
/
)
/
(
exp
)
р
/
1
(
(1.12)
Метрополис с соавторами [7] предложил в качестве
π
i
взять распределение
Больцмана
∑
−
−
=
i
i
i
i
kT
U
kT
U
)
/
exp(
/
)
/
exp(
р
(1.13)
В результате среднее значение любой физической величины M можно записать в виде
∑
=
>≈
<
m
i
i
M
m
M
1
)
/
1
(
(1.14)
Ансамбль, состоящий из m конфигураций, получают путем задания вероят- ностей перехода от одной конфигурации к другой. Вероятность перехода от i-й конфигурации к j-й p
ij
считают зависящей от энергий этих конфигураций, а точнее от величины
(
U
j
–
U
i
)/
kT
:
p
ji
=
p
ij
exp[ -(
U
j
–
U
i
)/
kT]
(1.15)
Таким образом, строят простые цепи Маркова, т.е. последовательности случай- ных событий, в которых вероятность определенного события зависит от исхода предыдущего испытания. В соответствии с условием микроскопической обрати- мости вероятности p
ij
должны удовлетворять условиям:
p
ij
exp( -
U
i
/
kT) = p
ji
exp
( -U
j
/
kT)
и
Σ p
ij
=
1
(1.16)
Обычно полагают, что при
U
j
≤ U
i
p
ij
=
w
i
(1.17) и при
U
j
> U
i
p
ij
= w
i
exp[ -(
U
j
–
U
i
)/
kT ],
где w
i
– вероятность появления некоторой конфигурации при случайном выборе
(с использованием последовательности равномерно распределенных случайных чисел). Используя центральную предельную теорему теории вероятностей может быть доказано [8], что рассматриваемая цепь Маркова задает распределение, асимптотически стремящееся к каноническому.
18
На практике алгоритм Метрополиса реализуют следующим образом. Пусть заданы потенциал взаимодействия, конфигурация системы (начальное располо- жение частиц в элементарной ячейке моделирования) и температура Т. Рассчиты- вают потенциальную энергию системы U
i
и вносят случайное изменение в кон- фигурацию (случайным образом выбирают k-ю частицу в ячейке и смещают ее).
При этом энергия системы становится равной U
j
. Если U
j
< U
i
, то считают, что система перешла в новое состояние. Если U
j
> U
i
, то сравнивают величину exp[-
(U
j
–U
i
)/kT ] со случайным числом
ξ ∈ (0,1). Если ξ ≤ exp[ -(U
j
–U
i
)/kT ], то счита- ют, что система перешла в j-е состояние. Если же
ξ > exp[-(U
j
–U
i
)/kT ], то переход в новое состояние не происходит, k-я частица сохраняет свои прежние координа- ты. При этом j-ую конфигурацию в цепи не учитывают, а рассматривают прежнее расположение частиц, соответствующее энергии U
i
. Таким образом, чем больше значение энергии имеет система при случайном изменении конфигурации, тем с меньшей вероятностью она переходит в это состояние.
Максимальную величину сдвига и поворота молекулы выбирают так, чтобы количество принятых и отвергнутых конфигураций было приблизительно одина- ковым. При других значениях отношения принятых конфигураций к их полному числу, как правило, наблюдается более медленная сходимость результатов к ка- ноническому значению.
В результате генерирования цепей Маркова длиной в несколько миллионов конфигураций отбрасывают начальный неравновесный участок цепи, а на равно- весном участке отбирают m статистически независимых молекулярных конфигу- раций, по которым рассчитывают средние значения физических величин. Оче- видно, что расчет средних значений сопряжен с большим количеством вычисле- ний, выполнение которых стало возможным только при появлении быстродейст- вующих ЭВМ.
Одно из преимуществ метода Монте-Карло состоит в том, что алгоритм легко адаптировать к любому статистическому ансамблю. Например, при моделирова- нии системы в NPT-ансамбле необходимо периодически изменять объем ячейки, а в уравнении (1.17) вместо разности потенциальных энергий использовать раз- ность энтальпий:
∆H = ∆U + P∆V – kT ln(1+∆V/V)
N
,
(1.18)
19 где P – давление, V – объем ячейки,
∆V – изменение объема ячейки.
Усредняя термодинамические функции по конфигурациям на равновесном участке цепи Маркова по уравнению (1.14), легко рассчитать конфигурационную энтальпию и мольный объем :
Н
конфиг
= <
U> + P <V>, V
m
= <
V> N
A
/
N
(1.19)
Для расчета других термодинамических функций – изобарической теплоем- кости С
P
, изотермической сжимаемости k
T
и коэффициента температурного рас- ширения
α, необходимо проводить моделирования системы при различающихся параметрах состояния, а затем находить конечные разности. Оценку указанных термодинамических функций можно сделать и по результатам одного моделиро- вания, используя флуктуационные формулы:
>
<
>
><
<
−
>
<
=
∂
∂
=
>
<
>
<
−
>
<
=
∂
∂
−
=
>
<
−
>
<
=
∂
∂
=
V
NkT
V
H
HV
T
V
V
V
NkT
V
V
P
V
V
k
NkT
H
H
T
H
C
P
T
T
P
P
2 2
2 2
2 2
1
,
1
,
α
(1.20)
Для решения большинства задач вполне достаточно проведения вычислений в каноническом ансамбле. Однако если требуется, то можно использовать и боль- шой канонический ансамбль [11, 12].
1.3. Периодические граничные условия
Очевидно, что какими бы мощными не были компьютеры, невозможно ре- шать уравнения для макро объемов жидкости, содержащих порядка 10 23
молекул.
Свойства системы, состоящей из сотен - тысяч молекул существенно отличаются от макро свойств жидкости. Во-первых, энергетические, динамические характе- ристики молекул, находящихся вблизи поверхности и внутри объема микрокапли различны. Влияние поверхностных эффектов тем значительнее, чем меньше раз- мер исследуемого объекта. Во-вторых, при наличии границы раздела фаз с тече- нием времени меняется число частиц и объем рассматриваемой системы. Следо- вательно, невозможно провести расчеты такого объекта в относительно простом каноническом ансамбле (NVT или NPT). Для преодоления этих затруднений ис-
20 пользуют специально разработанные методики расчета. Остановимся на двух наиболее часто применяемых методах.
Для минимизации влияния указанных эффектов используют периодические граничные условия. В элементарную ячейку моделирования, которую выбирают чаще всего в форме куба, помещают N частиц. Длину ребра ячейки L рассчиты- вают по экспериментальному зна- чению плотности жидкости. Все бесконечное пространство запол- няют аналогичными ячейками – об- разами основной ячейки. Так, на плоскости основную ячейку окру- жают восемь образов (см. рис. 1.1.), а в трехмерном пространстве – два- дцать шесть. Процедура моделиро- вания по методу Монте-Карло учи- тывает независимые смещения час- тиц в основной ячейке; при этом те же смещения одновременно испытывают частицы во всех образах. Если в резуль- тате смещения частица из основной ячейки выйдет за ее пределы, то через проти- воположную грань ячейки входит новая частица, идентичная уходящей. При вы- числении полной энергии учитывают взаимодействия частиц ячейки не только между собой, но и с частицами в ячейках-образах. Таким способом удается со- хранить постоянной среднюю численную плотность и минимизировать влияние поверхностных эффектов.
Поскольку каждая молекула находится на конечном расстоянии от своих об- разов, возникает задача корректного расчета потенциальной энергии. Количество соседних молекул окружающих выделенную молекулу возрастает пропорцио- нально третей степени расстояния, при этом энергия парного взаимодействия, как правило, убывает более быстро. Используют различные способы преодоления данного затруднения. В простейшем случае ограничивают область действия по- тенциалов межмолекулярного взаимодействия, т.е. начиная с некоторого расстоя-
• •
• •
•
• •
• •
•5`
• •
• •
•
• •
• •
•
•1 •2 3
• 4•
•5
• •
• •
•
• •
• •
•
• •
• •
•
• •
• •
•
1 2 3 4 5 6 7 8 9 ... 33