Моделирование микроканонического ансамбля методом Монте Карло в пакете MATLAB (стр. 5 )

Результаты выполнения описанной последовательности команд представлены на рис. 6-8.

Рис. 6. Зависимость мгновенных значений полной энергии системы от времени

Рис. 7. Зависимость мгновенных значений энергии демона от времени

Рис. 8. Функция распределения энергии демона

Задача 6.
Так как изначально неизвестны конфигурации спинов, обеспечивающие искомое значение энергии, в качестве начальной конфигурации мы выбираем конфигурацию с минимальной энергией, поэтому для достижения равновесного состояния с заданной полной энергией системе требуется некоторое конечное время (время релаксации). Это хорошо видно из рис. 11, 12. Аппроксимируйте зависимости Esystem(t), Edemon(t), используя метод наименьших квадратов. Оцените время релаксации полной энергии системы system и время релаксации энергии демона demon. Как зависят system, demon от числа спинов Nspin, образующих систему? Сравните полученные результаты, с результатами полученными при моделировании случайных блужданий.
Задача 7.
Создайте собственные функции, позволяющие исключать переходный период и вычислять равновесные значения средней полной энергии системы <Esystem>, среднюю энергию демона <Edemon>, среднюю намагниченность системы <M>, средний квадрат намагниченности <M2>. Сравните значения равновесной температуры системы с соответствующим значением температуры, полученной без учета процесса релаксации системы к равновесному состоянию.
Задача 8.
После достижения равновесного состояния значения полной энергии системы и энергии демона, как видно из рис. 11, 12, испытывают случайные флуктуации около средних значений. Исследуйте как зависит величина флуктуаций средней полной энергии системы <Esystem>, средней энергии демона <Edemon>, средней намагниченности системы <M>, среднего квадрата намагниченности <M2> от Nspin. Сравните полученные результаты с результатами, полученными при моделировании случайных блужданий.
Задача 9.
Вычислите <Esystem> для систем, состоящих из различного числа спинов N. Сравните полученные результаты с точным ответом для бесконечной одномерной решетки

(17)

Как полученные результаты зависят от числа спинов N и количества шагов метода Монте-Карло?
Задача 10.
Проведите вычисления и постройте зависимость <M2> от температуры T для системы, состоящей из N спинов (N = 200). Аппроксимируйте данную зависимость методом наименьших квадратов. Зависит ли вид аппроксимирующей от числа спинов.
Задача 11.
Введите ненулевое магнитное поле h и вычислите <Edemon>, <M>, <M2>, как функции от h для заданной полной энергии. Определите связь <Edemon> с температурой для h0. Сравните значения равновесных температур для случаев h=0, h0 при одинаковых значениях средней полной энергии.

Обобщим, использованный выше метод микроканонического ансамбля, использованный нами для одномерной системы спинов, на случай двумерной системы (рис. 9).
Как и в одномерном случае выбираем периодические (тороидальные) краевые условия, т. е. будем представлять решетку кольцом, в котором спины sij, находящиеся в узлах (i=1, j=1,...,N) взаимодействуют со спинами, находящимися в узлах (i=N, j=1,...,N), а спины, расположенные в узлах (i=1,...,N, j=1), взаимодействуют со спинами, расположенными в узлах (i=1,...,N, j=N) (рис. 10). Это обеспечивает равное число взаимодействий для всех спинов в изучаемой системе.

Рис. 9. Квадратная двумерная решетка, состоящая из N = 16 спинов

Рис. 10. Задание граничных условий в двумерном случае

Взаимодействие выбранного спина с соседями можно рассматривать, как взаимодействие с одним спином, величина которого равна алгебраической сумме величин четырех соседних спинов, которая, как очевидно, может принимать три значения: 0, 2, 4. Следовательно, минимально возможное значение изменение энергии при опрокидывании одного спина составляет 2J-(-2J)=4J, в отличие от одномерного случая в котором квант изменения энергии равнялся 2J. Отмеченное обстоятельство определяет необходимость замены в (16) 2J 4J:

(18)

Алгоритм моделирования двумерной системы, состоящей из N спинов, методом Монте-Карло реализуется следующей последовательностью действий:
1. Задание числа спинов решетки Nspin.
2. Задание числа шагов метода Монте-Карло на спин.
3. Задание ориентации спинов в узлах квадратной решетки в момент времени t=1 (начальной конфигурации системы).
4. Выбор случайным образом одного из спинов системы.
5. Вычисление пробного изменения энергии.
6. Если пробное изменение приводит к уменьшению энергии системы, то система отдает энергию демону и новая конфигурация принимается.
7. Если пробное изменение увеличивает энергию системы, то новая конфигурация принимается в том случае, если демон имеет достаточную энергию для передачи ее системе.
8. Если пробное изменение не меняет энергию системы, то принимается новая конфигурация.
9. Повторение пп. 4-8 (число повторений равно числу спинов в системе).
10. Повторение пп. 4-9 (число повторений равно числу шагов метода Монте-Карло на спин).
Для реализации описанного алгоритма моделирования двумерной системы спинов методом Монте-Карло в пакете MATLAB создадим файл Ising2.m, содержащий описание функции, возвращающей мгновенные значения: энергии системы, энергии демона, намагниченности, числа принятия решений, а также мгновенные конфигурации спинов.

% листинг файла Ising2.m
function [Es,Ed,SpM,A,S] = Ising2(Nspin,J,h,Esi,NTrial)
% функция, возвращающая мгновенные значения: энергии системы,
% энергии демона, намагниченности, числа принятия решений, а также
% мгновенные конфигурации спинов
% Nspin число спинов решетки
% J константа обменного взаимодействия
% h напряженность внешнего магнитного поля
% Esi конечная энергия системы
% NTrial число испытаний
Ns=Nspin.^0.5; % число спинов вдоль одной стороны
s=ones(Ns,Ns); % начальная конфигурация спинов
Esystem=-(J+h)*Nspin; % начальная энергия системы
Edemon=4*J*floor((Esi-Esystem)/(4*J)); % начальная энергия демона
Es(1)=Esystem; Ed(1)=Edemon; S=s; k=1; for i=1:NTrial  Accept=0;  for j=1:Nspin
 % случайный выбор узла сетки
 Ix=floor(Ns*rand(1)+1);  Iy=floor(Ns*rand(1)+1);
 % граничные условия
 if Ix==1  Left=Ns;  else  Left=Ix-1;  end;  if Ix==Ns  Right=1;  else;  Right=Ix+1;  end;  if Iy==1  Down=Ns;  else;  Down=Iy-1;  end;  if Iy==Ns  Up=1;  else  Up=Iy+1;  end;
 % пробное изменение энергии
 de=2*s(Iy,Ix)*(-h+J*(s(Iy,Left)+s(Iy,Right)+s(Down,Ix)+s(Up,Ix)));  if de<=Edemon % принятие пробного изменения энергии
 s(Iy,Ix)=-s(Iy,Ix);  Accept=Accept+1;  Edemon=Edemon-de;  Esystem=Esystem+de;  end;  k=k+1;  Es(k)=Esystem;  Ed(k)=Edemon;  A(k-1)=Accept;  s1=sum(s);  SpM(k)=sum(s1);  S=cat(3,S,s);  end; end; A=A/NTrial;