Диссертация (1150258), страница 10
Текст из файла (страница 10)
Основы метода Монте-КарлоВ истории вычислительной химии метод Монте-Карло (МК) занимает особоеместо. Из всех основных расчетных методов, применяемых в этой области химии,именно метод Монте-Карло был впервые реализован на программируемой электронновычислительной машине (1953 г. [150]). Два других важнейших метода – молекулярнаядинамика и квантовохимические расчеты – появились на несколько лет позже (1959 г.[151] и 1956 г. [152], соответственно). При решении физико-химических задач спомощьюметодаМонте-Карлоосуществляетсяпостроениевыборкиизконфигурационной части фазового пространства исследуемой системы (импульсноеподпространство не рассматривается) в условиях заданного статистического ансамбля.Вычисление математических ожиданий доступных для расчета параметров системыпроизводится усреднением этих параметров по построенной последовательностиконфигураций.42Использование аппарата ансамблей Гиббса позволяет определять различныетермодинамические свойства системы A(r), если известна плотность распределения дляданного ансамбля ρ(r):A(r ) A(r ) (r ) d r ,(1)где <A(r)> – среднее по ансамблю значение термодинамической величины.
Например,дляканоническогоансамбляплотностьраспределениявконфигурационномподпространстве имеет вид: (r ) гдеU(r)–потенциалexp U (r ) k B T ,Z NVTвзаимодействия,kB–(2)константаБольцмана,ZNVT–конфигурационный интеграл:Z NVT exp U (r ) k B T dr(3)Таким образом, для расчета любой характеристики необходимо определить значениеконфигурационного интеграла.За исключением ограниченного набора систем с предельно упрощенной формойзависимости U(r) такие интегралы не могут быть рассчитаны аналитически, и для ихвычисления используются приближенные численные методы.
Одним из таких методов иявляется метод Монте-Карло. Для вычисления интегралов вида (3) метод Монте-Карло вобщем случае предполагает построение огромного количества различных конфигурацийсистемы, многие из которых будут иметь очень высокую энергию, а значит, крайненизкую вероятность. Вклад таких конфигураций в средние значения рассчитываемыхсвойств системы будет очень малым.Метод, предложенный Николасом Метрополисом и сотр. [150] позволилзначительно повысить статистическую эффективность выбора конфигураций.
Основнаяидея этого метода заключалась в использовании статистического веса для выбораконфигураций, имеющих высокую вероятность. Таким образом, в выборку не попадаликонфигурации с высокой энергией и малой вероятностью, и число конфигураций,необходимых для исследования свойств системы, резко сокращалось.Реализация метода Метрополиса основывается на применении аппарата цепейМаркова.
Цепью Маркова называется последовательность статистических испытаний,которая удовлетворяет следующим условиям:431. Исход каждого испытания зависит только от предшествующего испытания и независит от более ранних испытаний.2. Каждое испытание принадлежит к конечному множеству возможных исходов.Каждому состоянию i, в котором может находиться система, присваиваетсяопределенная вероятность ρi. Кроме того, задаются вероятности перехода πmn из одногосостояния m в другое состояние n. Полный набор вероятностей перехода πmn может бытьпредставлен в виде матрицы π размером N×N, где N – количество всех возможныхсостояний. Такую матрицу называют матрицей перехода. Каждый ряд матрицыперехода составлен из вероятностей того, что система из текущего состояния mперейдет в какое-нибудь состояние n.
По свойствам вероятностей сумма элементовкаждого ряда этой матрицы будет равна единице:mn1(4)nРаспределение вероятностей состояний {ρi}, применение к которому матрицыперехода π позволяет с ненулевой вероятностью попасть из любого состояния m влюбоесостояниераспределением.nзаДляконечноечисловероятностейпереходов,предельногоназываетсяраспределенияпредельным(предельныхвероятностей) справедливо выражение:m mn n(5)mДля физических объектов, находящихся в равновесии, определено предельноераспределениезаданноговероятностейсостояний,термодинамическогоансамбля.котороеОднакосоответствуетвероятностираспределениюпереходовπmnнеизвестны.Для поиска элементов матрицы перехода π, которые должны удовлетворитьусловиям (4) и (5), дополнительно используют принцип микроскопической обратимости m mn n nm(6)Фактически, принцип микроскопической обратимости соответствует условию (5):mm mn n nm n nm nm(7)nПотому для построения фазовых траекторий в выбранном ансамбле используюталгоритмы, основанные на построении матриц перехода π, удовлетворяющих условиям44(4) и (6).
Первый такой алгоритм был предложен Метрополисом и сотр. и носитназвание «несимметричного решения». Рассматривается два случая: mn mn mn mn n m n m , m nn m , m n(8)где αmn – произвольная симметричная матрица. Симметричность матрицы αобеспечивает соблюдение принципа микроскопической обратимости.
Несимметричноерешение не исключает вероятности, что система останется в исходном состоянии: mm 1 mn(9)nmГлавным достоинство алгоритма, предложенного Метрополисом, является то, что дляпостроения фазовой траектории необходимо знать только соотношение вероятностейдвух состояний ρn/ρm. Таким образом, пропадает необходимость рассчитыватьнепосредственно конфигурационные интегралы.Вероятности состояний и, соответственно, переходов из одного состояния вдругое, как было упомянуто выше, зависят от термодинамического ансамбля, в которомрассматривается система.2.2. Метод Монте-Карло в каноническом ансамблеКанонический ансамбль это ансамбль закрытых систем, находящихся в условияхпостоянства температуры T, объема V и числа частиц N.
Вероятности состояний зависятот их потенциальной энергии, а условная вероятность перехода из одного состояния вдругое равна отношению вероятностей этих состояний и, следовательно, зависит отизменения энергии при таком переходе.Переход из одного состояния в другое в другое происходит путем измененияконфигурации системы. В каноническом ансамбле переход из текущего состояние вследующее может быть осуществлен либо смещением молекул относительноначального расположения, либо вращением молекул относительно их центров масс(когда это возможно). Формальных ограничений на число одновременно смещаемыхили вращаемых молекул не существует, однако практика использования метода МонтеКарло показывает, что более эффективным является подход, при котором смещают иливращают не все молекулы одновременно, а лишь одну из них.
Ниже описана45последовательность действий, необходимых для построения новой конфигурациисистемы.1. Выбор случайной молекулы.2. Изменение положения выбранной молекулы путем её смещения иливращения.Для смещения молекулы координаты её центра масс изменяют на случайнуювеличину вдоль всех координатных осей. В плотных системах при смещении молекулына большое расстояние (порядка нескольких ангстрем или более) резко возрастаетвероятность получения энергетически невыгодной конфигурации. Поэтому припрактическомиспользованииметодаМонте-Карловканоническомансамблерасстояние, на которое может быть перемещена молекула, часто ограничиваютвеличиной ∆rmax, которая задается до начала расчета. Смещения вдоль каждого изнаправлений будут определяться какxin xim 21 1xmaxyin yim 2 2 1y max(10)z in z im 2 3 1z maxгде x max , y max и z max – максимально возможные смещения в направлениях x, y и z, ξ1,ξ2, ξ3 – случайные числа от 0 до 1.Ориентацию молекулы при расчетах методом Монте-Карло удобно задавать спомощью углов Эйлера.
Тогда изменение ориентации может быть достигнуто путемнебольшого изменения каждого из углов случайным образом.in im 21 1 max in im 2 2 1 max in im 2 3 1 max(11)где max , max и max – максимально возможные изменения углов Эйлера, ξ1, ξ2, ξ3 –случайные числа от 0 до 1.3. Оценка вероятности перехода.Как было упомянуто выше, вероятность перехода из состояния m в состояние nопределяется изменением энергии системы при этом переходе, ∆Umn:46U mn U n U m(12)Если новая конфигурация n имеет более низкую энергию (∆Umn ≤ 0), то вероятностьтакого состояния выше, чем исходного состояния m, и такая конфигурацияпринимается.
Если при переходе энергия системы повышается ((∆Umn > 0),тогда,согласно (8), новая конфигурация принимается с вероятностью ρn/ρm:1 n Z NVTexp( U n / k B T ) exp( U n / k B T ) exp( U mn / k B T ) 1 exp( U mn / k B T ) m Z NVT exp( U m / k B T )exp( U n / k B T )(13)Чтобы определить, принимается ли данная конфигурация, её вероятность сравниваетсясо случайным числом ξ, которое может иметь значение от 0 до 1.
Если вероятностьперехода больше случайного числа ξ, то система переходит в новое состояние n. Если,напротив, вероятность перехода меньше случайного числа ξ, то система остается впрежнем состоянии m.Описанная последовательность действий, позволяющая осуществлять переход отодного состояния к другому, называется шагом метода Монте-Карло. При расчетах вканоническом ансамбле производятся шаги двух типов: смещения и вращения. То, какоедействие будет производиться с молекулой на данном шаге, может быть определенолибо случайным образом, либо с помощью заранее заданной последовательности.Обычно шаги смещения и вращения совершают в приблизительно равном соотношении(то есть, при случайном выборе типа шага вероятности смещения и вращения должныбыть близки).Для получения достаточно полной выборки наиболее вероятных при заданныхусловиях состояний требуется проделать значительное количество таких шагов.
Какправило, длительность расчета методом МК измеряется в миллионах шагов (от десятковдо нескольких сотен, в зависимости от сложности системы).2.3. Метод Монте-Карло в большом каноническом и «гибридном» ансамбляхБольшой канонический μVT-ансамбль является ансамблем открытых систем, гдеобъем и температура системы закреплены, однако число частиц в системе можетизменяться. Реализация метода Монте-Карло для большого канонического ансамбля47была предложена Г.Э.Норманом и В.С.Филиновым в 1969 году [153,154].
Позжеобширная методическая работа по развитию метода была проведена Адамсом [155,156].С практической точки зрения большой канонический ансамбль в методе МонтеКарло можно рассматривать как расширение канонического ансамбля по числу молекул.В результате такого расширения в наборе закрепляемых переменных происходит заменаэкстенсивного термодинамического параметра (число частиц) на интенсивныйпараметр, контролирующий обмен веществом между исследуемой системой игипотетическимрезервуаром.Такимконтролирующимпараметромявляетсяхимический потенциал. Для реализации расчетов в таком ансамбле, к шагам смещенияи вращения молекулы, введенным в каноническом ансамбле, должен добавитьсяалгоритм, позволяющий менять число молекул в системе – шаги типа «вброс» и«выброс». Такие шаги позволяют регулировать число частиц в ячейке в соответствии сзаданным химическим потенциалом путем либо добавления молекулы в случайноеместо в ячейке моделирования, либо удаления из ячейки случайно выбранной частицы.Если частица удаляется из системы (шаг типа «выброс»), то отношениевероятностей двух состояний определяется по формулеnN exp U mn / k B T lnmzV (14)где z – так называемая «абсолютная активность»:z1exp / k BT 3(15)Λ – длина волны де Бройля: h2 2mk B T1/ 2(16)Вероятность принятия шага «вброса» может быть представлена какnzV exp U mn / k B T lnmN 1(17)В выражениях для вероятностей обоих типов шагов N обозначает число частиц висходном состоянии m.