ПРИМЕНЕНИЕ МЕТОДА МОНТЕ-КАРЛО В ЗАДАЧАХ МОДЕЛИРОВАНИЯ ФАЗОВЫХ ПЕРЕХОДОВ

УДК 519.688:004.94

И.Л. Детченков, М.И. Рвачева, А.Г. Масловская

ПРИМЕНЕНИЕ МЕТОДА МОНТЕ-КАРЛО В ЗАДАЧАХ МОДЕЛИРОВАНИЯ

ФАЗОВЫХ ПЕРЕХОДОВ

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

APPLICATION OF MONTE-CARLO METHOD IN PROBLEMS OF PHASE TRANSITION SIMULATION

The algorithmization and program implementation of two-dimensional Ising model were performed in the article. The stochastic model is based on application of the Metropolis algorithm realized by means of the Monte-Carlo method. Computational scheme was presented to calculate the physical quantities characterizing the ferromagnetic phase transition ofthe second order. Test problems are used to demonstrate the results of the computing experiment.

Введение

В настоящее время моделирование фазовых переходов и связанных с ними критических явлений - одна из важнейших прикладных задач статистической теории и привлекает внимание ученых в различных областях научного знания [1-6]. Модель для изучения магнитных фазовых переходов была впервые предложена в 1920 г. В. Ленцем и вошла в историю под именем его ученика Изинга. Можно заметить, что фундаментальные законы, описывающие динамику поведения объектов естественных и технических наук, нашли масштабное применения для формализации экономических и финансовых процессов. Так, в работе [5] представлена концепция финансовых рынков, обнаруживающих нелинейное статистическое поведение, а также показана возможность описания характеристик подобных систем на основе физической модели Изинга. В основе такого подхода лежит принцип аналогии. Представление финансового рынка в виде комплексной системы инвесторов-агентов, взаимодействующих между собой, оказалось достаточным для описания стохастической динамики такой системы с помощью модели Изинга.

Аналитическое изучение фазовых переходов в предложенной модели вызывает большие затруднения, при этом метод Монте-Карло [4, 7, 8] предоставляет возможности численного моделирования систем большой размерности, постановки и проведения вычислительного эксперимента, направленного на изучение закономерностей, характеризующих магнитные материалы и фазовые переходы в них. Метод Монте-Карло нашел приложение и в статистическом моделировании поведения характеристик макроскопических систем, поскольку предусматривает относительно несложную процедуру определения средних величин в каноническом ансамбле.

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

Постановка задачи моделирования

Рассмотрим ^-мерную кристаллическую решетку, содержащую N=1? узлов, где Ь - характерный размер решетки. Решетка расположена в плоскости XOY. В узлах находятся электроны, обладающие магнитным моментом (спином). Спин $г может принимать значение si — +1, если он сонаправлен с осью 02 и значение si — -1, если направлен противоположно, как показано на рис. 1. Любая конфигурация задается набором переменных s1,S2,...SN для всех узлов решетки. Энергия системы при наличии магнитного поля будет определяться по формуле [1, 2, 9]:

Рис. 1. Схема расположения спинов в узлах кристаллической решетки.

Е — -Е ^,^ - Мб Е ^ , эВ

(и) /=1

где Ji ] - потенциал взаимодействия спинов, эВ (или Дж,

ра, Дж/Тл.

Рис. 2 а и рис. 2 б демонстрируют ферромагнитное состояние системы. В этом случае потенциал взаимодействия для двух соседних спинов J > 0 и данное состояние оказывается энергетически выгоднее (Е — - J ). Если же потенциал взаимодействия J < 0 , то энергетически более выгодным становится реализация антиферромагнитного состояния системы (рис. 2 в) с Е — + J .

а б в

Рис. 2. Схематическое изображение конфигурации спинов в феромагнитном (а, б) и антиферромагнитном (в) состояниях системы.

Для описания состояний физической системы вводят равновесные статистические характеристики - среднюю энергию (Е) и среднюю намагниченность ^м) с использованием средних по ансамблю значений рассматриваемых физических величин [2, 3]:

Е = - X Е ехр

М = - X Мр ехр

^ • Т

где 7 = X ехр

( Е Л __Р_

^ ■ Т

- сумма по всем Р микросостояниям системы; Т - температура системы, К;

кт = 1.38-10-23 - постоянная Больцмана, Дж/К.

Макроскопические характеристики системы - теплоемкость С и магнитная восприимчивость % могут быть определены следующим образом [3, 6]:

кТ - Т и

АУ - кТ - Т

Е^-(Е)2

^т2^ - (ш)2 ^,

гДе Е = — X Е„ ехр 2s = 1 S

кТ • Т

; 2 = X ехр

кТ - Т

;(т) = — X ш ехр

кТ - Т

В простейшем случае вычисления можно проводить не по схеме Гиббса (2)-(3), а с использованием значений генеральных математических ожиданий:

X Е X М

Е=^, М=у=1

В математической модели суммарная намагниченность системы определяется как сумма: N

М = X у. и в размерных величинах выражается через единицы ¡и^ .

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

), индукция внешнего магнитного поля (В), потенциал взаимодействия спинов (У ), характерный размер решетки (Е), постоянная Больцмана

(кТ ). Переменные модели и величины, подлежащие определению: конфигурация спинов (у7), энергия

системы (Е), магнитный момент системы (М), теплоемкость системы (С), магнитная восприимчивость системы (х).

Алгоритм решения прикладной задачи и его программная реализация

Для решения сформулированной задачи используем метод статистических испытаний и алгоритм Метрополиса [2, 3]. Формализуем алгоритм с помощью блок-схемы, представленной на рис. 3.

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

Рис. 3. Блок-схема алгоритма реализации модели.

Рис. 4. Структура программного приложения, предназначенного для моделирования

характеристик магнитных систем.

Исследование характеристик магнитного фазового перехода II рода по данным вычислительного эксперимента

Целый ряд чистых материалов (например, железо, никель) и различных сплавов в отсутствие внешнего магнитного поля имеют спонтанную намагниченность [1]. Данные свойства проявляются только при температуре ниже температуры Кюри. При температурах выше температуры Кюри спонтанная намагниченность исчезает, т.е. происходит переход из упорядоченного режима, в котором спины расположены упорядоченно, в хаотический режим, в котором ориентация спинов случайна.

Моделирование проводилось с использованием нормированных параметров. Значения параметров кТ, ¡л0, ¡лБ, А V полагались равными единице; энергия измерялась в единицах кТ • Г . В качестве входных параметров требуется инициализировать число узлов решетки N константу обменного взаимодействия J, величину внешнего магнитного поля В, температуру системы Г (или массив значений температур) и количество испытаний Р для реализации метода Монте-Карло.

Вычислительный эксперимент № 1. Целью данного эксперимента является демонстрация результатов модельного расчета конфигурации спинов, динамики изменения энергии и намагниченности при фиксированной температуре и нулевом внешнем поле. Установим начальные параметры: число узлов решетки N=256; константа обменного взаимодействия J=1 отн. ед.; внешнее магнитное поле В=0 отн. ед.; температура системы Г=1.5 отн. ед.; количество испытаний при реализации метода Монте-Карло Р=1500. Роль такта времени играет номер итерации.

Модельные представления конфигурации системы в начальный момент времени и спустя продолжительное время продемонстрированы на рис. 5 (начальная конфигурация - все спины ориентированы вверх) и рис. 6 (соседние спины направлены противоположно). Незакрашенными ячейками показаны ориентированные вверх спины, закрашенными - вниз. Время выполнения расчетов - 3 сек.

а б В

Рис. 5. Конфигурации спинов системы в момент времени t = 0 - а, t = 2000 - б, t = 10000 отн. ед. - в (N = 256, J = 1 отн. ед., В = 0 отн. ед., Г = 1.5 отн. ед.).

а Б в

Рис. 6. Конфигурации спинов системы в момент времени t = 0 - а, t = 2000 - б, t = 10000 отн. ед. - в (N = 256, J = 1 отн. ед., В = 0 отн. ед., Г = 1.5 отн. ед.).

Представленные результаты свидетельствуют, что при J<0 независимо от первоначальной конфигурации спинов система релаксирует к ферромагнитному состоянию - все спины ориентированы вверх, как показано на рис. 5 в, или все спины ориентированы «вниз» (рис. 6 в). Таким образом, при J=1 энергетически более выгодным состоянием является то, в котором большинство спинов будет сонаправлены, что соответствует ферромагнитным состоянию. Значению J= -1 должно соответствовать антиферромагнитное состояние системы. На рис. 7 приведены модельные конфигурации системы в различные моменты времени при значении параметра моделирования 3= -1. В этом случае независимо от начальной конфигурации система релаксирует к антиферромагнитному состоянию.

а б в

Рис. 7. Конфигурации спинов системы в момент времени t = 0 - а, t = 2000 - б, t = 10000 отн. ед. - в (N = 256, 3 = -1 отн. ед., В = 0 отн. ед., Т = 1.5 отн. ед.).

На рис. 8 приведены результаты расчета мгновенных значений энергии и намагниченности от времени для параметров моделирования, соответствующих представленным на рис. 7 результатам.

Рис. 8. Зависимость мгновенных значений энергии (а) и намагниченности (б) системы от времени при

начальной ориентации спинов «вверх».

Результат расчета намагниченности в данном случае свидетельствует о релаксации к антиферромагнитному состоянию.

Вычислительный эксперимент № 2. Цель второго эксперимента - моделирование зависимости макроскопических характеристик от температуры. Расчету подлежали температурные зависимости магнитной восприимчивости (5), теплоемкости (4), энергии (2) и намагниченности системы (3). Для получения этих зависимостей была проведена серия вычислений при следующих параметрах модели: число узлов решетки N=256; константа обменного взаимодействия 3={1; -1}; внешнее магнитное поле В=0 отн. ед.; температура системы изменялась от 0.1 отн. ед. до 4.1 отн. ед., с шагом 0.2 отн. ед.; количество испытаний в методе Монте-Карло Р=1500. Первоначальная конфигурация спинов соответствовала ферромагнитному состоянию - все спины ориентированы «вверх». Расчет при данных параметрах составил примерно 60 сек. На рис. 9 а изображен график зависимости теплоемкости от температуры, на рис. 9 б - график зависимости магнитной восприимчивости от температуры при 3=1 отн. ед.

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

Т, отн ед. Т, отн ед.

Рис. 9. График зависимости теплоемкости от температуры - а, график температурной зависимости

магнитной восприимчивости - б при J = 1 отн. ед.

Вычисленные зависимости средней энергии и средней намагниченности системы от температуры при J=1 отн. ед. позволяют отметить, что при низких температурах Т<ТС решетка будет спонтанно намагничиваться (свойство класса ферромагнетиков), в этом интервале намагниченность принимает наибольшие значения, в то время как энергия соответствует минимальным значениям. По достижению температуры Кюри происходит резкое падение намагниченности (и возрастание энергии), что соответствует переходу в парамагнитное состояние. Результаты моделирования зависимости средней энергии и средней намагниченности системы от температуры при J = —1 отн. ед. показывают, что вещество в антиферромагнитном состоянии слабо намагничивается, при любой температуре его средняя намагниченность приблизительно равна нулю, поэтому далее будем рассматривать только ферромагнитное состояние.

Вычислительный эксперимент № 3. Целью данного эксперимента является определение температуры фазового перехода. Начальные параметры моделирования: число узлов решетки N=256; константа обменного взаимодействия J=1 отн. ед.; внешнее магнитное поле В=0 отн. ед.; температура системы изменялась от 0.1 до 4.1 отн. ед., с шагом 0.2 отн. ед.; количество испытаний ^=1500. Расчет при данных параметрах занял 60 сек.

Представляет определенный интерес установление значения температуры Кюри по графику зависимости намагниченности от температуры. Используя известные сведения из математического анализа, рассмотрим температурные зависимости первой и второй производных от (М). Из полученных графиков видно наличие критического поведения производных намагниченности в окрестности точки фазового перехода. Температура Кюри определяется точкой минимума графика функции первой производной и нулем графика функции второй производной. Численная оценка этого значения дает величину Тс = 2.39 отн. ед., что с точностью до численной погрешности в расчете производной соответствует значению 2.27 отн. ед., полученному для аналитического решения задачи моделирования фазовых переходов [2].

Вычислительный эксперимент № 4. Целью эксперимента являлось изучение влияния размера решетки на значения магнитной восприимчивости и теплоемкости. Для получения этой зависимости была проведена серия вычислений при начальных параметрах: число узлов решетки N = {100, 144, 196, 256}; константа обменного взаимодействия J=1 отн. ед.; внешнее магнитное поле В = 0 отн. ед.; температура системы изменялась от 0.1 до 4.1 отн. ед., с шагом 0.2 отн. ед.; количество испытаний ^=1500. Моделирование проводилось для ферромагнитного состояния системы. На рис. 10 а и 10 б

представлены зависимости теплоемкости и магнитнои восприимчивости от температуры, при разных параметрах решетки.

Рис. 10. Зависимость теплоемкости от температуры - а, зависимость магнитной восприимчивости от температуры - б при различных значениях параметра решетки.

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

Вычислительный эксперимент № 5. Цель эксперимента заключалась в изучении влияния внешнего магнитного поля на поведение модели. Для исследования этой зависимости была проведена серия вычислений при начальных параметрах: число узлов решетки N=256; константа обменного взаимодействия J=1 отн. ед.; внешнее магнитное поле 5={0, 0.01, 0.02} отн. ед.; температура системы изменялась от 0.1 до 4.1 отн. ед., с шагом 0.2 отн. ед.; количество испытаний ^=1500. Моделирование проводилось для ферромагнитного состояния системы. На рис. 11 а и 11 б представлены зависимости теплоемкости и магнитной восприимчивости от температуры при разных значениях внешнего магнитного поля В.

Рис. 11. Зависимость теплоемкости от температуры - а, зависимость магнитной восприимчивости от температуры - б при различных значениях внешнего магнитного поля.

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

Вычислительный эксперимент № 6. Целью эксперимента являлось получение зависимости намагниченности от внешнего магнитного поля при разных температурах. Для исследования этой зависимости была проведена серия вычислений при начальных параметрах: число узлов решетки N=256; константа обменного взаимодействия J=1 отн. ед.; внешнее магнитное поле B изменялось от -2 до 2 отн. ед., с шагом 0.1 отн. ед.; температура системы Т={1.5; 3; 10} отн. ед.; количество испытаний -Р=1500. Моделирование проводилось для ферромагнитного состояния системы. На рис. 12 представлены зависимости, полученные в результате расчетов.

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

В, отн. ед. В, отн. ед.

.00 С-I-1-1-1-1-I-1-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

Д отн. ед. в

Рис. 12. Зависимость намагниченности от внешнего магнитного поля при Т = 1.5 отн. ед. - а,

при Т = 3 отн. ед. - б, при Т = 10 отн. ед. - в.

Заключение

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

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