Критические явления в термодинамике спиновых систем тема диссертации и автореферата по ВАК РФ 01.04.02, кандидат наук Солдатов Константин Сергеевич

• Солдатов Константин Сергеевич
• кандидат наук
• 2019

Оглавление диссертации кандидат наук Солдатов Константин Сергеевич

Введение

Глава 1. Критические явления и фазовые переходы в векторных

моделях

1.1 Термодинамические состояния и средние величины

1.2 Двумерная модель Изинга с взаимодействием между ближайшими соседями

1.3 Расчет плотности вероятности состояний в модели Изинга на квадратной решетке

1.4 Расчет термодинамических и магнитных характеристик из статистической суммы

1.5 Краткие выводы по главе

Глава 2. Конечно-размерный скейлинг в ферромагнитных спиновых

системах на решетке пирохлора

2.1 Критические индексы

2.2 Масштабная инвариантность и конечно-размерный скейлинг

2.3 Структура решетки пирохлора и вычислительный алгоритм

2.4 Критические индексы ферромагнитной модели Изинга на

решетке пирохлора

2.5 Критические индексы ферромагнитной XY модели на решетке пирохлора

2.6 Критические индексы ферромагнитной модели Гейзенберга на решетке пирохлора

2.7 Универсальность конечно-размерного скейлинга

2.8 Краткие выводы по главе

Стр.

Глава 3. Критическое поведение намагниченности

фрустрированных решеток Изинга с антиферромагнитным взаимодействием во внешнем поле

3.1 Модель и вычислительный алгоритм

3.2 Критическое поведение намагниченности антиферромагнитной модели Изинга на треугольной решетке

3.3 Критическое поведение намагниченности антиферромагнитной модели Изинга на решетке кагоме

3.4 Краткие выводы по главе

Глава 4. Критическое поведение остаточной энтропии

антиферромагнитных фрустрированных систем во

внешнем поле

4.1 Критическое поведение в моделях без разбавлений

4.2 Критические явления в антиферромагнитной модели Изинга с разбавлением на решетке пирохлора

4.3 Алгоритм Ванга-Ландау и расчёт остаточной энтропии

4.4 Критическое поведение намагниченности и остаточной энтропии неразбавленных антиферромагнитных фрустрированных систем

во внешнем магнитном поле

4.4.1 Добавочные критические явления в разбавленной антиферромагнитной модели Изинга на решетке пирохлора

4.4.2 Добавочные критические явления в раззбавленной антиферромагнитной модели Изинга на решетке кагоме

4.4.3 Добавочные критические явления в разбавленной антиферромагнитной модели Изинга на треугольной решетке

4.4.4 Природа больших пиков энтропии и плато намагниченности на примере решетки пирохлора

4.5 Краткие выводы по главе

Заключение

Список литературы

Стр.

Приложение А. Высокопроизводительные вычисления и алгоритмы

Приложение Б. Параллельный алгоритм точного расчета коэффициентов статистической суммы на синтаксисе языка C++

Приложение В. Параллельный GPU алгоритм для модели Изинга на

решетке пирохлора на синтаксисе языка C++

Список сокращений и условных обозначений

ФМ, FM ферромагнетик, ferromagnetic

АФМ, AFM антиферромагнетик, antiferromagnetic

ПМ, PM парамагнетик, paramagnetic

СС, SG спиновое стекло, spin glass

ИСЛ, ASI искусственный спиновый лед, artificial spin ice

ПГУ, PBC периодические граничные условия, periodic boundary conditions

ПП, OP параметр порядка, order parameter

ОС, GS основное состояние, ground state

ПВС, DOS плотность вероятности состояний, density of states

МК, MK алгоритм Монте-Карло, Monte Carlo algorithm

ВЛ, WL алгоритм Ванга-Ландау, Wang and Landau algorithm

ФП, PT фазовый переход, phase transition

ГПУ, GPU графическое процессорное устройство, graphics processing unit

ЦПУ, CPU центральное процессорное устройство, central processing unit

Введение диссертации (часть автореферата) на тему «Критические явления в термодинамике спиновых систем»

Введение

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

Универсальность и масштабируемость являются двумя важнейшими свойствами при изучении критических явлений [1; 2]. Критические явления - это собирательное название класса эффектов, которые наблюдаются при особых условиях, например, при температуре Кюри в магнетиках, или при достижении критических значений внешних магнитных полей. Критические явления, связанные с фазовыми переходами второго рода, подразделяются на ограниченное число классов универсальности в зависимости от специфических параметров материала, фундаментальной симметрии системы, её пространственной размерности и числа компонент параметра порядка.

К основным методам статистической физики для изучения фазовых переходов и критических явлений относятся аналитические решения, численные расчеты Монте-Карло [3], а также теория ренормализационной группы [4].

Актуальность исследования критических явлений в физике магнетизма обусловлена необходимостью развития модельных представлений о магнитных явлениях и природе ферромагнетизма, что предоставляет возможности для практического использования этого явления. Знание свойств высоко-коррелированных систем открывает возможности для решения широкого класса фундаментальных задач в области физики, математики, информатики, химии и т.д. Большое внимание исследователей к критическим явлениям в термодинамике спиновых систем обусловлено наличием новых интересных свойств, которыми обладают магнитные материалы. Открытие соединений спинового льда, таких как Бу2Т1207 и Ио2Т1207 на решетке пирохлора ускорило изучение фрустраций [5; 6]. Существование остаточной макроскопической

энтропии является открытым вопросом, который представляет интерес для развития физики критических явлений фрустрированных систем. Впервые макроскопическое вырождение обсуждалось Паулингом для водяного льда [7]. В материалах спинового льда магнитные ионы (Бу3+ или Но3+) занимают узлы решетки пирохлора, которая состоит из тетраэдров с общими вершинами. Поле локального кристаллического окружения выравнивает магнитные моменты в направлениях, соединяющих центры двух тетраэдров при низких температурах [8; 9]. В низкотемпературном состоянии спинового льда магнитные моменты сильно ограничены локально и подчиняются, так называемым «правилам льда», как и в водяном льду, т. е. в основном состоянии два спина направлены внутрь и два наружу в каждом тетраэдре решетки пирохлора.

Эффекты критического поведения намагниченности и остаточной энтропии под действием магнитного поля, и в частности, наличие нескольких плато намагниченности, были изучены теоретически [10-12] и экспериментально [13-16]. В работе [12] Исаков и др. наблюдали критические явления, сопровождающиеся большим пиком энтропии между двумя плато.

Взаимовлияние эффектов разбавления и фрустраций является также актуальной темой исследования критического поведения спинового льда. Ке и др. [17] изучали эффекты разбавления путем замещения магнитных ионов Бу3+ или Но3+ немагнитными ионами Y3+. Экспериментально обнаружено немонотонное поведение «остаточной» энтропии как функции концентрации разбавления, также было рассмотрено обобщение теории остаточной энтропии Паулинга [17]. Результаты экспериментальных исследований в сочетании с Монте-Карло моделированием были приведены в работах [18; 19].

Целью диссертационной работы является фундаментальное изучение термодинамики фрустрированных спиновых систем в критических точках с помощью численных расчетов. Для достижения поставленной цели были сформулированы и решены следующие задачи:

1. Произвести высокопроизводительные вычисления для ферромагнитной модели Изинга на решетке пирохлора. Определить критическую температуру методом конечно-размерного скейлинга. Вычислить критические температуры для классической XY модели и классической модели Гей-зенберга на решетке пирохлора.

2. Исследовать критическое поведение ферромагнитной модели Изинга на решетке пирохлора. Произвести расчет критических индексов корреля-

ционной длины и намагниченности. Сравнить полученные результаты с данными для простой кубической решетки. Исследовать универсальность скейлинговых функций. Выполнить аналогичные расчеты для классической ХУ модели и классической модели Гейзенберга на решетке пирохлора.

3. Исследовать критическое поведение намагниченности систем «кагоме-льда» во внешнем магнитном поле, произвести сравнительный анализ с антиферромагнитной моделью Изинга на двумерной решетке кагоме и треугольной решетке.

4. Исследовать природу происхождения больших пиков энтропии и плато намагниченности при критических значениях внешнего магнитного поля в модели Изинга на решетке пирохлора. Выполнить аналогичный анализ для модели Изинга на треугольной и кагоме решетках.

Научная новизна:

1. Методом конечно-размерного скейлинга получены значения критических индексов и критических температур для классической ферромагнитной модели Изинга, ХУ модели и классической модели Гейзенберга на решетке пирохлора.

2. Вычислена количественная мера универсальности скейлинговых функций и вычислены метрические коэффициенты в ферромагнитной модели Изинга на решетке пирохлора. Выполнены аналогичные расчеты для классической ХУ модели и классической модели Гейзенберга на решетке пирохлора.

3. Исследовано поведение намагниченности АФМ модели Изинга на двумерной решетке кагоме и треугольной решетке. Наблюдается скачкообразный переход в плато намагниченности при критических значениях внешнего поля. Этот эффект возникает вследствие конкуренции между обменной энергией и энергией Зеемана, которые имеют различные значения для каждой решетки.

4. Установлена природа критических явлений, сопровождающихся большими пиками остаточной энтропии и плато намагниченности АФМ систем спинового льда во внешнем поле.

Фундаментальная значимость работы состоит в том, что разработанный универсальный метод конечно-размерного скейлинга для численного расчета критических индексов и критических температур в ферромагнитных спино-

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

Практическая значимость подтверждается набором разработанных и зарегистрированных в Роспатенте РФ высокопроизводительных программных комплексов и пакетов программ для ЭВМ [20-24], которые могут быть использованы при исследовании численными методами различных моделей спиновых систем, например, систем с конкурирующими взаимодействиями или большим координационным числом решетки.

Основные положения, выносимые на защиту:

1. В классических ферромагнитных моделях на решетке пирохлора: критическая температура Тс = 4.21394, критический индекс корреляционной длины V = 0.629 и критический индекс намагниченности в = 0.327 в модели Изинга, Тс = 2.02850, V = 0.672 и в = 0.346 в модели ХУ, Тс = 1.31695, V = 0.711 и в = 0.366 в модели Гейзенберга.

2. Значения метрических коэффициентов 0\ и решетки пирохлора равны 1.0 и 0.91 в классической ферромагнитной модели Изинга, 0.96 и 0.91 в классической ферромагнитной ХУ модели, 0.92 и 0.91 в классической ферромагнитной модели Гейзенберга, соответственно.

3. Скачкообразный переход между плато намагниченности в антиферромагнитной модели Изинга на треугольной решетке происходит при критическом значении внешнего поля к/З = 6 и при критическом значении внешнего поля к/З = 4 в антиферромагнитной модели Изинга на решетке кагоме. Разбавление этих структур вакансиями приводит к появлению добавочных скачков намагниченности.

4. Критическое поведение остаточной энтропии и намагниченности в разбавленной антиферромагнитной модели Изинга на решетке пирохлора во внешнем поле имеет место из-за дискретной структуры плотности состояний, которая трансформируется, когда внешнее магнитное поле воздействует на систему. Разбавление системы вакансиями приводит к добавочным пикам энтропии и скачкам намагниченности при значениях внешнего поля к/З = 3,6,9,12.

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

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

1. Конференция Far Eastern Natural Sciences 2015, Дальневосточный Федеральный Университет, Владивосток, 2015.

2. Конференция SPB OPEN 2016, Санкт-Петербургский Академический Университет РАН, Санкт-Петербург, 2016.

3. Семинар Tokyo Metropolitan University, Токио, 2017.

4. Конференция DPG Spring Meeting 2018, Technical University Berlin, Берлин, 2018.

5. Семинар стипендиатов программ «Михаил Ломоносов» и «Иммануил Кант», центр «Digital October», Москва, 2018.

6. Конференция ASCO-NANOMAT 2018, Дальневосточный Федеральный Университет, Владивосток, 2018.

Личный вклад. Содержание диссертации и основные положения, выносимые на защиту, отражают персональный вклад автора в опубликованные работы. Задачи, представленные в диссертации, были решены автором. Вклад автора в работы, выполненные в соавторстве, считается равнозначным.

Публикации. Основные результаты по теме диссертации изложены в 5 печатных изданиях, из них 4 изданы в журналах, индексируемых Web of Science, 1 статья проиндексирована в ВАК. Получены 5 свидетельств на программы ЭВМ. Общий список работ, в которых опубликованы материалы диссертации, содержит 10 работ.

Объем и структура работы. Диссертация состоит из введения, четырёх глав, заключения и трёх приложений. Полный объём диссертации составляет 118 страниц, включая 42 рисунка и 3 таблицы. Список литературы содержит 69 наименований.

Глава 1. Критические явления и фазовые переходы в векторных моделях

В данной главе приводятся основные понятия статистической физики систем взаимодействующих частиц, а также представлены результаты исследования плотности состояний простой квадратной решетки Изинга с взаимодействием между ближайшими соседями. Для численных расчетов был использован параллельный алгоритм, осуществляющий полный перебор всех возможных конфигураций системы. Были исследованы решетки размерами 4 х 4, 5 х 5 и 6 х 6 с учетом периодических граничных условий. Основные результаты, полученные автором, опубликованы в работе [25].

1.1 Термодинамические состояния и средние величины

Сочетание значений температуры Тс и давления Рс, при которых исчезает различие в свойствах жидкой и газообразной фаз вещества называется критической точкой. В критической точке плотность жидкости и её насыщенного пара становятся равны, а поверхностное натяжение жидкости падает до нуля, поэтому граница разделения жидкость-пар исчезает. Подобные критические явления наблюдаются во многих других системах: магнит теряет свои ферромагнит-

и Т /* V-/ V-/ и

ные свойства в точке Кюри, жидкий гелий теряет свойство сверхтекучести в Л-точке. Если определенное характерное свойство системы полностью исчезает при некоторой температуре, то говорят о фазовом переходе второго рода. Если же превращение захватывает одну область образца за другой, то такой переход называют фазовым переходом первого рода. Было замечено, что различные фазовые переходы второго рода имеют много общих особенностей в поведении [26]. Статистические свойства различных веществ вблизи критической точки подчиняются поразительно простым общим законам, однако, статистическая теория критических явлений полностью до сих пор не разработана.

Впервые фазовый переход в твердом теле наблюдал Гопкинсон в 1889 году [27]. Позднее, в 1895 г. Пьер Кюри подтвердил существование критической точки и экспериментальным путем определил температурную зависимость маг-

нитной восприимчивости. С тех пор для критической температуры магнитных и диэлектрических веществ используют название «точка Кюри».

В ферро - и антиферромагнетиках фазовый переход представляет собой коллективное упорядочение направлений элементарных магнитных моментов (спинов) [28]. Появление такого упорядочения может вызываться как обменным взаимодействием частиц, так и внешними источниками. Различные геометрические структуры магнетиков имеют общую природу, определяемую статистическими свойствами системы частиц в целом. Возникает естественное предположение, что значение температуры фазового перехода определяется спецификой самой системы, в то время как поведение различных термодинамических (а быть может и кинетических) величин вблизи критической точки является общим свойством всех систем многих частиц [29]. Для проверки правильности этого предположения следовало бы, в различных моделях систем многих частиц, отличающихся законами взаимодействия, и пользуясь методами статистической физики и численного моделирования обнаружить фазовый переход и сравнить термодинамическое поведение различных систем вблизи температуры фазового перехода. Статистическая физика имеет дело с усредненными характеристиками физической системы. Очевидными примерами являются атмосфера внутри комнаты, вода в чайнике, или атомы в постоянном магните. Такие системы составлены из огромного числа индивидуальных компонент (атомов или молекул). Наблюдатель способен контролировать состояние компонент системы лишь в незначительной степени, и то не всегда. Он может лишь измерить небольшое число усредненных величин, характеризующих систему - температуру, плотность или намагниченность. Задача статистической физики состоит в том, чтобы предсказывать соотношения между наблюдаемыми макроскопическими величинами, располагая лишь данными о микроскопических силах, действующих между компонентами системы [30]. Для вычисления термодинамических функций, исходя из миксроскопических сил, дейстующих между компонентами системи, рассмотрим систему с набором состояний 8 и гамильтонианом Еа. Рассмотрим статистическую сумму

где к - постоянная Больцмана; Суммирование выполняется по всем допустимым состояниям системы. Тогда выражение для свободной энергии ^ имеет вид

(1.1)

в

^ = -кТ 1п

(1.2)

а вероятность найти систему в состоянии в равна

г-1ехр-^-. (О)

Если X - некоторая измеряемая характеристика системы, такая как полная энергия или намагниченность, имеющая значение для состояния б, то ее измеренное среднее термодинамическое значение равно

ж—ч —Е (X) = ЯХаехр . (1.4)

в

Таким образом, главная задача равновесной статистической механики состоит в вычислении суммы по состояниям. Для систем с непрерывным спектром эта сумма превращается в интеграл, а для квантово-механических систем — в сумму диагональных элементов матрицы плотности. Результат такого вычисления дает Z и Е как функции и любых других переменных, входящих в Е3, например магнитного поля. Термодинамические и магнитные характеристики затем можно получить посредством дифференцирования. К сожалению, для любых систем макроскопического размера вычисление Z представляет собой невероятно трудную задачу, которая ввиду своей математической сложности может быть решена лишь в некоторых простейших случаях. При решении данной задачи желательно выбрать как можно более простую модель, чтобы получить точное решение. Вместе с тем, модель после статистического расчета должна привести к появлению фазового перехода. Суть подхода к решению состоит в том, чтобы попытаться найти модели, для которых Еа достаточно простая функция, чтобы статистическая сумма (1.1) могла быть вычислена точно. Возможно, что это не даст нам информации о значениях термодинамических функций реальных систем, но мы получим качественное представление о том, как могут вести себя системы, в особенности вблизи Тс.

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

1.2 Двумерная модель Изинга с взаимодействием между ближайшими

соседями

Представим себе магнетик, состоящий из молекул, расположенных в узлах регулярной решетки. Предположим, что имеется N таких узлов, пронумерованных индексами. Представим теперь каждую молекулу, как микроскопическую магнитную стрелку, которая направлена всегда либо вдоль некоторой заданной оси, либо точно в противоположном направлении. Таким образом, каждая молекула имеет две возможные конфигурации, которые можно описывать с помощью «спиновой» переменной 0{, принимающей значения +1 (параллельно оси) или -1 (антипараллельно оси) [30].

Таким образом, модель Изинга содержит четыре основных предположения: 1) единичные векторы фиксированы в узлах кристаллической решетки; 2) каждый вектор может находиться лишь в двух состояниях, скажем, для ¿-го узла решетки эти состояния характеризуются переменной 0{, принимающей два значения: +1 или -1 (для ферромагнетика — два направления элементарного магнитного момента); 3) учитывается взаимодействие лишь между расположенными по соседству единичными векторами, т. е. энергия всей системы

Е = -3 ^ ого3, (1.5)

м

где г, ] - соседние узлы решетки. 4) Пусть о = [о"1...,о>] обозначает набор значений всех N спинов. Тогда имеется 2м различных наборов 0{, и каждый из них описывает некоторое состояние системы.

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

тогда Гамильтониан системы N частиц во внешнем поле равен

3 м м

Н = -- ^ аг,3ъъ - цН ^ а,, (1.6)

" г,3=1 3 =1

где а^ = 1 если узлы г и у ближайшие соседи, и а= 0 во всех остальных случаях.

Гамильтониан представляет собой функцию, состоящую из двух частей. Первая часть включает вклад межмолекулярных сил внутри магнетика, а вторая — вклад от взаимодействия каждого спина с внешним магнитным полем. Ферромагнитному упорядочению соответствует 3 > 0, антиферромагнитному 3 < 0. Разумеется, термодинамические и магнитные свойства системы можно определить из статистической суммы, которая в данном случае принимает вид

3 М Н м

2 = X) ... X) ехр(кТ2^+ а^ (1Л)

ст»=±1 ах=±1 3=1

Например, внутренняя энергия и и среднеквадратичная намагниченность М, отнесенные к одной частице, равны, соответственно,

" = £ ^ С«)

« = ¡1 ^ ('■»

Теплоемкость в расчете на один спин для системы, состоящей из пары спинов Изинга вычисляется как

С = 1 (1.10)

2 6 Т ^

В любой физической системе мы предполагаем все взаимодействия инвариантными по отношению к обращению времени.

Трудность вычисления статистической суммы для решетки состоящей из конечного числа спинов Изинга в общем случае

3 М Н м

я = ... X) ехр(а^а'аз + а^ (1.11)

Сг=±1 ам=±1 г,] 3=1

заключается в расчете кратности вырождения спинового избытка М^ (разницы между числом спинов вниз и числом спинов вверх) по энергии Е^ (сумма энергий всех парных взаимодействий). Стоит также отметить, что сумма кратностей

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

- Задается число битовых строк, столбцов, общее число спинов.

- Инициализируется векторный массив для хранения информации о кратности вырождения спинового избытка по энергии, а также векторный массив для хранения сортированной информации.

- Инициализируется массив, содержащий целочисленные переменные di, набор двоичных представлений которых соответствует столбцам решетки. -Инициализируется массив bi, необходимый для хранения числа единиц, содержащихся в битовом столбце, а также массив десятичных значений cj, соответствующих сдвинутым кольцевым образом битовым массивам.

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

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

- Расчет энергии обменных парных взаимодействий между спинами в узлах столбцов решетки выполняется следующим образом: применяется логическое «ИЛИ» между всеми парами di и щ+\, включая периодические граничные условия. Затем, к предыдущему результату применяется логическое «НЕ». Битовые значения, которые не несут в себе информацию о системе обнуляются путем вычитания из вектора результата 232 — 2П (степень 32 для «int» типа); Число единиц в результирующем битовом векторе, т.е. число пар спинов, обменная энергия которых отрицательна, сохраняется в переменной Еп.

- Подсчет энергии в строках ai(di;j) может быть выполнен следующим образом: каждое десятичное число «j» - индекс элемента массива целочисленных переменных Cj - ставится в соответствие десятичному числу «Cj», бинарное представление которого есть циклический сдвиг значащих бит числа «j». Применяется логическое «ИЛИ» между всеми парами целых чисел di и Cj, где j = di, к предыдущему результату применяется «НЕ», незначащие битовые значения обнуляются. Ссылка на элемент массива bi используется для определения количества единиц и суммирования с Еп.

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

Etot = 2(-Еп + N). (1.12)

Правильность полученных результатов расчета можно проверить по биномиальным коэффициентам. Известно, что сумма кратностей вырождения спинового избытка Mi фиксированного значения по всем разрешенным для него значениям энергии Ei есть биномиальный коэффициент, соответствующий определенному разбиению номера строки в треугольнике Паскаля. С помощью данного метода были вычислены статистические суммы для квадратичных решеток, состоящих из 16 и 25 и 36 спинов Изинга. При дальнейшем увеличении размеров системы, время расчета растет экспоненциально. Например, для решетки 10 х 10, необходимо осуществить полный перебор 2100 конфигураций, на что даже для экзафлопсной вычислительной системы (1018 Mips) потребуется около пенталиона (1015) секунд, что составляет 31,7 миллионов лет [32].

1.3 Расчет плотности вероятности состояний в модели Изинга на

квадратной решетке

Чтобы получить плотность вероятности состояний феррмагнитной модели Изинга на квадратной решетке, был произведен полный перебор всех возможных конфигураций. Были получены графики ПВС относительно средней намагниченности Mi и относительно значений энергии Ei. Отношение числа вырождения Ni какого-либо состояния к общему числу 2N конфигураций - есть вероятность появления данного состояния Р при Т ^ ж. Данные характеристики в наглядной форме представлены на графиках ПВС (рис. 1.1, 1.2, 1.3).

Большая часть состояний сосредоточена в близких к нулевому значению энергетических областях, и вблизи значений средней намагниченности Mi = 0. Такие конфигурации имеют наибольшую вероятность быть полученными с помощью приближенных методов, к примеру, методом случайных блужданий или методом Монте-Карло. Самым низковерятным при Т ^ ж является основное энергетическое состояние, т.е. когда все спины имеют одинаковую ориентацию.

ф +

—+-

4-тг-

-30 -20

о. 10"3

+ * * :

+ +

* +

: * Ф ф

+ + ——Ф-

-15 -10 -5

0 М,

Рисунок 1.1 — Плотность вероятности состояний модели Изинга размером Ь = 4 х 4 на квадратной решетке по энергии (слева) и намагниченности (справа)

+ +-Н-+ -Н-

* : + + + + + t + + +

1 1 - + + + + + + + + +

$ ф + + + +

+ + + + ф +

+ +

-10

Е,-

+ ++ Г+ + Ф + $ + + фф++Ф + + +

ф + + Ф*Ф + : 4- + + + + + + + + + ф

т т + + + +Ф + + Т ф+ + т .•V

+ + + + + + + + + + + + + + + + + + + + + +

+ + + + + +

+ +

0 М,

Рисунок 1.2 — Плотность вероятности состояний модели Изинга размером Ь = 5 х 5 на квадратной решетке по энергии (слева) и намагниченности (справа)

Стоит также отметить, что для решетки 5x5, отсутствуют значения средней намагниченности М^ = 0 и энергии Е^ = 0. Это свойство связано с нечетностью длины стороны решетки, и оно будет выполняться для любых таких квадратичных решеток (например 7 х 7, 9 х 9 и т.д.).

Список литературы диссертационного исследования кандидат наук Солдатов Константин Сергеевич, 2019 год

Список литературы

1. Stanley H.E. Introduction to Phase Transitions and Critical Phenomena. — New York: Oxford University Press, 1971. — 336 pp.

2. Hu Chin-Kun. Historical Review on Analytic, Monte Carlo, and Renormalization Group Approaches to Critical Phenomena of Some Lattice Models // Chinese J. Phys. — 2014. — no. 52. — P. 76.

3. Metropolis N., Rosenbluth A.W. et al. Equation of State Calculations by Fast Computing Machines // J. Chem. Phys. — 1953. — no. 21. — P. 6.

4. Wilson K. G. Renormalization Group and Critical Phenomena. I. Renormalization Group and the Kadanoff Scaling Picture // Phys. Rev. B. — 1971. — no. 4. — P. 9.

5. Harris M. J., Bramwell S. T. et al. Geometrical Frustration in the Ferromagnetic Pyrochlore Ho2Ti2C>7 // Phys. Rev. Lett. — 1997. — Vol. 79. — Pp. 2554-2557.

6. Ramirez A. P., Hayashi A. et al. Zero-point entropy in 'spin ice' // Nature. — 1999. — Vol. 399. — P. 333-335.

7. Pauling L. The Structure and Entropy of Ice and of Other Crystals with Some Randomness of Atomic Arrangement // Journal of the American Chemical Society. — 1935. — Vol. 57, no. 12. — Pp. 2680-2684.

8. Bramwell S. T., Gingras M. J. P. Spin Ice State in Frustrated Magnetic Pyrochlore Materials // Science. — 2001. — Vol. 294, no. 5546. — Pp. 1495-1501.

9. Diep H. T. Frustrated Spin Systems. — WORLD SCIENTIFIC, 2013.

10. Harris M. J., et al. Bramwell. Liquid-Gas Critical Behavior in a Frustrated Pyrochlore Ferromagnet // Phys. Rev. Lett. — 1998. — Vol. 81. — Pp. 4496-4499.

11. Rashba E. I. Spin currents in thermodynamic equilibrium: The challenge of discerning transport currents // Phys. Rev. B. — 2003. — Vol. 68. — P. 241315.

12. Maultzsch J., Reich S., Thomsen C. Double-resonant Raman scattering in graphite: Interference effects, selection rules, and phonon dispersion // Phys. Rev. B. — 2004. — Vol. 70. — P. 155403.

13. Hiroi Z., Matsuhira K. et al. Specific Heat of Kagome Ice in the Pyrochlore Oxide Dy2Ti2O7 // Journal of the Physical Society of Japan. — 2003. — Vol. 72, no. 2.

- Pp. 411-418.

14. Sakakibara T., Tayama T. et al. Observation of a Liquid-Gas-Type Transition in the Pyrochlore Spin Ice Compound Dy2Ti2O7 in a Magnetic Field // Phys. Rev. Lett. — 2003. — May. — Vol. 90. — P. 207205.

15. Fukazawa H., Melko R. G., Higashinaka et al. Magnetic anisotropy of the spin-ice compound Dy2Ti2O7 // Phys. Rev. B. — 2002. — Vol. 65. — P. 054410.

16. How H., Vittoria C. Surface retarded modes in multilayered structures: Parallel magnetization // Phys. Rev. B. — 1989. — Vol. 39. — Pp. 6823-6830.

17. Ke X, Freitas R. S. et al. Nonmonotonic Zero-Point Entropy in Diluted Spin Ice // Phys. Rev. Lett. — 2007. — Vol. 99. — P. 137203.

18. Lin T., Ke X. et al. Nonmonotonic residual entropy in diluted spin ice: A comparison between Monte Carlo simulations of diluted dipolar spin ice models and experimental results // Phys. Rev. B. — 2014. — Vol. 90. — P. 214433.

19. Scharffe S., Breunig O. et al. Suppression of Pauling's residual entropy in the dilute spin ice Dy2Ti2O7 // Phys. Rev. B. — 2015. — Vol. 92. — P. 180405.

20. Солдатов К.С., Нефедев К.В., Капитан В.Ю. Расчет разрешенных значений энергии и спинового избытка в 2-D модели Изинга для простой квадратной решетки, Св. №2015660080 от 21.09.2015г // Оф. бюллетень «Программы для ЭВМ. Базы данных. Топологии интегральных микросхем». — 2015.

21. Солдатов К.С., Нефедев К.В., Капитан В.Ю. Параллельный расчет коэффициентов статистической суммы 2-D модели Изинга с конечным числом частиц, Св. №2015660128 от 22.09.2015 // Оф. бюллетень «Программы для ЭВМ. Базы данных. Топологии интегральных микросхем». — 2015.

22. Солдатов К.С., Нефедев К.В., Капитан В.Ю. Моделирование магнитных многослойных структур, Св. №2017611233 от 12.01.2017г // Оф. бюллетень «Программы для ЭВМ. Базы данных. Топологии интегральных микросхем».

— 2017.

23. Солдатов К.С., Нефедев К.В., др. В.Ю. Капитан и. Двумерная плотность состояний сложных спиновых систем, Св. №2017610581 от 12.01.2017г // Оф. бюллетень «Программы для ЭВМ. Базы данных. Топологии интегральных микросхем». — 2017.

24. Солдатов К.С., Нефедев К.В., др. В.Ю. Капитан и. Программа создания образцов диполей, Св. №2017612886 от 01.06.2017г // Оф. бюллетень «Программы для ЭВМ. Базы данных. Топологии интегральных микросхем». — 2017.

25. Soldatov K. et al. Approaches to numerical solution of 2D Ising model // Journal of Physics: Conference Series. — 2016. — Vol. 741.

26. Aharoni A. An Introduction to the Theory of Ferromagnetism. — 2000. — 11.

27. Hopkinson J. Magnetic and Other Physical Properties of Iron at a High Temperature // Philosophical Transactions of the Royal Society of London Series A. — 1889. — Vol. 180. — Pp. 443-465.

28. Isihara A. Statistical Physics. — Academic Press, 1971.

29. Кринчик Г. С. Физика магнитных явлений. — Издательство Московского университета, 1976.

30. Baxter R. J.Exactly Solved Models in Statistical Mechanics. — London: Academic Press, 1982.

31. Onsager L. Crystal Statistics. I. A Two-Dimensional Model with an Order-Disorder Transition // Phys. Rev. — 1944. — Vol. 65.

32. Nefedev K., Peretyatko A. Superlinear Speedup of Parallel Calculation of Finite Number Ising Spins Partition Function // HPC - UA. — 2013.

33. Soldatov K., Nefedev K., Okabe Y. et al. Large-scale calculation of ferromagnetic spin systems on the pyrochlore lattice // Physics Letters A. — 2017. — Vol. 381. — 6 pp.

34. Rushbrooke G. S. On the Thermodynamics of the Critical Region for the Ising Problem // The Journal of Chemical Physics. — 1963. — no. 39. — Pp. 842-843.

35. Griffiths Robert B. Thermodynamic Inequality Near the Critical Point for Ferro-magnets and Fluids // Phys. Rev. Lett. — 1965. - Vol. 14. - Pp. 623-624.

36. Griffiths R. B. Ferromagnets and Simple Fluids near the Critical Point: Some Thermodynamic Inequalities // The Journal of Chemical Physics. — 1965. — Vol. 43. — Pp. 1958-1968.

37. Fisher M. E. The theory of equilibrium critical phenomena // Reports on Progress in Physics. — 1967. — Vol. 30. — Pp. 615-730.

38. Fisher M. E. Magnetic Critical Point Exponents—Their Interrelations and Meaning // Journal of Applied Physics. — 1967. — Vol. 38. — Pp. 981-990.

39. Stauffer D. Scaling theory of percolation clusters // Physics Reports. — 1979. — Vol. 54. — Pp. 1 - 74.

40. Binder K. Finite size scaling analysis of ising model block distribution functions //

Zeitschrift für Physik B Condensed Matter. — 1981. — no. 43. — P. 22.

41. Domb C., Green M. S., Lebowitz J. L. Phase Transitions and Critical Phenomena, vol. 8. — London: Academic Pres, 1983. — 518 pp.

42. Blote H. W. J., Luijten E., Heringa J. R. Ising universality in three dimensions: a Monte Carlo study // Journal of Physics A: Mathematical and General. — 1995.

— no. 28. — P. 25.

43. Campostrini M., Hasenbusch M. et al. Critical behavior of the three-dimensional XY universality class // Phys. Rev. B. — 2001. — no. 63. — P. 28.

44. Gottlob A. P., Hasenbusch M. Critical behaviour of the 3D XY-model: a Monte Carlo study // Physica A: Statistical Mechanics and its Applications. — 1993. — no. 201. — P. 21.

45. Campostrini M., Hasenbusch M. et al. Critical exponents and equation of state of the three-dimensional Heisenberg universality class // Phys. Rev. B. — 2002. — no. 65. — P. 21.

46. Holm C., Janke W. Critical exponents of the classical three-dimensional Heisenberg model: A single-cluster Monte Carlo study // Phys. Rev. B. — 1993. — no. 48.

— P. 15.

47. Soldatov K., Peretyatko A. et al. Comparison of diluted antiferromagnetic Ising models on frustrated lattices in a magnetic field // Physics Letters A. — 2019. — Vol. 383, no. 12. — Pp. 1229 - 1234.

48. Peretyatko A., Nefedev K., Okabe Y. Interplay of dilution and magnetic field in the nearest-neighbor spin-ice model on the pyrochlore lattice // Phys. Rev. B. — 2017. — no. 95. — P. 5.

49. Hukushima K., Nemoto K. Exchange Monte Carlo Method and Application to Spin Glass Simulations // Journal of the Physical Society of Japan. — 1995. — no. 65. — P. 5.

50. Andriushchenko P., Soldatov K. et al. Large peaks in the entropy of the diluted nearest-neighbor spin-ice model on the pyrochlore lattice in a [111] magnetic field // Phys. Rev. E. — 2019. — Vol. 99. — 12 pp.

51. Wannier G. H. Antiferromagnetism. The Triangular Ising Net // Phys. Rev. — 1950. — Vol. 79. — Pp. 357-364.

52. Kano K., Naya S. Antiferromagnetism. The Kagome Ising Net // Progress of Theoretical Physics. — 1953. — Vol. 10, no. 2. — Pp. 158-172.

53. Udagawa M., Ogata M., Hiroi Z. Exact Result of Ground-State Entropy for Ising Pyrochlore Magnets under a Magnetic Field along [111] Axis // Journal of the Physical Society of Japan. — 2002. — Vol. 71, no. 10. — Pp. 2365-2368.

54. Metcalf B. D., Yang C. P. Degeneracy of antiferromagnetic Ising lattices at critical magnetic field and zero temperature // Phys. Rev. B. — 1978. — Vol. 18. — Pp. 2304-2307.

55. Baxter R. J., Tsang S.K. Entropy of hard hexagons // Journal of Physics A: Mathematical and General. — 1980. — Vol. 13, no. 3. — Pp. 1023-1030.

56. Nagle John F. New Series-Expansion Method for the Dimer Problem // Phys. Rev. — 1966. — Vol. 152. — Pp. 190-197.

57. Nagle J. F. Lattice Statistics of Hydrogen Bonded Crystals. I. The Residual Entropy of Ice // Journal of Mathematical Physics. — 1966. — Vol. 7, no. 8. — Pp. 1484-1491.

58. Wang F., Landau D. P. Determining the density of states for classical statistical models: A random walk algorithm to produce a flat histogram // Phys. Rev. E. — 2001. — Vol. 64. — P. 056101.

59. Shevchenko Y., Nefedev K., Okabe Y. Entropy of diluted antiferromagnetic Ising models on frustrated lattices using the Wang-Landau method // Phys. Rev. E. — 2017. — Vol. 95. — P. 052132.

60. Isakov S. V., Gregor K., Moessner et al. Dipolar Spin Correlations in Classical Pyrochlore Magnets // Phys. Rev. Lett. — 2004. — Vol. 93. — P. 167204.

61. Ferreyra M. V., Giordano G. et al. Thermodynamics of the classical spin-ice model with nearest neighbour interactions using the Wang-Landau algorithm // The European Physical Journal B. — 2016. — Vol. 89, no. 2. — P. 51.

62. Schulz B. J., Binder K. et al. Avoiding boundary effects in Wang-Landau sampling // Phys. Rev. E. — 2003. — Vol. 67. — P. 067102.

63. Berg B. A., Hansmann U., Neuhaus T. Simulation of an ensemble with varying magnetic field: A numerical determination of the order-order interface tension in the D=2 Ising model // Phys. Rev. B. — 1993. — Vol. 47. — Pp. 497-500.

64. Zhou C., Bhatt R. N. Understanding and improving the Wang-Landau algorithm // Phys. Rev. E. — 2005. — Vol. 72. — P. 025701.

65. Belardinelli R. E., Manzi S., Pereyra V. D. Analysis of the convergence of the Wang-Landau algorithms in the calculation of multidimensional integrals // Phys. Rev. E. — 2008. — Vol. 78. — P. 067701.

66. Okabe Y., Tomita Y, Yamaguchi C. Application of new Monte Carlo algorithms to random spin systems // Computer Physics Communications. — 2002". — Vol. 146, no. 1. — Pp. 63 - 68.

67. Zukovic M., Borovsky M., Bobak A. Phase diagram of a diluted triangular lattice Ising antiferromagnet in a field // Physics Letters A. — 2010. — Vol. 374, no. 41.

— Pp. 4260 - 4264.

68. Федотов И. А. Модели параллельного программирования. — Москва: Салон

- Пресс, 2013.

69. Amdahl Gene M. Validity of the single processor approach to achieving large scale computing capabilities // A FIPS computer conference. — 1967. — Pp. 3-4.

Приложение А Высокопроизводительные вычисления и алгоритмы

За последние годы компьютеры и вычислительные мощности все быстрее приближаются к своему пределу тактовых частот. Производительность процессоров (точнее их отдельных ядер) перестала нарастать огромными темпами, и оптимальной возможностью для её повышения является непосредственное увеличение количества процессоров. В связи с этим, параллельное программирование стало высоко востребованной дисциплиной и с каждым днем приобретает все большую актуальность. Суть параллельного программирования заключается в построении схемы параллельных вычислений, что нередко представляет собой весьма нетривиальную задачу. Большинство параллельных вычислений реализуется в виде параллельных программ. При этом в существующем последовательном алгоритме выделяются независимые друг от друга последовательные ветви, которые могут быть выполнены параллельно, и с этим учетом пишется программа. Нередко берется уже существующая последовательная программа, которая путем добавления неких конструкций распараллеливания преобразуется в параллельно-последовательную. Чем сложнее исходная последовательная программа, чем хуже она структурирована, чем запутаннее зависимости между отдельными ее ветвями, тем сложнее оказывается выполнить ее распараллеливание.

С одной стороны, для описания параллельных программ возможно уйти от популярных на сегодняшний день императивных языков программирования (C++, Java, Pascal), описывающих конкретный порядок выполняемых действий. В данном случае требуется использование декларативных языков, т.е. таких, которые описывают непосредственно ожидаемый результат, оставляя на долю компилятора выбор порядка действий и инструкции по их выполнению. Примерами таких языков являются язык логического программирования Prolog и языки функционального программирования Haskell, Erlang, Lisp. Несмотря на потенциальную мощь, декларативные языки не приобрели широкой популярности, так как для эффективной компиляции кода в любом случае приходится использовать императивные языки. Различие между императивным и декларативным

программированием в чем-то напоминает различие между прямым вычислением функции и решением уравнения, т.е. обращения функции [68].

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

Естественно, что, используя параллельную систему с р вычислительными устройствами, пользователь ожидает получить ускорение своей программы в р раз по сравнению с последовательным вариантом. Но действительность практически всегда оказывается далека от идеала. Предположим, что структура информационных зависимостей программы определена и доля операций, которые нужно выполнять последовательно, равна £ Крайние случаи в значениях Г соответствуют полностью параллельным (Г = 0) и полностью последовательным (Г = 1) программам. Тогда для того, чтобы оценить, какое ускорение 8 может быть получено на компьютере из р процессоров при данном значении Г, можно воспользоваться законом Амдала [69]

^ < , (АЛ)

р

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

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

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

1. Написание последовательной программы;

2. Полноценная отладка последовательной программы;

3. Распараллеливание программы с использованием выбранного средства;

4. Отладка параллельной программы.

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

Стандарт обмена сообщениями Massage Passing Interface (MPI) расширяет возможности вычислителей и делает возможным параллельную алгоритмизацию, что позволяет решать задачи, которые не могли быть решены ранее. В частности, при помощи MPI был произведен точный расчет статистической суммы для решетки, состоящей из 36 спинов Изинга. В последовательном коде число вложенных циклов эквивалентно длине вектора щ, т.е. числу битовых векторов. При распараллеливании, полный перебор всех возможных конфигураций и расчет величины энергии и спинового избытка распределен между всеми задействованными процессами за исключением нулевого, в котором происходит прием результатов вычислений в виде сообщений от остальных процессов и их сборка в единый массив. Итерации программы раздаются процессам последовательными интервалами значений их номеров, при этом длина интервала равна ближайшему сверху целому от деления количества итераций на количество процессов, из за чего последний процесс получает в общем случае меньшее количество подзадач. В силу специфики последовательного алгоритма, количество итераций не кратно количеству процессов, что приводит к неравномерной загрузке вычислительного кластера.

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

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

п = аЖ + 6,0 < Ь < Ж (А.2)

Ширина интервала итераций для каждого процесса Пк, к = 0,...,Ж — 1 может быть вычислена следующим образом

I а + 1, к < &; пк = { (А.3)

I а, к ^ 6.

Первые Ь интервалов оказываются при этом шире остальных на единицу. Левая граница каждого интервала п1к, вычисляется на основе тех же соотношений:

к

п{ = ^2 Пт — Пк = +

то=0

Таким образом, на каждый процесс распределяется интервал итераций

К; пк + пк).

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

к, к < 6; 6, к ^ Ь.

(А.4)

Приложение Б

Параллельный алгоритм точного расчета коэффициентов статистической

суммы на синтаксисе языка C++

10

15

20

25

#include<iostream> #include<vector> #include<math.h> #include <fstream> #include<mpi.h> #include <algorithm> #include<time.h> #include<stdio.h> #include<stdlib.h> using namespace std; int main(int argc, char *argv[]) { int rank,size;

long long size1,size2; double starttime, endtime; MPI_Init(&argc, &argv); starttime = MPI_Wtime();

MPI_Comm_rank (MPI_COMM_WORLD, &rank); /* get current process id * /

MPI_Comm_size (MPI_COMM_WORLD, &size); /* get number of processes */

MPI_Status status;

ofstream out("output6x6.txt");

unsigned const short int nstr=6; //число битовых строк unsigned const short int ncol=6; //число битовых столбцов unsigned const long long N=pow(2.,(double)nstr);//максимальное дес

ятичное значение переменной 2Astr в битовом столбце unsigned const short int Nspin=nstr*ncol; //общее число спинов long long eones=0л0-N;//(int)(pow(2.,32.)-1)-N;//десятичное предст авление, максималь ного числа единичных битов не участвующих в расчетах, но хранящихся в переменных long long roled=N-1; //переменная кольцевого сдвига long long n=N/2; //максимальное десятичное число, для двоичного пр едставления которо го кольцевой побитовый сдвиг осуществляется простым побитовым сдвигом влево

5

35

40

45

50

55

int flag=0;

unsigned long long ones=0; //переменная для хранения числа единиц,

используется для подсчета спинового избытка unsigned long long xstr=0; //переменная для хранения числа положит ельных парных взаи модействий, используется для подсчета энерг ии

long long in[ncol];//десятичные значения переменных (битовых столб цов)

long long inn^colHN^Z/десятичные значения, соответствующие сдви

нутым кольцевым об разом битовым массивам long long d=(long long)pow(2.,(double)(nstr))-1,//маска, максималь

ное число единиц в битовом столбце long long Nconf=(long long int)pow(2.,(double)(nstr*ncol)); //по

лное число конфигу раций long long a = Nconf / (size-1), b = Nconf % (size-1); long long nloc = a +((rank-1<b)?1:0);

long long nleft = a * (rank-1) +(((rank)<b)?rank-1:b); for(long long i=0;i<ncol; i++) {in[i]=0;

for(long long j=0; j<N; j++) {

if(j<n) inn[i][j]=(j<<1); else inn[i][j]=(j<<1)-roled;

}}

long long k=0;

long long ns=0, i=0, j=0;

if(rank!=0){ //только для рангов не равных 0

vector<long long>imas; //массив для хранения информации о кратнос ти вырождении спи нового избытка по энергии, данные организован ы триадами

vector<long long> imasr,'//массив для хранения сортированной инфор мации

for(long long int nk=nleft; nk<nleft+nloc;++nk) { ns=nk;

for (long long i=0;i<ncol;i++){ in[i] = ns&d; ns = ns>>ncol;

//подсчет числа спинов "вверх" ones=0;

for(j=0; j<ncol; j++) { k=in[j];

70

75

80

85

90

95

while(k)

{k&=k-1; ++ones;

// подсчет энергии парных взаимодействий в строках xstr=0;

k=~(in[0]Ain[ncol-1]);

k-=eones;

while(k)

{ k&=k-1; ++xstr;

for(j=0; j<ncol-1; j++) { k=~(in[j]Ain[j+1]); k-=eones; while(k) { k&=k-1;

++xstr;}

// подсчет энергии парных взаимодействий в столбцах for(j=0; j<ncol; j++) { k=~(in[j]Ainn[j][in[j]]); k-=eones; while(k) { k&=k-1; ++xstr;

ones=Nspin-2*ones; xstr=(-2)*xstr+2*Nspin; flag=0;

for(long long i=1; i<imas.size(); i+=3) if(imas[i]==ones && imas[i+1]==xstr)

{imas[i-1]++; flag++; break;

}

if(flag==0) {imas.push_back(1); imas.push_back(ones);

imas.push_back(xstr); }

//сортировка

110

115

120

125

130

135

for(long long j=Nspin; j>=-Nspin; j--) for(long long i=0; i<imas.size(); i+=3) if(j==imas[i+1]) {imasr.push_back(imas[i]);

imasr.push_back(imas[i+1]); imasr.push_back(imas[i+2]);

}

/*

long long sizev=imasr.size();

MPI_Send(&sizev,1, MPI_LONG_LONG, 0, rank, MPI_COMM_WORLD); //cout<<"rank="<<rank<<" sizev="<<sizev<<endl;

MPI_Send(&imasr.front(),sizev, MPI_LONG_LONG, 0, rank, MPI_COMM_ WORLD);

} //конец условия "только для рангов не равных 0" if(rank==0){ //только для рут

vector<long long> imasr1,'//массив для хранения результатов вычисле ний

MPI_Recv(&size1, 1, MPI_LONG_LONG, 1, 1, MPI_COMM_WORLD, &status); imasr1.resize(size1);

MPI_Recv(&imasr1.front(), size1, MPI_LONG_LONG, 1, 1, MPI_COMM_ WORLD, &status);

vector<long long>imas2;// массив для хранения сообщений от процесс ов

for(long long rank1=2; rank1<size; rank1++) // прием сообщен

ий от остальных процессов

{

MPI_Recv(&size1, 1, MPI_LONG_LONG, rank1, rank1, MPI_COMM_WORLD,

&status); imas2.resize(size1);

MPI_Recv(&imas2.front(), size1, MPI_LONG_LONG, rank1, rank1, MPI _COMM_WORLD, &status); int flag=0;

long long imrsize=imasr1.size(); long long it=1; for(it=1; it<imas2.size(); it+=3) {flag=0;

for(long long ki=1; ki<imrsize; ki+=3) {

if(imas2[it]==imasr1[ki] && imas2[it+1]==imasr1[ki+1]) {

imasr1[ki-1]+=imas2[it-1]; flag++;

150

155

160

165

170

}

}

if(flag==0)

{imasr1.push_back(imas2[it-1]); imasr1.push_back(imas2[it]); imasr1.push_back(imas2[it+1]);

}

imas2.erase(imas2.begin(),imas2.end());

long long jj=0;

for(long long i=0; i<imasr1.size(); i++) {//cout<<imasr1[i]<<" ";

out<<imasr1[i]<<" "; jj++;

if(jj==3){

//cout<<endl; out<<endl; jj=0;}

}

for(long long i=0; i<imasr1.size(); i+=3) { k+=imasr1[i];} //cout<<k<<endl; out<<k<<endl; //system("pause");

}//конец условия "только для рут" endtime = MPI_Wtime();

printf("That took %f seconds\n",endtime-starttime);

MPI_Finalize(); return 0;}

Приложение В

Параллельный GPU алгоритм для модели Изинга на решетке пирохлора на

синтаксисе языка C++

include <stdlib.h> #include <stdio.h> #include <string.h> #include <math.h> 5 #include <assert.h>

#include <curand_kernel.h>

_global_ void device_function_spinset (int*,unsigned int

);

_global_ void set_rand_init (curandState_t*,int,

unsigned int );

10 _global_ void device_function_init_K (double, int*,

curandState_t*,

unsigned int*, unsigned int*,

unsigned int*,

unsigned int*,

unsigned int*,int*,

unsigned int

);

15 _global_ void device_function_analysis_K (unsigned int*,

unsigned int );

_global_ void device_ReduceLabels (unsigned int*,

unsigned int*,

unsigned int*,

unsigned int );

_global_ void device_function_spin_flip (int*, int*, unsigned

int*,

unsigned int

);

20 _global_ void device_function_sum (int*, int*, int*,

unsigned int*,

unsigned int*,

unsigned int );

device_ unsigned int root_find

);

(unsigned int*, int

#define rand_block 16384 // number of total thread number in random number

25

#define nx 16 #define ny nx #define nz nx #define nxy (nx*ny)

30 #define nxyz (nx*ny*nz)

int nla = (16*nxyz);

void period_sc(); void period();

35 unsigned int mem_sc = sizeof(int)*(6*nxyz); int* nn_sc = (int*) malloc(mem_sc); unsigned int mem_nn = sizeof(int)*(6*16*nxyz); int* nn = (int*) malloc(mem_nn);

40 // parameters of pyrochlore lattice int main(int argc,char** argv) { /*------------Device Init

*/

cudaSetDevice(0);

/*

Generation of Nearest Neighbor List -----*/

unsigned int N_site[6];

unsigned int mem_h_N = sizeof(int)*(16*nxyz); int* h_N_NNList = (int*) malloc(mem_h_N);

unsigned int count, Nmax_NNList, tmp;

period_sc(); period();

// Reduction of Nearest Neighbor List Nmax_NNList=0; 60 for(int i=0; i< nla; i++) {

70

75

80

85

90

95

count=0;

for(int j=0; j< 6; j++) { if(i>nn[i*6+j]){

N_site[count]=nn[i*6+j]; count++;

}

}

// Sort of the number of Nearest Neighbor for(int j=1; j< count; j++) { for(int k=0; k<j; k++) { if(N_site[j]<N_site[k]){ tmp = N_site[j]; N_site[j] = N_site[k]; N_site[k] = tmp;

}

}

for(int j=0; j< count; j++) { nn[i*6+j]=N_site[j];

}

h_N_NNList[i]=count; Nmax_NNList=max(count,Nmax_NNList); // end of for(int i=0; i< nla; i++)

// Output

for(int i=0; i< nla; i++) { printf("%d ", i);

for(int j=0; j< h_N_NNList[i]; j++) { printf("%d ", nn[i*6+j]);

}

printf("\n");

}

----------- Variables for GPU command

------------------------------*/

unsigned int threads, grid;

threads = 512;

grid = (nla+512-1)/512;

}

/

*

*

/

/

*

120

125

130

unsigned int threads_rand, grid_rand; threads_rand = 128;

grid_rand = (rand_block+128-1)/threads_rand;

float flip_timer, total_timer; cudaEvent_t start_timer, stop_timer; cudaEventCreate( &start_timer ); cudaEventCreate( &stop_timer );

/*------------Size of array

------------------------------------------*/

= sizeof(int)*(nla); = sizeof(int)*(nla*Nmax_

= sizeof(int)*(1); = sizeof(curandState_t)*(

unsigned int mem_N

unsigned int mem_NNList NNList);

unsigned int mem_1

unsigned int mem_rand rand_block);

unsigned int mem_measured_quantity = sizeof(int)*(grid);

unsigned int mem_measured_magnet = sizeof(int)*(grid);

------------ Assignment of memory for host

--------------------------*/

int* int* int* int*

h_NNList h_flag

h_out_energy h_out_magnet

(int*) malloc(mem_NNList);

(int*) malloc(mem_1);

(int*) malloc(mem_measured_quantity);

(int*) malloc(mem_measured_magnet );

----------- Copy the Nearest Neighbor List to devie memory

-----------*/

for(int i=0; i< nla; i++) {

for(int j=0; j< h_N_NNList[i]; j++) { h_NNList[i+j *nla] =nn[i*6+j];

}

---------- Assignment of memory for device

------------------------*/

int* d_spin; int* d_spin_new; unsigned int* d_label; unsigned int* NNList;

ConeectList; N_NNList; N_ConeectList; curandState_t* d_random_data; int* d_flag; int* d_out_energy;

unsigned int* unsigned int* unsigned int*

/

*

/

*

/

*

140

145

150

155

160

165

170

int* d_out_magnet;

cudaMalloc( void**) &d_spin mem_ _N );

cudaMalloc( void**) &d_spin_new mem_ _N );

cudaMalloc( void**) &d_label mem_ _N );

cudaMalloc( void**) &NNList mem_ _NNList );

cudaMalloc( void**) &ConeectList mem_ _NNList );

cudaMalloc( void**) &N_NNList mem_ _N );

cudaMalloc( void**) &N_ConeectList mem_ _N );

cudaMalloc( void**) &d_random_data mem_ rand);

cudaMalloc( void**) &d_flag mem_ _1 );

cudaMalloc( void**) &d_out_energy mem_ _measured_ _quantity);

cudaMalloc( void**) &d_out_magnet mem_ _measured_ magnet );

// Transfer the Nearest Neighbor List cudaMemcpy(NNList, h_NNList, mem_NNList, cudaMemcpyHostToDevice);

cudaMemcpy(N_NNList, h_N_NNList, mem_N, cudaMemcpyHostToDevice);

Set of initial random numbers ----------*/

int iri1=12345;

set_rand_init<<<grid_rand,threads_rand>>>(d_random_data,iri1, nla);

cudaDeviceSynchronize();

/*------------ Temperature and Monte Carlo steps

----------------------*/

double t, beta; double boltz;

const double tstart = 3.8, tunit = 0.05;

const int tnum = 17; // const double tstart = 4.2075, tunit // const int tnum = 2 6; // const double tstart = 4.2105, tunit // const int tnum = 14;

0.0005; // nx=32 0.0005; // nx=4 8

185

190

195

200

205

// const double tstart = 4.212, tunit = 0.00025; // nx=64

// const int tnum = 17;

// const double tstart = 4.213, tunit = 0.00025; // nx=96

// const int tnum = 9;

const int mcsl = 1000, mcs2 = 20000; // const int mcs1 = 10000, mcs2 = 200000;

/*------------ Variables for measured quantities

----------------------*/

double H, M, f2order, fm[5];

double fm2, fm4, fe1, fe2, cv;

double fnla2 = (double)nla*(double)nla;

/*------------ Start of GPU calculation

-------------------------------*/

printf("#\n");