Алгоритм Метрополиса - Метод Монте Карло в химическом моделировании

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

Wi = (1.10)

Если в некоторый фиксированный объем помещать случайным образом молекулы, энергия взаимодействия между которыми задается набором потенциальных функций, то в зависимости от конфигурации системы больцмановский множитель exp (- Ui /кТ) может принимать различные значения[11-17]. Некоторые конфигурации дают значительный вклад в канонические средние, а некоторые - практически нулевой (например, когда две частицы сближены настолько, что между ними имеется сильное отталкивание).

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

Среднее по ансамблю от любой физической величины M рассчитывают по формуле

<M>= iexp(-)/= (1/11)

Где i - номер конфигурации (среднее берется по всем рассмотренным конфигурациям системы)[16]. Поскольку усреднение (1.11.) проводят по конечному числу конфигураций m со смещенной выборкой, для исключения влияния смещения каждую конфигурацию необходимо брать с весом 1/рi :

<M>? exp(-)/ (3)

Метрополис с соавторами [15-16] предложил в качестве рi взять распределение Больцмана

Рi = exp()/ (1.13)

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

<M>? (1/m) (1.14)

Ансамбль, состоящий из m конфигураций, получают путем задания вероятностей перехода от одной конфигурации к другой. Вероятность перехода от i-й конфигурации к j-й pij считают зависящей от энергий этих конфигураций, а точнее от величины (Uj - Ui)/kT:

Рij = Pji exp[-(Uj - Ui)/kT] (1.15)

Таким образом, строят простые цепи Маркова, т. е. последовательности случайных событий, в которых вероятность определенного события зависит от исхода предыдущего испытания. В соответствии с условием микроскопической обратимости вероятности pij должны удовлетворять условиям:

Рij=exp() Pji=exp() и (1.16)

Обычно полагают, что при

Uj ? Ui Pji = wi (1.17)

И при Uj >Ui Рij = wi exp[-(Uj - Ui)/kT]

Где wi - вероятность появления некоторой конфигурации при случайном выборе (с использованием последовательности равномерно распределенных случайных чисел). Используя центральную предельную теорему теории вероятностей может быть доказано [17], что рассматриваемая цепь Маркова задает распределение, асимптотически стремящееся к каноническому.

На практике алгоритм Метрополиса реализуют следующим образом. Пусть заданы потенциал взаимодействия, конфигурация системы (начальное расположение частиц в элементарной ячейке моделирования) и температура Т. Рассчитывают потенциальную энергию системы Ui и вносят случайное изменение в конфигурацию (случайным образом выбирают к - ю частицу в ячейке и смещают ее). При этом энергия системы становится равной Uj. Если Uj < Ui, то считают, что система перешла в новое состояние. Если Uj < Ui, то сравнивают величину exp() со случайным числом оє(0,1). Если о? exp[-(Uj - Ui)/kT], то считают, что система перешла в j-е состояние. Если же о> exp[-(Uj - Ui)/kТ ], то переход в новое состояние не происходит, к - я частица сохраняет свои прежние координаты [18]. При этом j-ую конфигурацию в цепи не учитывают, а рассматривают прежнее расположение частиц, соответствующее энергии U Таким образом, чем больше значение энергии имеет система при случайном изменении конфигурации, тем с меньшей вероятностью она переходит в это состояние.

Очевидно, что расчет средних значений сопряжен с большим количеством вычислений, выполнение которых стало возможным только при появлении быстродействующих ЭВМ[19].

Одно из преимуществ метода Монте-Карло состоит в том, что алгоритм легко адаптировать к любому статистическому ансамблю. Например, при моделировании системы в NPT-ансамбле необходимо периодически изменять объем ячейки, а в уравнении (1.17) вместо разности потенциальных энергий использовать разность энтальпий:

ДH = ДU + PДV - kТ ln(1+ДV/V)N, (1.18)

Где P - давление, V - объем ячейки, ДV - изменение объема ячейки.

Усредняя термодинамические функции по конфигурациям на равновесном участке цепи Маркова по уравнению (1.14), легко рассчитать конфигурационную энтальпию и мольный объем :

Нконф = <U> + P<V>, Vm = <V>NA/N (1.19)

Для расчета других термодинамических функций - изобарической теплоемкости СР, изотермической сжимаемости kT и коэффициента температурного расширения б, необходимо проводить моделирования системы при различающихся параметрах состояния, а затем находить конечные разности[20]. Оценку указанных термодинамических функций можно сделать и по результатам одного моделирования, используя флуктуационные формулы:

KT==),

==),

Для решения большинства задач вполне достаточно проведения вычислений в каноническом ансамбле. Однако если требуется, то можно использовать и большой канонический ансамбль [21, 22].

Похожие статьи




Алгоритм Метрополиса - Метод Монте Карло в химическом моделировании

Предыдущая | Следующая