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

• Название: Практикум вычислительные методы исследования молекулярной динамики аксенова е. В., Кшевецкий м. С. Учебнометодическое пособие (для студентов физического факультета)
• Дата: 03.09.2020
• Размер: 466.62 Kb.
• Имя файла: VychMethodMolDyn.pdf
• Тип: Практикум; #136653
• страница: 1 из 2

САНКТ-ПЕТЕРБУРГСКИЙ
ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
ФИЗИЧЕСКИЙ ПРАКТИКУМ
ВЫЧИСЛИТЕЛЬНЫЕ МЕТОДЫ ИССЛЕДОВАНИЯ
МОЛЕКУЛЯРНОЙ ДИНАМИКИ
АКСЕНОВА Е.В., КШЕВЕЦКИЙ М.С.
Учебно-методическое пособие
(для студентов физического факультета)
Санкт-Петербург
2009

Вычислительные методы исследования молекулярной динамики. – СПб.:
СПбГУ, 2009. – 50 с.
Авторы-составители: к.ф.-м.н. Аксенова Е.В., к.ф.-м.н. Кшевецкий М.С.
Рецензенты:
заведующий кафедрой молекулярной биофизики проф., д.ф.-м.н. Трусов А.А.,
заведующий кафедрой статистической физики проф., д.ф.-м.н. Щёкин А.К.
(Санкт-Петербургский государственный университет)
Печатается по решению методической комиссии физического факультета
СПбГУ.
Рекомендовано Ученым советом физического факультета СПбГУ.
Настоящее пособие содержит вспомогательный материал в виде описа- ния лабораторных работ по модулю физического практикума “Вычисли- тельные методы исследования молекулярной динамики”, проводимого на кафедре статистической физики СПбГУ для студентов физического фа- культета.
Пособие позволяет студентам выполнять подготовку к занятиям и мо- жет быть полезным для всех изучающих статистическую физику и термо- динамику и желающих освоить численные методы моделирования физиче- ских процессов.

Содержание
Введение
5 1 Метод молекулярной динамики
6 1.1
Классическая молекулярная динамика . . . . . . . . . . . . .
6 1.2
Потенциал межмолекулярного взаимодействия . . . . . . . .
7 1.3
Уравнения движения . . . . . . . . . . . . . . . . . . . . . . .
8 1.4
Граничные условия . . . . . . . . . . . . . . . . . . . . . . . .
10 1.5
Вычисление средних . . . . . . . . . . . . . . . . . . . . . . .
14 1.6
Начальные условия . . . . . . . . . . . . . . . . . . . . . . . .
15 1.7
Обезразмеривание . . . . . . . . . . . . . . . . . . . . . . . .
15
Лабораторная работа 1.1. Моделирование динамики замкнутой системы многих частиц . . . . . . . . . . . . . . . . . . . . .
17
Лабораторная работа 1.2. Моделирование динамики системы многих частиц при постоянной температуре или давлении
19
Приложение. Комментарии к коду программ . . . . . . . . . . . .
20 2 Метод Монте-Карло
21 2.1
Вычисление интегралов простейшим методом Монте-Карло
21 2.2
Анализ погрешности метода Монте-Карло
22 2.3
Метод Монте-Карло и канонический ансамбль . . . . . . . .
24 2.4
Двумерная модель Изинга . . . . . . . . . . . . . . . . . . . .
29 2.5
Микроканонический ансамбль
33
Лабораторная работа 2.1. Моделирование системы методом
Метрополиса в каноническом ансамбле . . . . . . . . . . . .
36
Лабораторная работа 2.2. Моделирование системы методом
Монте-Карло в микроканоническом ансамбле . . . . . . . . .
37
Лабораторная работа 2.3. Двумерная модель Изинга
38
Приложение. Комментарии к коду программ . . . . . . . . . . . .
39 3 Броуновская динамика
41 3.1
Уравнение Ланжевена . . . . . . . . . . . . . . . . . . . . . .
41 3.2
Решение уравнения Ланжевена . . . . . . . . . . . . . . . . .
42 3

3.3
Численное решение уравнения Ланжевена . . . . . . . . . . .
44 3.4
Анализ особенностей реализации случайной силы в численной схеме
45 3.5
Сила, действующая продолжительное время . . . . . . . . .
47
Лабораторная работа 3. Броуновская динамика
48
Приложение. Комментарии к коду программы . . . . . . . . . . .
49
Список литературы
49 4

Введение
Настоящее пособие представляет собой описания лабораторных работ по моделированию систем многих частиц. В частности, описываются та- кие известные методы, как метод молекулярной динамики и метод Монте-
Карло. Наряду с самими методами в пособии обсуждается их реализация для описания конкретных систем. Пособие ориентировано на студентов третьего курса физического факультета. Для понимания излагаемого ма- териала необходимы базовые знания по термодинамике и статистической физике.
Пособие состоит из трех разделов. Первый раздел посвящен методу мо- лекулярной динамики. Рассматривается классическая молекулярная дина- мика для систем частиц, взаимодействие между которыми описывается с помощью потенциала Леннард-Джонса. Исследуются следующие системы:
замкнутая система, система, находящаяся в термостате, и система, нахо- дящаяся под постоянным внешним давлением. Особое внимание уделено заданию граничных условий в зависимости от того, какие параметры си- стемы являются постоянными. Также рассмотрен вопрос о роли обезраз- меривания в задачах численного моделирования.
Во втором разделе изучается моделирование с помощью метода Монте-
Карло. Сначала приводится способ расчета одномерного интеграла мето- дом Монте-Карло и выполняется анализ погрешности этого расчета. За- тем дается краткое описание микроканонического ансамбля и описывается моделирование замкнутой системы методом Монте-Карло. Также в этом разделе обсуждается канонический ансамбль и моделирование систем с по- мощью алгоритма Метрополиса. В заключении этого раздела приводится описание одномерной и двумерной моделей Изинга.
В третьем разделе рассматривается движение броуновской частицы. Это движение описывается уравнением Ланжевена со случайной силой. Приво- дится решение этого уравнения и численный алгоритм для описания дви- жения броуновской частицы. Также обсуждается модификация численного алгоритма в случае, если случайная сила действует конечное время.
Каждый раздел заканчивается заданиями к лабораторным работам и
5

комментариями к коду программ. Лабораторные работы выполняются сту- дентами самостоятельно. Прилагаемая к каждой лабораторной работе про- грамма на языке C снабжена подробными комментариями. Задача студен- тов — разобраться в предлагаемом коде и на его основе выполнить задания.
Общие рекомендации к выполнению лабораторных работ
• Необходимо до начала моделирования провести обезразмеривание опи- сывающих задачу уравнений, выявить безразмерные комбинации па- раметров модели и дальнейшие действия производить в безразмерных величинах.
• Необходим контроль точности результатов и устойчивости применяе- мого численного метода.
• При выводе результатов в графической форме графики должны быть построены так, как это принято в научной литературе (с указанием того, какие величины отложены по осям, масштабами и т. д.).
• Отчет по лабораторной работе должен включать: название работы, ис- полнитель, группа; результаты тестирования программы; результаты,
полученные в ходе выполнения задания (графики, расчеты); анализ результатов.
1
Метод молекулярной динамики
1.1
Классическая молекулярная динамика
Метод молекулярной динамики (МД) — это метод, в котором временн´
ая эволюция системы взаимодействующих атомов или частиц отслеживается интегрированием их уравнений движения [1–5]. Для описания движения атомов или частиц применяется классическая механика.
При моделировании силы межатомного взаимодействия представляют в форме классических потенциальных сил (как градиент потенциальной энергии системы). Наборы конфигураций, получаемые в ходе расчетов ме- тодом молекулярной динамики, распределены в соответствии с некоторой
6

статистической функцией распределения, например, отвечающей микро- каноническому распределению. Заметим, что точное знание траекторий движения частиц системы на больших промежутках времени не является необходимым для получения результатов макроскопического (термодина- мического) характера.
Метод классической МД имеет границы применимости. Он не приме- ним, когда начинают играть роль квантовые эффекты (легкие атомы, низ- кие температуры и т. д.) Также необходимо, чтобы времена на которых рассматривается поведение системы, были больше, чем время релаксации исследуемых физических величин. Кроме того его не удается использовать для систем с очень большим количеством частиц. Это связано с тем, что время вычислений растет пропорционально квадрату количества частиц.
Метод молекулярной динамики, изначально разработанный в теорети- ческой физике, получил большое распространение в науке о веществе, а также в биохимии и биофизике. Он играет важную роль в определении структуры белка и уточнении его свойств.
1.2
Потенциал межмолекулярного взаимодействия
Определим модель рассматриваемой нами системы. Будем считать ди- намику системы классической, а молекулы — химически инертными ша- риками. Сила взаимодействия любых двух молекул зависит только от рас- стояния между ними. Тогда полная потенциальная энергия системы U ,
состоящей из N частиц, определяется суммой энергий двухчастичных вза- имодействий
U = u(r
12
) + u(r
13
) + u(r
23
) + . . . =
N
X
iu(r ij
),
(1)
где u(r ij
) — энергия взаимодействия двух частиц с номерами i и j, r ij
=
|r i
− r j
| — расстояние между этими частицами. Взаимодействие вида (1)
соответствует газам и некоторым жидкостям.
Выражение для u(r) можно построить на основе квантовомеханического расчета. Однако нам будет достаточно выбрать в качестве u(r) простую феноменологическую формулу. Будем исходить из того, что на малых r
7

частицы испытывают сильное отталкивание, а на больших расстояниях —
слабое притяжение. Одной из наиболее употребительных формул для u(r)
является потенциал Леннард-Джонса [2, 6]
u(r) = 4ε


σ
r

12
−

σ
r

6

(2)
Здесь первое слагаемое в скобках описывает отталкивание, а второе — при- тяжение. Слагаемое, описывающее притяжение, может быть получено в результате теоретических расчетов, а вклад, описывающий отталкивание,
выбран из соображений удобства численных расчетов. График этого по- тенциала показан на Рис. 1.
r/
u/
Рис. 1: Потенциал Леннард-Джонса. Здесь r приведено в единицах σ, а u(r) — в единицах ε.
Заметим, что потенциал Леннард-Джонса параметризуется “длиной” σ и
“энергией” ε. При r = σ потенциал u(r) = 0. Параметр ε представляет собой глубину потенциальной ямы в точке минимума u(r), достигаемого при r =
2 1/6
σ. Заметим, что данный потенциал является короткодействующим и u(r) практически равно 0 при r > 2.5σ.
1.3
Уравнения движения
В рассматриваемой нами системе частиц гамильтониан имеет вид
H =
N
X
i=1
mv
2
i
2
+ U,
(3)
8

где m — масса частицы, v i
— модуль вектора скорости v i
частицы с номером i, U определено в уравнении (1). Уравнения движения можно получить из гамильтониана (3). Они имеют вид dr i
dt
= v i
,
m dv i
dt
= f i
=
N
X
j=1,j6=i
F(r ij
).
(4)
Здесь t — время, f i
— полная сила, действующая на частицу с номером i со стороны других частиц системы, F(r ij
) — сила, действующая на частицу с номером i со стороны частицы с номером j
F(r ij
) = −∇
i u(r ij
) = −
∂u(r ij
)
∂r ij r
i
− r j
r ij
(5)
1.3.1
Численное решение уравнений движения
Численное решение системы уравнений (4) несложно построить, если известны координаты и скорости всех частиц в начальный момент времени.
Для этого достаточно написать разложение координат и скоростей в ряд
Тейлора и использовать уравнения (4)
r i
(t + ∆t) = r i
(t) + v i
(t)∆t +
1 2
a i
(t)∆t
2
v i
(t + ∆t) = v i
(t) + a i
(t)∆t,
(6)
где a i
= f i
/m — ускорение частицы с номером i. Подобные формулы мож- но использовать, если все отброшенные члены ряда малы по сравнению с учтенными.
Однако не все численные методы МД будут одинаково хорошо работать на большим временах. Внешне это будет проявляться как несохранение од- ной или нескольких физических величин (например, энергии или импульса системы). Классическим примером плохого метода с быстрорастущей по- грешностью является метод Эйлера [2]. Мы воспользуемся для численного решения системы (4) алгоритмом Верле в скоростной форме [2]
r i
(t + ∆t) = r i
(t) + v i
(t)∆t +
1 2
a i
(t)∆t
2
v i
(t + ∆t) = v i
(t) +
1 2
[a i
(t) + a i
(t + ∆t)]∆t.
(7)
9

То есть, зная значения координат, скоростей и ускорений частиц в момент времени t мы сначала находим координаты в момент времени t+∆t (первое уравнение (7)), затем с помощью выражения (5) рассчитываем ускорения в момент времени t + ∆t по найденным уже координатам и, наконец, с помощью второго уравнения (7) мы находим скорости в момент времени t + ∆t.
1.4
Граничные условия
Будем считать, что система представляет собой параллелепипед со сто- ронами L
x
, L
y и L
z и ограничена плоскостями α = 0 и α = L
α
, где
α = x, y, z. Такой выбор объема системы удобен для численных расчетов и не накладывает никаких дополнительных ограничений в случае жидкостей и газов. Мы рассмотрим три типа систем с постоянным числом частиц N :
• замкнутая система объема V , с заданным значением энергии E (мик- роканонический ансамбль [7]);
• система объема V , находящаяся в термостате с температурой T
w
(ка- нонический ансамбль [7]);
• система находящаяся под действием внешнего давления P
ext
(канони- ческий ансамбль);
и для каждой из них опишем соответствующие граничные условия.
1.4.1
Граничные условия для замкнутой системы с заданным значением энергии
Для такого типа систем наиболее типично использование периодических граничных условий (ПГУ, [1]). Эти условия, фактически, устремляют объ- ем системы к бесконечности, что приводит к уменьшению влияния ограни- чивающих поверхностей на поведение системы.
В этом разделе для простоты будем считать, что L
x
= L
y
= L
z
≡ L.
Таким образом, система представляет из себя кубическую ячейку. Пред- ставим, что эта кубическая ячейка вместе со всеми находящимися в ней
10

частицами повторяется в пространстве бесконечное количество раз. Мате- матически это означает, что для любой наблюдаемой величины Ψ выпол- няется соотношение
Ψ(r) = Ψ(r + Ln),
n = (n x
, n y
, n z
),
(8)
где n x
, n y
и n z
— любые целые числа. Фактически, это означает, что если частица вылетела за границу исходной ячейки, то такая же частица вле- тела через противоположную стенку с той же скоростью. Заметим, что в случае ПГУ изменится выражение (1) для потенциальной энергии. Теперь мы должны учесть взаимодействия частиц исходной клетки с частицами во всех остальных клетках
U =
N
X
iu(r ij
) +
X
n6=0
N
X
i,j=1
u(|r i
− (r j
+ Ln)|).
(9)
Используемый нами потенциал (2) короткодействующий, поэтому в сум- ме (9) мы ограничимся только ближайшими к данной частице соседями,
вне зависимости от того, лежат они в исходной клетке или нет. То есть, под расстоянием r ij между частицами с номерами i и j мы будем понимать r
ij
= min n
|r i
− r j
+ Ln|
и пользоваться прежней формулой (1). Тогда частица в исходной клетке взаимодействует только с N −1 частицами, причем ровно по одной частице для каждого j, лежащими в исходной клетке или же в ближайших соседях исходной клетки. Фактически, мы обрезали потенциал U на расстояниях r
c
>
√
3L/2.
При выполнении численного моделирования необходимо проверять вы- лет частицы за пределы ячейки и возвращать частицу назад в ячейку в соответствии с ПГУ.
1.4.2
Граничные условия для системы с постоянной температурой
Рассмотрим систему объема V , состоящую из N частиц, которая нахо- дится в контакте с термостатом. Будем предполагать что термостат поддер- живает постоянную температуру T
w стенок системы. Периодические гра- ничные условия, изложенные выше, не учитывают взаимодействие систе- мы с термостатом, поэтому мы будем использовать другой вид граничных
11

условий. Будем считать, что частица при ударе о стенку прилипает к ней,
приобретает температуру стенки, а затем отлетает от стенки в произволь- ном направлении. В этом случае модуль вектора скорости отлетающей ча- стицы может быть найден из условия mv
2
i
2
=
3 2
k
B
T
w
,
(10)
При выполнении численного моделирования необходимо выявлять фак- ты столкновения частиц со стенками термостата, вычислять новые скоро- сти столкнувшихся со стенкой частиц в соответствии с формулой (10), а сами частицы размещать на поверхности той стенки, с которой произошло столкновение.
1.4.3
Граничные условия для системы, находящейся под постоянным давлением
Теперь рассмотрим случай системы находящейся под постоянным дав- лением. Для простоты рассмотрим систему, одна из стенок которой по- движна (поршень двигающийся вдоль оси z). Будем считать, что поршень находится под действием постоянного внешнего давления P
ext
(Рис. 2). В
этом случае при столкновениях частиц с поршнем будет происходить об- мен энергией, и это будет вызывать изменение объема системы. Опишем взаимодействие частиц с поршнем следующим образом:
∆V =L
x
L
y
∆L
z
,
∆L
z
=v p
∆t,
∆v iz
= − 2(v iz
− v p
),
(11)
где ∆V — изменение объема системы, ∆L
z
— изменение линейного размера системы вдоль оси z, ∆v iz
— изменение проекции v iz скорости частицы при столкновении с поршнем, а v p
— скорость движения поршня. Очевидно, что скорость движения поршня будет зависеть от разности давлений внутри, P ,
и снаружи, P
ext
, системы, а также от массы m p
поршня. Несложно написать следующую модельную формулу для скорости движения поршня:
v p
= −(P
ext
− P )
L
x
L
y
∆t m
p
(12)
12

Здесь масса m p
является подгоночным параметром, который следует выби- рать так, чтобы v p
∆t  σ. С другой стороны, слишком большие значения величины m p
невыгодны с вычислительной точки зрения.
Рис. 2: Система, находящаяся под постоянным внешним давлением P
ext
. Верхняя стенка си- стемы представляет собой поршень.
Нам осталось найти давление P внутри системы. Для этого мы будем использовать теорему вириала [2]
P V =
1 3
N
X
i=1
(mv
2
i
+ r i
· f i
),
(13)
где f i
— полная сила, действующая на частицу с номером i со стороны всех остальных частиц.
При выполнении численного моделирования необходимо выявлять фак- ты столкновения частиц со стенками и поршнем. В случае столкновения частиц со стенками могут быть выбраны условия упругого отражения ча- стиц от стенки, а при столкновении частиц с поршнем следует использовать выражения (11).
Для определенности будем вычислять новое положение поршня перед выполнением каждого шага по времени. Для этого нам нужно найти те- кущее значение давления P , скорости поршня v p
и координату нового по- ложения поршня L
z
. Тогда при выполнении шага по времени мы можем просчитать столкновения со стенками и поршнем и выполнить необходи-
13

мые изменения координат и скоростей частиц. В случае столкновения ча- стиц с поршнем обычно достаточно разместить столкнувшиеся частицы на поверхности поршня.
1.5
Вычисление средних
Метод молекулярной динамики позволяет найти нам мгновенные значе- ния микроскопических величин системы, таких как координаты и скорости частиц. С помощью микроскопических величин можно вычислить мгно- венные значения макроскопических величин системы, например: энергии,
полного импульса, момента инерции. Однако мгновенные значения редко представляют практический интерес, поскольку в экспериментах наблю- даются средние значения физических величин. В статистической физике различают два вида средних: среднее по ансамблю и среднее по времени [7].
В соответствии с эргодической гипотезой [7] среднее по времени совпадает со средним по ансамблю.
Важной характеристикой системы является температура. Она может быть выражена через среднюю кинетическую энергию частиц системы по формуле
*
N
X
i=1
mv
2
i
2
+
=
s
2
k
B
T,
(14)
где h· · · i обозначают усреднение по времени или по ансамблю, а s — число степеней свободы системы. В случае замкнутой системы s = 3(N − 1), а в случае систем с постоянной температурой или давлением s = 3N .
В молекулярной динамике мы изучаем эволюцию во времени одной и той же системы, поэтому для вычисления средних значений физических величин воспользуемся усреднением по времени. Так, среднее физической величины Ψ за k шагов моделирования будем искать по формуле hΨi k
=
1
k∆t
Z
k∆t
0
Ψ(t)dt ≈
1
k k−1
X
j=0
Ψ(j∆t) =
1
k k−1
X
j=0
Ψ
j
(15)
14

1.6
Начальные условия
Нам осталось задать координаты и скорости частиц в начальный мо- мент времени. Будем считать, что в начальный момент времени частицы распределены по кубической клетке равномерно, то есть, система разбита на кубики с объемом V /N , и в центре каждого кубика находится частица.
Скорости частиц определим следующим образом. Зададим максималь- ное значение проекции скорости частиц v max в начальный момент времени.
Определим сначала проекции скоростей всех частиц случайным образом из интервала [−v max
, v max
]. При этом распределение скоростей в началь- ный момент времени будет равномерным. Затем сосчитаем полный импульс системы по формуле p =
N
X
i=1
mv i
Поскольку суммарный импульс системы в начальный момент времени дол- жен обращаться в ноль, вычтем из скорости каждой частицы в начальный момент времени p/N . Полученные скорости мы и будем использовать в качестве начальных. Заметим, что в сложных задачах, когда, например,
система очень медленно приходит в состояние равновесия, используют “сту- пенчатое” задание начальных условий. Сначала задают начальные условия относительно произвольно. Затем запускают моделирование. В какой-то момент делают снимок координат и скоростей всех частиц системы. Этот снимок берут в качестве начальных условий для следующего этапа моде- лирования. Эти шаги выполняют несколько раз.
1.7
Обезразмеривание
В задачах моделирования крайне важно предварительно провести про- цедуру обезразмеривания, то есть, переписать все уравнения, описывающие эволюцию системы в безразмерном виде. Кроме того процедура обезразме- ривания часто позволяет выявить характерные масштабы в исследуемых системах, определить границы применимости тех или иных методов. Также это крайне удобно с точки зрения вычислений: не приходится задумывать- ся, с величинами какого порядка мы работаем, что с чем сравниваем, и мал
15

ли конкретный параметр или велик.
Обезразмеривание преследует две цели: (1) упростить уравнения, ис- ключив из них “лишние” константы и величины; (2) определить фактиче- ское количество и набор параметров, определяющих динамику системы.
Проведем обезразмеривание в нашей модели.
В нашей задаче размерными величинами являются параметры потенци- ала (2): σ — длина и ε — энергия. Еще одной размерной величиной является масса частицы m. Объем системы V имеет такую же размерность, как и σ
3
,
и нет необходимости включать его в список величин, с помощью которых мы будем обезразмеривать задачу. Аналогично для микроканонического ансамбля суммарная энергия системы имеет такую же размерность, что и
ε. Таким образом, нам достаточно использовать три величины σ, ε и m,
чтобы обезразмерить уравнения движения. На основе σ, ε и m мы можем строить любые размерности, встречающиеся в нашей задаче.
Если все величины, имеющие размерность длины, поделить на σ, то мы получим безразмерные расстояния. Аналогичным образом энергетические характеристики системы можно поделить на ε и получить безразмерные значения. Из выражения для кинетической энергии несложно получить,
что размерность времени будет иметь величина τ = σ(m/ε)
1/2
. Тогда без- размерное время будет иметь вид t
∗
= t/τ . Здесь и далее безразмерные величины мы будем обозначать звездочкой. Уравнения движения (4) в без- размерной форме примут вид dr
∗
i dt
∗
= v
∗
i
,
dv
∗
i dt
∗
= a
∗
i
(16)
Аналогичным образом можно переписать уравнения для моделирования
(7) и обезразмерить начальные условия. Обратим внимание, что наша зада- ча стала полностью безразмерной, в ней не осталось ни одного размерного параметра. Именно обезразмеренные уравнения мы и будем использовать при моделировании.
16

Лабораторная работа 1.1. Моделирование динамики замкнутой системы многих частиц
Целью настоящей работы является исследование поведения замкнутой системы многих частиц методом молекулярной динамики.
Содержание работы
1. Выпишите математическую модель для моделирования динамики си- стемы. Проведите обезразмеривание системы уравнений. Определите какие безразмерные величины являются начальными данными, а ка- кие параметрами математической модели.
2. Оцените при каких значениях начальных скоростей частиц в резуль- тате моделирования получится газ, а при каких значениях могут об- разовываться капли. Оценку выполните из тех соображений, что в газообразном состоянии кинетическая энергия, приходящаяся на одну частицу, должна быть больше, чем минимальное значение потенциаль- ной энергии u(r ij
), взятое по модулю. Полученную оценку используй- те для выбора значения параметра v max
. Задайте v max много меньше полученной оценки и сравните значения потенциальной и кинетиче- ской энергий системы по прошествии некоторого времени, сделайте выводы. Какие еще условия должны соблюдаться, чтобы в результате моделирования получался газ?
3. Изучите реализацию метода молекулярной динамики (код программы прилагается). Пояснения к коду программы смотрите в Приложении к разделу 1.
4. Изучите распределение скоростей частиц в газообразной фазе. Напи- шите функцию для расчета средних моментов скоростей. Средний мо- мент n-того порядка проекции скорости на ось x вычисляется по фор- муле hM
(n)
x i
k
=
*
N
X
i=1
v n
ix
+
k
,
17

Проведите усреднение при k = 100. Сравните полученные результаты с первыми пятью моментами для распределения Максвелла.
5. Исследуйте, как быстро система приходит в состояние равновесия в зависимости от значений N , v max и ∆t. Критерием прихода системы к равновесию считайте установление постоянной температуры, которая связана с кинетической энергией системы соотношением (14). Объяс- ните, почему в (14) для рассматриваемой системы s = 3(N − 1).
6. Убедитесь, что полная энергия системы при моделировании сохраня- ется. Для этого изучите отклонение текущего значения полной энергии от его начального значения. Будем считать, что энергия сохраняет- ся, если модуль отклонения текущего значения энергии от начального много меньше начального значения энергии.
7. Напишите функцию для расчета радиальной функции распределения g(r). Радиальная функция распределения позволяет определить веро- ятность нахождения частицы на расстоянии r от заданной частицы.
Радиальная функция распределения имеет вид g(r) =
V
N
n(r)
4πr
2
,
где n(r)∆r — среднее число частиц, расположенных на расстояниях от r до r + ∆r от заданной частицы. Радиальную функцию несложно получить, поскольку в каждый момент времени известны координа- ты всех частиц. Для нахождения n(r) разобьем диапазон возможных расстояний между частицами r на интервалы длиной ∆r и посчитаем среднее число частиц, попадающих в каждый интервал. Это число ча- стиц и будет как раз n(r)∆r. Усреднение n(r) выполняйте по k = 100
шагам по времени.
8. Получите радиальную функцию в газообразной фазе и в случае нали- чия капель. Сделайте выводы. На каких расстояниях r имеет смысл рассчитывать радиальную функцию распределения в нашей модели?
18

Лабораторная работа 1.2. Моделирование динамики системы многих частиц при постоянной температуре или давлении
Целью настоящей работы является исследование поведения системы многих частиц при постоянной температуре или давлении методом моле- кулярной динамики.
Содержание работы
1. Выпишите математическую модель для моделирования динамики си- стемы. Проведите обезразмеривание системы уравнений. Определите какие безразмерные величины являются начальными данными, а ка- кие параметрами математической модели.
2. Оцените при каких значениях начальных скоростей частиц в резуль- тате моделирования получится газ, а при каких значениях могут об- разовываться капли. Оценку выполните из тех соображений, что в газообразном состоянии кинетическая энергия, приходящаяся на одну частицу, должна быть больше, чем минимальное значение потенциаль- ной энергии u(r ij
), взятое по модулю. Полученную оценку используй- те для выбора значения температуры стенки T
w
, исходя из выраже- ния (14). Задайте T
w много меньше полученной оценки и сравните зна- чения потенциальной и кинетической энергий системы по прошествии некоторого времени, сделайте выводы. Какие еще условия должны со- блюдаться, чтобы в результате моделирования получался газ?
3. Изучите реализацию метода молекулярной динамики (код программы прилагается). Пояснения к коду программы смотрите в Приложении к разделу 1.
4. Изучите распределение скоростей частиц в газообразной фазе. Напи- шите функцию для расчета средних моментов скоростей. Средний мо- мент n-того порядка проекции скорости на ось x вычисляется по фор- муле hM
(n)
x i
k
=
*
N
X
i=1
v n
ix
+
k
,
19

Проведите усреднение при k = 100. Сравните полученные результаты с первыми пятью моментами распределения Максвелла.
5. Исследуйте, как быстро система приходит в состояние равновесия в зависимости от значений N , v max
, T
w и ∆t. Критерием прихода си- стемы к равновесию считайте установление постоянной температуры,
которая связана с кинетической энергией системы соотношением (14).
Объясните, почему в (14) для рассматриваемой системы s = 3N .
6. Напишите функции для расчета средней энергии системы и средне- квадратичного отклонения энергии системы. Убедитесь, что средняя энергия системы сохраняется после прихода системы к равновесию.
Усреднение проводите по 100 шагам моделирования.
7. Модифицируйте программу для случая системы при постоянном внешнем давлении. Используйте рекомендации из раздела 1.4.3. Убе- дитесь в работоспособности алгоритма. Проверьте, что система при- ходит в равновесие.
Приложение. Комментарии к коду программ
Прилагаемая программа на языке С снабжена подробными коммента- риями. Программа реализует моделирование методом молекулярной дина- мики в безразмерных переменных.
Сначала приведены функции, вычисляющие силу, действующую на ча- стицу, кинетическую и потенциальную энергию системы. Особое внимание следует обратить на функцию time_step, реализующую расчет координат,
скоростей и ускорений для одного шага по времени. Далее следуют функ- ции, которые собирают данные для построения распределения по скоро- стям. Эти функции можно не изучать подробно.
Важная функция initialize задает начальные координаты и скорости частиц. В этой функции следует разобраться подробно.
Далее идет вспомогательный блок, необходимый для графического отоб- ражения распределения по координатам.
20

Наконец, основной является функция main, где собственно и происходит выполнение моделирования.
При запуске программы следует ввести число частиц N , максимальное значение проекции скорости v max и шаг по времени ∆t. В случае работы 1.2.
также потребуется ввести температуру T
w
. Начнется процесс моделирова- ния. При нажатии комбинации клавиш Ctrl+C, вычисления приостановят- ся. Далее следуйте инструкциям на экране.
2
Метод Монте-Карло
2.1
Вычисление интегралов простейшим методом Монте-Карло
Представим себе прямоугольник высотой h и длиной b − a такой, что функция f (x) лежит внутри него (Рис. 3). Мы хотим сосчитать
F =
Z
b a
f (x)dx.
(17)
a
b
h
f(x)
Рис. 3:
Генерируем n пар случайных равномерно распределенных чисел x l
и y
l
, удовлетворяющих условиям a ≤ x l
≤ b и 0 ≤ y l
≤ h. Доля точек
(x l
, y l
), которые удовлетворяют условию y l
≤ f (x l
), представляет собой оценку отношения интеграла от функции f (x) к площади прямоугольника.
Отсюда оценка F
n определяется выражением
F
n
= A
n s
n
(18)
21

где n s
— число точек, лежащих под кривой, n — общее количество точек,
а A — площадь прямоугольника.
Другая разновидность метода Монте-Карло основывается на теореме,
гласящей, что интеграл (18) определяется средним значением подынте- гральной функции f (x) на отрезке a ≤ x ≤ b. Для вычисления этого сред- него возьмем x l
случайным образом и произведем выборку значений f (x).
Оценка F
n одномерного интеграла (18) методом “выборочного среднего”
выражается формулой
F
n
= (b − a)hf i =
b − a n
n
X
l=1
f (x l
),
(19)
где x l
случайные числа, равномерно распределенные на отрезке a ≤ x l
≤
b, а n — количество испытаний. Оказалось, что формула (19) особенно эффективна при расчете многомерных интегралов.
Может показаться странным, что с помощью детерминированного ком- пьютера удается генерировать последовательности “случайных” чисел. Со- здание качественного генератора случайных чисел представляет собой до- статочно сложную задачу. Пока нас вполне устроит использование “псев- дослучайных” последовательностей, генераторы которых предусмотрены в языках программирования.
Годом рождения метода Монте-Карло считается 1949 год, когда в свет выходит статья Метрополиса и Улама “Метод Монте-Карло” [8]. Название метода происходит от названия города в княжестве Монако, широко извест- ного своими казино, поскольку именно рулетка является одним из самых широко известных генераторов случайных чисел. Название было предло- жено Метрополисом в честь его дяди, который был азартным игроком.
2.2
Анализ погрешности метода Монте-Карло
Одной возможной мерой погрешности является дисперсия ∆
2
, опреде- ляемая выражением
∆
2
= hf
2
i − hf i
2
,
(20)
22

где hf i =
1
n n
X
l=1
f (x l
)
(21)
и hf
2
i =
1
n n
X
l=1
f (x l
)
2
(22)
Величина ∆
2
называется стандартным отклонением. Если бы функция f не зависела от x, то ∆ равнялось бы нулю. Одним из условий, которому должна удовлетворять мера погрешности, является ее убывание с ростом n. Один из способов уменьшения погрешности заключается в проведении серии дополнительных расчетов по n испытаний в каждом. Каждый такой расчет дает среднее значение, или измерение, которое мы обозначим Ψ
α
Эти измерения, вообще говоря, не будут одинаковыми, поскольку каждое измерение производится со своей последовательностью случайных чисел.
Качественно величина разброса измерений служит мерой погрешности од- ного измерения. Предположим, имеется набор из j измерений {Ψ
α
}, со- стоящих из одинакового числа испытаний. Удобной мерой разброса этих измерений является стандартное отклонение средних ∆
j
, которое мы опре- делим как
∆
2
j
= hΨ
2
i − hΨi
2
,
(23)
где hΨi =
1
j j
X
α=1
Ψ
α
=
1
jn j
X
α=1
n
X
l=1
Ψ(x lα
)
(24)
и hΨ
2
i =
1
j j
X
α=1
Ψ
2
α
(25)
Можно аналитически получить оценку
∆
j
= ∆/
√
n.
(26)
Таким образом, использование статистически независимых серий позволя- ет уменьшить погрешность при выполнении расчетов.
23

2.3
Метод Монте-Карло и канонический ансамбль
Используем теперь метод Монте-Карло для моделирования равновесных свойств систем с многими степенями свободы. Рассмотрим сначала систему,
у которой число частиц N и объем V фиксированы. Будем считать, что си- стема находится в тепловом контакте с окружающей средой. За счет этого теплового контакта осуществляется обмен энергией между лабораторной системой и окружающей средой в виде тепла. Обычно исследуемая систе- ма мала по сравнению с окружающей средой. Большая система, имеющая намного больше степеней свободы, называется резервуаром или термоста- том. Предположим, кроме того, что влиянием внешних параметров, таких,
как гравитационные и магнитные поля, можно пренебречь. Заметим, что условие постоянства полной энергии здесь относится к составной системе,
состоящей из исследуемой системы и резервуара, а энергия исследуемой системы может меняться. Представим себе большое число воображаемых копий исследуемой системы и резервуара. Рассматриваемые как единое це- лое, исследуемая система вместе с резервуаром являются изолированными и могут быть описаны микроканоническим ансамблем. Поскольку нас ин- тересуют равновесные значения физических величин, описывающих иссле- дуемую систему, желательно знать вероятность P
s
, с которой исследуемая система обнаруживается в состоянии s с энергией E
s
. Ансамбль, который описывает распределение вероятностей состояний исследуемой системы,
находящейся в тепловом равновесии с резервуаром, называется канони- ческим. Вообще, в качестве исследуемой системы может выступать любая макроскопическая система, которая гораздо меньше резервуара. Исследу- емой системой может быть столь малый объект, как отдельная частица,
если ее можно надежно выявить среди частиц резервуара.
Вероятность P
s в каноническом ансамбле описывается распределением
Гиббса [7]
P
s
=
1
Z
exp (−βE
s
) ,
(27)
где β = 1/(k
B
T ), T — температура, Z — нормировочный множитель (ста- тистическая сумма), выбираемый из условия равенства единице суммы по всем состояниям системы. Распределение вероятностей (27) также назы-
24

вается каноническим распределением. Отметим, что распределение Гибб- са характеризуется температурой. Если предположить, что формула (27)
применима для любой исследуемой системы, находящейся в тепловом рав- новесии с резервуаром, то мы увидим, что в каноническом ансамбле вся- кое макросостояние определяется параметрами T , N и V . По сравнению с этим в микроканоническом ансамбле макросостояние характеризуется ве- личинами E, N и V .
Для классических систем выражение, в точности соответствующее фор- муле (27), получается для функции распределения в фазовом пространстве
ρ(p , q) =
1
Z
exp (−βE(p , q)) ,
(28)
где q и p обозначают совокупность всех координат и импульсов частиц системы, нормировочный множитель Z определяется из условия
1
Z
Z
exp (−βE(p , q)) dpdq = 1.
Среднее значение любой физической величины может быть вычислено с помощью распределения Гиббса по формуле hΨi =
1
Z
Z
Ψ(p , q) exp (−βE(p , q)) dpdq.
(29)
Задачу вычисления средних можно упростить, если учесть, что потенциал взаимодействия между частицами не зависит от импульсов частиц. Поэто- му здесь можно выделить интегралы по импульсам и сосчитать их с по- мощью распределения Максвелла. Заметим, что для определения средних значений в системе из N частиц требуется вычислять 6N мерный инте- грал или, после выделения интегралов по импульсам, 3N мерный. Прямой численный расчет такого объекта может быть весьма затруднительным.
Мы воспользуемся для расчета методом Монте-Карло, который, как было рассмотрено выше, состоит в интегрировании по случайной выборке точек,
вместо интегрирования по регулярной решетке.
Как и прежде будем считать, что система представляет собой парал- лелепипед со сторонами L
x
, L
y и L
z
, взаимодействие между частицами описывается с помощью потенциала Леннард-Джонса (2) и на систему на- ложены периодические граничные условия. В этом случае скорости частиц
25

непрерывны. Энергия любой конфигурации будет зависеть от положений частиц, и полная энергия есть сумма кинетических и потенциальных энер- гий отдельных частиц.
Казалось бы, наиболее простой способ выполнения интегрирования со- стоит в том, чтобы поместить каждую из N частиц в случайное положение внутри системы. Затем можно сосчитать энергию получившейся конфигу- рации и придать ей вес exp (−βE(p , q)). Задавая таким образом конфигу- рации многократно, мы сможем сосчитать интеграл с помощью формулы вида (19). Однако такой способ оказывается крайне невыгодным с вычис- лительной точки зрения, поскольку в подавляющем большинстве случаев вес exp (−βE(p , q)) случайной конфигурации будет очень мал. Рассмотрим модификацию метода Монте-Карло, в которой вместо случайного выбора конфигураций и придания им веса exp (−βE(p , q)), мы будем выбирать конфигурации с вероятностью exp (−βE(p , q)) и взвешивать их равномер- но (с одинаковым весом).
Поместим N частиц в любую конфигурацию, например, на регулярную решетку. Затем выберем случайным образом одну из частиц и изменим ее координаты и скорость следующим образом x
i
→ x i
+ σδ,
y i
→ y i
+ σδ,
z i
→ z i
+ σδ,
v ix
→ v ix
+ v m
δ,
v iy
→ v iy
+ v m
δ,
v iz
→ v iz
+ v m
δ,
(30)
где v m
— характерная скорость частиц, например, среднеквадратичная ско- рость, δ — параметр корректировки координат и скоростей, лежащий в пределах от 0 до 1. Чем меньше этот параметр, тем слабее будет изменение текущей конфигурации системы.
Затем сосчитаем изменение энергии системы ∆E, вызванное смещением частицы. Если энергия системы уменьшилась или не изменилась, то есть,
∆E ≤ 0, то мы принимаем изменение конфигурации системы, присвоив ча- стице ее новое положение и скорость. Если же ∆E > 0, то мы принимаем изменение конфигурации с вероятностью W = exp (−β∆E). Это означает,
что если мы возьмем случайное число ξ, лежащее в пределах от 0 до 1, и если ξ ≤ W , то мы принимаем новые положение и скорость частицы. Если же ξ > W , то мы возвращаем частицу в ее исходное положение. Затем, в
26

независимости от того, приняли мы или нет новые положение и скорость частицы, мы считаем, что система перешла в новую конфигурацию и мож- но выполнять следующий шаг моделирования.
После выполнения достаточно большого числа n таких шагов мы можем вычислить средние значения интересующих нас величин по формуле hΨi =
1
n n
X
j=1
Ψ
j
,
(31)
где Ψ
j
— значение величины Ψ после выполнения j шагов.
Заметим, что при смещении частиц согласно (30) с отличной от нуля вероятностью, любая частица сможет достичь любой точки системы за до- статочно большое число шагов. А значит, мы сможем побывать во всех воз- можных конфигурациях системы. Конфигурации генерируются с вероят- ностью, пропорциональной требуемой вероятности, и поэтому все средние превращаются в арифметические средние, как в формуле (31). Оценить этим способом статистическую сумму Z невозможно. Роль переменной ξ
можно пояснить следующим образом. Величина W принимает значения от
0 до 1. Представим себе отрезок [0, 1], который разбили на 2 отрезка с дли- нами W и 1 − W . Затем на отрезок [0, 1] бросается случайное число ξ. В
зависимости от того, на какой из двух отрезков W и 1 − W попало число
ξ, и принимается решение о переходе системы в новое состояние.
Описанный выше алгоритм был предложен Метрополисом и др. [9]. В
задачах статистической механики выражения “метод Монте-Карло” и “ме- тод Метрополиса” — почти синонимы.
Сформулируем теперь по шагам алгоритм Метрополиса на примере си- стемы частиц.
2.3.1
Алгоритм Метрополиса для канонического ансамбля
1. Формируем начальную конфигурацию.
2. Производим случайное пробное изменение в конфигурации. Напри- мер, выбираем случайную частицу и пробуем переместить ее на слу- чайное расстояние и случайным образом изменить ее скорость.
27

3. Вычисляем ∆E, то есть, изменение энергии системы, обусловленное произведенным пробным изменением конфигурации.
4. Если ∆E меньше или равно нулю, то принимаем новую конфигурацию и переходим к шагу 8.
5. Если ∆E положительно, вычисляем “вероятность перехода”
W = exp(−β∆E).
6. Генерируем случайное число ξ в интервале (0, 1).
7. Если ξ ≤ W , то новую конфигурацию принимаем, в противном случае сохраняем предыдущую конфигурацию.
8. Определяем значения требуемых физических величин в текущей кон- фигурации.
9. Повторяем шаги 2–8 для получения достаточного числа n конфигура- ций или “испытаний”.
10. Вычисляем средние по конфигурациям по формуле (31).
11. Выполняем шаги 1–10 для получения статистически независимых се- рий измерений. Вычисляем средние значения по средним значениям серий по формуле (24).
2.3.2
Выборка по значимости
Алгоритм Метрополиса основана на так называемом методе выборки по значимости. Поясним этот метод на простом примере. Как мы видели ранее, оценка погрешности метода Монте-Карло пропорциональна ∆. По- этому желательно использовать метод вычисления, который уменьшает ∆
и повышает эффективность каждого испытания. В качестве примера этого метода введем положительную функцию p(x) такую, что
Z
b a
p(x)dx = 1.
(32)
28

Тогда интеграл (17) можно переписать в виде
F =
Z
b a
 f (x)
p(x)

p(x)dx.
(33)
Мы можем вычислить интеграл (33), производя выборку в соответствии с
“распределением вероятности” p(x) и конструируя сумму
F
n
=
1
n n
X
i=1
f (x i
)
p(x i
)
(34)
Фактически, это выражение означает, что мы считаем вклад в сумму слага- емых вида f (x i
)/p(x i
), но каждое слагаемое берется с весом p(x i
). Заметим,
что для однородного случая p(x) = 1/(b − a) и (34) переходит в (19). Мы хотим выбрать функцию p(x) так, чтобы дисперсия подынтегрального вы- ражения f (x)/p(x) была минимальна. Выбор функции p(x) представляет сложную проблему. Однако, в реальных физических задачах ее удается по- добрать, поскольку зачастую известны статистические свойства изучаемых систем. Удобно выбирать функцию p(x), которая ведет себя как f (x) там,
где f (x) велика. Если мы сумеем подобрать подходящую p(x), то подын- тегральное выражение будет медленно меняющейся функцией и, следова- тельно, дисперсия уменьшится.
2.4
Двумерная модель Изинга
Одной из простейших моделей, используемых в статистической физике для моделирования фазовых переходов в магнитных веществах или би- нарных составах, является модель Изинга. Данная модель относится к широкому классу решеточных моделей, в которых рассматриваются ло- кальные взаимодействия. В простейшем случае это взаимодействия между ближайшими узлами решетки. В магнитных системах локальные взаимо- действия обусловлены спинами, расположенными в узлах решетки. Спины могут представлять собой, например, магнитные моменты атомов в твер- дом теле, взаимодействующие друг с другом и внешним магнитным полем.
Модель была предложена Ленцем и исследована его дипломником Изин- гом с целью изучения фазового перехода из парамагнитного состояния в ферромагнитное. Изинг рассчитал термодинамические свойства модели в
29

одномерной постановке и нашел, что в ней фазовый переход отсутствует.
Однако в двумерном и трехмерном случаях модель Изинга действительно обнаруживает переход.
Рассмотрим решетку, содержащую N узлов, и предположим, что с каж- дым узлом решетки i связано число s i
, где s i
= +1, если спин ориентирован
“вверх” (или вдоль внешнего поля), и s i
= −1, если он ориентирован “вниз”.
Любая конкретная конфигурация, то есть, микросостояние решетки, зада- ется набором переменных {s
1
, s
2
, . . . , s
N
} для всех узлов решетки.
Полная энергия E при наличии магнитного поля H в модели Изинга равняется
E = −J
N
X
hi,ji,is i
s j
− H
N
X
i=1
s i
,
(35)
где первая сумма в (35) берется по всем ближайшим соседним парам спинов hi, ji, а вторая сумма — по всем спинам решетки. Константа обменного вза- имодействия J является мерой силы взаимодействия между ближайшими соседними спинами. Если J > 0, то состояния ↑↑ и ↓↓. которые характери- зуются одинаковой ориентацией спинов ближайших соседей, энергетически выгоднее по сравнению с состояниями ↑↓ и ↓↑, у которых соседние спины ориентированы в противоположные стороны.
E = + J
E = - J
Рис. 4: Энергия взаимодействия между ближайшими соседними спинами в отсутствии внеш- него магнитного поля
Следовательно, можно ожидать, что для J > 0 состояние с наименьшей полной энергией является ферромагнитным, то есть, в среднем суммарное число спинов, ориентированных в одном направлении, не равно нулю. Если
J < 0, предпочтительнее состояния ↑↓ и ↓↑, для которых соседние спины антипараллельны, и можно ожидать, что состояние с наименьшей энергией является антиферромагнитным, то есть, спины упорядочены через один.
Если наложить внешнее магнитное поле H, направленное вверх, то спины
↑ и ↓ приобретают дополнительную внутреннюю энергию, равную −H и
30

+H соответственно. Отметим, что H измеряется в таких единицах, что магнитный момент на спин равен 1.
В основу модели Изинга положен ряд упрощающих предположений.
Пренебрегается кинетической энергией атомов, связанных с узлами решет- ки. В энергии взаимодействия учитывается вклад только ближайших со- седей и предусматривается только два дискретных состояния для спинов.
Наибольшее распространение для спиновых систем Изинга получила ди- намика “опрокидывания спина”. В этой динамике спин выбирается случай- ным образом и пробное изменение (испытание Монте-Карло) соответствует опрокидыванию спина из состояния ↑ в ↓ или из ↓ в ↑.
В этом разделе мы рассмотрим моделирование ферромагнетиков. Фер- ромагнетики проявляют спонтанную намагниченность, то есть, намагни- ченность при отсутствии внешнего магнитного поля. Эта ненулевая на- магниченность присутствует только при температуре ниже вполне опреде- ленной температуры T
c
, называемой температурой Кюри или критической температурой. При температурах T > T
c намагниченность пропадает. Тем самым T
c отделяет хаотическую фазу при T > T
c от ферромагнитной фазы при T < T
c
Механизм ферромагнетизма имеет квантовомеханическую природу. Не- смотря на это, изучение классической модели Изинга в двумерной и трех- мерной постановках полезно для понимания свойств магнитных систем в окрестности фазового перехода. Однако в силу своего классического харак- тера и пренебрежения другими спиновыми компонентами модель Изинга не дает полного описания ферромагнетизма, особенно при температурах значительно ниже T
c
. В частности, в модели Изинга предполагается, что отдельные моменты локализованы, и поэтому модель Изинга не применима к металлам, таким, как железо и никель. Для исследования свойств моде- ли Изинга нам необходимо конкретизировать, какие физические свойства представляют интерес, и разработать программу их вычисления. К рас- сматриваемым равновесным характеристикам относятся средняя энергия hEi, средняя намагниченность hM i, теплоемкость C и магнитная воспри- имчивость χ.
Суммарный магнитный момент, или намагниченность M , определяется
31

формулой
M =
N
X
i=1
s i
(36)
Интерес представляют средние значения величины hM i и среднеквадра- тичное отклонение hM
2
i − hM i
2
как функции температуры системы и на- ложенного магнитного поля.
Один из методов измерения C при постоянном внешнем магнитном поле вытекает из определения
C =
∂hEi
∂T
(37)
Другой метод измерения С основан на использовании связи теплоемкости со статистическими флуктуациями полной энергии в каноническом ансам- бле
C =
1
k
B
T
2
(hE
2
i − hEi
2
).
(38)
Магнитная восприимчивость χ является другим примером “функции от- клика”, поскольку она характеризует способность системы откликаться или опрокидываться за счет изменения во внешнем магнитном поле. Изотер- мическая магнитная восприимчивость при нулевом поле определяется по- средством термодинамической производной
χ = lim
H→0
∂hM i
∂H
(39)
Восприимчивость при нулевом поле можно связать с флуктуациями намаг- ниченности в системе:
χ =
1
k
B
T
(hM
2
i − hM i
2
),
(40)
где hM i и hM
2
i отвечают нулевому полю. Соотношения (38) и (40) явля- ются примерами общего соотношения между функциями отклика и равно- весными флуктуациями.
Поскольку нас интересуют свойства бесконечной системы, нужно учесть краевые условия. В качестве краевого условия мы выберем периодические граничные условия. Фактически, это означает, что в энергию (35) требуется включить дополнительные члены, описывающие взаимодействие узлов с номерами 1 и N .
32

Теперь, когда мы конкретизировали некоторые интересующие нас рав- новесные величины, нам надо реализовать алгоритм Метрополиса. В от- личие от предыдущего раздела, описываемая система не является класси- ческой. Поэтому здесь мы будем пользоваться распределением (27), а не
(28). Сам же алгоритм Метрополиса совпадает с приведенным в предыду- щем разделе. Отличие состоит только в том, что на втором шаге вместо изменения координат и скорости случайно выбранной частицы мы будем изменять ее спин на противоположный.
2.5
Микроканонический ансамбль
Рассмотрим теперь замкнутую систему, у которой число частиц N , объ- ем V и полная энергия E фиксированы. Предположим, кроме того, что система изолированная, то есть, влиянием внешних параметров можно пре- небречь. В общем случае замкнутая макроскопическая система стремится перейти в стационарное равновесное состояние с максимальным беспоряд- ком, или энтропией. Макросостояние системы характеризуется величина- ми N , V и E. На микроскопическом уровне существует в общем случае громадное число различных способов, иначе говоря конфигураций, в ко- торых может реализовываться данное макросостояние (N, V, E). Каждая конкретная конфигурация, или микросостояние, является достижимой, ес- ли его характеристики соответствуют данному макросостоянию.
Раз у нас нет никаких причин предпочесть одно микросостояние дру- гому, разумно постулировать, что в любой данный момент времени систе- ма с равной вероятностью может оказаться в любом из своих достижимых микросостояний. Чтобы точнее выразить этот постулат равных априорных вероятностей, представим себе изолированную систему с Ω достижимыми состояниями. Вероятность найти систему в микросостоянии s равна
P
s
=
(
1/Ω,
если s достижимо,
0
в противном случае.
(41)
Сумма P
s по всем Ω равна единице. Средние от физических величин можно определять двумя способами. В обычном лабораторном эксперименте из- мерение физических величин производят в течение достаточно большого
33

промежутка времени, чтобы в системе успело реализоваться большое число ее достижимых микросостояний. С точки зрения такого усреднения по вре- мени смысл вероятностей в формуле (41) заключается в том, что величина
P
s дает долю времени, которую одна система за время последовательности наблюдений находится в данном микроскопическом состоянии s.
Несмотря на простой смысл средних по времени, удобно сформулиро- вать статистические средние в данный момент времени. Вместо проведе- ния измерений на одной системе представим себе совокупность, или ан- самбль, систем, которая составлена из идентичных воображаемых копий,
характеризующихся одним и тем же макросостоянием. Количество систем в ансамбле равно числу возможных микросостояний. Ансамбль систем, ха- рактеризуемый величинами (N, V, E) и описываемый распределением ве- роятностей вида (41), называется микроканоническим ансамблем. Предпо- ложим, что некоторая физическая величина Ψ имеет значение Ψ
s
, когда система находится в состоянии s. Тогда среднее от Ψ по ансамблю дается выражением hΨi =
X
s
Ψ
s
P
s
,
(42)
где P
s определено в (41).
2.5.1
Моделирование методом Монте-Карло
Зададимся вопросом, как можно выполнить усреднение по ансамблю при заданных N , V и E. Поскольку в системе фиксирована энергия, то алгоритм Метрополиса здесь не применим. В качестве одного из способов можно было бы перебрать все микросостояния c заданной энергией и вы- числить средние по ансамблю от требуемых физических величин, как мы это делали в наших примерах. Очевидная процедура — это зафиксировать
N и V , изменить случайным образом координаты и скорости отдельных ча- стиц и принять конфигурацию, если она имеет требуемую полную энергию
E. Данная процедура, однако, весьма неэффективна, поскольку большин- ство конфигураций не будут иметь требуемую полную энергию и должны быть отвергнуты. Следуя духу метода Монте-Карло, желательно разрабо- тать практический метод получения репрезентативной выборки из полного
34

числа микросостояний.
Эффективная процедура Монте-Карло состоит в следующем. Предста- вим себе макроскопическую систему, которая составлена из двух “подси- стем”, а именно исходной рассматриваемой системы, называемой далее си- стемой, и подсистемы, состоящей из одного элемента. Эта дополнительная степень свободы называется “демон”. Демон путешествует по системе и из- меняет ее конфигурацию. Если требуемое изменение уменьшает энергию системы, то демон забирает себе избыток энергии. Если же требуемое из- менение увеличивает энергию системы, то демон наоборот отдает энергию системе. Единственное ограничение состоит в том, что демон не может от- дать больше энергии, чем у него есть в данный момент, поэтому, если у демона недостаточно энергии, изменение конфигурации системы не проис- ходит. В итоге алгоритм принимает следующий вид:
1. Формируем начальную конфигурацию. Начальная энергия демона равна 0.
2. Выбираем случайным образом частицу и производим пробное измене- ние ее координат и скорости согласно выражению (30).
3. Вычисляем изменение энергии системы ∆E, обусловленное изменени- ем координат и скорости частицы.
4. Если пробное изменение уменьшает или не изменяет энергию системы
(∆E ≤ 0), то новая конфигурация принимается, и энергия демона увеличивается на ∆E.
5. Если пробное изменение увеличивает энергию системы (∆E > 0), и энергия демона больше, чем ∆E, то новая конфигурация принимает- ся, и энергия демона уменьшается на ∆E. В противном случае новая конфигурация отбрасывается, и частица сохраняет свои старые коор- динаты и скорость.
6. Повторяем шаги 2–5 для получения достаточного числа n конфигура- ций или “испытаний”.
7. Вычисляем средние по конфигурациям.
35

8. Выполняем шаги 1–7 для получения статистически независимых серий измерений. Вычисляем средние значения по средним значениям серий по формуле (24).
Через достаточно большое количество шагов, необходимых для установ- ления равновесия, демон и система достигнут компромисса и “договорятся”
на какую-то среднюю для каждого энергию. Полная энергия остается по- стоянной, и, поскольку демон представляет только одну степень свободы по сравнению с многими степенями свободы системы, можно предположить,
что флуктуации энергии системы будут малы.
Заметим, что, с другой стороны, мы можем изучать распределение энер- гии демона E
d
, а всю окружающую систему рассматривать как резервуар.
Демон является исследуемой системой, микроканоническое состояние кото- рой определяется только ее энергией. Демон находится с окружающей его системой в тепловом равновесии и его энергию можно описывать в рам- ках канонического ансамбля. Следовательно, вероятность P (E
d
) того, что демон имеет энергию E
d имеет вид (27), где теперь E
s
= E
d
Распределение (27) обеспечивает простой способ вычисления темпера- туры T по средней энергии hE
d i демона. Поскольку hE
d i равняется hE
d i =
R E exp (−βE) dE
R exp (−βE) dE
= β
−1
= k
B
T ,
(43)
то T есть средняя энергия демона, поделенная на k
B
. Заметим, что ре- зультат hE
d i = k
B
T в (43) справедлив, только если энергия демона может принимать непрерывные значения.
Лабораторная работа 2.1. Моделирование системы методом Метрополиса в каноническом ансамбле
Целью настоящей работы является исследование канонического ансам- бля методом Метрополиса.
Содержание работы
1. Выпишите математическую модель (выражения для энергии соот- ветствующего ансамбля) для моделирования методом Монте-Карло.
36

Определите набор входных параметров.
2. Выполните обезразмеривание уравнений и начальных параметров.
3. Выполните программу. Найдите средние значения кинетической и по- тенциальной энергий путем усреднения по последним k шагам метода
Монте-Карло (самостоятельно модифицируйте программу). Значения k выбирайте самостоятельно таким образом, чтобы средние значения от серии к серии изменялись несильно.
4. Найдите распределение квадрата модуля скорости, усредняя по по- следним k шагам метода Монте-Карло.
5. Найдите наиболее вероятное значение квадрата модуля скорости ча- стицы, исходя из найденного распределения.
6. Вычислите среднюю энергию отдельно взятой частицы и среднюю энергию системы на частицу.
7. Моделирование канонического ансамбля методом Монте-Карло про- водится при фиксированной температуре. По аналогии с демоном для микроканонического ансамбля, найдите температуру газа, используя в качестве демона одну из частиц газа. Сравните найденную темпера- туру с температурой термостата.
Лабораторная работа 2.2. Моделирование системы методом Монте-Карло в микроканоническом ансамбле
Целью настоящей работы является исследование микроканонического ансамбля методом Монте-Карло.
Содержание работы
1. – 5. См. шаги 1.–5. из лабораторной работы 2.1.
6. Вычислите среднюю энергию демона и среднюю энергию системы на частицу. Зависят ли результаты, полученные в п. п. 3 и 4, от того,