CC BY
10
Использование графических процессоров для моделирования дифракционных характеристик наноразмерных структур
В.С. Неверов, НИЦ «Курчатовский институт», аспирант,
vs-never@hotmail. com
Доклад посвящен общим принципам и различным аспектам реализации вычислений интенсивности излучения, рассеянного наноразмерными структурами, на графических процессорах с использованием технологии CUDA. Показано, что перенос расчетов на графические процессоры ускоряет вычисления в десятки раз по сравнению с расчетами на центральных процессорах.
Физическая задача
Рентгеновская, нейтронная и электронная дифрактометрии являются наиболее эффективными инструментами исследования атомной структуры вещества. Вышеуказанные диагностики заключаются в рассеянии пучка фотонов, нейтронов или электронов на исследуемом образце и последующей обработке получаемой дифракционной картины (зависимости интенсивности рассеянного излучения от углов рассеяния). В случае кристаллических образцов (монокристаллов или поликристаллов с размерами областей кристалличности > 0.1 мкм) дифрактограммы имеют резкие пики, именуемые брегговскими. В этом случае восстановление структуры атомной ячейки кристалла возможно с применением стандартных методов кристаллографии (в том числе методами Фурье-анализа). Однако анализ результатов рассеяния на квазиаморфных материалах, для которых пространственное упорядочение структуры сохраняется лишь на нескольких нанометрах или нескольких десятках нанометров, имеет ряд сложностей, связанных с существенным отклонением от поликристаллической структуризации. Широкие пики на получаемых в таких случаях дифрактограммах нельзя трактовать как брегговские, поэтому дифракционный профиль содержит недостаточно информации для восстановления структуры сложного материала. Стандартные методы кристаллографии не дают исчерпывающей информации об объекте. Проблемы, возникающие при анализе квазиаморфных материалов, описаны подробно в [1]. В то же время информация об исследуемом веществе от других диагностик, а также информация об условиях получения вещества, помогут
существенно сузить класс возможных решении - наноразмерных структур, отвечающих за широкие дифракционные пики. В этом случае наноструктуру образца в дипазоне размеров порядка 1-10 нм возможно восстановить доступными вычислительными средствами путем сопоставления экспериментальной дифрактограммы с референсными. Получить референсные дифрактограммы на наноразмерных структурах можно и экспериментально, однако значительно проще рассчитать их численно.
Хотя рентгеновской, нейтронной и электронной дифракции отвечают разные физические процессы - в первом случае рассеяние фотонов происходит на электронной оболочке атома, во втором, рассеяние нейтронов происходит на нуклонах ядра, а в третьем, электроны взаимодействуют и с электронной оболочкой атома, и с протонами в ядре - все три типа дифракции в широком диапазоне параметров математически описываются одинаково.
Амплитуда рассеянного излучения на ансамбле из N атомов с координатами г имеет вид:
А(ч (я Уч\ (1)
где ч - вектор рассеяния, а
Г, (я) - атомный форм фактор, в котором учтен характер рассеяния излучения (рентгеновского, нейтронного или электронного) на отдельном атоме. Интенсивность вычисляется, как квадрат амплитуды:
I (ч )= А (ч )А(ч ). (2)
В случае рассеяния на порошке структура образца изотропна, усреднение интенсивности по пространственным углам приводит к формуле Дебая:
-Д Д вт( агн)
1 ( я ) = Ц ( я ) (я , (3)
г= 0 ] = 0 Ч'у
где г]к=|гргк| - расстояние между атомами ) и к. Здесь и далее будем считать, что форм-факторы атомов действительны (например, для широко используемой дифракции холодных нейтронов это не выполняется только для нескольких химических элементов).
Хотя при исследовании квазиаморфных материалов наибольший интерес представляет именно порошковая дифракция, созданный численный код производит вычисления и двумерных дифракционных картин, получающихся при рассеянии на монокристаллах.
Решение на графических процессорах
Автору известно, по крайней мере, о двух кодах, моделирующих интенсивность рассеянного излучения на наноструктурах с использованием графических процессоров (далее GPU) [2,3]. Ни один из них пока не доступен для использования сторонними исследователями. Кроме этих кодов существует код, моделирующий малоугловое рентгеновское рассеяние, доступный в исходном коде (https://simtk.org/home/simsaxs). Однако задача малоуглового рассеяния отлична от рассматриваемой. В то же время автором ранее был создан численный код, вычисляющий интенсивности рассеянного излучения на произвольно упорядоченных ансамблях наноструктур с использованием центральных процессоров (далее CPU) и распараллеленный при помощи MPI и OpenMP [4]. Этот код применялся в задаче восстановления наноструктурного состава углеродистых пленок, извлеченных из вакуумной камеры токамака Т-10 [5,6,7]. Поскольку задача не требует двойной точности вычислений (подробно анализ точности см. в [2]) и является хорошо распараллеливаемой, было решено перенести вычисления на GPU. Для реализации использовалась архитектура CUDA (далее будем придерживаться терминологии, используемой в [8]).
Задача вычисления интенсивности рассеянного излучения разбивается на две подзадачи:
• вычисление абсолютных координат атомов всего ансамбля по их относительным координатам (координатам в ячейке, отдельной структуре и ее пространственных копиях),
• вычисление интенсивности рассеянного излучения на ансамбле атомов по формуле Дебая (порошковая дифракция, одномерный профиль), или по общей формуле (дифракция на монокристалле, двумерная картина).
Вычисление абсолютных координат атомов
В качестве начальных данных для кода задаются координаты атомов в ячейке кристалла или структуре, симметрично эквивалентные позиции и номера ячейки вдоль векторов трансляции (для кристаллов), а также положения структуры и ее копий во всем рассматриваемом ансамбле наноструктур. Перед расчетом интенсивности рассеянного излучения удобно рассчитать абсолютные координаты атомов в ансамбле. Поскольку положение каждого атома рассчитывается независимо от остальных, эти вычисления можно произвести на GPU. Вычисления представляют собой линейные операции по преобразованию координат и ничего более, поэтому приводить их здесь мы не будем. Помимо координат каждый атом характеризуется целочисленным идентификатором, характеризующим тип атома. Таким образом, после
расчета абсолютных координат в глобальной памяти GPU размещается массив 16-байтных структур (float3 r, int id), полностью описывающих атомы в ансамбле. Отсюда следует ограничение на число атомов в ансамбле. Для ансамблей из 107 атомов массив будет занимать ~ 153 МБ, при доступных на настоящий момент 1-2 ГБ памяти на одну видеокарту. Это и делает возможным уже отмеченное ранее восстановление наноструктурного состава широкого класса аморфных сред в диапазоне размеров порядка 1-10 нм.
Вычисление интенсивности рассеянного излучения на ансамбле атомов по формуле Дебая
Перед выполнением расчета интенсивности в глобальную память GPU должны быть загружены все форм-факторы атомов и массив значений модуля вектора рассеяния q.
Формула Дебая (3) для интенсивности рассеянного излучения симметрична относительно замены i на j. Поэтому полное число значений интенсивности, которые нужно рассчитать, равно mN(N+1)/2, где m - число точек сетки модуля вектора рассеяния, а N - полное число атомов в ансамбле. С одной стороны, нужно добиться максимальной параллельности в вычислении этих значений, а с другой - избежать конфликтов записи этих значений в глобальную память. Очевидно, что разместить в глобальной памяти массив значений интенсивности размером m-N(N+1)/2, что позволило бы избежать конфликтов записи, при больших N - невозможно. Достичь компромисса можно, разбив атомные данные между Ng сетками, состоящими из Nb х Nb блоков. На рисунке 1 приведена схема распределения атомов между блоками тредов, принадлежащими разным сеткам. Частичные суммы из формулы (3) будут вычисляться на сетках последовательно с синхронизацией: dim3 grid(Nb,Nb); for (i=0;i<Ng;i++){
//вызов функции-ядра, вычисляющей интенсивность //рассеяния для атомов, попавших в i-ую сетку DebyeKernel<<<grid,Nt>>> (...) ; cudaThreadSynchronize() ;
}
Синхронизация позволит сократить размер массива значений интенсивности в глобальной памяти до Nb х Nb х m. Писать значения в этот массив внутри функции-ядра можно без синхронизации. Задача двойного суммирования массива интенсивности по Nb будет возложена на CPU после выполнения всех вычислений. Число сеток Ng должно быть таким, чтобы с одной стороны массив значений интенсивности Nb х Nb х m помещался в глобальной памяти, а с другой - время суммирования этого массива на CPU было много меньше времени его
вычисления на GPU. Кроме этого, существует дополнительный фактор - ограничение времени исполнения функции-ядра. Если такое ограничение включено, функция-ядро должна успевать выполняться до того, как будет прервана. Таким образом, число сеток Ng, или их размер Nb х Nb - варьируемые параметры, зависящие от конкретных GPU.
атомы i
Рис. 1. Распределение атомов между блоками тредов, принадлежащих разным сеткам. Штрихами закрашены участки, которые не нужно рассчитывать в силу симметричности формулы (3) относительно замены i-j. Фиолетовая область - сетка из Nb х Nb блоков. Голубая область - блок из Nt тредов
Каждый блок, состоящий из Nt независимых тредов, отвечает за вычисление парциального вклада в интенсивность от Na i-ых и j-ых
атомов. Если Na = л/N t , каждое значение rij = lri ~ r j| будет рассчитано отдельным тредом. Разместив полученные значения rij- в общей памяти на чипе и синхронизовав треды, можно приступить к расчету интенсивности. Каждый тред при этом рассчитает со всеми rij-только свои значения интенсивности (ср. с алгоритмом в [2], где каждый тред считает все значения интенсивности со своим значением rj). Добиться этого можно, когда m кратно Nt, а поскольку m -регулируемый параметр, можно считать, что это выполняется всегда. Число тредов в блоке - тоже регулируемый параметр, но для современных GPU оптимальным значением для данной задачи является Nt =256 (Na = 16). Ниже приведен фрагмент функции-ядра, выполняющий расчет интенсивности по формуле Дебая.
//в общей памяти размешаются массив rij
shared float rij[Na][Na]; //и массив идентификаторов i-ых и j-ых атомов _shared_ short id[2*Na];
//копирование идентификаторов из глобальной памяти //в общую
if (threadIdx.x<iEnd) {
id[threadIdx.x]=at[iBegin+threadIdx.x].id;
}
if (threadIdx.x<jEnd){
id[Na+threadIdx.x]=at[j Begin+threadIdx.x] .id;
}
//вычисление значений rij и запись их в общую память, //каждый тред считает только свое значение short i=threadIdx.x/Na; short j=threadIdx.x%Na; if ( (i<iEnd)&&(j<jEnd)) {
rij [i] [j]=length(at[iBegin+i] .r-at[jBegin+j] .r);
}
//обязательная синхронизация перед расчетом //интенсивности, чтобы удостовериться, //что все rij заполнены
syncthreads(); //Nq=m/Nt - количество точек по q, //которые должен рассчитать один тред for (short iq=0;iq<Nq;iq++){ float lI=0, lfj[Na];
//копирование значения q в локальную память треда float lq=q[iq*Nb+threadIdx.x];
//копирование значений форм-факторов j-ых атомов //при данном q в локальную память for (j=0;j<j End;j ++) {
lfj[j]=f[id[Na+j]][iq*Nb+threadIdx.x];
}
for (i=0;i<iEnd;i++){
float lfi=f[id[i]][iq*Nb+threadIdx.x]; for (j=0;j<jEnd;j++){
float qrij=lq*rij[i][j];
//вычисление интенсивности для данного rij //и добавление к знач. в локальной памяти треда if (qrij <0.0001f) lI+=lfi*lfj[j]; else lI+=(lfi*lfj[j]*_sinf(qrij))/qrij;
}
}
//запись значений интенсивности в глобальную память I[gridDim.y*Nq*Nb*blockIdx.x+Nq*Nb*blockIdx.y+iq*Nb+ threadIdx.x]+=lI*c;
}
Для интенсивности нейтронного рассеяния функция-ядро будет немного отличаться от модуля вектора рассеяния q в силу независимости форм-факторов, но принцип вычислений сохранится. Слабым местом вышеприведенной функции-ядра является копирование форм-факторов j-ых атомов в локальную память. Это приведет к нехватке регистров и, как следствие, к неполной занятости GPU при расчетах. Копирования форм-факторов в локальную память можно избежать, если, как в [2], предварительно отсортировать атомы по типам. В этом случае все j-ые атомы внутри одного блока будут иметь один и тот же тип, и вместо Na значений форм-фактора в локальную память надо будет скопировать только одно значение. Однако в нашем коде на данный момент такая предварительная сортировка не реализована.
Вычисление двумерной картины рассеянного
излучения на ансамбле атомов
Двумерная картина рассчитывается в полярных координатах q,q по формулам (1)-(2).
Расстояний между атомами вычислять не нужно, поэтому этот случай проще, чем предыдущий. Сетка (q,q>) разбивается на равные квадраты. Каждый блок считает только свой квадрат, а каждый тред -только одну точку. Если время исполнения функции-ядра лимитировано драйвером, ее можно вызывать в цикле, разбив массив атомов между вызовами.
dim3 grid(m/Nt,n/Nt); dim3 block(Nt,Nt); int Nlaunch=RoundUp(N/NpLaunch); for (int i=0;i<Nlaunch;i++) { int iBegin=i*NpLaunch; int iEnd=min(iBegin+NpLaunch,N); ExactKernel<<<grid, block>>> (...) ; cudaThreadSynchronize();
}
Синхронизация между вызовами ядра нужна, потому что каждое ядро записывает свою часть суммы по атомам в одни и те же массивы действительной и мнимой частей амплитуды. Ниже приведен фрагмент функции-ядра, выполняющий расчет двумерной дифракционной картины.
//каждый тред вычисляет свои точки на ceTKe(q,fi)
int iq=Nb*blockIdx.x+threadIdx.x;
int ifi=Nb*blockIdx.y+threadIdx.y;
//значение модуля вектора рассеяния q копируется
//из глобальной памяти в локальную
float lq=q[iq];
//создаются локальные значения действительной и
//мнимой частей амплитуды рассеяния в точке
float lA_real=0, lA_image=0;
//по значениям модуля вектора рассеяния q
//и полярного угла fi производится расчет вектора
//рассеяния q_v
float sintheta=lq/q[m-1];
float costheta=sqrt(1-sintheta*sintheta); float cosfi=0, sinfi=0;
_sincosf(ifi*2.f*PI/n,&sinfi, &cosfi) ;
float3 q_v=make_float3(costheta*cosfi,costheta* sinfi,sintheta)*qi;
for (int i=iBegin;i<iEnd;i++){
//копирование структуры, описывающей i-ый атом //в локальную память atomdata lat=at[i]; float qr=dot(q_v,lat.r); _sincosf(qr,&sinfi,&cosfi);
//копирование значения форм-фактора в локальную память float lf=f[lat.id][iq]; lA_real+=lf*cosfi ; lA_image+=lf*sinfi;
}
//запись в глобальную память значений действительной //и мнимой частей амплитуды рассеяния A_real[iq*n+ifi]+=lA_real; A_image[iq*n+ifi]+=lA_image;
Интенсивность по известным действительной и мнимой частям амплитуды вычисляется также на GPU.
Результаты
Тестирование производительности проводилось на GPU - Nvidia GeForce GTX285 и CPU - Intel Core i7-970 3.2 ГГц (6 ядер). На CPU запускалась MPI версия кода, использующая все 6 ядер и вычисляющая с двойной точностью. График зависимости ускорений от числа атомов в ансамбле при использовании GPU приведен на рисунке 2. Для 64 103 атомов расчет интенсивности рассеянного излучения по формуле Дебая на GPU выполнился в 55 раз быстрее при 75% занятости GPU. Ограничением занятости явилась нехватка регистров. Время вычисления на GPU составило 2 мин. 33 с. Расчет двумерной дифракционной картины для того же числа атомов занял на GPU в 39 раз меньше времени, чем на CPU при занятости 100%. Время расчета двумерной дифракционной картины на GPU составило 18 с.
Компиляторы CPU: MS Visual C+ GPU: NVCC 4.0 ь 2010 у
-♦-ID проф -■-2D карп шь (формула Дебая на
GPU: Nvidia GeFo CPU: Intel Core i7 rce GTX285 970, б-core, 3.2GHz
1.0Е+03 8.0E+03 6.4E+04
Число атомов в структуре
Рис. 2. Зависимость ускорения вычислений интенсивности рассеянного излучения при использовании GPU от числа атомов в структуре
Выводы
Задача моделирования дифракционных характеристик наноразмерных структур является идеальной для реализации на GPU. Скорость вычислений за счет использования GPU увеличивается в десятки раз даже в сравнении с новейшими шестиядерными процессорами. Созданный код может использоваться в автоматизированных системах обработки экспериментальных данных рентгеновской, электронной и нейтронной дифракций, в том числе построенных по принципу настольного грида.
Благодарности
Автор благодарен А.Б. Кукушкину и В.А. Вознесенскому за полезные обсуждения, а также всем соавторам статей [4-7] - за стимулирующее сотрудничество.
Литература
1. Billinge S.J.L. The nanostructure problem. //Physics 3, 25 (2010), http ://physics. aps. org/ articles/v3/25
2. L0nning Reiten A. X-ray scattering simulations using GPU-enabled algorithms. //Norwegian University of Science and Technology, preprint http://www.nanowiki.no/images/0/0ffGpu_scattering_v1.pdf
3. Gelisio L., Azanza Ricardo C. L., Leoni M., Scardi P. Real-space calculation of powder diffraction patterns on graphics processing units. // J. Appl. Cryst. (2010). 43, pp. 647-653.
4. Neverov V.S., Voloshinov V.V., Afanasiev A.P., Kukushkin A.B., Marusov N.L., Semenov I.B., Stankevich V.G., Svechnikov N.Yu., Tarasov A.S.,
Veligzhanin А.А., Zubavichus Ya.V. Distributed computing and optimization algorithms for interpretation of x-ray scattering by carbon nanostructures in the deposited films from tokamak T-10 // Proceedings of the 4th Intern. Conf. "Distributed Computing and Grid-Technologies in Science and Education" (Dubna, June 28-July 3, 2010). - Dubna: JINR, Д-11-2010-140, 2010.-p.452.
5. Kukushkin A.B., Neverov V.S., Marusov N.L., Semenov I.B., Kolbasov B.N., Voloshinov V.V., Afanasiev A.P., Tarasov A.S., Stankevich V.G., Svechnikov N.Yu., Veligzhanin A.A., Zubavichus Ya.V., Chernozatonskii L.A., Few-nanometer-wide carbon toroids in the hydrocarbon films deposited in tokamak T-10.// Chem. Phys. Lett., 506, 2011, pp. 265-268.
6. Неверов В.С., Кукушкин А.Б., Марусов Н.Л., Семёнов И.Б., Волошинов В.В., Афанасьев А.П., Тарасов А.С., Велигжанин А.А., Зубавичус Я.В., Свечников Н.Ю., Станкевич В.Г., Моделирование рентгеновской дифракции на углеродных наноструктурах и определение их возможного топологического состава в осажденных пленках из токамака Т-10.// Вопросы атомной науки и техники, Сер. Термоядерный синтез, 2010, том. 1, стр. 7-21.
7. Неверов В.С., Кукушкин А.Б., Марусов Н.Л., Семёнов И.Б., Волошинов В.В., Афанасьев А.П., Тарасов А.С., Велигжанин А.А., Зубавичус Я.В., Свечников Н.Ю., Станкевич В.Г., Численное моделирование эффектов интерференции рентгеновского рассеяния углеродными наноструктурами в осаждённых плёнках из токамака Т-10.// Вопросы атомной науки и техники, Сер. Термоядерный синтез, 2011, том. 1, стр. 1323.
8. NVIDIA Corporation, CUDA C Programming Guide, Version 4.0, 5/6/2011, http://developer.download.nvidia.com/compute/cuda/4_0/toolkit/docs/CUDA_C_ Programming_Guide.pdf
