Метод Монте Карло. Некоторые сведения из теории вероятности и матстатистики, статистической физики и термодинамики. Цепи Маркова. Реализа

8. Средние от функций  случайных величин. Если  F (x1,x2,…,xn) - это функция непрерывных случайных величин, а w (x1, x2, …, xn) - плотность вероятности для событий x1, x2, …, xn, то статистическое среднее величины  F определяется как:

или можно записать:                                                                    , где W(F)- плотность вероятности заданного значения  F.

Теперь вспомним некоторые представления из термодинамики и статфизики.

Описание физической системы в терминах функций макроскопических состояний:

Макроскопическое состояние задается параметрами, измеряемыми в эксперименте: объем V, давление P, температура T.

 Одно макроскопическое  состояние ® множество микроскопических состояний.

Система N частиц в конфигурационном пространстве W, занимающая объем V при температуре T (N,V,T - ансамбль) с потенциальной энергией UW (XN), где XN характеризует совокупность всех координат xi. Весь объем V ® большое число пронумерованных ячеек s: система в состоянии Аi, если изображающая ее точка находится в i-й ячейке; среднее значение любой наблюдаемой физической величины <F> записывается в виде (распределение Гиббса или каноническое распределение):

 

Здесь N - число частиц, V - объем, Т - температура, b = (kT)-1; k =1,372´10-16 эрг/град -                        константа Больцмана, UN(Ai) - энергия взаимодействия частиц, x - их координаты. Задачи т/д и статфизики сводятся к вычислению фазового, или конфигурационного интеграла Q-1(N,V,T). С помощью такого интеграла можно найти средине значения т/д величин, например, энергию:                                                                               

 

 

Вероятность нахождения системы в состоянии Аi определяется как

 

w(i) = Q-1(N,V,T)exp[-bUN(Ai)]

 

Как это делается и вообще, что это такое? Пусть у нас есть система из N материальных точек. Тогда состояние системы можно изобразить с помощью 6N канонических переменных: 3N координат и 3N скоростей. Или с помощью математической абстракции - воображаемого 6N-мерного пространства, координатами которого являются 6N канонических переменных. Такое пространство называется фазовым или конфигурационным пространством. Точка в фазовом пространстве изображает состояние системы в данный момент времени. С течением времени точка перемещается в фазовом пространстве - так появляется фазовая траектория. С классической точки зрения состояние системы полностью определяется определением ее канонических переменных. Такое состояние системы называется микроскопическим. Макроскопическое состояние задается параметрами, измеряемыми в реальном эксперименте - макроскопическом опыте. Объем, давление, температура и т.д. Одному макроскопическому состоянию может соответствовать множество микроскопических состояний. Макроскопические параметры являются функциями канонических переменных Fk(x), где x - канонические переменные. Заданием всех Fk(x) не определяются все хk. Поэтому  из макроскопических измерений можно сделать лишь статистические выводы о микроскопических переменных х, то есть через задание плотности вероятности этих переменных. Это фазовая плотность вероятности или фазовое распределение. Зная фазовое распределение можно вычислить вероятность обнаружить систему в заданном фазовом объеме.

Фазовый ансамбль - это совокупность фазовых точек, т.е. микроскопических состояний системы. Каждому состоянию приписывается определенная плотность вероятности, тогда отдельная точка фазового ансамбля является случайной величиной. И такой фазовый ансамбль называется статистическим фазовым ансамблем. Нас интересуют системы, находящиеся в термодинамическом равновесии, т.е. когда фазовая плотность вероятности не зависит явно от времени. Система может быть разбита на подсистемы, и каждой будет соответствовать свой интеграл энергии, и больше никаких интегралов (могут быть интегралы импульсов, момента количества движения и т.д.). Системы, находящиеся во внешнем потенциальном поле, где интегралы импульса и момента количества движения не сохраняются, и все поведение системы определяется только интегралом энергии, называются эргодическими системами. С математической точки зрения эргодические системы это такие системы, фазовое пространство которых нельзя разделить на две области так, чтобы фазовая точка из одной области никогда бы не попала во вторую область. Фазовая точка как бы оббегает все подпространство заданной энергии, т.е. может проходить сколь угодно близко к любой из точек подпространства. Для эргодической системы среднее из точек по времени от любой величины F(x) - функции только энергии и не зависит от других интегралов движения. Но это не доказано строго, т.е. то, что реальные физические макроскопические системы являются эргодическими, теоретически не доказано. Это принимается как гипотеза. То есть, не доказано, что системы, моделирующие реальные макроскопические физические системы, являются эргодическими. Тогда можно считать, что равновесная фазовая плотность вероятности является функцией только энергии системы. Доказать это невозможно.

Итак, основное положение термодинамики гласит: макроскопическое состояние системы, т.е. значение ее внутренних параметров при термодинамическом равновесии, зависит от внешних параметров и энергии (или температуры, если температура определяется как функция энергии).

В термодинамике выделяют два типа систем: адиабатическая система - это изолированная от внешних тел и имеющая определенную, строго заданную энергию Е; и изотермическая система - система, находящаяся во взаимодействии с термостатом, имеющим заданную температуру Т: kT = q, q - т/д температура, Т- в градусах.

Свободная энергия тела - это работа, производимая над телом при обратимом процессе с постоянными температурой и объемом. Основное распределение термодинамики для канонического NVT-ансамбля - это распределение Гиббса или каноническое распределение.

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

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

В методе МК независимыми переменными, сохраняющими постоянные значения при моделировании, выбираются N, V и Т (N,V,Т-ансамбль). Молекулы двигаются случайным образом в соответствии с генератором случайных чисел, и каждое новое расположение либо принимается, либо отбрасывается с вероятностью, определенной по закону р » е-U/kT. После достаточно большого числа шагов и последующего усреднения можно получить оценки равновесных свойств системы.

Такие вычисления возможно только численными методами. Поэтому от многократных интегралов переходим к суммированию - интегральному суммированию. При расчетах таких интегральных сумм выявляются две основные сложности: малость элементарного объема dq и необходимость учета большого числа частиц N. Как преодолевают эти трудности? Первую - малость объема - с помощью построения марковских цепей. В общих чертах смысл такого подхода в возможности учета не всех слагаемых в интегральной сумме, а в выборе лишь основных, наиболее типичных слагаемых, определяющих значения интегральной суммы, т.е. проведение так называемой существенной выборки. Трудность, связанная с учетом большого числа частиц в исследуемой системе решается так. Свойства макросистем вполне проявляются уже на свойствах микрообъема. Фактически, как показывают прямые расчеты, достаточно рассматривать столько частиц, чтобы занимаемый ими объем лишь в несколько раз превышал объем, в котором проявляются корреляции между частицами. Т.о., за исключением некоторых особых точек можно ограничиться несколькими сотнями частиц. Т.к. нас интересует состояние системы вблизи положения равновесия, то для наших расчетов такое предположение вполне оправдано. Для устранения влияния флуктуаций искомые величины находятся путем усреднения по большому числу существенных состояний системы с помощью цепей Маркова. Погрешность, обусловленная поверхностными эффектами при ограничении размеров системы, может быть снижена путем расчетов для нескольких значений числа частиц в системе и экстраполяцией на бесконечно большое число частиц. Кроме того, можно выбрать определенные граничные условия для изучаемой системы. Конечность числа рассмотренных состояний системы приводит к появлению статистической ошибки, а рассмотрение относительно малого числа частиц дает систематическую погрешность.

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

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

Что такое упомянутые выше марковские цепи?

Под состоянием Аi понимается набор значений микроскопических (координаты, скорости и др.) и макроскопических (число частиц, объем, температура и др.) переменных, характеризующих систему. В зависимости от того, каким было исходное состояние, отличаются друг от друга и наборы переменных, входящих в Аi. Пусть F(Ai) - известная функция состояния Ai, вид этой функции различен для различных термодинамических величин - средних величин. Плотность вероятности - из распределения Гиббса, явный вид этой функции известен и является одинаковым для всех вычисляемых величин - средних значений различных термодинамических функций одной и той же системы. Поскольку выражения для вычисления среднего и распределения Гиббса общие, то можно разработать формальную схему для вычисления средних значений, не конкретизируя, что именно понимается под Ai. Совокупность всех возможных Ai рассматривается как набор значений случайной дискретной величины, распределенной с вероятностью wi. Производя достаточно большое число испытаний этой случайной величины можно образовать последовательность A1, A2, … AM и найти F(A1), …, F(AM)/

Из закона больших чисел получим:

 

если М достаточно велико. Т.е. наш алгоритм разбивается на две части: моделирующую, т.е. генерирующую последовательность Ai, и измерительную - вычисляющую средние значения термодинамических величин. Получение последовательности A1, A2, … AM, отвечающей заданному распределению вероятности, достигается с помощью аппарата цепей Маркова. Что такое - марковские цепи? Это случайные процессы, или последовательность случайных состояний, у которых вероятность каждого последующего состояния связана с предыдущим состоянием системы. То есть мы будем говорить, что совокупность всех возможных состояний Ai образует цепи Маркова с постоянными вероятностями переходов Ai ® Аj , равными

и удовлетворяющим условию нормировки:

 

Каковы свойства цепей Маркова? Переход в последовательность из состояния Ai в состояние Аj за n шагов может осуществляться по различным промежуточным состояниям. Обозначим рij - как суммарную вероятность реализации перехода Ai ® Аj за n шагов по всевозможным путям. Если все Ai образуют один эргодический класс (рij ³ 0 – условие эргодичности), т.е. из любого Ai достижимо любое Аj при конечном числе переходов. То есть мы получаем переход системы к стационарному распределению состояний независимо от выбора начального состояния.

В пространстве Ai существует большое число марковских цепей, реализующих искомое распределение. Конкретный выбор рij производится по соображениям простоты, а также по тому, какая из реализуемых с их помощью последовательностей A1, A2, … AM позволяет определить средние значения при меньшем числе шагов М и какая последовательность быстрее выходит на стационарный участок.

Конфигурации, соответствующие распределению w(i), получают следующим образом:

• генерируется последовательность конфигураций Аi , ..., АM, которая соответствует цепи Маркова с одношаговой  вероятностью перехода pij системы из состояния Аi в Аj;

При таком выборе pij, начиная с достаточно большого n, состояния Аi будут появляться с вероятностью, пропорциональной w(in), т.е. генерируется ансамбль с распределением w(i). Средние значения величины <F> будут вычисляться как средние арифметические на стационарном участке марковской цепи:

В нашем случае известны w(i), нужно определить pij .

Фактически испытание проводится следующим образом. Датчик случайных чисел выдает случайное число x, равномерно распределенной на интервале между 0 и 1. Если оказывается, что это число больше вероятности перехода рij то переход Аi®Аj считается осуществившимся и  в качестве Аi2 в последовательности A1, A2, … AM выбирается Aj1, если же случайное число меньше вероятности перехода, то в качестве Аi2 выбирается Аi1. Начиная с выбранного Аi2 вся процедура повторяется сначала. После проведения большого числа М таких последовательных испытаний определяется требующаяся последовательность и выполняя вычисления средних для различных т/д величин, находим искомые величины. Построенная таким образом последовательность A1, A2, … AM называется траекторией цепи Маркова. Пояснить построение такой последовательности можно на простом примере. В последовательности состояний Аi, реализуемых при решении конкретных задач статфизики, некоторые состояния могут повторяться за счет неосуществимых переходов. Но подавляющее большинство состояний Аi из всей возможной совокупности состояний вообще не встречается в этой последовательности. Но при этом полученная последовательность оказывается достаточно представительной, чтобы вычисления средних для т/д величин давали результаты с точностью порядка 1% и выше. Именно малость требуемой величины числа испытаний М по сравнению с общим числом возможных состояний Аi и объясняет эффективность существенной выборки метода МК, по сравнению с другими численными методами при вычислении многократных интегралов и сумм, встречающихся в выражении для т/д величин. Но М все же достаточно велико - порядка 1-10 млн. испытаний.