Кондратьев А.С., Романов В.П. Задачи по статистической физике (Кондратьев А.С., Романов В.П. Задачи по статистической физике.djvu), страница 15

После выполнения интегрирования по и' и т1" выражение (5.49) принимает Ч Ч вид а = йт 71пйат+Ь~-)-. (5.50) Ч Суммирование по и следует проводить, введя так называемый параметр обрезания — некоторое значение до, малое по сравнению с обратным атомным размером, поскольку говорить о параметре порядка имеет смысл только на расстояниях, где еще применимо понятие о макроскопических свойствах вещества.

Переходя от суммирования к интегрированию в (5.50) и вычисляя интеграл в сферической системе координат, находим йт2пМ ~'„з,„~~ат+5~') йуУ 1 з,„(~5,1~ з,„~а~, з~ 2 з 2атчо з/з — ~д + — г — — 2Я-] а с!д(дд~~]]) Обратим внимание, что интегрирование производится только пополовине п-пространства, поскольку т1 и т1 независимы. Ч -Ч 13. Исследовать методом Монте-Карло температурную зависимость намагниченности, восприимчивости и теплоемкости в двумерной квадратной решетке Изинга в окрестности точки фазового перехода В рода. В этой модели гамнльтониан системы записывается в виде И = ~12 И,.и., (5.Ы) с /' где р.

и р. принимают только два значения: +1 и -1, а суммн- Е / рование проводится по ближайшим соседям выделенного узла. Р е ш е н и е. Описанная модель может быть достаточно просто рассчитана и проанализирована на ЭВМ. Задается квадратная плоская решетка с числом узлов Ф = Е х Е. Можно рассматривать либо конечную решетку, либо решетку с периодическими граничными условиями, чтобы устранить влияние границ. 103 Второй способ предпочтительнее, поскольку при этом результаты значительно точнее.

При периодических граничных условиях считается, что решетка окружена точно такими же копиями системы. Поэтому, если каждый спин задавать парой индексов х и у(р ), где1 х«,(., 1 у Г., то х,у И„+ „=И„„. Н„„+ =Н „И„+ у+ =Н„„. Спины, находящиеся на границе, взаимодействуют с соседями основной решетки и с окружающими копиями. Процедура, проводимая на машине, состоит в следующем. Произвольно задается начальная конфигурация спинов.

Можно либо считать все ц равными +1 или — 1, либо использовать подпрограмму случайных чисел (для машины серии ЕС это подпрограмма г<А(чВЩ, когда величина е меняется случайным образом на промежутке (О, Ц. Нетрудно для каждого узла реализовать программу типа бросания монеты: при е ЯОБ ц =1, при е»05 ц =-1. Х,у х,у Далее полезно напечатать начальное состояние, чтобы иметь возможность следить за его эволюцией. Параметрами, которые вычисляются для такой системы, будут: намагниченность М = = Яц.

(полезно также иметь модуль этой величины ~М~ г — ~Яр.(), энергия системы, рассчитываемая по формуле (5.29) ! (ее удобно вычислять в единицах яТ). Для этого вводится параметр К = 1/(лТ) или К вЂ” «безразмерная температура». При вычислениях следует иметь в виду, что при Т 0 система ферромагнитная, а при Т «О антиферромагнитная. Восприимчивость ~ определяется из формулы для флуктуации намагниченности 1( = К(<М > — <М> ).

И наконец, зная энергию как функцию температуры, можно вычислить теплоемкость С = йЕ/П. Обратим внимание, что в точке фазового перехода при вычислении на ЭВМ у и С должны иметь максимумы, а <М> — обращаться в нуль, если двигаться со стороны низких температур. Вычисление средних значений проводится методом Монте-Карло. Он позволяет построить марковскую цепь так, чтобы конфи- -и!!ит) гурации появлялись с весом е ' .

Как альтернативный ваант можно было бы осуществить перебор всех возможных конфи- 104 гураций (их будет 2 ). Но при больших У такая процедура нерациональна. Аппарат марковской цепи позволяет организовать последовательную процедуру перебора конфигураций с выбранными специальным образом переходными вероятностями от одной конфигурации к другой так, чтобы конфигурации после некоторого па-нl(ит) чальиого участка уже имели вероятности е Итак, у нас есть основная плоская ячейка размером т. х Е. Прежде всего, с вероятностью 1/Е для каждой из двух координат (х,у) выбираем узел. Координаты каждого узла выбираем по формуле р = ~1+1, где ~Е округляется до меньшего из двух целых чисел А « ~Е « «А + 1.

Пусть номер этого узла равен У (х,у). Совершаем пробный переворот спина этого узла, т.е. заменяем р на о -р . Подсчитываем изменение энергии при таком перевороте: 0 би(ц, -Н ) = и(-ц„)-йр„). о "о ~о "о При вычислении изменения энергии, естественно, принимаются во внимание только ближайшие соседи. Если И/ — О, то переворот осуществляется, и мы приходим к новой конфигурации.

Если Иl » О. то переворот осуществляется с вероятностью -Ьи/(ит е, которая лежит в интервале О е ) 1. Для этого берется случайное число О «г. «1 и сравнивается с -биl(ьт) — йиl(ит) е . Если е «е, то переворот осуществляется; если ~ » е ), то перехода нет, и система остается в прежнем состоянии. Но это состояние учитывается фактически еще раз в счетчике числа конфигураций, и оно принимается во внимание при вычислении средних значений 6 <Х> = ЯХ,/6, (=1 где 6 — полное число конфигураций.

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

В качестве опоркой точки можно использовать критическую температуру Т, которая изс' вестна из точного решения Онзагера 21п~о = 1. с Даже для систем порядка 6 х 6 илн 10 х 10 «машинное» решение уже близко к точному, особенно если применять описанные выше периодические граничные условия. Составьте программу в соответствии с описанной процедурой расчета н проведите вычисления на ЭВМ. 14. Методом молекулярной динамики исследовать свойства системы из Ф частиц, взаимодействующих между собой посредством потенциала Леннарда — Джонса, помещенных в куб с ребром Е, проверить установление в системе максвелловского распределения по скоростям, вычислить корреляционную функцию скорости <ч(0)ч(1)>, рассчитать радиальную функцию распределения, вычислить энергию и теплоемкость С и исследовать поведение теплоемкостн вблизи критической точки.

Р е ш е н и е. Для системы из Ф частиц записываются уравнения движения глг,= Г,, ~ =1,2,...,У. (5. 52) Г, = — Я ЧФ(~г.— г.~) ! есть сила, действующая на ~'-частицу со стороны остальных частиц. Для устранения влияния стенок на систему накладываются периодические граничные условия, согласно которым считается, что куб находится в центре системы, состоящей из точно таких же кубов с тем же распределением частиц. Прн этом, если частица уходит из исследуемого куба, необходимо учесть, что из соседнего куба в него войдет частица. При вычислении сил Г.

необходимо учитывать только несколько блнг жайших соседей, т.е учитывать взаимодействие в области с фиксированным радиусом Р. Потенциал межмолекулярного взаимодействия в форме Леннарда-Джонса имеет вид Ф~г) = 4Я ((Г/Г) — (Г/Г) ] . 106 В начальный момент положения частиц задаются произвольно. Скорости выбираются таким образом, чтобы средняя кинетиче- ская энергия соответствовала определенной температуре и Ф Я ч.

= О, т.е. отсутствовали направленные потоки. Параметры г=! ! интересующей системы можно найти по справочникам. Например, для жидкого аргона при температуре Т = 95 К и плотности р = = 1,374 г/см они равны — — 1,26, ~~У = 0,934, — = 4,65 10!! с Численное интегрирование системы (5.52) проводится следующим образом. После 1-го шага по времени известны совокупности координат (г. 1, скоростей (ч. ) и ускорений (а; 1.

(!) (!) (!) ! ! Тогда скорости и координаты, спустя промежуток времени „(!+!) „(!) (!) !О ! ! ,<~), ~ „ („(~.~)„и) Этот результат обычно уточняется следующим образом. По значениям (г. ) находятся новые ускорения !'а. 1. Затем ! !+!) (!+!) )о ! скорости и координаты снова пересчитываются в предположении, что ускорение изменилось на этом шаге линейно: (Я+ ) (!) 1 !а(!) ~ (!+!Яй! !.(!+!) ( !+!) + 1 ) а(!+!) (!Ц ((!!)2 Число частиц выбирается равным 10 — 20, интервал времени -!3 й! 1О с. Максвелловское распределение по скоростям, которое должно установиться через некоторое число шагов, проверяется по известным значениям моментов максвелловской функции распределения <п> <и> ... <и "> и = 4 5.

Если значения, полученные на 1-м шаге, совпадают с точно вычисляемыми с помощью функции распределения Максвелла с погрешностью порядка 10 Я, то распределение по скоростям можно считать установившимся. Необходимо также проверить, выполняется ли закон сохранения энергии. Если он начинает заметно нарушаться, то надо уменьшить временной шаг й!. Автокорреляционная функция скоростей <ч(0)ч(!)> вычисляется путем усреднения по всем частицам с начальным моментом 107 с последующим усреднением по разным начальным моментам ОГ 02' '"' Ор' ж .(0).(1» = ~К К.(1„).(1„1).

р,"=1,, О! 0~ Радиальная функция распределения п(г), которая имеет физический смысл плотности вероятности нахождения частицы на расстоянии г от данной, вычисляется следующим образом. Вокруг каждой частицы проводятся шаровые слои толщиной примерно 0,1-0,2 диаметра частицы у и. Измеряется среднее число частиц в каждом слое. Функция й(г) определяется как отношение этого числа частиц к числу частиц при средней плотности Ф/(~, Таким образом, прн больших г величина д(г) стремится к единице.

Энергия системы вычисляется по формуле Е = ~ЯФ(~г.-г.~) + 2 ~;гпо.. е,1 с Для вычисления теплоемкостн необходимо выполнить расчеты энергии при разных температурах. Для исследования окрестности критической точки необходимо задать в системе критическую плотность и выполнить расчеты прн температурах, близких к критической, но несколько превосходящих ее. Составьте программу в соответствии с описанной процедурой и проведите расчеты на ЭВМ.