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

Российская академия наук Институт вычислительной математики

На правах рукописи ИБРАГИМОВ Ильгиз Владимирович

УДК 519.6

АЛГЕБРАИЧЕСКИЕ СТРУКТУРЫ И ПАРАЛЛЕЛЬНЫЕ ВЫЧИСЛЕНИЯ В ЗАДАЧАХ КВАНТОВОЙ ХИМИИ

05.13.16 — Применение вычислительной техники, математического моделирования, и математических методов в научных исследованиях

ДИССЕРТАЦИЯ на соискание ученой степени кандидата физико-математических наук

Научный руководитель д.ф.-м.н Тыртышников Е.Е.

Москва 1998

Содержание

Введение 4

1 Задачи и методы нахождения электронного строения молекул 8

1.1 Постановки химических задач................................8

1.2 Вычислительные методы нахождения электронного строения молекул ..................................................14

1.3 Метод Фока решения многоэлектронного уравнения Шре-дингера в одноконфигурационном приближении..........24

1.4 Условие минимума энергии..................................26

1.5 Условие единственности минимума энергии................29

1.6 Численное решение уравнений Хартри-Фока..............33

1.6.1 Метод теории возмущений для решения уравнений Хартри-Фока......................................33

1.6.2 Метод тестовых функций для решения уравнений Хартри-Фока............................................36

1.7 Актуальные проблемы........................................39

2 Базисы на основе регулярных конечных элементов 41

2.1 Вычисление интегралов ......................................46

2.2 Свойства дискретизованного функционала Хартри-Фока 54 2.2.1 Вычисление значения дискретизованного функционала энергии Хартри-Фока..........................54

2.2.2 Вычисление первых производных дискретизован-ного функционала энергии Хартри-Фока..... 56

2.2.3 Вычисление вторых производных дискретизован-ного функционала энергии Хартри-Фока..... 57

2.2.4 О точных преобразованиях дискретизованного функционала энергии Хартри-Фока........ 58

2.2.5 Приближенное вычисление значения дискретизованного функционала энергии Хартри-Фока ... 60

2.3 Решение итерационными методами дискретных уравнений Хартри-Фока....................... 63

2.4 Построение предобуславливателей для метода Давидсона 67

2.4.1 Метод решения систем линейных уравнений с матрицей вида М + И.................. 68

2.4.2 Аппроксимация задачи С + И задачей М. + И . . 72

3 Численные эксперименты 75

4 Принципы распараллеливания программного комплекса по решению уравнений Хартри-Фока 83

4.1 Параллельные алгоритмы БПФ и скалярного произведения векторов.......................... 84

4.2 Теоретические исследования параллельных алгоритмов БПФ и вычисления суммы векторов............ 86

4.3 Численные эксперименты.................. 95

4.4 Параллельная реализация программного комплекса по решению уравнений Хартри-Фока............. 102

Заключение 109

Литература 110

Введение

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

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

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

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

Оператор Хартри-Фока задан в трехмерном пространстве и состоит из суммы трех операторов:

Т — кинетической энергии электронов;

У\ — потенциальной энергии электронов в поле ядер атомов;

V2 — потенциальной энергии электронов в поле, которое создается всеми остальными электронами.

По порядку линейный оператор Т больше двух остальных операторов У\ и V-?.

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

Для решения такой задачи методом дискретизации в этой работе предлагается использовать в качестве базисных функций конечные элементы на регулярной сетке. В этом случае дискретный оператор Хартри-Фока записывается в виде:

I -ь + бе, - - Е с;тсг + ¿(ад (т £ с;сЛ I Сз = О, V; - 1,..., м,

где Ь и 5 — ленточные трехуровневые [2] симметричные матрицы, причем £ может быть единичной; Т — циркулянтная трехуровневая положительно определенная симметричная матрица; С{ — ленточные трехуровневые симметричные матрицы, составленные из элементов г-ого собственного вектора С^; Е^ — ^'-ый собственный вектор и соответствующее ему собственное значение; М — число электронных пар в системе.

Такой подход позволяет применять итерационные методы типа Ланцоша [3] и Давидсона [4] так, что арифметические затраты умножения вектора на дискретный оператор имеют порядок N ./V против ТУ2 для базисных функций неструктурированного вида и Д^4 для базисных функций с бесконечным носителем. Также разработаны эффективные предобуславливатели для такой задачи с тем же порядком арифметической сложности.

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

При разработке этих коммуникационных примитивов была применена модель, в которой время передачи данных между процессорами выражается в виде сг + где а — время задержки для подготовительной работы (измеряется в количестве операций с плавающей точкой, которые можно сделать за время этой задержки), т — скорость передачи (слово/время выполнения одной операции с плавающей точкой), N — число слов в передаваемом пакете. Как оказалось, такая сравнительно простая модель позволяет хорошо описать реально наблюдаемые величины на компьютерах типа МВС-100/1000 [5].

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

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

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

Для двух других коммуникационных примитивов дело обстоит намного сложнее. В зависимости.от того, каковы величины <т и т параллельного компьютера, каков размер задачи (N) и число используемых процессоров (Р), меняется параллельный алгоритм, который при данных параметрах будет более эффективен. Поэтому предлагается следующий подход для нахождения наиболее эффективного алгоритма. Построим класс алгоритмов, которые будут зависеть от одного параметра, причем этот параметр должен, по возможности хорошо описывать все возможные параллельные алгоритмы. Далее в терминах модели а и т выпишем время работы алгоритма в зависимости о этого параметра, и минимизируем время работы алгоритма по этому параметру при заданных сг, г, N и Р. Для операции вычисления скалярного произведения такой подход однажды уже был применен для эффективной реализации этого алгоритма на компьютере CRAY — T3D [5].

Был выписан класс параллельных алгоритмов, зависящих от параметра (допустимые значения этого параметра — целые числа на отрезке [l,log2P]) и реализован на компьютере МВС-100. В результате удалось получить, что эти алгоритмы хорошо масштабируемы для реального числа процессоров.

На основе этих коммуникационных примитивов был разработан комплекс программ по решению уравнений Хартри-Фока. Было получено, что с помощью этого комплекса можно решить очень большие задачи (область 643, число электронных пар до 19, это почти 10 миллионов неизвестных) за реальное время, так как наблюдаемая производительность достигала 70 MFlop/s.

Глава 1

Задачи и методы нахождения электронного строения молекул

1.1 Постановки химических задач

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

На основе величины Е делаются выводы о химической устойчивости и реакционной способности молекул. Пусть дана химическая реакция

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

Я

Я

Я

Я

Рисунок 1. Вращение двух групп СЯ3— молекулы этана друг

здесь А — одна или несколько исходных молекул, В — один или несколько продуктов реакции. Величина Ев — Ед характеризует скорость и выход реакции.

Большой интерес вызывает значение Е некоторых переходных состояний молекул, называемых интермедиатами. Последние не устойчивы и не могут соответствовать какому-то реально существующему в природе веществу. Они образуются в процессе химических реакций. Часто время их существования (до превращения в другое вещество) меньше Ю-12 секунды. Так, если реакция (1.1) протекает через образование промежуточного продукта С, то есть А —> С —> Я, то величина Ее — Ед соответствует энергии активации, которая также значительно влияет на скорость реакции.

Вторая производная от Е по расстоянию между двумя атомами в молекуле соответствует частоте колебания молекулы по связи между этими атомами. Эта частота напрямую связана с наблюдаемой физической характеристикой молекулы — ее инфракрасным (ИК) спектром. Более того, на ИК-спектр может влиять не только зависимость Е от расстояния, но и зависимость Е от угла поворота двух частей молекулы относительно связи между двумя атомами (рис. 1).

Еще одна очень важная характеристика молекул — их электронная структура. Чтобы описать электронную структуру, надо учесть фун-

относительно друга.

А^гВ

(1.1)

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

Уравнение, которое служит математической моделью электронов, известно как уравнение Шредингера [8]. Для одноэлектронной системы оно имеет вид

д2ф д2ф д2ф 8т

где т — масса электрона, Е — его полная энергия, V — его потенциальная энергия в поле внешних сил, 7г — постоянная Планка. В физическом смысле ¡^ф*ф(10, — вероятность нахождения электрона в объеме О, а ф — принято называть орбиталями, или электронными облаками. Если V — это потенциальная энергия электрона в поле одного ядра атома, то известен аналитический вид функций ф и соответствующих значений Е:

mqj п 9)

(L2)

(1.3)

здесь qe — заряд электрона,

г,0,ф — сферические координаты,

L^+i1 — обобщенные полиномы Лагерра [8],

Yim — некоторая полиномиальная функция от cos в, sin в, cos ф и sin ф,

п,1,т — целые числа, такие, что для некоторого п — 0, ...,оо: I — 0,..., п — 1 и m = — I,..., 0,..., I:

Как видно, каждое собственное значение Еп вырождено п2 раз. Минимальное значение Е соответствует устойчивому в природе состоянию и наблюдается при п = I = т = 0.

я

Ц

/

/ я

С-

с

\

я

я

1

2

Рисунок 2. 1 циклопропан, 2 циклобутан.

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

Итак, очерчены числовые характеристики молекулы, которые связаны с их измеряемыми физическими величинами.

Далее определим "объект" и "инструмент": класс молекул, их конкретные числовые характеристики, которыми будем интересоваться, и, только потом, методы, с помощью которых будем вычислять эти величины (главы 1.2-1.6).

Будем изучать свойства веществ, содержащих малые циклы, то есть циклы, состоящие из 3-х или 4-х атомов, самыми простыми представителями которых являются циклопропан и циклобутан (рис. 2).

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

Рисунок 3. Малые циклы с насыщенными группами.

необычной для циклических молекул является способность трехчленного цикла к раскрытию и расширению цикла [7].

Несомненно, причину специфического поведения соединений этого класса следует искать в особенностях их молекулярного строения. Так угол связи между углеродами в циклопропане 60°, а в циклобута-не < 90° (из-за не плоского строения), в то время, как для предельных углеводородов соответствующий угол близок к тетраэдрическому, то есть агссоэ (—§) — 109°. Это приводит к аномальному строению электронной структуры и высоким значениям полной энергии.

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

Рассмотрим несколько подклассов этих веществ и акцентируем внимание на задачах, которые необходимо решить для таких веществ.

Малые циклы с насыщенными группами. В первую очередь очень интересным является зависимость Е молекулы при колебании и вращении молекул по связи В,2 с циклом (рис. 3). К сожалению не всегда из экспериментальных данных можно получить эти величины, так как такие вещества очень реакционно способны и не выделяются в чистом виде. В то же время такая зависимость очень информативная: по ней можно получить длины волн поглощения в ИК спектре молекулы. В данной задаче для практических нужд достаточно знать Е с точностью до 3 — 4-х знаков после запятой.

Малые циклы с ненасыщенными группами. Прежде всего интересна задача вращения молекулы по связи цикл —

сн2

н

3

4

5

Рисунок 4. Малые циклы с ненасыщенными группами.

атом при двойной связи (рис. 4). Сложность данной задачи заключается в том, что электроны двойной связи существенно взаимодействуют с электронами цикла. Так для очень простого соединения 3 (рис. 4), где Яг = Д3 = = Н, В,2 — СНз нет как достоверных экспериментальных данных, так и достоверных данных расчетов [9]. А для такого рода соединений очень важно знать электронную структуру (которую можно получить только из расчетных данных) для того, чтобы на основе ее предсказывать реакционную способность молекулы.

Еще менее изученными являются соединения с двойной связью в цикле или соединения, у которых один из углеродных атомов цикла замещен на О, N или 5. Для некоторых из таких соединений, например, 5 (рис. 4), работ по ядерной или электронной структуре практически нет, в то же время существуют работы по синтезу и использованию в химическом синтезе [9].

Все эти молекулы объединяют несколько свойств.

1) При значительном числе (до 20 атомов в молекуле) линейные размеры молекул очень незначительны. Почти для всех этих молекул2 расстояния между самыми да�