Метод Монте-Карло в статистической физике

Размещено на http://www.allbest.ru/

Введение

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

Метод Монте-Карло был разработан в 1949 году в работах Дж. Неймана (J. Neumann), С. Улама (S. Ulam), H. Метрополиса (N. Metropolis). Метод Монте-Карло основан на методе статистического моделирования, который был известен ещё в XVIII веке. Классическим примером такого моделирования является «игла Бюффона»: получение числа ? путём случайного бросания иглы на горизонтальную поверхность, расчерченную сеткой равноотстоящих параллельных линий [13]. С появлением ЭВМ метод был усовершенствован и получил в 1949 году название «метод Монте-Карло». В настоящее время основные усилия исследователей направлены на создание эффективных Монте-Карло алгоритмов различных физических, химических и социальных процессов для параллельных вычислительных систем. Кроме того, различные варианты метода Монте-Карло можно использовать для решения задач оптимизации.

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

Целью настоящей работы является изучение метода Монте-Карло (алгоритм Метрополиса-Гастингса) на примере жидкости находящейся в объеме V, парные взаимодействия между частицами которой описываются потенциалом Леннарда-Джонса.

1. Метод Монте-Карло в статистической физике

1.1 Работа алгоритмов Монте-Карло

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

Существенным для понимания метода Монте-Карло является то, что вместо решения аналитической задачи для приближенного решения можно моделировать случайный процесс и использовать такие статистические оценки, как вероятности, математические ожидания.

Оказывается, что очень часто вместо вычисления ряда сложных аналитических выражений можно экспериментально определить соответствующие статистические оценки.

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

Статистический метод имеет ряд преимуществ, таких как точность и не сложность математических расчетов, возможность получения динамических характеристик, отсутствие дополнительной избыточности, вводимой для целей синхронизации [14].

случайный число жидкость идеальный



1.2 Основные свойства стандартного случайного числа

Основным инструментом для моделирования случайных величин на ЭВМ является подходящий генератор стандартных случайных чисел, дающий выборочные значения аi случайной величины а, равномерно распределенной в интервале .

Перечислим вероятностные характеристики случайной величины а. Распределение этой величины является абсолютно непрерывным с плотностью.

, . (1)

Функция распределения равна

(2)

Несложно вычислить математическое ожидание и дисперсию

, (3)

(4)

Следующее свойство равномерно распределенных точек является важным для обоснования алгоритмов метода Монте-Карло. Если одномерная точка равномерно распределена в области конечного объема , то она также равномерно распределена в произвольной подобласти объема при условии попадания в эту подобласть; при этом [12].

1.3 Генераторы случайных чисел и псевдослучайные числа

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

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

Другим способом генерации хорошей выборки псевослучайных чисел является создание набора из большого количества случайных чисел. Хотя следует помнить, что свойства источника псевдослучайных чисел будут распространяться и на такие наборы.

Рассмотрим последовательность чисел: . Возникает вопрос, является ли она случайной? По определению случайной величиной называется величина, которая может принимать одно из множества значений в результате опыта, причём появление того или иного значения этой величины до её измерения нельзя предсказать с абсолютной точностью. Если между символами нет зависимости, то последовательность данных чисел будет случайной.

В качестве количественной меры неопределенности случайной величины обычно вводят энтропию [10]. Генераторам случайных чисел (ГСЧ) для формирования случайных чисел необходимо начальное значение, которое обычно получают с помощью источников энтропии. Генератор псевдослучайных чисел (ГПСЧ) использует единственное начальное значение, откуда и следует его псевдослучайность, а ГСЧ всегда формирует случайное число, имея высококачественную случайную величину, предоставленную различными источниками энтропии.

Следует отметить, что в настоящее время не существует такого детерминированный алгоритма, который может генерировать полностью случайные числа. Как правило, ГПСЧ может только аппроксимировать некоторые свойства случайных чисел. Любой ГПСЧ с ограниченными ресурсами рано или поздно зацикливается - начинает повторять одну и ту же последовательность чисел. Длина циклов ГПСЧ зависит от самого генератора и составляет около , где - размер внутреннего состояния в битах. Следует отметить, что линейные конгруэнтные Конгруэнтное число - это натуральное число, равное площади прямоугольного треугольника со сторонами, длины которых выражаются рациональными числами. генераторы и LFSR LFSR (linear feedback shift register) - регистр сдвига битовых слов, у которого значение входного (вдвигаемого) бита равно линейной булевой функции от значений остальных битов регистра до сдвига.-генераторы обладают максимальными циклами порядка . Если порождаемая последовательность ГПСЧ сходится к слишком коротким циклам, то такой ГПСЧ становится предсказуемым и непригодным для практических приложений.

1.4 Метод Метрополиса-Гастингса

Алгоритм Метрополиса-Гастингса используется для сложных функций распределения. В нем вспомогательная функция распределения меняется со временем.

Алгоритм Метрополиса-Гастингса позволяет производить выборку из любой функции распределения. Он основан на создании цепи Маркова, то есть на каждом шаге алгоритма новое выбранное значение зависит только от предыдущего . Данный алгоритм использует вспомогательную функцию распределения , зависящую от . На каждом шаге для этой функции генерируется случайное значение . Затем с вероятностью

(5)

(или с вероятностью , если ), выбранное значение принимается как новое: , а иначе оставляется старое: .

Возьмем, например, в качестве вспомогательной функции нормальную функцию распределения:

. (6)

Для выбранного случайного значения , следующее выбираемое значение случайной величины будет для функции .

Теперь найдем произведение , где

(7)

является отношением вероятностей между промежуточным значением и предыдущим, а

(8)

_ это отношение между вероятностями пойти из в или обратно. Если симметрична, то второй множитель равен . Случайное значение на новом шаге выбирается по следующему правилу:

Если , то ,

Иначе если , то

Алгоритм стартует из случайного значения , и некоторое количество шагов работает «вхолостую», чтобы «забыть» о начальном значении.

Лучше всего алгоритм работает в случае, когда форма вспомогательной функции близка к форме целевой функции , но добиться этого зачастую невозможно. Для решения этой проблемы вспомогательную функцию настраивают в ходе подготовительной стадии работы алгоритма. Например, для нормального распределения настраивают его параметр так, чтобы доля «принятых» случайных значений (то есть тех, для которых ) была близка к . Если слишком мала, то значения будут получаться слишком близкими и доля принятых будет высока. Если слишком велика, то с большой вероятностью новые значения будут выскакивать в зоны малой вероятности , отчего доля принятых значений окажется низкой.

Для численной реализации этот метод можно представить в следующем виде [3]:

1) выбираем начальное состояние системы;

2) вычисляем энергию исходной системы ;

3) выбираем (случайным образом) частицу с номером ;

4) помещаем выбранную частицу в точку , где - исходные координаты i частицы, величины выбираются из интервала ;

5) вычисляем энергию системы для новой конфигурации системы;

6) если новая энергия системы оказывается меньше исходной , то мы сохраняем новое состояние системы и ее энергию переходим на шаг 3);

7) в случае если , то мы генерируем случайное число , где ;

8) если , то мы оставляем новую конфигурацию системы и переходим на шаг 3).

Рисунок 1 иллюстрирует схему работы алгоритма Метрополиса-Гастингса.

Рисунок 1 - Допущение перемещений частиц с увеличением энергии системы в алгоритме Метрополиса-Гастингса

1.5 Другие методы исследования свойств жидкостей

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

Еще один из способов моделирования жидкостей - решение уравнений гидродинамики, которое по своей сути является уравнением движения жидкости, внутри которой выделяется некоторый объем. В случае идеальной жидкости сила из второго закона Ньютона сводится к давлению окружающей объём жидкости и воздействию внешних силовых полей. Если поле представляет собой силы инерции или гравитации, так что эта сила пропорциональна напряжённости поля и массе элемента объёма. Тогда

,

где S _ поверхность выделенного объёма, _ напряжённость поля. Переходя, согласно формуле Гаусса _ Остроградского, от поверхностного интеграла к объёмному и учитывая, что , где _ плотность жидкости в данной точке, получим:

В силу произвольности объёма подынтегральные функции должны быть равны в любой точке:

Выражая полную производную через конвективную производную и частную производную

получаем уравнение Эйлера для движения идеальной жидкости в поле тяжести:

где _ плотность жидкости, _ давление в жидкости, _ вектор скорости жидкости, _ вектор напряжённости силового поля, _ оператор набла для трёхмерного пространства [11].

1.6 Граничные условия

Решить проблемы, связанные с наличием у системы границ, обычно удается с помощью периодических граничных условий. Рассмотрим, например, трехмерную решетку. Самый верхний и нижний слои трехмерной решетки, так же как и другие грани куба, в этом случае рассматриваются как соседи. На рисунке 2 это показано на примере квадратной решетки размером .



Рисунок 2 - Два варианта задания периодических граничных условий для квадратной решетки размером

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

В ряде случаев при изучении объемных характеристик модели используются граничные условия, отличающиеся от строго периодических. Например, для изучения свободной поверхности системы применяют граничные условия со свободной границей. Тогда «отсутствующие спины» за свободной поверхностью просто полагаются равными нулю. Этот случай показан на рисунке 3. При условии свободной границы рекомендуется делать выборку на границе (и в прилегающих слоях) больше, чем внутри объема [8].



Рисунок 3 - Квадратная решетка размером (узлы решетки обозначены точками, а связи между нами прямыми линиями) с граничными условиями со свободной границей по всем направлениям решетки (слева) или только в одном направлении (справа)

1.7 Потенциал парного взаимодействия частиц. Потенциал Леннарда-Джонса

Потенциал Леннарда-Джонса позволяет моделировать парное взаимодействие неполярных молекул [7]. Он описывает зависимость энергии взаимодействия двух частиц от расстояния между ними (см. рисунок 4). Этот потенциал позволяет передать свойства реального взаимодействия сферических неполярных молекул, поэтому широко используется в расчётах и при компьютерном моделировании.

Потенциал Леннарда-Джонса записывается в следующем виде:

(9)

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

Рисунок 4 - Характерный вид потенциала Леннарда-Джонса (синяя линия), и модельного ступенчатого потенциала (красная линия)

Второй член в формуле (9) отвечает за взаимное притяжение на расстояниях порядка , в то время как на расстояниях меньших работвет первый член , который отвечает за взаимное отталкивание частиц.

На практике также могут использоваться модельные ступенчатые потенциалы, которые позволяют упрощать расчеты многократно (он показан на рисунке 4 красными линиями):

где коэффициент определяет ширину прямоугольной потенциальной ямы.



1.8 Радиальная функция распределения частиц

Ближний порядок - упорядоченность во взаимном расположении атомов или молекул в веществе, которая повторяется лишь на расстояниях, соизмеримых с расстояниями между атомами, то есть ближний порядок - это наличие закономерности в расположении соседних атомов или молекул (см. рисунок 5).

Понятие ближнего порядка вводится через парную функцию распределения . Для этого представляется в виде:

где - одночастичная функция распределения, а - расстояние между двумя молекулами.

Рисунок 5 - Ближний порядок в жидкостях

Функция носит название радиальной функции распределения. В основе такого представления парной функции распределения лежит предположение об однородности жидкости и изотропности потенциала взаимодействия.

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

Однако для реального вещества ситуация иная. На рисунке 6 показана характерная радиальная функция распределения для жидкости Леннард-Джонса вблизи тройной точки. Она имеет осцилляции, затухающие с ростом r. Таким образом, вероятность найти молекулы на расстояниях, соответствующих локальным максимумам больше, нежели на расстояниях, соответствующих локальным минимумам в жидкости присутствует ближний порядок.

При увеличении температуры или уменьшении плотности ближний порядок становится всё менее отчётливым. Для разреженного реального газа

,

где - потенциал парного взаимодействия частиц, k - постоянная Больцмана, T - температура системы.

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

Рисунок 6 - Радиальная функция распределения для жидкости Леннарда - Джонса вблизи тройной точки ( )

Следует также отметить, что рассчитываемая теоретически радиальная функция распределения напрямую связана с наблюдаемой на практике величиной , описывающей структурные свойства вещества. Величина , также называемая структурным фактором связана с через трехмерное преобразование Фурье [1]:

(10)

здесь - волновой вектор ( _ волновое число), а интегрирование ведется в области пространства от до .

Получаемый экспериментально спектральный фактор также может быть преобразован в радиальную функцию распределения с помощью обратного преобразования Фурье [1]:

(11)



2 Практическая часть

2.1 Генерация начального состояния

В качестве начального распределения выбираем случайное распределение частиц в объеме (рисунок 7). Для заданного случайного распределения частиц рассчитывается начальная энергия системы - .

Основные параметры жидкости, используемые для генерации начального распределения следующие:

1) жидкость состоит из частиц жидкого аргона (Ar) c плотностью при температуре ;

2) объем системы рассчитывается по общим формулам , где - число плотности жидкого аргона, - контактное расстояние между частицами в безразмерных единицах.

Рисунок 7 - Начальное распределение частиц в объеме

Для генерации начального состояния был написан программный код на языке Си (приложение A).



Заключение

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

Для исследования структуры жидкого аргона был написан программный код, позволяющий генерировать случайное начальное распределение N частиц с заданной плотностью ?. В следующей работе будет реализован алгоритм Метрополиса-Гастингса, позволяющий приводить данную систему к состоянию термодинамического равновесия и рассчитывать радиальную функцию распределения. Также будет проведено численное исследование рассмотренной здесь Леннард-Джонсовской жидкости и реализовать расчет структурных факторов для такой системы при разных начальных параметрах.



Библиографический список

1. Allen M. P. and Tildesley D. J Computer simulation of liquids [Книга]. - Oxford : Clarendon Press, 1987. - стр. 402.

2. Frenkel D. and Smit B. Understanding molecular simulation: from algorithms to applications [Книга]. - Cornwall : Academic Press, 1996. - стр. 23-51.

3. Landau D. P. and Binder K. A Guide to Monte Carlo simulation in statistical physics [Книга]. - New York : Cambridge University Press, 2009. - стр. 490.

4. Metropolis N. and Ulam S. The Monte Carlo Method [Журнал] // Journal of the American Statistical Association. - 1949 r.. - 247 : Т. 44. - стр. 335--341.

5. Monte Carlo and Molecular Dynamics of Condensed Matter Physics [Конференция] / ред. Binder K. Ciccotti G.. - Como : [б.н.], 1996.

6. S. Fishman G. Monte Carlo : concepts, algorithms, and applications [Книга]. - Berlin : Springer, 1996. - стр. 110.

7. Srinivasan S. G. and Baskes M. I. On the Lennard-Jones EAM potential [Журнал] // Proceedings of the Royal Society A . - 2004 r.. - Т. 460. - стр. 1649-1672.

8. Биндер К. и Хеегман Д. В. Моделирование методом Монте-Карло в статистической физике: Введение: Пер. с англ. В. Н. Задкова. [Книга]. - Москва : Наука. Физматлит., 1995. - стр. 144.

9. Гулд Х. и Тобочник Я. Компьютерное моделирование в физике [Книга]. - Москва : Мир, 1990. - Т. 1 : стр. 28.

10. Духин А. А. Теория информации [Книга]. - Москва : Гелиос АРВ, 2007. стр. 248.

11. Ландау Л. Д., Лифшиц, Е. М. Гидродинамика [Книга]. - Москва : [б.н.], 1986. - Т. VI.

12. Михайлов Г. А. и Войтишек А. В. Численное статистическое моделирование. Методы Монте-Карло : учеб. пособие для студ. вузов Книга Москва : Издательский центр «Академия», 2006.

13. Федотов Н. Г. Методы стохастической геометрии в распознавании образов [Книга]. - Москва : Радио и связь, 1990. - стр. 142.

14. Цыпкин Я. З. и Попков Ю. С. Теория нелинейных импульсных систем Книга - Москва : Наука, 1973. - стр. 414.



Приложение A

Заполнение объема частицами жидкого аргона

/*Программа для заполнения объема V частицами жидкого аргона (N = 100) при плотности rho = 0.84. Заполнение происходит случайным образом, с условием неперекрывания частиц между собой. */

//====================Версия 05.2016=================

#include <stdlib.h>

#include <cstdlib>

#include <iostream>

#include <fstream>

#include <ctime>

#include <iomanip>

#include <cmath>

using namespace std;

const int N=1000; // число частиц в кубическом объеме

const double sigma = 1.0; // контактное расстояние между частицами

const double n=0.84/(sigma*sigma*sigma); // число плотности

const double V=N/n; // объем

const double L=pow(V,1.0/3.0); // длина стороны

double sqr(double x) // функция возведения в квадрат

{

return x*x;

}

double func(double dx,double dy,double dz)

{

return sqrt(sqr(dx)+sqr(dy)+sqr(dz));

}

int main(void)

{

setlocale(LC_ALL, "rus");

bool flag;

double *x = new double[N];

double *y = new double[N];

double *z = new double[N];

for( int i=0;i<=N;i++)

{

if (i>=1)

{

flag = true; //

while(flag) //

{

x[i]=L*rand()/RAND_MAX;

y[i]=L*rand()/RAND_MAX;

z[i]=L*rand()/RAND_MAX;

int counter=0;

for (int j=0; j<i; j++)

{

double dr=func(x[i]-x[j],y[i]-y[j],z[i]-z[j]);

if (dr<=sigma)

counter++;

if (counter>0)

flag=true;

else

flag=false;

}

}

}

else

{

x[i]=L*rand()/RAND_MAX;

y[i]=L*rand()/RAND_MAX;

z[i]=L*rand()/RAND_MAX;

}

}

ofstream fout;

fout.open("coords.txt",ios_base::out);

for (int i=0;i<N;i++)

fout<<i<<"\t"<<x[i]<<"\t"<<y[i]<<"\t"<<z[i]<<endl;

fout.close();

free(x); // освобождение памяти

free(y);

free(z);

return 0;}

Размещено на Allbest.ur