КАРДИНАЛЬНОЕ УСКОРЕНИЕ РАСЧЕТОВ ГИГАНТСКИХ БИОМОЛЕКУЛ МЕТОДАМИ КВАНТОВОЙ ХИМИИ, ТРЕБУЮЩИМИ ПРИМЕНЕНИЯ СУПЕРЭВМ И/ИЛИ GRID-СИСТЕМ

ББК Г511.4:В192.252 ГРНТИ 31.15.15, 50.07.05 УДК 544.18:004.67

Н. А. Аникин, А. Ю. Мускатин, М. Б. Кузьминский. С. Н. Леднев, А. В. Смирнов, А. И. Русаков

Кардинальное ускорение расчетов гигантских биомолекул методами квантовой химии, требующими применения суперЭВМ и/или ОИ.ГО-систем

Аннотация. Расчеты электронной структуры молекул квантовохимическими методами давно проводятся с использованием суперЭВМ. Сегодня они проводятся на лидере суперкомпьютерного списка TOP500 и будут осуществляться на первом в США экзафлопсном суперкомпьютере.

Краткий обзор современных методов квантовой химии и их применения на суперЭВМ для расчетов в первую очередь больших молекул показывает необходимость применения ускоренных аппроксимационных методик для реализации возможностей проведения таких расчетов. Это особенно актуально для массовых расчетов таких гигантских биомолекул, как докинг-комплексы белок-лиганд.

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

Изложенная методика позволила достигнуть ускорения на порядки и предполагается к использованию в будущих неэмпирических расчетах. Описанные методики и программы для необходимых массовых расчетов докинг-комплексов естественно вписываются в пакетную систему обработки заданий и могут использоваться в GRID-среде. Такая GRID-система создается на вычислительных ресурсах ЯрГУ и ИОХ РАН на базе стандартных в рамках EGI программных средств UMD 4).

Работа выполнена при финансовой поддержке РФФИ, проект №18-07-00657.

© Н. А. Аникин<, А. Ю. Мускатин<, М. Б. Кузьминский!, С. Н. Леднев< , А. В. Смирнов!, А. И. Русаков! , 2020

© Институт органической химии им. Н.Д. Зелинского РАН!1&2, , 2020 © Ярославский государственный университет им. П.Г. Демидова!4-5 , 2020

© Программные системы: теория и приложения (дизайн), Ц^-цл1

Введение

Методы квантовой химии, применяемые для расчета молекул, являются традиционной областью использования суперкомпьютеров (см., например, [1]). Квантовохимические расчеты осуществляются на первом в списке TOP500 (на ноябрь 2019) суперкомпьютере Summit [2],[3]. Знаменитый планируемый первым в США экзафлопс-суперкомпьютер Аргоннской национальной лаборатории [4],[5] будет использоваться для решения задач квантовой химии [6].

В неэмпирических и полуэмпирических квантовохимических методах для решения исходных интегродифференциальных уравнений используется аппроксимация волновых функций с применением базисных наборов функций. Далее мы имеем в виду только традиционные для молекул, экспоненциально убывающие с расстоянием от ядер атомов базисные функции—атомные орбитали, (AO) общего функционального вида в сферической системе координат

(1) ФМ,ф) = Ylm(9,<f)rl R(r),

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

вкe-a"r .

С ростом числа N используемых АО растет точность расчета, однако сильно возрастает и его ресурсоемкость, особенно время расчета. После расчета необходимых базирующихся на АО молекулярных интегралов, число которых есть O(N4), решение сводится обычно к некоторым задачам линейной алгебры. В разных типах этих интегралов используется до четырех центрированных на разных атомах AO, причем для больших молекул более всего лимитируют время расчета именно четырехцентровые интегралы.

Основные применяемые сегодня методы квантовой химии требуют для расчета электронного строения, включающего полную энергию E молекулярной системы, процессорного времени от O(N3) для традиционных быстрых и менее точных полуэмпирических методов до O(N!) для полного метода конфигурационного взаимодействия. Широко используемые методы теории функционала плотности (DFT) требуют O(N4). Число N линейно пропорциональна числу атомов, и поэтому суперЭВМ требуются для расчета E либо из-за необходимости

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

Часто используемые полуэмпирические методы и методы ХФ/БЕТ традиционно требуют нахождения собственных значений и векторов получаемой матрицы одноэлектронного гамильтониана (фокиана Г) размерности N х N:

(3) ГС = еБС,

где Б — матрица интегралов перекрывания, а выше речь шла об интегралах, необходимых для расчета матрицы Г. Решение этого уравнения требует O(N3) процессорного времени, дает молекулярные орбитали (МО) как линейную комбинацию АО с коэффициентами из матрицы С, и выполняется итерационно. Для очень больших молекул возможно отбрасывать малые интегралы (их доля существенно возрастает с ростом расстояний между удаленными «концами» молекулы), что дает уменьшение времени расчета нужных интегралов до О(№).

Для расчетов наиболее быстрыми полуэмпирическими методами и БЕТ таких больших молекул появились ускоренные вычислительные методики с отказом от традиционной диагонализации Г из уравнения (3) и с линейным О^) масштабированием времени расчета для гигантских молекул (более тысячи атомов), но с очень большим коэффициентом перед О^). Для гигантских молекул из многих тысяч атомов для достижения линейного масштабирования необходимо еще применять так называемый метод быстрых мультиполей ЕММ [7].

Полуэмпирические методы гораздо быстрее (но и менее точны), чем БЕТ, т.к. не требуют расчета такого большого числа интегралов. Имеется сообщение о проведении полуэмпирического расчета молекулы из 100 тысяч атомов с применением гибридного распараллеливания ОрепМР+МР1 на 1024 процессорах [8]. Использование суперЭВМ Ломоносов сделало возможным квантовые расчеты полных докинг-комплексов в РФ за разумное время с применением только полуэмпирических методов [9].

Более точные неэмпирические расчеты, в том числе по наиболее быстрому и широко используемому БЕТ, могут быть в 1000 раз медленнее, чем полуэмпирические [10]. В [11] приведены результаты расчетов методом БЕТ, адаптированным в программном комплексе

ONETEP на O(N) — расчеты больших молекул. Эти расчеты производились на входившем в TOP500 суперкомпьютере Iridis 4 на базе двухпроцессорных узлов с восьмиядерными Intel Xeon E5-2670/2.6 ГГц и межсоединения Infiniband FDR c использованием гибридного OpenMP + MPI распараллеливания.

Расчет такого вытянутого фибриллярного белка (реальные глобулярные белки будут считаться значительно медленнее) из 13696 атомов занял 7,3 час. При этом DFT-расчет белка на суперЭВМ Cray XC30 [12] проводился вообще не с применением гауссовых AO (соответствующих формулам (1) и (2)), а в базисе плоских волн, используемом обычно для расчетов твердых тел, а не молекул.

Все сказанное относится к так называемым одноточечным расчетам с фиксированной геометрией расположения атомов молекулы, но в ряде случаев требуется еще оптимизация геометрии (декартовых или внутренних координат) с поиском минимума E, что существенно увеличивает время расчета из-за необходимости одноточечных расчетов при разных геометриях, особенно для белков с их тысячами атомов (и соответствующих степеней свободы оптимизации геометрии).

Если же нужен расчет не только E, а и других термодинамических величин, то нужны другие классические для применения суперкомпьютеров методы вычислительной химии — молекулярной динамики, которые и на суперкомпьютерах сейчас основываются обычно на простейших неквантовых расчетах E, а на вышеуказанных квантовохимических уровнях в ближайшем будущем нереальны.

Поскольку полная оптимизации геометрии требует на порядки более высоких времен расчета, когда речь заходит о возможностях таких вычислений по DFT для нескольких сотен (менее тысячи атомов) с прицелом на расчет фрагментов белков [1З]. В этой работе оптимизация геометрии по программе Jaguar проводилась для молекулярных систем, содержащих до 600 атомов (до 10 тысяч базисных функций).

Расчеты, проведенные на известной Техасской суперЭВМ Stampede на базе связанных Infiniband FDR-межсоединением двухпроцессорных узлов с восьмиядерными Intel Xeon E5-2680/2,7 ГГц с использованием гибридного распараллеливания OpenMP + MPI, на 256 ядрах требуют около 76 минут, при этом достигается ускорение примерно в 114 раз. На рисунке 1 представлена зависимость времени расчета от числа используемых процессорных ядер, показывающая существенное снижение эффективности распараллеливания на большем их числе.

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

^ 100000

■3 10000

к Е ш О.

Рисунок 1. Зависимость времени расчета от числа использованных процессорных ядер [13]

для исследования фрагмента. Информация про подобные подходы, в том числе для актуальных докинг-комплексов содержится в [14]. Известной квантовохимической реализацией такого подхода является метод фрагментных молекулярных орбиталей ЕМО.

Такие расчеты выполняются и на суперкомпьютерах. Например, в расчете [15] на суперЭВМ Ломоносов-2 использовался вариант ЕМО-ТББЕТ. Однако мы такое направление в статье не рассматриваем.

Некоторые квантовохимические программные комплексы специально ориентируются на эффективные расчеты больших молекул на суперкомпьютерах, например, МТСИеш дал возможность считать по БЕТ молекулы с сотнями атомов [16] на японском К-компьютере, первым в мире достигшем производительности 10 РЕЬОРЯ. Едва ли не самыми востребованными сегодня являются расчеты огромных органических молекул, докинг-комплексов белков (протеинов, содержащих обычно больше 1000 атомов, но бывает немало и гигантских белков, содержащих порядка 100 тысяч и более атомов) с лигандами, что актуально, в частности, для задач конструирования лекарств.

При этом в каждом исследовании необходимы расчеты тысяч докинг-комплексов. Даже сейчас на фоне активного развития быстрых полуэмпирических методы квантовой химии и БЕТ для таких больших молекул некоторые авторы считают квантовохимические расчеты таких молекулярных структур целиком невыполнимыми в течение разумного времени [14].

Нами предложены новые аппроксимационные усовершенствования методов квантовой химии и разработаны соответствующие программные средства на Еог1гап-95, которые ориентированы исключительно на очень большие молекулы, главным образом массовые расчеты актуальных докинг-комплексов.

Один из современных способов ускорения расчета кулоновских двухэлектронных интегралов больших молекулах (одна из лимитирующих стадий) основан на аппроксимации произведения пар базисных функций АО от координат одного электрона в виде линейной комбинации гауссовых функций, называемых вспомогательными функциями плотности (сокращённо ВФП)

(4) ^¿И)^(г!) = /п(г1) ~ зЛГО =53 ВрХр(гч),

где ф —базисные (атомные) функции (орбитали, АО) отдельных атомов, а хр — ВФП. Для двухэлектронных интегралов используются укороченные обозначения:

(5) (// = Ц ^(г1)ф&(г1)

При этом различных ВФП для хорошей точности аппроксимации надо брать примерно в два раза больше, чем АО, что во много раз меньше квадрата числа произведений пар АО.

Для повышения точности используется прием сведения всех четырехцентровых интегралов от АО (требует О^4) вычислений) к трехцентровым интегралам от пар АО, которых О^3) с точностью до малых второго порядка по отклонению / от д:

(/п|/т) = (Зи + (/и - 3п)|3т + (/т - Ят)) =

= (Зи|Зт) + (/и - Зи|Зт) + (Зи 1 /т - Зт) + (/и - Зи^т - Зт) ~ ~ (/и|Зт) + (Зи|/т) - (Зи|Зт)

Здесь дт — линейная комбинация ВФП хр, а (/и|дт) — трехцентровые интегралы.

Для ускорения расчета двухэлектронных трехцентровых кулоновских интегралов (/и|дт), которые замедляют расчеты больших молекул, мы предлагаем следующее:

• Каждая функция х для исходного базиса с высокой точностью аппроксимируется комбинацией наших универсальных базовых гауссовых орбиталей (БГО) через аппроксимацию каждой входящей в х гауссовой экспонентом через комбинацию четырех БГО с показателями экспонент, близких к показателю гауссовой экспоненты из аппроксимируемой функции

(6) в-Ьг&2 « £ спв-апГ&2.

• Показатели экспонент БГО образуют специально подобранную геометрическую прогрессию, достаточную для хорошей аппроксимации всех актуальных ВФП х. Для повышения точности этой аппроксимации необходим набор показателей экспонент БГО, близких друг к другу, но при этом стандартная 64-разрядная точность становится недостаточной, что на порядки повышает погрешности округления. Мы преодолели это использованием вместо отдельных БГО наборов их линейных комбинаций с коэффициентами из матрицы Б-1/2, где Б — близкая к вырожденности матрица 4 х 4 интегралов перекрывания между БГО с разными, но близкими показателями экспонент.

Все это используется для последующего быстрого вычисления всех интегралов в виде сплайнов от Д. Интегралы для любых разнообразных базисов быстро рассчитываются в виде линейный комбинации заранее затабулированных сплайнов от Д для интегралов от универсальных БГО из нашей создаваемой БД.

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

Проведенная нами систематическая сплайн-аппроксимация интегралов перекрывания и кинетической энергии в зависимости от межъядерного расстояния для всех пар БГО с отобранными нами для аппроксимаций показателями экспонент показала вполне достаточный уровень точности расчета: точность интегралов порядка от 10-10 до 10-8 (чаще всего порядка 10-9, и лишь редко доходит до 4 • 10-8). Этого достаточно для вполне приемлемой точности квантовохимических расчетов при аппроксимации интегралов перекрывания и кинетической энергии от всех пар АО для широко распространенных расширенных базисов 6-3Ю* и 6-31++С** для всех атомов 1-3 периодов, образующих ковалентные связи.

Апробация проведена на расчете порядка миллиона интегралов перекрывания и кинетической энергии (исключая пренебрежимо малые, что быстро оценивается еще до этого расчета) для молекулы белка IMMUNOPHILIN FKBP-12 из 1663 атомов.

Расчет всех интегралов перекрывания и кинетической энергии белка FKBP-12 в базисе 6-31G* ускорен более чем в 20 раз. Сейчас нами разрабатывается более сложная аппроксимация для сверхбыстрого вычисления вклада кулоновских трех- и четырехцентровых интегралов в фокиан F — одну из лимитирующих стадий неэмпирческих (в том числе DFT) расчетов гигантских молекул. Метод открывает новые перспективы ускорения вычислений других типов молекулярных интегралов в произвольных центрированных на атомах базисах.

На рисунке 2 приведено изображение всех приблизительно 1700 атомов рассчитанного нами относительно небольшого белка — комплекса человеческого иммуноглобулина с иммунным супрессором IMMUNOPHILIN FKBP-12 из Protein Data Bank [17].

Выше уже отмечено, что многие другие белки еще крупнее. Но даже этот белок значительно крупнее большинства обычных органических молекул. Для наглядности мы привели на рисунке ориентировочное число вышеуказанных двухэлектронных кулоновских интегралов для использованного нами и очень часто применяемого в БГТ базиса 6-ЗЮ*. Для многих других белков число интегралов на порядки больше. Поэтому для весьма актуальных массовых (от

Рисунок 2. Типичный относительно небольшой протеин 1700 атомов)

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

По сравнению с эффективным линейно масштабируемым (заменяющим диагонализацию матриц) методом CG-DMS время расчета уменьшилось от 5 до 30 раз, по сравнению с другим линейно масштабируемым методом DivCon—от 18 до 120; выигрыш в быстродействии увеличивается с ростом размера белка.

При этом отличие наших рассчитанных ДЕ от результатов стандартного расчета всего докинг-комплекса с диагонализацией полной матрицы F составляет лишь около 0,4 ккал/моль (коэффициент линейной корреляции равен 0,997), т.е. все рассчитанные значения ДЕ весьма близки к результатам стандартной диагонализации полного фокиана.

В этом нашем базовом вышеописанном методе активная часть формировалась как универсальная (единая) для всех разнообразных лигандов из полного набора, содержащего типично более тысячи лигандов. Из-за этого для лигандов, близких к атомам из, например, «левой» части полости без нужды учитывался вклад многих атомов из других частей полости белка, удаленных от данного лиганда, т.е. N было завышенным.

Индивидуальное выделение активной (и соответственно замороженной) части белка для каждого лиганда приводило бы к чрезмерному росту расхода времени расчета на стадии расчета вклада замороженной части белка в F. Поэтому мы усовершенствовали этот базовый метод разбиением всего набора лигандов на группы лигандов [21], где каждая группа лигандов и соответствующая ей активная часть белка локализованы в своей отдельной части полости белка, что уменьшает размер необходимой его активной части.

Это ускоряет расчет, поскольку у каждой группы лигандов теперь имеется своя активная часть белка меньшего размера, чем была общая для всех лигандов, что отображено на рисунке 3.

Так, для вышеупомянутого выше белка FKBP-12 набор из 1308 лигандов был разбит на 17 групп, что уменьшило время расчета в 2,4 раза на одном процессоре AMD Athlon XP 3000+/2,2 ГГц до немного более суток. При этом среднеквадратичное отличие энергий образования докинг-комплексов от величин, полученных по нашему

базовому методу без образования групп, равнялось всего лишь 0,15 ккал/моль, т.е. точность расчета не ухудшилась.

Но мы и далее усовершенствовали (по сравнению с [21]) наш алгоритм образования групп, что ускорило расчет еще в 1,65 раза до 17,5 час, а среднеквадратичное отклонение энергии от не сгруппированного расчета стало 0,08 против 0,15 ккал/моль в предыдущем варианте алгоритма, т.е. расчет стал и быстрее, и точнее.

Указанное время расчетов полуэмпирическим методом AM1 не включает оптимизацию геометрии, и было бы желательно применять более точный, но требующий на несколько порядков больше времени расчета метод DFT. Поэтому наш метод массовых расчетов докинг-комплексов и его усовершенствования далее мы планируем перенести и программно реализовать для расчетов более точным методом DFT. Там ресурсоемкий расчет вклада двухэлектроных кулоновских интегралов предполагается сильно ускорить на основе нашей описанного в предыдущем разделе статьи аппроксимационного метода.

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

Проект РФФИ 18-07-00657 создает такую GRID-систему, включаю-щуую вычислительные ресурсы небольшого гетерогенного кластера в ИОХ РАН и часто (в том числе и самой Nvidia) относимого к суперкомпьютерам сервера DGX-1 c 8 GPU V100 (c пиковой DP-производи-тельностью 7,8 Терафопс каждый) и кластера в ЯрГУ.

Кластеры базируются на CentOS 7, а GRID-среда развивается на поддерживаемых EGI программных средствах репозитория UMD 4. GRID-сервер в ИОХ РАН базируется на промежуточном программном обеспечении ARC 15.3.19 из UMD-4.8.2, которое поддерживает все основные языки описания задач в GRID-системах. Служба ARC Resource-coupled EXecution service (A-REX) обеспечивает интерфейс с локальными системами пакетной обработки, в качестве которых у нас применяется Torque 4.2.10. Для обеспечения оптимальной работы с большими наборами данных используются программные средства dCache 3.2.21 из UMD-4.7.0. Данные результатов квантовохимических расчетов хранятся в формате CML версии 3 [22].

Распараллеливание программных комплексов, реализующих наши методические разработки, описанные выше, осуществляется c применением OpenMP, что позволяет не только эффективнее использовать аппаратные средства гетерогенного кластера в ИОХ РАН, но и задействовать в будущем все GPU V100 из DGX-1 в ЯрГУ, поскольку последний стандарт OpenMP 5.0 позволит очень существенно улучшить достигаемую производительность, в том числе за счет эффективной поддержки работы с иерархией памяти в GPU, а реальные реализации OpenMP в компиляторах дают возможность применения нескольких ускорителей на сервере.

Заключение

• Дан краткий обзор современных методов квантовой химии с ориентацией их применения для расчетов больших и гигантских (из тысяч атомов) молекул. Показано, что для актуальных массовых расчетов докинг-комплексов применение даже наиболее быстрых полуэмпирических методов и DFT требует не только применения суперЭВМ, но и разработки и программной реализации дополнительных вычислительных методов для кардинального ускорения расчетов.

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

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

• В этих целях на вычислительных ресурсах ИОХ РАН и Яр-ГУ создается GRID-система на основе ставшего современным европейским EGI-стандартом репозитория UMD-4.