16.7.   МОДЕЛИРОВАНИЕ КЛАССИЧЕСКИХ ЖИДКОСТЕЙ

Из повседневного опыта мы хорошо знаем о существовании трех фаз вещества: твердой, жидкой и газообразной (рис. 16.5). Исходя из этого опыта мы в состоянии различить эти три фазы по их свойствам. Например, твердые тела обладают жесткостью в противоположность газам и жидкостям, которые текут под действием касательного напряжения. Нашей задачей в этом разделе является получение с помощью метода Монте-Карло более детального представления о качественных различиях всех трех фаз.

Моделирование классических жидкостей методом Монте-Карло существенно упрощается за счет того обстоятельства, что координаты скорости (импульса) ии на что не влияют. Например, вклад координат скорости в среднюю энергию составляет ^ на степень свободы. Следовательно, нам нужно производить выборку только положений молекул, т.е. «конфигурационных» степеней свободы. Возможно ли такое упрощение для квантовых систем?

К числу характеристик жидкости, представляющих интерес с точки зрения физики, относятся средняя энергия, удельная теплоемкость и уравнение состояния. Еще одной важной величиной является парная корреляционная функция g(r), которую мы определим ниже. Предположим, что в области объема V содержится N частиц, так что средняя плотность равна р = N/V. (В двумерном и одномерном случаях объем V заменяется соответственно на площадь и длину.) Выберите одну из частиц, которую назовем опорной частицей, в качестве начала системы координат. Тогда вероятность обнаружить вторую частицу в интервале между г и г + Лг равна pg(r)dt, где элемент «объема» dx = 4nr2dr (d = 3), 2mdr (d = 2) и 2dr (d = 1). Мы предполагаем, что g(r) —> —* 0 при г —» 0, поскольку частицы не могут проникать друг в друга. Кроме того, мы предполагаем, что g(r) —> 1 при г —» оо, поскольку неоднородность в окрестности частицы ограничена конечным размером области. Заметим, что если проинтегрировать pg(r)dr по dr, то получится N-1, т.е. полное число частиц в системе за вычетом опорной частицы. В задачах 16.16—16.18 мы найдем, что g(r) является мерой флуктуации плотности и тем самым локального «порядка» системы.

Исходя из функции g(r) можно также получить несколько термодинамических характеристик. Например, если в системе присутствуют только двухчастичные взаимодействия, то можно показать, что средняя потенциальная энергия на частицу может быть представлена в виде

Можно также воспользоваться вириалом силы (см. гл. 6) и записать уравнение состояния в виде

Твердые сферы в d-мерном случае. С целью разделить эффекты короткодействующего отталкивания и длиннодействующего притяжения исследуем сначала модель твердых сфер с потенциалом межчастичиого взаимодействия

Такая модель широко изучалась в одномерной (твердые стержни), двумерной (твердые диски) и трехмерной постановках как с помощью метода Монте-Карло, так и метода молекулярной динамики.

Поскольку взаимодействие (16.30) вообще не содержит притягивающей части, то нет и никакого фазового перехода газ — жидкость. Имеется ли фазовый переход между жидкими фазами прн низких плотностях и твердыми фазами при высоких плотностях? Может ли твердая фаза образоваться при отсутствии притягивающего взаимодействия? В дальнейшем мы получим ответы на эти вопросы.

Какие физические величины имеет смысл рассматривать для системы твердых сфер? Никаких тепловых характеристик, таких, как средняя потенциальная энергия, нет, поскольку для твердых сфер эта величина всегда равна нулю. Наибольший интерес представляет g(r), поскольку она дает информацию о корреляциях частиц и уравнении состояния. Если потенциал имеет вид (16.30), то можно показать, что выражение (16.29) приводится к виду

Мы будем вычислять g(r) для различных значений г и затем экстраполировать получаемые результаты на г = с.

Алгоритм Метрополиса для твердых сфер легко реализовать, выбирая случайным образом частицу и перемещая ее в новое пробное положение. Если новое положение перекрывает другую частицу, данное перемещение отвергается и сохраняется старая конфигурация; в противном случае перемещение принимается. Максимальное смещение б разумно, хотя и не обязательно оптимально, выбирать из тех соображений, чтобы при данном выборе S принималась приблизительно половина пробных состояний. Самую большую трудность в реализации этого алгоритма представляет выявление пересечения двух частиц. Если число частиц не слишком велико, достаточно вычислять расстояния между пробной частицей и всеми остальными частицами. Эта простая процедура используется в программе harddisk для системы твердых дисков. Для систем ббльшего числа частиц данная процедура занимает слишком много времени. Как говорилось в гл. 6, лучше применять другой метод: разбивать систему на «ячейки» и проверять частицы только в своей и соседних ячейках.

Выбрать начальные координаты дисков сложнее, чем может показаться на первый взгляд. Один метод заключается в размещении дисков в ящике случайным образом. Если какой-нибудь диск накладывается на уже имеющийся, то его бросают куда-то в другое место. Таким способом удается довольно быстро вычислить возможную начальную конфигурацию в случае низкой плотности. Что касается достижения более высоких плотностей, то можно поступать следующим образом. Представьте себе, что стенки сдвигаются по направлению к центру ящика до тех пор,  пока одна из стенок едва не коснется одного из дисков. Затем диски перемещаются какое-то число шагов Монте-Карло и стенки вновь сдвигаются. Очевидно, что даииая процедура тем тяжелее, чем больше плотность.

Другой метод отправляется от упорядоченного простраиствеииого распределения частиц. Можно, например, начать с дисков, образующих треугольную решетку максимальной плотности, которая иас интересует. Далее можио расширить решетку, пересчитав пропорционально координаты дисков. В результате в новой конфигурации ие возникнет ии одного пересечения дисков. Даииая процедура осуществляется в подпрограмме initial. Заметим, что для выполнения условия, чтобы каждый диск имел одинаковое число ближайших соседей, необходимо рассматривать только четное число частиц. Систему единиц удобно выбрать такой, чтобы расстояния измерялись в радиусах твердого кбра о*. Приведенная плотность равна р* = р/о^, где d—размерность системы.

В подпрограмме move выбирается случайная частица, генерируются ее пробные координаты с поправкой иа периодические краевые условия и определяется ее пересечение с другими частицами.

Существуют по крайней мере два способа обработки периодических краевых условий. До сих пор для проверки положения частицы мы прибегали к инструкциям IF. В подпрограммах cell и separation примени-

ется другой метод. Он основан иа использовании свойств функции truncate, имеющейся в языке True BASIC. Данная функция, когда ее второй аргумент равен нулю, обрезает все знаки дробной части первого аргумента. Например, truncate^. 3,0) = 9, truncate(-9.3,0) = -9 и truncate(9.7,0) = 9. В качестве побочного упражнения читателю предоставляется возможность определить, какой метод иа вашем компьютере работает быстрее.

Заметим, что функция truncate(x.O) языка True BASIC эквивалеитиа функции int(x) в Фортране. В противоположность этому функция int языка True BASIC выдает наибольшее целое, не превосходящее значения ее аргумента. Следовательно, в языке True BASIC имеем ш/(9.3) = 9, mf(-9.3) = -10 и ш/(9.7) = 9.

Теперь, когда определено «правильное» положение пробного перемещения, мы обращаемся к подпрограмме overlap с целью определить, не перекрывается ли данным пробным перемещением положение какого-нибудь другого диска.

В подпрограмме correl вычисляется массив gcum, содержащий полное число частиц, отстоящих на расстояние от г до г + Lr от заданной опорной частицы. Обратите внимание на то, что массив gcum находится для g N опорных частиц. Несмотря на то что подпрограмма correl вызывается после каждого шага Монте-Карло, возможно, стоит вычислять корреляции пореже.

В подпрограмме output получается нормированная парная корреляционная функция g(r), для чего массив gcum делится на число выборок, плотность и площадь 2nrdr кольца, отстоящего от опорной частицы на

расстояние г. Поскольку gcum найдено для = N опорных частиц, мы де-

1

лим этот массив также иа g N, избегая тем самым двойного счета.

Последнюю конфигурацию просчитанного варианта следует сохранить в файле, с тем чтобы можно было ее использовать в качестве начальной конфигурации в следующем варианте.

В качестве проверки программы hard_disk рассмотрим сначала одномерную систему твердых стержней. В силу простоты этой системы уравнение состояния и g(r) могут быть вычислены точно. Например, уравнение состояния имеет вид

при этом плотность равна р = N/L. Поскольку твердые стержни не могут проходить друг сквозь друга, то легко видеть, что исключаемый объем составляет Mr и, значит, доступный объем равен L - Ncr. Заметим, что форма уравнения состояния (16.32) совпадает с уравнением Ван-дерВаальса (см. гл. 6) при ненулевом вкладе притягивающей части взаимодействия.