Параллельное решение задачи Хартри—Фока для молекулы графена: масштабируемость и гиперэффективность Текст научной статьи по специальности «Математика»

УДК 004.942

А. И. Свитенков, Д. М. Спельников, В. Г. Маслов, А. В. Бухановский

ПАРАЛЛЕЛЬНОЕ РЕШЕНИЕ ЗАДАЧИ ХАРТРИ—ФОКА

ДЛЯ МОЛЕКУЛЫ ГРАФЕНА: МАСШТАБИРУЕМОСТЬ И ГИПЕРЭФФЕКТИВНОСТЬ

Рассмотрен параллельный метод решения уравнения Хартри—Фока, основанный на алгоритме DC. Предложенный метод позволяет частично решать задачу самосогласования локально на вычислительных узлах, что служит сокращению накладных расходов.

Ключевые слова: задача самосогласования, графен, параллельный алгоритм, уравнение Хартри—Фока, накладные расходы, квантовая химия, масштабируемость.

Введение. Распространенным подходом к приближенному решению уравнения Шре-дингера является так называемое одноэлектронное приближение, при котором каждый электрон рассматривается как независимый, движущийся в усредненном поле ядер атомов и других электронов. Этот подход используется в методах Хартри—Фока, функционала плотности или приближении сильносвязанных электронов. Уточнение решения за счет учета корреляционной энергии электронов приводит к критическому росту вычислительной сложности и на практике возможно при моделировании молекулярных систем лишь из небольшого числа атомов [1].

Для моделирования свойств по возможности больших систем используются наиболее простые одноэлектронные ab initio или даже полуэмпирические приближения [2]. Однако традиционная постановка задачи и в этих случаях приводит к кубическому росту вычислительной сложности при увеличении числа атомов, что делает ее неприменимой для систем, содержащих 1000 атомов и более. В этой связи в настоящее время широко обсуждаются линейно масштабируемые методы, значительно расширяющие границы применения квантовой химии [3]. Однако увеличение размеров систем до десятков и сотен тысяч атомов дополнительно требует параллельной реализации указанных алгоритмов. В работе представлен опыт применения линейно масштабируемого метода "Divide-and-conquer" (DC) для квантово-химического уравнения Хартри—Фока. Моделированию подлежала электронная структура молекулярных соединений типа „графен" и „графан".

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

Параллельный алгоритм, основанный на пространственной декомпозиции молекулярной системы, аналогичной декомпозиции в алгоритме DC, приводит не только к снижению времени выполнения одной итерации, но и к уменьшению числа итераций самосогласования, требуемых для сходимости. Относительно полного времени решения задачи Хартри—Фока, таким образом, наблюдается гиперускорение. Для объяснения описанного эффекта необходимо рассмотреть свойства уравнения Хартри—Фока и его решений.

Постановка задачи. Уравнение Хартри—Фока. В методе Хартри—Фока гамильтониан молекулярной системы в уравнении Шредингера заменяется приближенным одноэлектронным аналогом — фокианом. В таком случае приближенное уравнение Шредингера приобретает вид [5]:

1 2

" ^ У к Vк О) + VVк =гкVк (г).

(1)

Оно должно быть решено относительно волновых функции ^к(г) в действительной области пространства в некоторой окрестности неподвижных центров атомных ядер. Здесь V — оператор эффективного потенциала, в котором движется электрон: V=V0+Vd+Vx; У0 описывает вклад взаимодействия электрона и атомных ядер; Vd и Vx описывают взаимодействие с остальными электронами системы:

(

^ V' (Г) = £ \

У (Г ')|2

Л

6IV'

г - г'

V /(г),

Т/ () (Г ')Ук (г О _, КУк (г) =Х .1 ■-;-^

' *к ^

г - Г'

V' (г).

(2)

При разложении по выбранному набору базисных функций к-я орбиталь представляет собой линейную комбинацию базисных функций фг с неизвестными коэффициентами Си [6]:

м

Vк = Е Скг Фг ,

ц=1

относительно которых уравнение Хартри—Фока может быть записано следующим образом:

N

Е=гкСкг , (3)

'=1

где Г — матрица фокиана в соответствующем базисном разложении. Выражение (3) представляет собой задачу на собственные числа и собственные векторы матрицы Г. Матрица плотности Р определяется выражением:

р у = Е СкгСк' .

к=косс

(4)

Как видно, матрица фокиана, в свою очередь, зависит от матрицы плотности. Таким образом, уравнения (2)—(4) формируют так называемую самосогласованную задачу, решение которой сводится к итерационному процессу с последовательным уточнением матрицы Р и соответствующей матрицы Г до достижения сходимости.

Из уравнения (3) видна кубическая сложность предлагаемого алгоритма относительно размера матриц Р и Г, однако при использовании базиса сильно локализованных функций внедиагональные элементы матрицы плотности довольно быстро затухают с расстоянием. Это приводит к тому, что реальное число отличных от нуля матричных элементов оказывается ~N, а не ~N [7]. То же относится и к матрице фокиана в этом представлении. Соответственно все действия с такими матрицами имеют трудоемкость ниже О^ ), что отражает локальный характер квантовой механики и так или иначе используется всеми линейно масштабируемыми алгоритмами решения задачи Хартри—Фока.

Алгоритм БС и его параллельная реализация. С помощью алгоритма БС общая матрица плотности строится на основе решения, полученного не для всей системы в целом, а для некоторых перекрывающихся фрагментов. Величина буферной зоны задается посредством величины отсечения интегралов перекрывания базисных функций, значения меньше которой считаются нулевыми. Атомы, базисные функции которых не перекрываются, считаются не взаимодействующими непосредственно [8].

Для каждой подобласти проводится отдельный расчет субматрицы плотности. Объединяющим условием для всей рассчитываемой системы является только энергия уровня Ферми. У полученной субматрицы плотности для подобласти исключаются из рассмотрения все „углы".

На рис. 1 приведена матрица плотности для центральной области. Из полной матрицы плотности (слева) в дальнейших вычислениях используется только „крестообразная" часть (справа) [6].

Рис. 1

Неперекрывающиеся области (центральные) суммируются с единичным весом, а перекрывающиеся — с весом 0,5.

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

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

На рис. 2, а представлена блок-схема параллельной версии БС алгоритма. Под блоком вычисления субматриц плотности на узлах понимается последовательность действий, приведенная на рис. 2, б. Если итерации внешнего цикла не выполняются, то блок-схемы описывают параллельный вариант исходного алгоритма БС, в противном случае — модифицированный вариант, а на рис. 2, б приведен итерационный процесс, выполняемый локально.

Рис. 2

Теоретически оценить параллельное ускорение модифицированного и исходного вариантов алгоритма БС позволяет выражение:

5 =_2 ^ + ^_-. (5)

<N+dfilC. IN+CV^Nn

Ш Пл xK у 2 где N — число атомов в системе, а и Р — коэффициенты трудоемкости, П — число узлов, п — число вычислительных ядер на узле. Операция диагонализации матрицы характеризуется линейностью алгоритма DC, в то время как обновление матрицы гамильтониана имеет квадратичную сложность с очень малым коэффициентом. Поскольку для систем размером ~105 атомов нелинейность проявляется, в модель введено квадратичное слагаемое. Третье слагаемое в знаменателе относится к накладным расходам: d — длина сообщения, приходящаяся на один атом системы, с — количество соседей для каждого атома, к — геометрический фактор, связывающий площадь фрагмента с длиной его границы, т — скорость MPI передачи данных; K — среднее число локальных итераций, приходящихся на одну глобальную.

Идея модификации алгоритма DC состоит именно в том, что локальные итерации позволяют существенно повысить скорость сходимости глобального процесса самосогласования, тогда K > 1, и накладные расходы снижаются в соответствующее количество раз. Случай K = 1 соответствует параллельному ускорению исходного алгоритма DC.

Измерение производительности. Параллельная производительность модифицированного алгоритма DC измерялась в ходе моделирования электронной плотности молекул графена различного размера. Запуски производились на суперкомпьютере МГУ „Ломоносов" на 128, 256 и 512 узлах (8 вычислительных ядер на узел). В таблице приведены результаты измерения полного времени решения задачи, времени вычисления одной глобальной итерации и индексы эффективности использования ресурсов. Отсутствие данных по ускорению объясняется невозможностью проведения последовательного расчета для систем таких размеров.

Результаты измерений времени решения задач Хартри—Фока на 128/256/512 узлах

Размер молекулы Полное время выполнения, с Эффективность Время выполнения одной итерации, с Эффективность

20802 2118/1469/1048 0,71/0,71 3,0/1,9/1,3 0,78/0,73

46202 4892/2839/1847 0,86/0,76 7,9/4,1/2,7 0,96/0,75

98562 32346/12835/7513 1,26/0,85 16,5/8,9/4,9 0,93/0,9

Из таблицы видно, что для системы максимального размера при переходе от 128 к 256 ядрам наблюдается гиперэффективность использования вычислительных ресурсов ~1,26. Эффект наблюдается для полного времени решения задачи, максимально достигаемая эффективность для одной итерации меньше единицы. Механизм возникновения обнаруженного эффекта в общих чертах следующий. Скорость сходимости процесса самосогласования для фрагментов молекулярной системы заметно зависит от числа атомов, поэтому при уменьшении размера фрагментов в определенных условиях выигрыш по времени от ускорения сходимости может превысить возрастающие затраты на пересылку данных. Однако такое качественное описание не позволяет ответить на важный вопрос об эффективности использования именно модифицированного алгоритма БС. Сходимость процесса самосогласования для молекулы целиком может оказаться лучшей, чем для фрагментов, это приведет к проигрышу относительно исходного алгоритма. Если рассматривать модель (5) как ускорение относительно одной глобальной итерации, полное ускорение запишется как 50=5п0(п/по)1, где п0 — число итераций до сходимости в немо-дифицированном алгоритме БС, П1 — число локальных итераций до сходимости,

по — число глобальных итераций. Коэффициент при £ в правой части этого выражения будем называть коэффициентом гиперэффективности. В выражении (5) К = п/.

Модель сходимости процесса Хартри—Фока. Для теоретической оценки ускорения предлагаемой параллельной модификации алгоритма БС необходимо рассмотреть зависимость числа необходимых итераций от размера молекулы. Этот непростой процесс приближенно можно описать из следующих общих соображений. Известно, что масштабируемость того или иного алгоритма определяется связанностью данных, характерной для решаемой задачи. Рассмотрим реакцию на точечное возмущение решения уравнения (1). В качестве такого практически точечного возмущения зафиксируем вариацию на сферической поверхности а некоторого малого радиуса. Решение будем искать во внешней области. Выражение (1) может быть представлено как уравнение Пуассона, если понимать его решение в виде итерационного процесса, в котором источник определяется видом ^-функции, найденным на предыдущей итерации, тогда:

Vк(Ч) = -\^nG(q, т)5^+ 2{от)(вк - 0.т, (6)

где т) — функция Грина уравнения Пуассона. Это уравнение может быть решено итерационно, в результате будет получен ряд п-кратных интегральных слагаемых и остаточный член, порождаемый первым и вторым слагаемыми в уравнении (6) соответственно. Сложность вычисления каждого из слагаемых пропорциональна размеру области интегрирования в степени кратности интеграла (интеграл по поверхности возмущения масштабируется как единица). Если область интегрирования всегда соответствует области решения задачи, конкретно — числу атомов в молекуле, то и решение задачи Хартри—Фока будет иметь экспоненциальную сложность, и таким же образом будет расти число итераций сходимости процесса самосогласования. Напротив, если отклик на точечное возмущение локален и можно определить область интегрирования, не зависящую от размера системы, то сложность будет постоянной. Более точное описание сходимости должно включать в себя учет убывания членов ряда и возможную локализацию областей интегрирования в них, этим объясняется отсутствие места для полиномиальной сложности.

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

Для измерения окрестности релаксации точечного возмущения варьировались диагональные элементы матрицы плотности, относящиеся к выделенному атому вдалеке от границ молекулы графена. Вариация производилась после достижения сходимости процесса самосогласования, она составляла 10 % от точной величины и удерживалась до повторного достижения самосогласования. Полученная матрица плотности по модулю вычиталась из точной. На рис. 3 приведен пример локального и нелокального отклика для молекулы графена, содержащей 572 атома, при пороге БС отсечения 10-3 (а) и 10-4 (б). Точечное возмущение находится в центре изображений.

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

Измерения скорости сходимости подтвердили экспоненциальное возрастание числа итераций для случая нелокального отклика и его постоянство с момента увеличения молекулы до размеров, превышающих область отклика. Для параметра БС-отсечения 10 описанное изменение поведения соответствует молекуле графена размером 10x10 бензольных колец (всего 282 атома), для параметров отсечения 210-4 и 110-4 экспоненциальный рост числа

итераций наблюдается при увеличении стороны квадратного листа графена до 30 и 40 бензольных колец соответственно.

а)_ б)

Рис. 3

В соответствии с такими предположениями модель зависимости числа итераций п0 от размера молекулы может быть записана так:

По( N) =

п, ае

N / N0

< п,

ае

^ / N0

N < N

0'

(7)

а, N > N

0

где N0 — размер молекулы, при котором наблюдается изменение поведения; он, как и коэффициент а, определялся экспериментально. Первое из выражений в фигурных скобках введено для корректного поведения рассматриваемой модели вблизи нуля. Модель коэффициента гиперэффективности требует также знания поведения коэффициента по, которое не исследовалось теоретически. Эксперимент показал слабую зависимость по от числа фрагментов П~100, которая аппроксимировалась линейно.

Итоговая модель ускорения имеет следующий вид:

^ = Я (N, П)

"Ь( N)

щ (N / П )по (П )

по (П) = Ь + сП.

(8)

На рис. 4, а представлены графики коэффициента гиперэффективности, построенные в соответствии с предлагаемой моделью, б — ускорения решения задачи самосогласования Хартри—Фока для квадратных листов графена 98562 (1) и 46202 атома (2) на базе моделей (5) и (8). Коэффициент гиперэффективности, т.е. ускорение, которое было бы достигнуто только за счет изменения числа итераций до сходимости, принимает значения больше единицы в обоих рассмотренных случаях. Это, однако, не гарантирует наличия гиперускорения, которое проявляется только для случая графенового листа максимального размера. Более того, измеренная в ходе экспериментальных исследований эффективность не отвечает максимальному ускорению, которое могло бы быть получено. На рис. 4, б угол наклона секущей, проведенной из нуля к точке кривой для большего листа, превысит биссектрису координатного угла, только если точка будет располагаться вблизи максимума кривой. Для листа графена меньших размеров такая точка вообще отсутствует. Однако для обеих кривых найдется секущая, проведенная через две точки и имеющая „гиперэффективный" угол наклона. Это означает, что проведенные измерения эффективности действительно ничего

не могут сообщить о реальном ускорении, а применение самого модифицированного алгоритма БС целесообразно только в ограниченном диапазоне числа пространственных фрагментов (область определения графика на рис. 4, а, соответствующая значениям, превышающим единицу).

а) nG

2

1

0

Заключение. Предложенная модификация алгоритма DC позволяет снизить накладные расходы, пересчитанные на одну итерацию самосогласования. Однако эффективность по отношению к полному времени решения задачи самосогласования демонстрирует более сложное поведение, обусловленное зависимостью числа итераций самосогласования от размера молекулы. Проведенные измерения показали высокую эффективность и гиперэффективность использования вычислительных ресурсов относительно базы, взятой при запуске на 128 узлах (8 ядер на узле), однако это не является гарантией столь же высокой эффективности относительно гипотетического последовательного исполнения.

Работа выполнена в рамках контракта 07.514.11.4146 ФЦП „Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2007—2013 годы".

СПИСОК ЛИТЕРАТУРЫ

1. Bernholdt D. E. Scalability of correlated electronic structure calculations on parallel computers: A case study of the RI-MP2 method // Parallel Computing. 2000. Vol. 26. P. 945—963.

2. Degoli E., Ossicini S. Engineering Quantum Confined Silicon Nanostructures: Ab-Initio Study of the Structural, Electronic and Optical Properties. 2009. Vol. 58. P. 203—279.

3. Nakai H., Kobayashi M. Linear-scaling electronic structure calculation program based on divide-and-conquer method. 2011. Vol. 4. P. 1145—1150.

4. Duchemina I., Gygi F. A scalable and accurate algorithm for the computation of Hartree-Fock exchange // Computer Physics Communications. 2010. Vol. 181, N 5. P. 855—860.

5. Alizadegan R., Hsia K. J., Martinez T. J. A divide and conquer real space finite-element Hartree-Fock method // J. of Chem. Phys. 2010. Vol. 132. P. 034101.

6. Goedecker S. Linear scaling electronic structure methods // Rev. Mod. Phys. 1999. Vol. 71, N 4. P. 1085-1123.

7. Bolliger C. Linear Scaling Electronic Structure Methods. July 2008 [Электронный ресурс]: <http://www.math.ethz.ch/~kressner/students/bolliger.pdf>.

8. Lin L., Lu J., Ying L., Car R. Fast algorithm for extracting the diagonal of the inverse matrix with application to the electronic structure of metallic systems // Commun. Math. Sci. 2009. Vol. 7, N 3. P. 755—777.

Сведения об авторах

Андрей Игоревич Свитенков — Санкт-Петербургский национальный исследовательский универси-

тет информационных технологий, механики и оптики; Научно-исследовательский институт Наукоемких компьютерных технологий; инженер; E-mail: [email protected]

Дмитрий Михайлович Спельников

Владимир Григорьевич Маслов

Александр Валерьевич Бухановский —

Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики; Научно-исследовательский институт Наукоемких компьютерных технологий; младший научный сотрудник; E-mail: [email protected] д-р физ.-мат. наук; Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики; Центр „Информационные оптические технологии"; ведущий научный сотрудник; E-mail: [email protected] д-р техн. наук; Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, Научно-исследовательский институт Наукоемких компьютерных технологий; директор НИИ НКТ; E-mail: [email protected]

Рекомендована кафедрой высокопроизводительных вычислений

Поступила в редакцию 18.06.13 г.