BC/NW 2008, №2 (13): 3.1

Характеристики трудоемкости и вычислительной эффективности прямых методов расчета электронной структуры больших молекул

Чернецов А.М., Шамаева О.Ю.

(Москва, Вычислительный центр им. А.А. Дородницына Российской академии наук, Московский энергетический институт (ТУ), Россия)

Расчеты электронной структуры гигантских молекул являются одними из самых сложных в современной науке и требуют использования высокопроизводительных вычислительных средств, таких как суперЭВМ и кластерные установки. Расчеты электронной структуры, в частности, биомолекул (белков, ДНК) и наночастиц, актуальны для ряда областей науки: химии, биохимии, физики конденсированного состояния вещества и др. В практическом плане эти расчеты важны для фармакологии, нанотехнологий, исследований явлений сверхпроводимости и др.

Основными методами расчетов электронной структуры молекул являются неэмпирические и полуэмпирические методы квантовой химии [1]. Неэмпирические методы расчета требуют слишком больших компьютерных ресурсов. Например, для расчета молекулы небольшого белка цитохром-С (1738 атомов) неэмпирическим методом DFT по специализированной программе ProteinDFT на кластере, содержащем 15 высокопроизводительных микропроцессоров HP Alpha 21264, требуется порядка 24 часов на выполнение одной итерации, а на весь расчет - больше одного месяца.

Для больших молекулярных систем, содержащих от 10³ до 10⁶ атомов, целесообразно применение полуэмпирических методов квантовой химии в так называемом приближении нулевого дифференциального перекрывания [1], в общем случае имеющих асимптотическую сложность расчета O(N³), где N – размерность задачи, пропорциональная числу атомов в молекуле. В основе математической модели для задачи расчета электронной структуры молекул полуэмпирическими методами - решение симметричной проблемы собственных значений. Поэтому асимптотическая сложность  этих методов определяется процессом диагонализации матриц и для сверхбольших молекул является неприемлемой даже при расчетах на суперЭВМ.

Для расчета физико-химических свойств нужны не сами собственные вектора, а матрица плотности P, являющаяся функцией от них.

При нахождении матрицы плотности P начальное приближение матрицы формируется из исходного фокиана F. Фокиан, в свою очередь, строится по определенным правилам из декартовых координат атомов, входящих в молекулу (см. [1, 2]). Сложность расчета фокиана для плотных матриц составляет O(N²).

Численный метод прямого нахождения матрицы плотности P –метод очистки (purification method) [3] - был разработан еще в 1960 г., однако в силу отсутствия в то время необходимых вычислительных ресурсов его применение было ограничено расчетом только небольших молекул. В 90‑х гг. ХХ века на основе этого метода были разработаны различные модификации, позволяющие ускорить процесс вычислений. Одной из таких модификаций является метод Пальцера-Манолополиса [4].

На основе метода Пальцера-Манолополиса в работе [5] предложен параллельный алгоритм расчета матрицы плотности для случая разреженных матриц с блочно - трехдиагональным портретом. На рис. 1 приведена структурная схема разреженных матриц с блочно-трехдиагональным портретом. Такую структуру матрицы имеют некоторые классы органических соединений, например полимеры, линейные несвёрнутые протеины и крупные фрагменты линейных цепей ДНК/РНК. Структурная схема параллельного алгоритма приведена на рис.2. Предложенный параллельный алгоритм реализован на языке Фортран 95 с использованием средств MPI и имеет линейную вычислительную сложность в зависимости от размерности задачи.

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

Рассмотрим оценки вычислительной эффективности процессов решения задачи расчета электронной структуры молекул указанными методами. Точнее, ускорение и эффективность, полученные за счет как параллельной реализации, так и применения свойств разреженности матрицы. Для этого приведем формулу итерационного расчета симметричной матрицы плотности P на (n+1)-ой итерации:

При этом начальная матрица плотности P формируется как фокиан F, нормализованный таким образом, что все его собственные значения лежат на отрезке [-1,1]. Для этого оцениваются минимальное Emin и максимальное Emax собственные значения. Останов итерационного процесса происходит при достижении заданной точности приближения. Рассмотрим более подробно использование свойств разреженности матриц с блочно-трехдиагональным портретом, представленных на рис. 1. Размеры каждого конкретного блока могут быть различны, поэтому введем параметр m=max(dim(Di)) - максимальный размер блока.

 

Рис. 1. Структурная схема разреженной матрицы с блочно-трехдиагональным портретом

 

Для такого вида матрицы целесообразно хранить массив блоков главной диагонали D1 и массивы блоков побочных диагоналей D2 и D3. При умножении матриц появляется ещё одна ненулевая побочная диагональ. Из квантово-химических свойств задачи мы можем ограничиться рассмотрением только 3-х блочных диагоналей – главной D1 и двух побочных D2 и D3. Так как P – симметричная матрица, то можно исключить массив D3 из дальнейшего рассмотрения.

Для формирования. компонент выражения (*) необходимо выполнить ряд матричных умножений. В результате формируются главная R1_i и побочная R2_j диагональ (где i=1..l, j=1..(l-1)):

Рис.2. Структура параллельного алгоритма расчета матрицы плотности методом Пальцера-Манолополиса.

Для приведенного на рис.2 алгоритма построим и исследуем теоретические временные оценки и оценки требуемого объема оперативной памяти при решении задачи расчета матрицы плотности P в зависимости от размерности, степени разреженности и числа вычислителей в кластере.

Рассмотрим оценки требуемого объёма оперативной памяти для одной итерации метода в случае реализации вычислений на одном вычислителе для плотных и разреженных (блочно-трехдиагональных) матриц.

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

Тогда, очевидно, число блоков главной диагонали nbl = [N/m]. Соответственно, требования к памяти для разреженных матриц имеет вид

V_разр(m) = [(2×m²×8×nbl + 2×m²×(nbl-1) ×8]×3 байт.

Общий объём памяти для плотных матриц можно оценить как:

V_плотн(N) = 3×N²×8 байт.

Следовательно, доля ненулевых блоков b определяется отношением:

b= V_плотн ¤V_разр.

Из анализа выведенной формулы получим, что при m=250 и N=7500 (nbl=30) доля ненулевых блоков составляет b=13.1%, а для m=250 и N = 7.5×10⁴ (nbl = 300) - b = 1.3%.

Перейдем к оценкам вычислительной сложности при следующих предположениях, которые могут быть получены из  выражения для P_n+1(*):

Пусть вероятности того, что коэффициент c_n £ ½ и c_n > ½

P{c_n £ ½} =q, P{c_n > ½} = 1-q.

Тогда трудоемкость вычисления одной итерации для матрицы плотности P зависит от:

Пусть u –время аддитивной операции, v – время мультипликативной операции.

Трудоемкости матричного умножения для плотных (T1) и разреженных (T2) матриц составляют соответственно:

Трудоемкости аддитивных операций для плотных (T3) и разреженных (T4) матриц составляют соответственно:

Тогда вычислительная сложность расчёта для плотных матриц есть

Для разреженных матриц трудоемкость вычисления составляет:

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

Рассмотрим случай использования модели распределенной памяти с множеством вычислителей p. Так как вычисления распределяются в среднем между вычислителями равномерно, то без ограничения общности положим, что каждый вычислитель обрабатывает одинаковое число блоков. В модели распределенной памяти необходимо учитывать время передач меду вычислителями, поэтому введем функцию F_swap(x) – взаимная передача x байт от одного процесса к другому (попарный вызов MPI_Send/MPI_Recv).

В этом случае для плотных матриц трудоемкость можно оценить как:

Для разреженных матриц получим:

Тогда можно оценить ускорение S(p) и эффективность Q(p) на p вычислителях как

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

Так, при реализации метода Пальцера-Манолополиса для плотных матриц максимально возможная размерность задачи составляла N= 8500 орбиталей (около 3500 атомов). При использовании свойств разреженности матрицы коэффициентов (случай блочно-трехдиагонального портрета) и при наличии оперативной памяти емкостью в 2 Гбайт можно производить вычисления матрицы плотности P существенно большей размерности, порядка N=150000 орбиталей (около 30000 атомов). Переход к гиперматричной схеме хранения позволяет работать со сверхбольшими молекулярными системами более общего вида.

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

Представленная работа проводилась в рамках РФФИ, проект 04-07-90220.

ЛИТЕРАТУРА

1.     Степанов Н.Ф. Квантовая механика и квантовая химия. – Изд-ва «Мир» и «МГУ», М., 2001, 518 с.

2.     Чернецов А.М. Распараллеливание процесса сборки фокиана. // Тринадцатая международная научно-техническая конференция студентов и аспирантов «Радиоэлектроника, электротехника и энергетика», в 3-х т., том 1. М. : МЭИ, 2007. – С. 379-380.

3.     McWeeny R., Rev. Mod. Phys. 32, 1960, p. 335.

4.     Anders M.N. Niklasson, C.J. Tumczak, Matt Challacombe. Trace resetting density matrix purification in O(N) self-consistent-field theory  J. of Chemical Physics, v. 118 N 15, 2003, pp 1-10.

5.     Чернецов А.М., Шамаева О.Ю. Расчеты разреженной матрицы плотности методом очистки на кластерных системах. Труды международной конференции «Информационные средства и технологии». в 3-х т.т. Т3. – М. : Янус-К, 2006. – С. 153-158.

6.     Писсанецки С. Технология разреженных матриц. – М. : Мир, 1988. – 412 с.