CC BY
29
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
УДК 004.932
РЕАЛИЗАЦИЯ ПАРАЛЛЕЛЬНЫХ АЛГОРИТМОВ ИССЛЕДОВАНИЯ ПСЕВДОСИММЕТРИИ АТОМНЫХ КЛАСТЕРОВ С ИСПОЛЬЗОВАНИЕМ ГРАФИЧЕСКИХ УСКОРИТЕЛЕЙ
© 2013 г. И.Н. Лозгачев, Н.В. Сомов
Нижегородский госуниверситет им. Н.И. Лобачевского ivan.lozgachev@itlab.unn.ru
Поступила в редакцию 14.09.2013
Разработаны параллельные алгоритмы исследования особенностей псевдосимметрии атомных кластеров. Создан вычислительный программный комплекс для исследования псевдосимметрических особенностей атомных кластеров и кристаллов.
Ключевые слова: кристаллическая структура, атомный кластер, псевдосимметрия, электронная плотность, инвариантность, высокопроизводительные вычисления, графический ускоритель.
Введение
Проблема установления взаимосвязи между химическим составом, атомной структурой и физическими свойствами кристаллов до сих пор является открытой. Квантово-механический подход к решению данной проблемы сопряжен с большими вычислительными затратами. В связи с этим большую пользу в исследовании проблемы «структура-свойство» несет выявление общих структурных характеристик, влияющих на физические свойства кристаллов. Одной из таких общих характеристик является симметрия.
Симметрия кристалла определяет набор физических свойств, которыми может обладать кристалл. В силу того, что кристаллическая решетка периодична в пространстве, характерными операциями для пространственной симметрии являются инверсия, повороты (на 180°, 90°, 60° и 120°), отражение в плоскости, инверсионные повороты и параллельный перенос, а также различные комбинации данных преобразований. Инвариантность атомной структуры кристалла относительно одной из 230 пространственных групп симметрии называется федоровской симметрией [1].
Однако для некоторых кристаллов характерно, что большая часть их атомной структуры инвариантна относительно надгруппы группы симметрии кристалла как целого. Такое свойство кристаллов называют псевдосимметрией, а кристаллы - псевдосимметричными [2, 3].
Псевдосимметрия влияет на многие физические свойства кристаллов: генерацию второй оп-
тической гармоники, пьезоэффект, пироэффект, сегнетоэлектрический эффект [4-6]. Псевдосим-метричное расположение атомов может приводить к закономерному погашению отдельных групп рефлексов на дифракционных картинах от псевдосимметричных кристаллов. Наличие псевдосимметрии в кристалле является признаком возможного фазового перехода второго рода и т.д. [7-10]. Вычисление количественных характеристик симметричности кристалла является затратным, т.к. элементарная ячейка кристалла может содержать несколько десятков тысяч атомов. Еще более затратной задачей является расчет данных характеристик для атомного кластера, состоящего из нескольких элементарных ячеек. Для эффективного решения такой задачи требуется применение параллельных алгоритмов. В данной работе предложены параллельные алгоритмы расчета степени инвариантности электронной плотности кристалла [11] и алгоритмы моделирования и исследования псевдосимметрических особенностей атомных кластеров кристаллов. Описаны результаты реализации данных алгоритмов для распределенных вычислительных сред с использованием технологий ОрепМР [12] и Ореп^ [13].
Вычисление степени инвариантности функции электронной плотности кристалла
Рассмотрим алгоритм расчета степени инвариантности электронной плотности относительно псевдосимметричного преобразования,
заданного матрицей обобщенного поворота и вектором параллельного переноса
Ч = {?, <*}, (1)
где ч - матрица обобщенного поворота, 1 -вектор трансляции [1]. В [11] был предложен функционал для количественной оценки степени инвариантности электронной плотности кристалла относительно заданного изометрического преобразования
Я < 2
ї [p(r )] =
Jp(r )p(qr ')dV Jp(r )2 dV
Т№Г
где
[p(r)] = K I G(Я )exp(- 2™(я , Г)),
Я
G(Я )= F (Я )F (- qT Я ),
K =
IIF(Я )
sin Qn
x
(8)
где 0max - максимальный угол дифракции, X -длина волны;
Я <-
1
d„
(9)
(2)
где р(г) - функция электронной плотности кристалла. Функционал (2) может принимать значения, лежащие в диапазоне (0,1]. Если исследуемый кристалл инвариантен относительно оператора Ч (Ч е G, где G - группа симметрии кристалла), то ^ [р(г )] = 1. Для кристаллов, асимметричных относительно преобразования Ч & G, Ъц[Р(г7)]-0.
Функцию электронной плотности можно разложить ряд Фурье. Тогда функционал (2) в обратном пространстве примет вид [10]
^ F (й У (- ЧТН )ехр(- 2то(й, 1))
КЮЬ —-------------^ I ----------------, (3)
где dmin - минимальное межплоскостное расстояние (разрешение).
Чем шире диапазон индексов суммирования h, к и l, тем с большей точностью рассчитывается электронная плотность кристалла. На практике было установлено, что оптимальным является диапазон индексов вектора обратной решетки
|Я| = 2.0 А" ', что соответствует dmin = 0.5 А и sin 0„
= 1.0 А-
где F(я) - структурная амплитуда, Й - вектор обратной решетки, ^ - транспонированная
матрица обобщенного поворота. Структурная амплитуда определяется следующим выражением
^й ЬУлИотехр^ • б)). (4)
і=1
где ^(Н) - атомный фактор, Оі - величина заселенности кристаллографической позиции, Т -
фактор Дебая-Валера [14]. Для удобства расчета выражение (3) перепишем в следующем виде
Па [Р(7)] = К У G(И )ехр(- 2%і(Й, 7)), (5)
(6) (7)
Диапазон суммирования по вектору обратной решетки й в (5) можно ограничить следующими величинами, принятыми в рентгеноструктурном анализе [15]
Вычисление функционала (5) выполняется в несколько этапов, алгоритм включает следующую последовательность действий:
1. Выполнение подготовительного этапа:
1.1. Размножение атомов структуры операторами группы симметрии кристалла;
1.2. Построение карты структурных амплитуд:
1.2.1. Определение размера карты структурных амплитуд (диапазонов суммирования h, к
и I);
1.2.2. Вычисление F(й) по формуле (4) для каждого узла карты.
2. Построение G-карты для скоростного расчета степени инвариантности электронной плотности относительно заданного преобразования :
2.1. Выполняется расчет карты G(й) по формуле (6);
2.2. По формуле (7) вычисляется нормировочный коэффициент К.
3. На завершающем этапе по G-карте вычисляется функционал (5).
Для распараллеливания каждого из этих этапов применяются различные подходы. Размножение атомов операторами симметрии не является трудоемкой задачей, поэтому выполняется последовательно в одном потоке. Алгоритм построения карты структурных амплитуд и G-карты реализован как для центрального процессора, так и для графического ускорителя. Разделение нагрузки между параллельными потоками происходит равномерно по узлам сетки так, что выполнение происходит независимо друг от друга: каждый поток вычисляет свое значение карты, поочередно выполняя приведенные действия. На последнем этапе вычислений запускается параллельный цикл, умножающий каждую ячейку G-карты на экспонентную составляю-
щую формулы (5), после чего с помощью редукции вычисляется значение степени инвариантности электронной плотности.
На центральном процессоре распараллеливание осуществляется с использованием стандартных директив OpenMP, реализующих параллельные циклы for и редукцию. Для выполнения на графическом ускорителе необходимо предварительно скопировать исходные данные в память устройства, что является трудоемкой операцией. Для достижения наибольшей эффективности обработки данных максимально используется локальное адресное пространство. За счет большого числа потоков на графических процессорах и отсутствия информационной зависимости между параллельными потоками предложенный алгоритм дает хорошее ускорение.
Полученная на первом этапе карта структурных амплитуд используется в других алгоритмах в качестве массива комплексных Фурье-коэффи-циентов разложения функции электронной плотности. Например, карту структурных амлитуд можно использовать как набор входных данных для решения задачи глобальной оптимизации функционала (2) [16, 17]. Однако для ее хранения требуется много памяти, и в случае больших объемов данных вычисление на графическом ускорителе становится невозможным. В данном случае вычисление функционала (5) производится прямым путем - итерации циклов разделяются между потоками, значения G(H) вычисляются непосредственно при суммировании.
Моделирование атомных кластеров
и исследование их псевдосимметрических особенностей
При исследовании симметрийных особенностей кристаллов разупорядоченных твердых растворов возникает задача построения атомного кластера кристалла, полученного в результате моделирования набора элементарных ячеек.
Атомный кластер кристалла представляет собой параллелепипед со сторонами aNx, bNy и cNz, где a, b и c - параметры элементарной ячейки исходного кристалла, Nx, Ny и Nz -число элементарных ячеек исходного кристалла, образующих атомный кластер в соответствующем направлении.
Рассмотрим алгоритм построения атомного кластера кристалла. Пусть нам известна атомная структура исходного кристалла. Имея информацию о независимой части атомной структуры и группе симметрии кристалла, несложно получить множество атомов S, заданное векто-
рами относительных координат г, характеризующими положения всех атомов в элементарной ячейке исходного кристалла.
Каждый атом исходной элементарной ячейки транслируется в соседние элементарные ячейки атомного кластера. Относительные координаты транслированного атома будут определяться следующим образом С (х,. + т)/Nx л
Г' = (У, + п)/NУ , (10)
(7,. + р)/ Nг
где X.
У,
V'
z. -
компоненты вектора r
т = 0, Nx -1, п = 0, Ny -1 и р = 0, Nz -1 - индексы элементарных ячеек исходного кристалла, образующих атомный кластер.
Наибольший интерес в настоящее время представляют атомные кластеры разупорядо-ченных кристаллов. Разупорядоченные кристаллы характеризуются наличием одной или нескольких смешанных кристаллографических позиций или позиций с дефицитной заселенностью - вакансиями. Разупорядоченность влияет на многие физические свойства кристаллов, в частности на оптические [18].
При моделировании кластера позиция Г' заселяется методом Монте-Карло с учетом вероятностей, заданных в модели исходного кристалла [10].
Для исследования псевдосимметрических особенностей атомных кластеров реализованы три типа численных экспериментов:
• исследование зависимости степени инвариантности электронной плотности атомного кластера от его линейных размеров;
• вычисление статистических характеристик (математического ожидания, дисперсии и т.п.) степени инвариантности электронной плотности для атомного кластера заданного размера;
• исследование зависимости степени инвариантности электронной плотности атомного кластера от его линейного размера и заселенности разупорядоченных позиций.
Параллельная реализация алгоритма генерации атомного кластера основана на распределении подъячеек между взаимодействующими потоками. Распараллеливание приведенных экспериментов организовано на более высоком уровне абстракции: каждый эксперимент дробится на подзадачи, отличающиеся размером кластера, количеством итераций или величиной заселенности разупорядоченной позиции. Время вычисления степени инвариантности электронной плотности атомного кластера растет пропорционально объему кластера, для равно-
и
мерного распределения нагрузки по вычислительным узлам необходима специальная логика. За единицу стоимости расчета возьмем время обработки атомного кластера размера 1^1x1, т.е. время обработки одной элементарной ячейки исходного кристалла. Тогда приближенно можно оценить стоимость расчета каждой подзадачи как произведение варьируемых параметров и всего эксперимента - как сумму стоимости всех задач следующим образом
Nx Ny Nz
C = О , I x^^ m x n x p , (11)
m=1 n=1 p=1
где C - стоимость эксперимента, О - число итераций по заселенности, I - число повторов. Так, например, стоимость расчета эксперимента с размерами кластера от 1x1x1 до 2x2x2, варьированием заселенности от 0.2 до 0.7 с шагом 0.1 и числом повторов 10 будет равна 27 x 6 x10 = 1620. Для наиболее эффективного распараллеливания нагрузки по узлам кластера остается лишь решить задачу равномерного распределения подзадач, учитывая их стоимость и стоимость эксперимента в целом. Данная задача является аналогом задачи об упаковке и принадлежит классу NP-полных алгоритмов. Для ее решения используется приближенный полиномиальный алгоритм First Fit Decreasing [19]. Реализация данного уровня параллелизма сделана с помощью технологии MPI, где каждый из процессов выполняет свою подзадачу на CPU или GPU вычислительной системы с разделенной памятью.
Вычислительный программный комплекс
На основе приведенных алгоритмов реализован программный комплекс, позволяющий исследовать псевдосимметрию атомных структур кристаллов и атомных кластеров. Информация о кристаллической структуре: параметры элементарной ячейки, координаты и заселенности атомов, тепловые параметры - извлекается из файла данных в формате CIF [20]. Данный формат принят как стандарт представления структурной информации Международным союзом кристаллографов [21].
Программный комплекс позволяет запускать вычислительные эксперименты как на системах с общей памятью, используя многопоточность центрального процессора и графического ускорителя, так и на системах с разделенной памятью. Использование технологии OpenCL расширяет список поддерживаемых графических ускорителей, поддерживая запуск на видеокартах NVidia, AMD и графических сопроцессорах Intel.
Результаты тестов производительности
Вычислительные эксперименты проводились c использованием следующей конфигурации:
• язык программирования C++;
• технологии распараллеливания: MPI,
OpenMP, OpenCL;
• операционная система Microsoft Windows Server 2008 HPC Edition SP2 x64;
• компилятор Microsoft Visual Studio 2010;
• 16 узлов по 2 CPU Intel Xeon L5630 2.13 ГГц, 4 ядра, 24 Гб оперативной памяти, NVidia Tesla M2050.
Полное время расчетов, выполняемых для атомного кластера, определяется как
ttotal = ^amp ^ ^calc , (12)
где t - время расчета узлов карты структурных амплитуд по формуле (4), величина tamp зависит от количества атомов в кластере и числа узлов карты структурных амплитуд; tcalc -время расчета G-карты по формуле (6) и функционала (5) на основе G-карты, величина tcalc зависит только от размеров G-карты.
Для выяснения истинного характера зависимости времени вычислений от различных параметров была проведена серия расчетов на модельных кристаллах. Модельные кристаллы для удобства расчетов имели кубическую элементарную ячейку, атомы генерировались случайным образом.
На рис. 1 представлены результаты измерений времени для параллельной CPU-версии и параллельной GPU-версии алгоритма расчета степени инвариантности электронной плотности в зависимости от количества атомов и объема ячейки.
Карта структурных амплитуд и G-карта всегда имеют одинаковое количество узлов. Пусть число узлов равно VM. Были установлены следующие зависимости времени расчетов от числа атомов и размера карт
tamp = кАШм , (13)
tcalc = 2kGVM , (14)
где N - число атомов, кА - линейный коэффициент, зависящий от производительности вычислительной системы и равный времени обработки одного атома для одного узла карты структурных амплитуд, kG - линейный коэффициент, зависящий от производительности вычислительной системы и равный времени обработки одного узла G-карты.
Количество атомов Объем ячейки, А3
♦ Время работы СРи ■ Время работы ОРИ ♦ Время работы СРи ■ Время работы ЄРи
Рис. 1. Время работы параллельных реализаций алгоритмов: а) зависимость времени расчетов от числа атомов в элементарной ячейке модельного кристалла (^яч= 45882 А3, эт0Д=1 А-1); б) зависимость времени расчетов от объема элементарной ячейки модельного кристалла (эт0/^=1 А- , число атомов в элементарной ячейке 4320)
Таблица
Результаты линейной аппроксимации зависимостей, приведенных на рис. 1
CPU GPU
Зависимость времени расчетов от числа атомов
а, с tcalc , с R а, с tcalc , с R
0.0325(1) 0.5(3) 0.9999 0.0021(1) 0.4(3) 0.9996
Зависимость времени расчетов от числа узлов G-карты
Р, с R Р, с R
0.0031(1) 0.9997 0.0002(1) 0.9999
Если число узлов карты постоянно, то зависимость времени расчетов от числа атомов аппроксимируется линейной функцией
ttotal = aN + ^calc , (15)
где a = kAVM = const - среднее время обработки
одного атома, t^ = const - время расчета G-карты и функционала (5).
При фиксированном числе атомов в элементарной ячейке кристалла зависимость времени расчетов от числа узлов в картах аппроксимируется следующей линейной функцией
ttotal =PVM , (16)
где р = kAN + 2kG = const - среднее время обра-
ботки одного узла обеих карт.
В таблице представлены результаты линейной аппроксимации методом наименьших квадратов, выполненные для зависимостей на рис. 1. Следует отметить, что величина VM, при фиксированном sin 0 / X , линейно связана с объемом элементарной ячейки кристалла соотношениями (8) или (9). Коэффициент корреляции в таблице обозначен R.
Рассмотрим результаты аппроксимации, полученные для зависимостей на рис. 1а. Время обработки одного атома а параллельной GPU-версией алгоритма в 15.5 раза меньше, чем время работы параллельной CPU-версии. Время
расчета /са1с с учетом погрешности отличается незначительно.
Время обработки одного узла обеих карт р для зависимостей, представленных на рис. 1б, также отличается в 15.5 раза для СРи- и ОРи-версий алгоритма.
Из данных, представленных в таблице, видно, что функции (15) и (16) хорошо аппроксимируют зависимости времени расчетов от числа атомов и количества узлов в картах, минимальный коэффициент корреляции составляет 0.9996. Для значений свободного члена ^с получены высокие погрешности (до 75%), это связано с тем, что суммарные погрешности измерения времени работы параллельного алгоритма превышают время вычисления и суммирования О-карты. Поэтому свободный член ^с невозможно достоверно измерить таким образом.
На рис. 2 показано ускорение алгоритма исследования псевдосимметрии атомных кластеров. В качестве объекта исследования был взят атомный кластер, построенный на основе элементарной ячейки кристалла К(Т1, 7г)0Р04 (ICSD-№ 66570) [22]. Эксперимент представлял собой построение атомных кластеров размеров от 1x1x1 до 3x3x3 с изменением заселенности атомов титана и циркония от 0 до 0.5 с шагом 0.1 и количеством повторов 5.
Число узлов вычислительного кластера
Рис. 2. Ускорение реализации алгоритма исследования псевдосимметрии атомных кластеров
Из графика (рис. 2) видно недостаточно эффективное ускорение (на 10 узлах кластера составляет приблизительно 5.5 раза), что может быть объяснено недостатками алгоритма First Fit Decreasing. В дальнейшем планируется заменить этот алгоритм на Best Fit [19].
Заключение
В работе описаны параллельный алгоритм расчета степени инвариантности электронной плотности кристалла и его реализация на вычислительных системах, использующих графические ускорители. Описан параллельный алгоритм построения атомных кластеров с разупоря-доченной структурой. Реализованы алгоритмы численных экспериментов, позволяющие про-водить статистические исследования зависимостей степени инвариантности электронной плотности от размера атомного кластера и его концентрационного состава.
Установлена зависимость времени работы параллельного CPU- и GPU-алгоритма от числа атомов и количества узлов карты структурных амплитуд и G-карты. Экспериментально доказано, что для расчетов эффективнее применять параллельную GPU-версию алгоритма.
Обнаружено падение ускорения алгоритма исследования псевдосимметрии атомных кластеров с ростом числа вычислительных узлов, вероятной причиной этого являются недостатки алгоритма First Fit Decreasing.
Работа выполнена при частичной финансовой поддержке программы «Научные и научно-педагогические кадры инновационной России на 2009-2013 годы», грант № 14.B37.21.1158.
Стсок лumepamуpы
1. Чупрунов Е.В., Хохлов А.Ф., Фаддеев М.А. Кристаллография: Учебник для вузов. М.: Изд-во физико-математической литературы, 2000. 496 c.
2. Сомов Н.В., Чупрунов Е.В. Трансляционная и инверсионная псевдосимметрия атомных структур кри-
сталлов органических и элементоорганических соединений // Кристаллография. 2009. № 4. С. 581-587.
3. Чупрунов Е.В. Федоровская псевдосимметрия кристаллов. Обзор // Кристаллография. 2007. №. 1. С. 5-16.
4. Ivanov V.A., Faddeev M.A., Chuprunov E.V. Pseudosymmetry and some characteristics of pyroelectric properties of crystals // Crystallography Reports. 2000. V. 45. № 5. P. 839.
5. Katkova M.R., Nosov S.S., Chuprunov E.V. Pseudosymmetry and ferroelectric phase transitions in the KTP structure type // Crystallography Reports. 2000. V. 45. № 4. P. 647.
6. Abrahams S.C., Kurtz S.K., Jamieson P.B. Atomic displacement relationship to Curie temperature and spontaneous polarization in displacive ferroelectric // Phys. Review. 1968. V. 172. № 2. P. 551-553.
7. Kroumova E., Aroyo M.I., Perez-Mato J.M. Prediction of new displacive ferroelectrics through systematic pseudosymmetry search. Results for materials with Pba2 and Pmc21 symmetry // Acta Crystallographica. 2002. V. 58. P. 921-933.
8. Igartua J.M., Aroyo M.I., Kroumova E. Search for Pnma materials with high-temperature structural phase transitions //Acta Crystallographica. 1999. V. 55. P. 177-185.
9. Иванов В.А. Дис... канд. физ.-мат. наук. Нижний Новгород, 2008.
10. Сомов Н.В. Дис. канд. физ.-мат. наук. Нижний Новгород, 2011.
11. Чупрунов Е.В., Солдатов Е.А., Тархова Т.Н. О количественных оценках симметричности кристаллических структур // Кристаллография. 1988. № 3. С. 759.
12. Официальный сайт спецификации OpenMP [Электронный ресурс]. Режим доступа: http://openmp.org/wp/.
13. Официальный сайт спецификации OpenCL [Электронный ресурс]. Режим доступа: http://www. khronos .org/opencl/.
14. Уманский Я.С. и др. Кристаллография, рентгенография и электронная микроскопия. М.: Металлургия, 1982. 632 с.
15. Порай-Кошиц М.А. Основы структурного анализа химических соединений: Учебное пособие. М.: Высшая школа, 1989. 192 с.
16. Gergel V.P., Strongin R.G. Parallel computing for globally optimal decision making on cluster systems // Future Generation Computer Systems. 2005. V. 21. № 5. P. 673-678.
17. Стронгин Р.Г., Гергель В.П., Баркалов К.А. Параллельные методы решения задач глобальной оптимизации // Известия высших учебных заведений. Приборостроение. 2009. Т. 52. № 10. С. 25-33.
18. Малов А.В., Марычев М.О., Рябочкина П.А. и др. Спектроскопические и структурные свойства кристаллов кальций-ниобий-галлиевого граната, активированных ионами Er3+ // Вестник Нижегородского университета им. Н.И. Лобачевского. 2008. № 6. С. 46-52.
19. Feitelson D., Rudolph L. Job scheduling strategies for parallel processing // Lecture Notes in Computer Science. 2005. V. 3834. 283 p.
20. Официальный сайт описания формата CIF [Электронный ресурс]. Режим доступа: http://it.iucr.org/G/.
21. Официальный сайт Международного союза 22. Hansen N.K., Protas J. Marnier G. The electron-
кристаллографов [Электронный ресурс]. Режим дос- density distribution in KTiOPO4 // Acta Crystallographi-
тупа: http://www.iucr.org/. ca. 1991. V. 47. P. 660-672.
IMPLEMENTATION OF PARALLEL ALGORITHMS FOR INVESTIGATION OF PSEUDOSYMMETRY OF ATOM CLUSTERS USING GRAPHICS ACCELERATORS
I.N. Lozgachev, N. V. Somov
Parallel algorithms for investigation of pseudosymmetry of atom clusters are developed. A software package to investigate pseudosymmetric features of atom clusters and crystals is created.
Keywords: crystal structure, atom cluster, pseudosymmetry, electron density, invariance, high performance computing, graphics accelerator.
