Расчет радиальной функции распределения в аргоне методом Монте-Карло

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

Количественной характеристикой упорядоченности в системе может служить так называемая радиальная (парная) функция распределения g(r), которая определяется как отношение вероятности нахождения произвольной пары частиц на расстоянии r друг от друга к вероятности их однородного некоррелированного распределения [1, 2]. Функция радиального распределения может быть рассчитана по формуле [1, 2]:

1 N g (r) = (— 2

\ PN £

n i ( r , r + A r ) \ 4 n r ¹ A r

где ni(r,r+Ar) - число атомов, находящихся в сферическом слое толщиной Ar на расстоянии r от i атома, р -плотность числа частиц системы. Величина р4nr2 Ar представляет собой число частиц в сферическом слое толщиной Ar, находящимся на расстоянии r от i-го атома, для некоррелированного распределения.

Если энергию U взаимодействия между N атомами, расположенными в точках 1,2...,N, можно представить в виде

U = t Ф ( r ,).                                    (2)

• 1 < i <  j <  N

где r i, - расстояние между i и j частицами. Для функции g ( r ) справедливы еще два уравнения, связывающие ее с термодинамическими свойствами системы:

E = N

³ kT + 2 р J ^ ( r ) dr

и

Р = p kT - ¹ р ² J 151 1 1 g ( r ) dr . ⁽⁴⁾ 6 ^j dr

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

Радиальная функция распределения g ( r ) зависит от р и Т, так что, строго говоря, следовало бы записывать ее в виде g ( r; р,T) . Пользуясь уравнениями (3) и (4), можно, зная g ( r ) во всем интервале значений р и Т, вычислить все термодинамические свойства жидкости, задавшись видом межатомного потенциала ф( r ) . Вот почему парная корреляционная функция важна в статистической термодинамике жидкостей, для которых справедливо предположение (2). Это предположение обычно считается верным, по крайней мере, как хорошее первое приближение для атомов простых жидкостей. Методом Монте-Карло принято называть такие численные методы, характерной особенностью которых является использование чисто стохастических элементов в отличие от чисто детерминистических уравнений метода молекулярной динамики. Частную форму этого подхода, применяемую в физике жидкого состояния, разработали Метрополис и др. [1-2].

Метод Монте-Карло состоит из построения последовательности молекулярных конфигураций путем случайных смещений частиц модельной системы. Каждая новая конфигурация принимается или отвергается; критерием служит вероятность конфигурации пропорциональная больцмановскому фактору данной конфигурации ехр(- вA E), в =1/кТ. Заметим, что предположение о парной аддитивности потенциала межмолекулярного взаимодействия здесь не существенно. Метод Монте-Карло пригоден при любом виде межчастичного взаимодействия, важно лишь, чтобы выполнялось условие эргодичности, что равносильно ограничениям на выбор периодических граничных условий.

Обычно матрица переходов молекулярных конфигураций строится следующим образом: в системе случайным образом или поочередно выбирается какая-либо одна частица и рассматривается ее случайное смещение. Если вследствие этого изменение полной конфигурационной энергии A E < 0, то переход считается приемлемым, и прежняя конфигурация заменяется новой. Если же A E > 0, то переход может произойти лишь с вероятностью р=ехр(- вA E). В этом случае машина случайным образом выбирает десятичную дробь в интервале между 0 и 1 и сравнивает ее с величиной ехр(- вA E). Если эта экспонента больше случайного числа, то переход частицы совершается, в противном случае этот переход отвергается.

Таким образом, вероятность переходов частиц в конфигурационном пространстве оказывается пропорциональной больцмановскому фактору ехр(- вA E). Поэтому полное среднее любой функции полученное в результате реализации вышеописанного процесса при п ^^ стремится к среднему по каноническому ансамблю.

Основная задача метода Монте-Карло в данном случае сводится к тому, чтобы получить радиальную функцию распределения (1), с помощью которой можно вычислить термодинамические параметры системы [5], в частности, изотермическую сжимаемость.

B _T - ¹ = (n ₀ ² kT ) ' <_пп >  ⁽⁵⁾

где n 0 – концентрация частиц, nn - длинноволновая асимптотика корреляционной функции [5]:

(nn ) = n 0 ( 1 + n 0 J ⁽ g ⁽ r ) - ¹ ) dr )

Для системы из N частиц радиальная функция распределения для трехмерного случая имеет вид:

g ( r ) =

V 1 dN ,

N 4 n r ² dr

а для двумерного случая:

g ( r ) =

S 1 dN

N 2 n r dr

где V, S – это соответственно объем или площадь системы, dr – шаг по расстоянию, с которым производится расчет, dN – число частиц, попадающих в слой толщины dr, подсчитываемые по расположению частиц в текущей конфигурации моделируемой системы.

Радиальная функция распределения g(r) после достаточно большого числа случайных смещений частиц усредняется. Далее получив из (5) сжимаемость, вычислим скорость звука по формуле us = Y_ BT i                                    (9)

n ₀

при заданной плотности и температуре аргона [3], где у - отношение теплоемкостей C p /C V .

Нами были проведены расчеты длинноволновой асимптотики корреляционной функции nn для аргона методом Монте-Карло для температур 100, 200 и 260˚С. В таблице приведены результаты расчета скорости звука u s в аргоне по вышеуказанной методике и экспериментальные данные u s взятые из справочника [4].

Таблица

Результаты расчетов скорости звука в аргоне

T, K	р , кг/м3	P, Па	u s , м/с (расчет)	U s’ , м/с [4]	A U sM, %
100	0,508	1,047 * 104	185,3	186,1	0,4
	5,084	1,023 * 105	183,2	184,1	0,5
200	0,339	1,395 * 104	262,1	263,5	0,5
	2,542	1,045 * 105	261,9	263,2	0,5
	24,742	1,002 * 106	259,9	261,7	0,6
	287,070	1,000 * 106	241,1	-
260	2,034	1,088 * 105	298,7	300,5	0,6
	18,810	1,002 * 106	298,0	300,5	0,8
	194,882	1,000 * 107	292,6	314,0	6,8

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