CC BY
47
Математическое моделирование. Оптимальное управление Вестник Нижегородского университета им. Н.И. Лобачевского, 2012, № 5 (2), с. 67-72
УДК 004.94
РАСЧЕТ ОПТИЧЕСКИХ СВОЙСТВ ПОЛУПРОВОДНИКОВЫХ НАНОКРИСТАЛЛОВ В РАМКАХ ТЕОРИИ ФУНКЦИОНАЛА ПЛОТНОСТИ С ИСПОЛЬЗОВАНИЕМ ПАРАЛЛЕЛЬНОГО ПРОГРАММИРОВАНИЯ НА ГРАФИЧЕСКИХ ПРОЦЕССОРАХ
© 2012 г. В.П. Гергель, А.В. Горшков, А.В. Линев, Д.В. Осокин, А.М. Сатанин, А.В. Швецов
Нижегородский госуниверситет им. Н.И. Лобачевского
alexshdze @mail. rn
Поступила в редакцию 10.09.2012
Предложен эффективный параллельный алгоритм расчета оптических свойств многоэлектронных систем на основе динамической версии метода функционала плотности. Алгоритм реализован на вычислительном кластере с графическими процессорами на узлах. Продемонстрировано, что алгоритм может быть использован для определения спектров как изолированных, так и взаимодействующих кремниевых нанокристаллов (квантовых точек).
Ключевые слова: параллельные вычисления, графические процессоры, теория функционала плотности, нанокристаллы, многоэлектронные эффекты.
Введение
Расчет оптических свойств многоэлектронных атомных и молекулярных систем представляет огромный интерес для спектроскопии, лазерной физики и оптоэлектронных приложений. В современной оптоэлектронике актуальной проблемой является исследование воздействия лазерного поля на отдельные кремниевые нанокристаллы (квантовые точки) и ансамбли взаимодействующих нанокристаллов, в которых состояния оккупируют большое число электронов. Интерес также представляют сложные гибридные структуры, в которых проявляется резонансная связь дискретных состояний экси-тонов в полупроводниковых квантовых точках с плазмонами в металлических наночастицах, приводящая к усилению фотолюминесценции полупроводниковых квантовых точек [1-4]. Для описания возбуждений в наноструктурах необходимо решить уравнение Шредингера для многоэлектронной волновой функции в многомерном пространстве и времени. В случае систем, состоящих из десятков и более частиц, подход, основанный на прямом решении уравнения Шредингера для многочастичной волновой функции, является крайне затратным. Упрощение возможно, если ограничиться кругом явлений, где проявляются свойства наблюдаемых, определяемых только электронной плотностью. В основе теории Кона - Шэма (теории функционала плотности) лежит тот факт, что наблюдаемые многоэлектронной системы могут
быть выражены через электронную плотность и для этого не требуется находить многоэлектронную волновую функцию. Теория Кона -Шэма расширенна Рунге и Гроссом на случай нестационарных систем [5] и позволяет рассчитывать динамический отклик многоэлектронных систем.
Однако расчеты даже в рамках теории функционала плотности в случае большого числа частиц для последовательной реализации на CPU требуют больших затрат времени. Одним из выходов является разработка параллельной реализации данного метода на графических процессорах.
Расчет спектров электронных возбуждений в нанокристаллах
Для определенности рассмотрим взаимодействие сильного лазерного поля с электронной подсистемой в кремниевой квантовой точке. Расчет отклика проведем в рамках нестационарной теории функционала плотности. Отправной точкой теории функционала плотности является уравнение Шредингера для многоэлектронной функции у(?) в присутствии нестационарного поля. Гамильтониан системы Н = Т + + Уея + представляет собой сумму кинетической энергии электронов, потенциала взаимодействия электронов с ядрами решетки и внешними полями, а также энергии электрон-электронного взаимодействия. Одноэлектронный гамильтониан в квантовой точке берется в
приближении эффективной массы [6-8]. С учетом сказанного, гамильтониан многоэлектронной системы записывается в виде
^2 Ыр Ыр_
д
гп—ф / (г, і )=
нчо=-—+£и(?, і)+
г=1
2т
г=1
(1)
і)=Е|фі(гі^,
і=1
ді
V2
2т,
■+и^5 [«; Ф о](г, і)
Ф і (г, і), (3)
где те - эффективная масса электронов, которая для конкретных материалов известна из экспериментальных работ. Далее выбрана атомная система единиц: е = П = т0 = 1, где т0 -масса свободного электрона. В выражении (1) фигурирует внешний потенциал и(г, і), который согласно теореме Рунге - Гросса однозначно связан с электронной плотностью п(г, і) для многоэлектронной системы, эволюционирующей из начального состояния у(і = 0). Следовательно, если известна лишь электронная плотность системы, эволюционирующей из заданного начального состояния у(і = 0), то она однозначно определяет внешний потенциал, который ее создает. Внешний потенциал полностью определяет гамильтониан системы, тогда нестационарное уравнение Шредингера может быть решено и получены наблюдаемые характеристики системы. Для этого выбирается система невзаимодействующих электронов, называемая КБ-системой, определенная таким образом, что точно воспроизводит электронную плотность истинной взаимодействующей системы. Тогда свойства истинной системы могут быть выражены через электронную плотность невзаимодействующей системы. Поскольку однозначное соответствие между зависящей от времени электронной плотностью и внешним потенциалом установлено для любого !¥ее, то оно справедливо и для \¥ее = 0, то есть и для КБ-сис-темы. Следовательно, внешний потенциал иКБ (г, і) = иКБ [п; Ф0](г, і) невзаимодействую-
щей системы, воспроизводящий электронную плотность п(г, і), эволюционирующую из состояния Ф0, однозначно определен. КБ-сос-тояние Ф0 невзаимодействующих электронов выбирается в виде детерминанта одночастичных орбиталей ф(г, 0). При этом КБ-состояние должно соответствовать заданной начальной электронной плотности, воспроизводить ее саму и ее первую производную по времени. Электронная плотность взаимодействующей системы может быть определена как
икБ [п; Ф0](г, і) = и[п; ^](г, і) + 3 , е2п(г', і)
(2)
где Фі (г, і) являются решениями уравнения Шредингера с потенциалом иКБ (г, і) :
- + ихс [и; ^0, Фо](г, ?). (4)
|г - г|
Здесь и(г, ?) - потенциал внешнего поля, второе слагаемое представляет собой потенциал Хаар-три ия(г, ?), третье слагаемое - обменно-корреляционный потенциал. Функциональная зависимость и(г, ?) определяется тем фактом, что внешний потенциал и(г, ?) является заданной функцией (складывается из потенциала, связанного с барьером на границе нанокристалла со средой, в которой он находится, и потенциальной энергии во внешнем электромагнитном поле волны). Для аппроксимации ихс( г, ?) используется адиабатическое приближение локальной плотности [9], которое хорошо зарекомендовало себя при расчетах спектров многоэлектронных систем. В этом случае ихс зависит лишь от значения электронной плотности в данной точке в данный момент времени. Таким образом, задача сводится к решению Ыр нестационарных уравнений Шредингера (3) с потенциалом (г, ?)
для орбиталей ф;- (г, ?).
Численный метод решения динамических уравнений Кона - Шэма
Для демонстрации общей методики исследуем динамику электронных состояний в сферической потенциальной яме с бесконечными стенками. Данная модельная система должна отражать свойства сферической квантовой точки с достаточно большим потенциальным барьером (это оправдано тем, что в реальных системах барьер составляет несколько электрон-вольт).
Согласно схеме решения нестационарных задач в методе функционала плотности необходимо найти точное основное состояние многоэлектронной системы с учетом электрон-электронного взаимодействия до включения возмущающего потенциала, зависящего от времени. Иными словами, используя самосогласованные уравнения Кона - Шэма необходимо решить стационарную задачу и найти уровни энергии в яме и соответствующие им орбитали. Основному состоянию будет соответствовать ситуация, когда заполнены нижние уровни энергии (по количеству электронов). К решению данной проблемы можно подойти и с другой стороны. Так как энергетические уровни и
соответствующие им орбитали в сферической яме без учета электрон-электронного взаимодействия известны, то можно расположить электроны на нижних уровнях и, медленно включая возмущение в виде потенциалов Хар-три и обменно-корреляционного потенциала, решать самосогласованные уравнения Кона -Шэма для нестационарной задачи. При этом система адиабатически выйдет на стационарные уровни энергии многоэлектронной системы в сферической яме. Если включение потенциалов достаточно медленное, то возбуждение электронной системы на высоколежащие уровни не происходит и электроны остаются на нижних уровнях, что и соответствует основному состоянию.
На первом шаге происходит расчет электронной плотности согласно формуле (2). Далее решением уравнения Пуассона определяется потенциал Хартри для данной и(г, ?):
ДиН (г, і) = -4п(п(г, і) - п+ (г)).
(5)
Здесь и+(г) - плотность положительных зарядов. Появление и+(г) в уравнении Пуассона вытекает из требования электронейтральности нанокристалла. Соответственно интегралы по объему от и+(г,) и и(г, ?) должны быть равны. Плотность и+(г) выбирается в виде равномерно размазанного по всему объему нанокристалла положительно заряженного желе. После этого определяется обменно-корреляционный потенциал по аппроксимирующей формуле [9]:
-'хс
(г, і) = -
1.222
(п)
(
- 0.0661п
1 +
11.4
(п)
(п)= 3 ----------.
у 4пп(г, і)
(6)
Далее решаются нестационарные уравнения Шредингера для Ыр орбиталей (для КБ-системы) на одном шаге по времени. В результате находятся орбитали КБ-системы на новом шаге по времени, по которым рассчитывается электронная плотность на новом шаге по времени и происходит переход к расчету и(г, ?) на следующем шаге по времени.
Основные вычислительные ресурсы тратятся на решение уравнения Пуассона и Ыр нестационарных уравнений Шредингера. Поэтому для ускорения расчетов необходимо выбрать наиболее эффективные параллельные алгоритмы для решения этих задач.
К числу наиболее быстрых последовательных методов решения уравнения Пуассона относятся итерационные методы (в частности метод сопряженных градиентов). Уравнение Пу-
ассона, записанное в конечных разностях, представляет собой систему линейных алгебраических уравнений (СЛАУ) с N неизвестными значениями потенциала Хартри в узлах сетки. Итерационные методы основаны на последовательном приближении к истинному решению СЛАУ из начального приближенного. Применительно к задаче решения уравнения Пуассона не требуется хранить матрицу системы, поскольку необходим лишь результат действия этой матрицы на вектор. Результатом является вектор, каждая из компонент которого выражается всего лишь через семь (для семиточечной аппроксимации лапласиана) элементов вектора, на который действовала матрица, а не через все N. Таким образом, сложность метода 0(^Щ, где N -число итераций, необходимых для решения задачи, не превосходит число N. При хорошем начальном приближении число N резко уменьшается. Применительно к методу Кона - Шэма стоит отметить, что электронная плотность, полученная на новом шаге по времени, не сильно отличается от полученной на предыдущем, поэтому в качестве приближенного решения уравнения Пуассона удобно брать решение ин(г, і), полученное на предыдущем шаге по времени. В этом случае число итераций N должно быть невелико и метод является одним из самых экономичных.
К числу наиболее быстрых последовательных методов решения нестационарного уравнения Шредингера стоит отнести метод Рунге -Кутта. Метод Рунге - Кутта требует 0(Щ действий. Недостатком метода является то, что он явно не сохраняет норму волновой функции. Однако метод Рунге - Кутта допускает блочное разбиение пространства. При этом вычисление волновой функции на новом шаге по времени в каждом блоке происходит независимо от остальных. После чего требуется лишь пересылка значений волновой функции в граничных узлах для систем с распределенной памятью и всего лишь расширение блоков на граничные узлы для системы с общей памятью. Таким образом, на основе метода Рунге - Кутта может быть построен эффективный параллельный алгоритм решения нестационарного уравнения Шредин-гера.
В настоящей работе для решения самосогласованных уравнений Кона - Шэма выбраны итерационный метод сопряженных градиентов и метод Рунге - Кутта, поскольку допускают блочное разбиение пространства, и, следовательно, будут хорошо сочетаться в параллельной реализации метода Кона - Шэма для расчета многоэлектронных задач.
Результаты расчетов
Оценка эффективности алгоритма
На вкладке рис. 1б приведено распределение электронной плотности в сферической квантовой точке. При расчете предполагалось, что в квантовой точке возбуждены три электрона, волновые функции которых в начальный момент времени соответствовали низшим собственным состояниям частицы в сферической яме радиуса a (0=36^^ a0 - радиус Бора). Время включения т потенциалов Хартри и обменнокорреляционного значительно большее всех обратных частот, существующих в нашей системе (характерные частоты оценивались по известным частотам переходов между уровнями сферической ямы). Поэтому можно считать, что включение потенциалов адиабатически медленное.
Для определения спектра достаточно подействовать на систему коротким импульсом q и наблюдать за колебаниями среднего значения
дипольного момента
(e • r) = J e • n(r, t) • rd3
электронной системы (рис. 1а). Воздействие на систему короткого электромагнитного импульса q приводит к возмущению волновой функции основного состояния и может быть описано как изменение фазы волновой функции на фактор - дг(у^у ехр( -щг)) [10]. При этом пики
спектральной функции дипольного момента соответствуют энергетическим уровням системы (рис. 1б).
а)
ю[10 3m0c2a2/ И] б)
Рис. 1
17
Предложенный алгоритм расчета оптических свойств полупроводниковых нанокристаллов в рамках теории функционала плотности на основе блочного разбиения пространства был реализован для гетерогенных систем с распределенной памятью (вычислительных кластеров с графическими процессорами на узлах).
Реализация метода Кона - Шэма предполагает решение двух типов уравнений на каждом шаге по времени - уравнения Шредингера и уравнения Пуассона. Каждое из этих уравнений в отдельности может быть решено как с использованием графических процессоров, так и с использованием систем с распределенной памятью [11-13]. В частности в процессе решения общей задачи был разработан алгоритм решения нестационарного уравнения Шредингера для GPU, в основе которого лежит метод Рун-ге - Кутта. Результаты проведенных экспериментов (рис. 2) демонстрируют эффективность применения графических процессоров для решения данной задачи. Видно, что GPU реализация быстрее однопоточной CPU версии более, чем в 30 раз при размерах трехмерной пространственной сетки от 60 узлов в каждом направлении. Эксперименты проводились с использованием центрального процессора Intel Xeon L5630 2.13 ГГц, 24 ГБ оперативной памяти и графического процессора NVidia Tesla X2070.
Число узлов сетки Рис. 2
Реализация метода Кона - Шэма предполагает следующую схему работы: на каждом шаге по времени необходимо решить уравнение Шредингера для каждой частицы, а затем выполнить пересчет потенциалов для всех частиц путем решения одного общего уравнения Пуассона. При работе с несколькими ОРИ, установленными в рамках нескольких узлов вычислительного кластера, актуальной является проблема передачи данных между разными графическими процессорами. Данные нужно передавать на двух этапах в рамках временного шага: после
r
3
0
решения уравнений Шредингера и после решения общего уравнения Пуассона. Причем необходимо сначала передать данные из памяти GPU в оперативную память узла кластера, потом по сети другому узлу кластера, а затем опять в память GPU. Объем вычислений на единицу передаваемых данных здесь невысок, что приводит к низкой эффективности распараллеливания метода Кона - Шэма в целом. Следует, однако, отметить, что для возможности выполнения расчетов в случае большого количества частиц при значительных размерах трехмерной пространственной сетки имеющегося объема памяти как в рамках графического процессора, так и в рамках вычислительного узла будет недостаточно. А значит использование систем с распределенной памятью все же актуально для решения поставленной задачи.
Рисунок 3 показывает зависимость ускорения параллельной GPU реализации метода Кона -Шэма для систем с распределенной памятью в зависимости от числа используемых графических процессоров. Эксперименты проводились на кластере из 16 узлов, на каждом узле установлен процессор Intel Xeon L5630 с частотой 2.13 Ггц, 24 ГБ оперативной памяти и графический процессор NVidia Tesla X2070. Используемая операционная система - MS Windows Server HPC Edition Service Pack 2. Результаты проведенных экспериментов показывают, что при увеличении количества используемых GPU время работы программы практически не меняется. Эффективность распараллеливания здесь при использовании 16 графических процессоров составляет примерно 12.5%.
2.5 ■
2.0
0.5
0 5 10 15 20
Кол-во ОРИ
Рис. 3
На рис. 4 приведены результаты следующего набора экспериментов: рассматривается параллельная ОРИ версия для систем с распределенной памятью, и производится ее запуск на кластере с графическими ускорителями на каждом из узлов. Осуществляются параллельные запуски алгоритма на разном числе узлов кластера, причем количество процессов (узлов) соответствует размеру пространственной сетки, ис-
пользуемой для расчетов. Например, при запуске в 4 процесса обсчитывается в 4 раза больше данных, чем при запуске в 1 процесс. Необходимость применения такого подхода обусловлена тем, что при возрастании размера сетки объем данных возрастает так, что не помещается в память одного графического процессора, и для ее обработки требуется несколько GPU. Эксперименты проводились на кластере из 68 узлов, на каждом узле установлен процессор Intel Xeon X5670 с частотой 2.93 ГГц, 24 ГБ оперативной памяти и графический процессор NVidia Tesla X2070. Используемая операционная система - Clustrx T-Platforms Edition.
„, 14
I 0.2
о
© ^_______,____,____,___,____,___,____,____
0 20 40 60 80
Кол-во ОРИ / Размер данных
Рис. 4
Рисунок 4 показывает отношение времени работы алгоритма в P процессов к времени работы алгоритма в 1 процесс, умноженное на P, так как при запуске в P процессов обсчитывается в P раз больше данных. Этот показатель назовем формальным ускорением. Анализ этого показателя можно осуществлять по аналогии с анализом реального ускорения. В отличие от предыдущего, данный набор экспериментов проводился при использовании большего количества ОРИ. Эффективность распараллеливания при использовании 68 ОРИ здесь составляет примерно 31.4%. Таким образом, использование дополнительных графических процессоров в данной задаче наиболее полезно при расчетах с дополнительным объемом данных.
Заключение
Разработанный алгоритм решения динамических уравнений Кона - Шэма на основе блочного разбиения пространства может давать существенное ускорение при реализации на гетерогенных системах с распределенной памятью. Показано, что использование дополнительных графических процессоров в данной задаче наиболее полезно при расчетах с дополнительным объемом данных. Продемонстрирована эффективность распараллеливания при использовании 68 ОРИ, которая составляет приблизительно
31.4%. На основе данного алгоритма разработана методика расчета оптических свойств многоэлектронных систем - полупроводниковых квантовых точек. Показано, что данный алгоритм может быть использован для определения спектров как изолированных, так и взаимодействующих кремниевых нанокристаллов (квантовых точек) различной формы.
Работа поддержана грантами РФФИ № 12-0731161-мола, № 12-07-00546-а и госконтрактом
№ 07.514.11.4147Минобрнауки.
Список литературы
1. Shimizu K.T., Woo W.K., Fisher B.R., Eisler H.J., Bawendi M.G. // Phys. Rev. Lett. 2002. V. 89. P. 117401.
2. Lee J., Govorov A.O., Dulka J., Kotov N.A. // Nano Lett. 2004. V. 4. P. 2323.
3. Wei Zhang, Govorov A.O., Bryant G.W. // Phys. Rev. Lett. 2006. V. 97. P. 146804.
4. Ratchford D., Shafiei F., Kim S., Gray S.K., Li X. // Nano Lett. 2011. V. 11. P. 1049-1054.
5. Runge E., Gross E.K.U. // Phys. Rev. Lett. 1984.
V. 52. P. 997-1000.
6. Ferconi M., Vignale G. // Phys. Rev. B. 1994. V. 50. P. 14722.
7. Hirose K., Wingreen N.S. // Phys. Rev. B. 1999. V. 59. P. 4604.
8. Koskinen M., Manninen M., Reimann S.M. // Phys. Rev. Lett. 1997. V. 79. P. 1389.
9. Лундквист С., Марч Н. Теория неоднородного электронного газа. М.: Мир, 1987.
10. Puente A., Serra L. // Phys. Rev. Lett. 1999. V. 83. P. 3266.
11. Broin C.O., Nicolopoulos L.A.A. An OpenCL implementation for the solution of TDSE on GPU and CPU architectures // Computer Physics Communications. 2012. 183. P. 2071-2080.
12. Маркус Е.Д. Исследование параллельной реализации метода расщепления для уравнения теплопроводности на кластерных вычислительных системах // Вестник Томск. гос. ун-та. 1(14), 2011. С. 73-78.
13. Игнатьев А. А., Затевахин М. А. Параллельный метод для решения уравнения Пуассона // Параллельные вычислительные технологии: Докл. Между-нар. науч. конф., Н. Новгород, Россия, 30 марта-3 апреля 2009. С. 491-495.
CALCULATION OF OPTICAL PROPERTIES OF SEMICONDUCTOR NANOCRYSTALS IN THE FRAMEWORK OF DENSITY FUNCTIONAL THEORY USING GPU PARALLEL PROGRAMMING
V.P. Gergel, A. V. Gorshkov, A. V. Linev, D. V. Osokin, A.M. Satanin, A. V. Shvetsov
An effective parallel algorithm is proposed for calculating optical properties of many-electron systems based on the dynamic density functional method. The algorithm has been realized on a GPU cluster. It has been demonstrated that the algorithm can be used to determine the spectra of both isolated and interacting silicon nanocrystals (quantum dots).
Keywords: parallel computing, graphics processing units, density functional theory, nanocrystals, many-electron effects.
