Новые методы решения электронных уравнений квантовой химии и их применение тема автореферата и диссертации по химии, 02.00.17 ВАК РФ

На правах рукописи

ж

Митин Александр Васильевич

Новые методы решения электронных уравнений квантовой химии и их применение

02.00.17 - Математическая и квантовая химия

АВТОРЕФЕРАТ диссертации на соискание ученой степени доктора физико-математических наук

005531794

Москва - 2013

005531794

Работа выполнена на кафедре физической химии Химического факультета Московского государственного университета имени М. В. Ломоносова и в Федеральном государственном бюджетном учреждение науки Объединённый институт высоких температур РАН

Официальные оппоненты:

доктор физико-математических наук, профессор Виницкий Сергей Ильич

Объединённый институт ядерных исследований, г. Дубна

доктор химических наук, профессор Дьячков Павел Николаевич

Институт общей и неорганической химии имени Н. С. Курнакова РАН

доктор физико-математических наук Николаев Александр Васильевич Научно-исследовательский институт ядерной физики имени Д. В. Скобельцына МГУ

Ведущая организация: Национальный Исследовательский Центр

"Курчатовский Институт"

Защита состоится 10 октября 2013 г. в 15:00 часов в аудитории 446 на заседании диссертационного совета Д 501.001.50 по химическим и физико-математическим наукам при Московском государственном университете имени М. В. Ломоносова по адресу: 119991 Москва, Ленинские горы, д. 1. стр. 3, МГУ, Химический факультет

С диссертацией можно ознакомиться в Научной библиотеке МГУ имени М. В. Ломоносова по адресу: г. Москва, Ломоносовский проспект, д. 27. Автореферат диссертации размещён на сайте ВАК: http://vak.ed.gov.ru/

Автореферат разослан « 3 » СЛ-^с-С^^ 2013 г.

Ученый секретарь диссертационного совета, Д 501.001.50, к.х.н.

Матушкина Н. Н.

Общая характеристика работы

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

Это хорошо видно на примере изменения представлений о размере молекулярной системы, которая может быть рассчитана методом Хартри-Фока, и прогрессе в релятивистских расчётах. Так, в начале 80-х годов, большим представлялся расчёт аминокислотной пары С-С в базисе из 315 функций [1]. В середине 90-х большими считались расчёты протеинов в базисе из ~ 3800 функций [2], а в настоящей работе, на основе предложенных новых методов и алгоритмов, были выполнены расчёты гем-цитохрома с в базисе из ~ 8900 функций. При этом аналогичные тестовые расчёты этого же цитохрома проводились в базисе из ~ 16000 функций.

В начале 70-х годов, релятивистские расчёты обычно выполнялись в од-ноцентровом приближение методом Дирака-Фока для молекул с одним тяжёлым атомом [3]. В настоящее время релятивистские расчёты сформировались как отдельное направление квантовой химии, с разными методами учёта релятивистских эффектов [4], включая двухкомпонентный релятивистский метод теории функционала плотности (ТФП), представленный в диссертации, позволяющий проводить расчёты соединений с тяжёлыми и сверхтяжёлыми атомами, в которых число таких атомов может достигать ~ 500 [5, 6].

Необходимо отметить, что важнейшей частью подобных работ является проведение расчётов экстремальных по своему характеру, нацеленных на достижение наиболее высокой точности рассчитываемых молекулярных параметров и выяснение пределов их применимости. Такие расчёты позволяют дать более полную оценку методам, их модификациям, выявить их сильные и слабые стороны и провести их объективное сравнение. При этом особенно важно сопоставление рассчитываемых теоретических параметров с имеющимися экспериментальными данными и их совместный анализ, что ведёт к наиболее корректному описанию исследуемых молекулярных систем. Это наглядно демонстрируют высокоточные расчёты молекул IIс.^ и Вег методом конфигурационного взаимодействия (КВ), представленные в диссертации, показывающие, что метод КВ позволяет проводить расчёты двухатомных молекул на пределе точности неэмпирических методов.

Настоящая диссертационная работа включает оба этих направления, что

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

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

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

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

Разработана эффективная программа двухкомпонентного релятивистского метода ТФП с использованием релятивистского эффективного остовно-го потенциала и проведён первый расчёт атома элемента (117) и его двухатомной молекулы. Разработан многоссылочный метод КВ с явной корреляцией большого числа электронов, позволяющий рассчитывать состояния высокой мультиплетности. Предложено более 20 новых методов решения обобщённой задачи на собственные значения, включая методы типа Ньютона и Ньютона-

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

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

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

Проведены высокоточные неэмпирические расчёты малых молекул, демонстрирующие расширенные возможности неэмпирических методов. Получено численно точное решение уравнения Хартри-Фока для Яг с предложенным модифицированным приближением ЛКАО, корректно описывающим корреляционные свойства МО. Методом КВ, используя оптимальную методологию, рассчитаны потенциальные кривые Х1Т>д состояний Нв2 и Ве^ на пределе точности неэмпирических расчётов: диссоционная энергия Вег отличается от экспериментальной менее чем 1.5 см-1, колебательные уровни энергии вычислены со среднеквадратичным отклонением ~ 4 см-1. Предсказано существование в основном состоянии В&2 колебательного уровня с и = 11 в двух вращательных состояниях. Найдена нелинейная зависимость релятивистских поправок от межъядерного расстояния, показывающая важность учёта релятивистских эффектов для точного расчёта колебательных уровней энергий.

Получена новая оценка энергии диссоциации Ваг на основе сравнительного анализа рассчитанных неэмпирических потенциальных кривых основных состояний Б г 2 и В (¡2 и имеющихся данных. Рассчитаны наиболее точные ИК и Рамановский спектры молекул Л/2О2 и А^О^.

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

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

более точные значения рассчитываемых величин по сравнению с использованием известных базисных функций. Метод G4 для термохимических расчётов соединений переходных металлов построен с использованием модифицированных базисов m6-31G для атомов переходных металлов, предложенных нами.

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

Апробация работы. Основные результаты диссертационной работы докладывались на конференциях: Конференция по квантовой химии, 20-22 сентября, 1983 г., Днепропетровск; III Всесоюзное совещание по изучению структуры молекул в газовой фазе, 4-7 июня 1984 г., Иваново; IX Всесоюзное совещание по квантовой химии. 14-15 июня 1985 г., Черноголовка; V Всесоюзное совещание по изучению структуры молекул в газовой фазе, 11-14 июня 1990 г., Иваново; X Всесоюзное совещание по квантовой химии. Координационное совещание по квантовой химии. 16-21 сентября 1991 г., Казань; 29 Symposium fur Theoretische Chemie, 27 Sept. - 1 Oct. 1993, Oberwiesenthal, Germany; 31 Symposium fiir Theoretische Chemie, 10-13 Oktober 1995, Loccum, Germany; 9th International Congress of Quantum Chemistry, June 9-14,1997, Emory University Conference Center, Atlanta, Georgia, USA; Fourth European Workshop, Quantum Systems in Chemistry and Physics, 22-27 April, 1999, INJEP, Marly-Le-Roi, Prance; International Congress on Computational and Applied Mathematics, 17-21 July, 2000, Katholieke Universiteit Leuven, Leuven, Belgium; X Conference on Current Trends in Computational Chemistry, 1-3 November, 2001, Jackson State University, Jackson, Mississippi, USA; Xlth International Congress of Quantum Chemistry 2003, 20-26 July, 2003, University of Bonn, Bonn, Germany; 40 Symposium for Theoretical Chemistry. Computational Chemistry, 19-23 September, 2004, Suhl, Germany; Svensk Teoretisk Kemi 2006, 4-5 May, 2006, Stockholm,Sweden; 47th Sanibel Symposium, 22-27 February, 2007, St. Simons Island, Georgia, USA; International Conference of Computational Methods in Sciences and Engineering 2007, 25-30 September 2007, Corfu, Greece; 48th Sanibel Symposium, 21-26 February, 2008, St. Simons Island, Georgia, USA; Thirteen International Workshop on Quantum Systems in Chemistry and Physics, 6-12 July, 2008, Michigan State University, East Lancing, MI 48824, USA; 12th V.A.Fock Meeting on Quantum and Computational Chemistry, 19-23 October, 2009, Kazan, Russia; Всероссийская конференция "Радиохимия - наука настоящего и будущего", 13-15 Апреля 2011, Химический Факультет, МГУ им. М. В. Ломоносова, Москва; 14 European Symposium on Gas Phase Electron

Diffraction, 24-28 June, 2011, Chemical Department, Moscow State University, Moscow, Russia.

Публикации. Материалы диссертации опубликованы в 41 печатной работе, из них 34 статьи в рецензируемых журналах, рекомендованных ВАК [1-34], две статьи в зарубежных сборниках трудов [35, 36] и 5 тезисов докладов на международных и российских конференциях [37 41]. Список этих публикаций приводится в конце автореферата.

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

Структура и объем диссертации. Диссертация состоит из Введения, 10 глав, разделённых на две части, Выводов и списка цитированной литературы. В работе 345 страницы, 83 таблицы и 32 рисунка. Список цитированной литературы включает 605 наименований.

Основное содержание работы

Часть I. Разработка эффективных методов решения уравнений теории электронной структуры молекул

Глава 1 диссертации начинается с анализа известной литературы по линейной масштабируемости матрицы плотности и оценке числа ненулевых двух-электронных интегралов. Проведённый анализ показал, что, наряду с матрицей плотности, ненулевые одно- и двух-электронные интегралы, главные элементы матрицы плотности и наибольшие JIKAO коэффициенты МО должны обладать свойствами линейной масштабируемости по числу базисных функций. Рассмотрение асмптотических свойств числа ненулевых двух-электронных интегралов (N2e(e)) позволило сформулировать и доказать теорему о линейной масштабируемости Л^е(е) по числу базисных функций (N) в больших молекулярных системах.

Линейная масштабируемость по N была также исследована чис-

ленно для трёх молекулярных систем. Наиболее полные исследования были выполнены для полимеров аланина, рассчитанных в базисе 6-31G. Число атомов в полимере менялось от 112 до 7912, а число сгруппированных базисных функций от 609 до 43509. Основное внимание в этих расчётах было уделено исследованию зависимости Л^е) от N для низкой точности вычисления интегралов, так как анализ метода вычисления матрицы Фока с использованием разности плотности для оценки вклада двух-электронных интегралов в

матрицу Фока показывает, что именно такие интегралы определяют масштабируемость этого метода по N. Действительно, двух-электронные интегралы вычисляются, когда выполняется неравенство тах\{1]\к1) * тахАОы\ > £, где тах\(ц\Ы) - оценка максимального интеграла в блоке интегралов, а тпахАЮы - соответствующий максимальный матричный элемент разности матриц плотности. Величина с обычно равна Ю-6, Ю-7 в расчётах больших молекулярных систем [2]. Если теперь учесть, что критерий сходимости по разнице матриц плотности имеет, как правило, тоже значение, то можно видеть, что число ненулевых двух-электронных интегралов, больших по абсолютной величине, определяет свойство масштабируемости такого метода вычисления матрицы Фока.

На рис. 1, а и 1, б представлена зависимость Л'ге(е) и её первых конечных разностей от N для точности вычисления интегралов 1.0Е-2 и 1.0Е-3, а на рис. 1, в и 1, г даны аналогичные зависимости для точности 5.0Е-4, З.ОЕ-4, 2.0Е-4, 1.5Е-4 и 1.0Е-4. Анализ полученных результатов показывает, что зависимость Лг21,(е) от N не является постоянной и уменьшается с увеличением N. Зависимость А2с(е) масштабируется линейно по N для низкой точности вычисления интегралов почти для всех N, тогда как для более высокой точности вычисления интегралов А72С(е) масштабируется быстрее, чем линейная. В частности, для больших по абсолютной величине интегралов (1.0Е-2 и больше) ]У2е(е) масштабируется линейно при любом числе базисных функций. Для интегралов больше 1.0Е-3 по абсолютной величине линейная масштабируемость наблюдается для числа базисных функций 20000 и более (см. рис. 1, а и, особенно, рис. 1, б для первых конечных разностей). Анализ зависимости Лг2с(£) от N для точности 5.0Е-4 и выше (рис. 1, в и 1, г) показывает, что, с практической точки зрения, Лг2е(е) линейно масштабируется начиная с точности 1.5Е-4. Для более высокой точности линейная мастабируе-мость не наблюдается для рассмотренного числа базисных функций. Однако, выпуклый характер зависимостей первых конечных разностей указывает, что степенная зависимость ЛГ2С(б) от N ниже второго порядка. Таким образом, результаты численных исследований находятся в полном соответствие с доказанной теоремой.

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

609 11334 22059 32784 43509 Число базисных функций

609 11334 22059 32784 Число базисных функций

(а)

(<0

5

^ 2Е+10

-•-5.0Е-04 -•-Э0Е-О4 -»-2.0Е-04 /

-«-1.5Е-04 — 1 0E-CJ4 / ^

11334 22059 32784 Число базисных функций

(в)

11334 22059 32764 Число базисных функций

(г)

Рис. 1. Зависимость Лг2р(с) (а,в) и первых конечных разностей (б,г) от N в расчётах полимеров алашта в базисе (¡-310 для различной точности вычисления интегралов. На рис. (а) и (б) левая и правая вертикальная шкала соответствует точности 1.0Е-2 и 1.0Е-3.

для различных размеров молекулярных систем. Эти числа были сосчитаны в Хартри-Фоковских и полуэмпирических расчётах полимеров аланина различной длины. Хартри-Фоковские расчёты были выполнены с базисом 6-31G, а в полуэмпирнческих расчётах использовался РМЗ гамильтониан [7]. Число базисных функций в этих расчётах изменялось от 609 до 3359 и от 527 до 7027 соответственно. Геометрия полимеров аланина была оптимизирована методом молекулярной механики, используя программу Tinker [8]. Полученные результаты представлены на рис. 2, а и 2, б. Зависимости, представленные на рисунках, показывают линейное масштабирование числа главных матричных элементов матрицы Фока в зависимости от числа базисных функций.

Было показано, что только наибольшие по абсолютной величине матричные элементы матрицы Фока важны в неэмпирических расчётах, а от-

-«-1.0Е-01

-*~1.0Е-02

-м-ЮЕ-ОЗ

-»-1.0Е-04

-•-1.0Е-05

-Ч-1.0Е-06

х 2 8Е*5 о

5

1158 1709 2259 2809 Число базисных функций

2152 3777 5402 Число базисных функций

(б)

Рис. 2. Зависимость главных матричных элементов матрицы Фока от N для различной точности их расчёта, полученной в расчёте полимеров аланина методом Хартри-Фока с базисом С-ЗЮ (а) и РМЗ полуэмпирическнм методом (б).

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

Установление линейной масштабируемости ЛГ2е(е) от /V, сделало актуальным проблему объяснения линейной масштабируемости вычисления матрицы Фока в методе быстрых мультиполей [9, 10] и квантово-химического древовидного кода [11], поскольку они давали объяснения основываясь на зависимости ЛГ2е(е) от N выше квадратичной [12].

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

2246 3896

Число базисных функций

2246 3346 4446 Число базисных функций

— -»-1.0Е-07 -♦-1.0Е-08 -Ф-1.0Е-09

2246 3696

Число базисных функций

(б)

2246 3346 Число базисных функций

Рис. 3. Зависимость времени вычисления матрицы Фока (а,б) и его первых конечных разностей (в,г) от N для разной точности вычисления двух-электронных интегралов, полученная в расчёте полимеров аланина в базисах 3-21ЭР (а,в) и 6-ЗШ (б,г).

тегралов в матрицу Фока.

Свойство масштабируемости обычного метода вычисления матрицы Фока с использованием разности матриц плотности для оценки вклада интегралов было численно исследовано в расчётах полимеров аланина методом Хартри-Фока с базисами ЭТО-ЗС, 3-218Р и б-ЗЮ при изменение числа блоков аланина от 10 до 100. Время вычисления матрицы Фока и его первых конечных разностей в зависимости от N для различной точности вычисления двух-электронных интегралов представлены на рис. 3, о, 3, б и рис. 3, в, 3, г для базисов 3-21БР и 6-ЗШ, соответственно.

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

Таблица 1. Время (в мин.) расчёта молекулы порфина (Дгл) методами Хартри-Фока и ТФП с различными базисами на ПК с Intel Pentium III 800 MHz CPU.

БФ Хартри-Фок ТФП/ВЗЬУР

Класс. Прямой Пр/Кл Класс. Прямой Пр/Кл

STO-3G 134 0.37 1.39 3.8 1.15 3.55 3.1

6-31G 244 1.93 6.17 3.2 3.59 10.5 2.9

6-311G 354 5.37 14.9 2.8 7.63 20.8 2.7

6-31G** 430 8.78 22.5 2.6 11.6 29.5 2.6

6-311G** 516 18.5 45.5 2.5 22.8 56.4 2.5

cc-pVDZ 406 12.2 43.9 3.6 15.3 50.2 3.3

cc-pVTZ 916 185.6 551.1 3.0 204.5 620.7 3.0

Таблица 2. Расчёты больших молекул методом Хартри-Фока. Требуемое дисковое пространство (в GB), время (в сек.) пользователя (Т(п)) и системы (Т(с)) при вычисление

матрицы Фока.

Молек. Базис AT. БФ Диск CPU T(n) T(c)

Indinavir cc-pVTZ 92 2280 519.6 Alpha, 1 GHz 25919 2544

(Gly)7o 6-31G 493 2953 17.6 Р III, 600 MGz 1094 160

6-31G 900 3900 57.5 P III, 600 MGz 6053 655

ДНА PDB 123D 3-21SP 405 3717 30.8 P III, 600 MGz 1485 229

ДНА PDB 330D 3-21SP 485 4453 45.4 P III, 600 MGz 2936 471

Ac-(Ala)ioo-NH2 6-31G 1009 5546 477.9 Alpha, 1 GHz 25504 2242

ДНА PDB 433D 6-31G 664 6080 363.5 Alpha, 1 GHz 16851 1474

Фока начинается примерно с 2000 базисных функций для базиса ЯТО-ЗИ и примерно с 3000 для валентно-расщеплённых базисов 3-213Р и 6-ЗЮ. С увеличением точности двух-электронных интегралов, начало линейного масштабирования вычисления матрицы Фока сдвигается в стороны большего числа базисных функций. Особенно хорошо это видно на рис. 3, в, 3, г, где представлены первые конечные разности.

Использование сжатия двух-электронных интегралов и их индексов позволяет не только придать методу вычисления матрицы Фока с использованием сохраняемых интегралов свойство линейной масштабируемости, но также является одним из путей повышения эффективности классического метода вычисления матрицы Фока. В связи с этим нами были разработаны наиболее эффективные алгоритмы блочного и индивидуального сжатия интегралов и их индексов. В результате был создан модифицированный обычный метод решения уравнения Хартри-Фока и ТФП со сжатием данных и линейно масштабируемым вычислением матрицы Фока. Сравнение этого метода с прямым методом на примере расчёта молекулы порфина, представленное в Табл. 1, показывает, что новый метод существенно превосходит последний.

Для демонстрации расширенных возможностей обычного метода ССП со сжатием данных нами были проведены расчёты методом Хартри-Фока больших молекулярных систем, детали которых представлены в Табл. 2. Эти данные показывают, что квантово-механнческие неэмпирические расчёты больших молекулярных систем возможны и могут проводиться на современных

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

Однако, не все шаги нового метода линейно масштабируются по N. Исследование нового алгоритма в расчётах с числом базисных функций до 15871 показало, что он будет эффективно работать до примерно 30000 функций для базисов типа 6-ЗШ. После этого шаги с кубической зависимостью от N (обобщённая задача на собственные значения, вычисление вектора ошибки) будут составлять заметную часть общего времени расчёта, приводя тем самым к существенно нелинейной масштабируемости этого алгоритма по N.

Число атомов

(а) (б)

(в) (г)

Рис. 4. Теоретическая зависимость времени выполнения векторных вычислений от длины вектора (а), время умножения матрицы Фока на орбитали (б), время ортогонализации орбиталей (в) и время вычисления матрицы плотности (г). Все времена даны в сек.

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

расчётах полимеров аланина с гамильтонианом РМЗ. Времена выполнения основных шагов этого метода в зависимости от числа атомов, а также теоретическая зависимость времени выполнения векторных вычислений от длины вектора, представлены на рис. 4. Полученные зависимости похожи друг на друга, также как и на теоретическую зависимость времени выполнения векторных операций от длины вектора. Это показывает, что техника разреженных матриц трансформируют основные кубически зависимые шаги метода ИСО в линейно зависимые.

В Главе 2 представлены результаты исследования проблемы сходимости ССП итераций по двум направлениям. Первое связано с нарушением принципа заполнения орбиталей или перекрыванием спектров орбитальных энергий заполненных и виртуальных орбиталей, проявляющееся в виде ос-цилляций полной энергии в ССП итерациях. Оно может возникать как для неточных орбиталей в ССП итерациях, так и быть обусловленным свойствами оператора Фока и Кона-Шама.

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

В качестве модифицированного оператора Фока нами было предложено использовать оператор

Fз = ^ - DтJDт ,

где .Р - оператор Фока, .7 - оператор кулоновского взаимодействия и От -полная матрица плотности. Этот оператор сдвигает орбитальные энергии на величины определяемые оператором 3. Тестирование этого оператора в расчёте состояний 2ПГ и 2П,; молекулы ЫО показало, что итерации с оператором Фока -Р не сходились, тогда как с оператором /з решения были получены за 26 итераций для 2ПГ и за 67 итераций для 2П; состояния.

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

(а) (б)

Рис. 5. Изменение полной энергии молекулы ВС в состояние 2Пг на начальных итерациях (а) и в конце итерационного процесса (б) при различной величине сдвига уровней.

-60.00

-60.50

-61.00 -61,40

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

В связи с этим нами были рассмотрены итерации с модифицированным оператором Фока вблизи самосогласованного решения (линейное приближение) и получено неравенство, определяющее выбор параметров динамического сдвига уровней. Оно позволило дать формальное описание метода. Метод был протестирован в расчётах состояний с открытыми оболочками молекул ВвЫ и ВС. Так, в расчёте состояния 2П, молекулы ВС решение было получено за 35 итераций при использование метода динамического сдвига уровней, тогда как наименьшее число итераций в методе постоянного сдвига уровней было 50 при оптимальном выборе величины сдвига.

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

Однако, все эти выводы были проверены для молекулярных систем, образованных атомами первых двух периодов, в которых плотность спектра орбитальных энергий не является высокой. В наших исследованиях сложного молекулярного комплекса: железо - порфириновое кольцо аксиально координированного имидазольными группами (\FePIrri2]), с высокой плотностью спектра орбитальных энергий, было обнаружено необычно большое изменение орбитальных энергий ¿-орбиталей переходного металла в нейтральном комплексе и его одно- и двух-зарядных катионах, приводящее к принципиальному изменению порядка следования орбиталей по энергии. На рис. 6 видно,

к *

(а) (б) (в) (а)

Рис. 6. Высшие заполненные орбитали [РрР/гл.г] комплекса (а,в) и его катиона (б.г) рассчитанные методом теории функционала плотности с ВЗЬУР функционалом (а,б) и методом Хартри-Фока (в,г) в базисе 6-ЗЮ**.

что при ионизации основного синглетного состояния комплекса \FePIrii2] с замкнутыми оболочками, метод ТФП показывает значительную перестройку верхних орбиталей, тогда как в методе Хартри-Фока подобная перестройка не наблюдается. При этом Малликеновский и КВО анализ заселённости ТФП орбиталей не согласуется с наблюдаемой перестройкой ТФП орбиталей, в то время как для Хартри-Фоковских орбиталей наблюдается полное согласие.

Наше объяснение наблюдаемого нарушения исходит из того, что в методе Кона-Шама одночастичные ортогональные функции фг(г) - орбитали вводятся для представления электронной плотности в виде

N

п{г) = ^Шг)\2 ,

г=1

где N - число электронов, дающее возможность вывода уравнения Кона-Шама уже для орбиталей. Такое представление полной электронной плотности очень похоже на аналогичное представление плотности в методе Хартри-Фока. Однако, метод Хартри-Фока строится на основе представления об ор-биталях, позволяющего определить открытую оболочку. В методе ТФП орбитали вводятся как некоторые ортогональные объекты для представления

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

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

= (Ь)"Ак,

где Дк = Ск — Ск-1, а Ск - вектор, полученный на к-ой итерации. Отсюда, в пределе бесконечного числа итераций, получается уравнение

= ск + [/ - (ьуг1 (Ск+Р - Ск), (1)

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

Лх0 + оА) = а/(х0 + к) + (1 - а)/(®0) + 0{Ь?) ,

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

Всего было рассмотрено три класса экстраполяционных методов: линейные экстраполяционные методы, линейные проекционные экстраполяцион-ные методы, экстраполяциооные методы на основе специальных функций и предложено много новых экстраполяционных методов. Среди них выделим четырёх-точечный линейный экстраполяционных метод

С* = — Ск+р,]Ск+(ц р = 1 2 ¿=12

ск+С1+р,] — Ск+р,] — ск+,и + Си

и линейный проекционный экстрополяционный метод

РГ

С*=Ск +

£г

и-1- (А<)Р

(Ск+Р — Ск) ■

В Главе 3 представлены результаты полученные при разработке нового алгоритма оптимизации геометрии многоатомных молекул. Обычно оптимизационные алгоритмы строятся на основе квази-Ньютоновских методов [15] с использованием первых производных и некоторого метода обновления матрицы Гесса. Скорость сходимости таких методов определяется главным образом

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

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

В связи с этим нами были проанализированы методы обновления матрицы гессиана Давидона-Флетчера-Пауэла (DFP), симметричная формула ранга один (SRI) и формула Бройдена-Флетчера-Голдфарба-Шано (BFGS). Опыт применения этих формул показал, что в общем случае оптимизационные алгоритмы с формулой BFGS, имеют наибольшую скорость сходимости. Однако, нами было замечено, что для положительно определённых квадратичных функций доказана теорема, утверждающая, что оптимизационный алгоритм с SRI формулой и линейным методом определения оптимальной длины шага достигает минимума, по крайней мере, на. п + 1 шаге и, если требуется n + 1 шаг, то полученный гессиан является действительно матрицей вторых производных в точке минимума [16]. Эта теорема не верна в общем случае для алгоритма с BFGS формулой, если не используется точный линейный метод определения оптимальной длины шага. Таким образом, эта теорема указывает не некоторое преимущество SRI формулы перед BFGS формулой для квадратичных функций. Поэтому можно полагать, что использование SRI формулы в процедуре оптимизации геометрии молекул будет более успешным, чем использование BFGS формулы.

В результате нами был разработан метод оптимизации геометрии молекул основанный на использование формулы SRI для обновления матрицы гессиана. Его другой отличительной чертой явилось использование модифицированного критерия окончания оптимизационного процесса.

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

1. малые изменения координат атомов вызывают малое изменение полной энергии молекулы,

2. малые изменения координат атомов вызывают большое изменение полной энергии молекулы,

3. большие изменения координат атомов вызывают малое изменение пол-

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

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

Первоначально, для тестовых расчётов был выбран набор из 30 молекул [17]. Однако, оптимизация геометрии этих молекул показала, что найденные нами минимумы семи молекул находятся ниже по энергии, чем минимумы указанные в [17]. Для шести молекул оптимизационный процесс прямо сходился к точкам минимума начиная с исходной геометрии. Сравнение геометрий молекул в новых минимумах с указанными в [17] показало, что различие связано с вращениями одной группы атомов относительно других. Ещё для одной молекулы оптимизационный процесс сходился к известной точке минимума. Однако, характер сходимости позволил предположить, что между начальной и конечной геометриями имеется критическая точка и, следовательно, другой конформер этой молекулы может иметь близкую энергию. В этой связи начальная геометрия молекулы была изменена и новая точка минимума с более низкой энергией была обнаружена. По этим причинам только 23 молекулы были использованы в наших тестовых расчётах.

Известно, что математические свойства формул обновления матрицы Гесса не зависят от координатной системы в которой они записаны. В тоже

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

Результаты расчётов показывают, что в среднем по всем расчётам, оптимизационный метод с BFGS формулой требует на 25.8 % и 33.2 % больше вычислений полных энергий и градиентов полных энергий по геометрии, соответственно, чем метод с использованием SRI формулы. При этом в декартовых координатах оптимизационный метод с BFGS формулой требует большего числа вычислений энергии и градиентов на 35.5 % и 45.8 %, соответственно, а во внутренних координатах оптимизационный метод с BFGS формулой требует большего числа вычислений энергии и градиентов на 26.2 % и 34.1 %, соответственно.

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

В Главе 4 представлен эффективный алгоритм параллелизации по узлам компьютерного кластера методов Хартри-Фока и ТФП. Важным преимуществом такой параллелизации для классического метода ССП с построением матрицы Фока из сохраняемых интегралов состоит в том, что ввод интегралов с внешних устройств на всех узлах может совершаться параллельно, суммируя тем самым скорость передачи данных между внешними устройствами и оперативной памятью. При этом отметим, что создание в последнее время твёрдотельных дисковых накопителей, очевидно даёт новые преимущества классическому методу ССП решения уравнений Хартри-Фока и Кона-Шама, поскольку новый тип дисков обеспечивает существенно более высокую скоростью передачи данных по сравнению с традиционными дисками.

Известные алгоритмы параллелизации ССП методов решения уравне-

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

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

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

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

балансировки узлов, было выбрано число двух-электронных интегралов. С такой целевой функцией потеря времени, из-за несбалансированности узлов, составляла около 2% в тестовых расчётах.

Таблица 3. Время (в мин.) расчёта молекул параллельными прямым (Пр.) и полупрямым (ППр.) методами Хартри-Фока и ТФП в зависимости от числа узлов па кластере персональных компьютеров оснащённых Pentium III 800 MHz CPU и Ultra ATA 100 дисками.

Число узлов

Мол. Мет. Базис БФ 1 2 3 4 5

SÍ17H20 HF cc-pVDZ 406 Пр. 81.2 42.4 29.9 22.3 18.5

Кл. 28.3 15.3 11.1 8.8 7.0

Пр./К л. 2.9 2.8 2.7 2.5 2.6

Taxol SVWN 6-3 IG 660 Пр. 395.9 202.8 137.3 104.9 85.4

ППр. 262.5 145.4 98.0 75.4 60.8

Пр./ППр. 1.5 1.4 1.4 1.4 1.4

Yohimbine BLYP cc-pVDZ 494 Пр. 289.8 145.9 99.1 76.6 60.6

ППр. 163.9 88.3 59.6 45.6 36.8

Пр./ППр. 1.8 1.7 1.7 1.7 1.6

Porp hiñe B3LYP cc-pVTZ 916 Пр. 614.7 319.1 209.6 159.3 130.4

ППр. 305.1 182.0 124.8 92.9 77.7

Пр./ППр. 2.0 1.8 1.7 1.7 1.7

Время расчёта молекул параллельными методами Хартри-Фока и ТФП в зависимости от числа узлов на кластере персональных компьютеров представлено в Табл. 3. Эти результаты показывают, что разработанный классический и полупрямой метода ССП заметно опережают прямой метод ССП. Достигаемое улучшение здесь несколько меньше, представленного ранее. Главная причина заключается в меньшей скорости дисков на узловых ПК. В вычислениях представленных ранее использовался ПК с RAID массивом, тогда как узлы кластера имели одиночный Ultra ATA 100 диск.

В Главе 5 представлена эффективная реализация релятивистского двух-компонентного метода ТФП с использованием релятивистского эффективного остовного псевдопотенциала (RECP), а также расчёт потенциалов ионизации и энергий сродства к электрону атомов галогенидов от Вг до At и их двухатомных молекул, включая первые расчёты атома элемента (117) и его двухатомной молекулы, выполненные с разработанной программой.

Хороши известно, что квантово-механические расчёты молекулярных систем, содержащие атомы тяжёлых и сверх тяжёлых элементов, основанные на решение уравнения Шрёдингера в различных приближениях, имеют точность которая в значительной степени определяется тем, что уравнение Шрёдингера является нерелятивистским [18]. В нерелятивистских расчётах релятивистские эффекты можно учесть путём использования RECP [19, 20]. В этом подходе внутренние электроны представляются некоторым полулокальным потенциалом. Первоначально он был введён для рассмотрения только

валентных электронов в нерелятивистском случае [21, 22]. Для молекулярных систем образованных лёгкими атомами, где спин-орбитальным взаимодействием можно часто пренебречь, скалярно-релятивистские расчёты можно проводить в рамках обычных нерелятивистских вычислений, поскольку релятивистские эффекты частично включены в 11ЕСР. Для молекулярных систем, включающих атомы тяжелее талия, спин-орбитальным взаимодействием нельзя пренебрегать. Поэтому скалярно-релятивистские расчёты в этом случае становятся слишком неточными. Для включения спин-орбитального взаимодействия в И.ЕСР гамильтониан необходимо иметь спин зависимый ИЕСР, означающий что действия ИЕСР на орбитали правильной симметрии, например р\ и рз, является различным.

Для заданного спин-зависимого ИЕСР гамильтониана имеется ряд подходов включения спин-орбитального взаимодействия в квантово-механиче-ские расчёты. В случае много-конфигурационных волновых функций можно вводить спин-орбитальное взаимодействие при учёте конфигурационного взаимодействия с использованием орбиталей, полученных в скалярно-реля-тивистском расчёте. Этот подход успешно работает, если орбитальная релаксация, вызванная спин-орбитальным взаимодействием, не является критической. В противном случае возникают проблемы со сходимостью КВ разложения. Для одночастичных волновых функций методов Хартри-Фока и ТФП наиболее прямой подход к учёту спин-орбитального взаимодействия заключается в использование двухкомпонентного гамильтониана и спинорных ор-

где фш и ф{[] - компоненты </>; со спином а и /3 соответственно. На основе этого подхода в нашей работе была реализована программа для неэмпирических расчётов атомов и молекул релятивистским двухкомпонентньш методом Хартри-Фока и методами ТФП с использованием 11ЕСР.

Релятивистские 11ЕСР расчёты требуют вычисления спин-орбитальных интегралов и интегралов с эффективным псевдопотенциалом. Программа вычисляющая эти интегралы была разработана на основе алгоритмов, представленных в [23, 24]. Она вычисляет интегралы от декартовых гауссовых функций без ограничений на угловой момент и без ограничений на угловой момент операторов в полулокальных потенциалах. При этом используется эффективная предварительная оценка значений интегралов, позволяющая опускать вычисление интегралов у которых радиальная часть меньше заданной точности. В результате, программа демонстрирует высокую эффективность. Например, расчёт ИЕСР интегралов без учёта симметрии для 2*2*2 единичной ячейки

биталей:

РЬ, состоящей из 63 атомов, с LANL2DZ RECP и 504 базисными функциями

потребовал около 83 сек. на ПК с Pentium-4, 3.0 GHz CPU.

Скалярно-релятивистские и двухкомпонентные

расчёты, выполненные для ОСНОВНЫХ состояний ато- Таблица 4. Потенциал

MOB Вт, I, At, (117) ИХ катионов И анионов, ПОЗВО- ионизации (IP) и энергия

сродства к электрону лили оценить их потенциалы ионизации и энергии ^^ ^ ^ атома

сродства к электрону. Для атома элемента (117) они рассчиташ,ые двухкомпо-

представлены в Табл. 4. Расчёты с использованием нентными методами.

RECP отличаются лишь незначительно от реляти- Method

зультатов достигает порядка 1 эВ. Двухкомпонент-

Method IP EA

HF RECP 6.42 0.44

HF DK6 6.44 0.37

HF ZORA 6.45 0.37

BP86 RECP 7.56 1.65

BP86 DK8 7.62 1.59

BP86 ZORA 7.62 1.61

B3LYP RECP 7.52 1.53

B3LYP DK6 7.57 1.47

B3LYP ZORA 7.58 1.47

соответствие с экспериментальными значениями для атомов Вг, I, АЬ. Максимальное отклонение результатов, полученных с использованием RECP, составляет только 0.08 эВ.

Влияние спин-орбитальных эффектов на энергию сродства к электрону заметно уже для атома Вг. Они быстро растут с увеличением заряда ядра: для атомов Вг, I, Л< и (117) они составляют 0.17, 0.33, 0.77 и 1.44 эВ соответственно. Это увеличение пропорционально Я2, в соответствие с общим правилом для релятивистских эффектов.

Спектроскопические константы Яе, Ше, ШеХе И Д, были рассчитаны ДЛЯ ОС- Таблица 5. Спектроскопические постоян-новных состояние молекул Вг2, /2, и "ые: (в А);<".*« (в "') и (в эВ)

анализ показывает, что наблюдается хо-

(117)2 это согласие лучше для результатов, полученных двухкомпонентными, чем скалярно-релятивистскими методами. При этом для этой молекулы наблюдается полное согласие между результатами, полученными двухкомпонентным '

методом с использованием RECP, и результатами с учётом всех электронов,

Метод Яе

Скалярно-релятивистские расчёты

HF RECP 3.017 134 0.148 0.36

HF DI<6 3.058 130 0.130 0.39

BP86 RECP 3.048 122 0.180 1.83

BP86 DK6 3.082 119 0.203 1.82

B3LYP RECP 3.053 122 0.173 1.62

B3LYP DK6 3.088 120 0.206 1.62

2х-компоментные расчёты

HF RECP 3.487 78 0.092 -0.74

HF DK6 3.502 77 0.092 -0.73

BP86 RECP 3.526 68 0.121 0.72

BP86 DK6 3.523 68 0.134 0.70

B3LYP RECP 3.533 69 0.114 0.51

B3LYP DK6 3.533 69 0.131 0.49

¿^-компонентные расчёты (все электр.)

HF 3.500 77 0.089

ВР86 3.522 67 0.092

B3LYP 3.532 68 0.085

полученными с использованием гамильтониана Дугласа-Крола. Сравнение двухкомпонентных и четырёхкомпонентных методов показывает, что Ле и ше, рассчитанные этими методами, являются практически одинаковыми.

Анализ результатов, полученных скалярно-релятивистским и двухком-понентным методами, показывает, что учёт спин-орбитального взаимодействия ведёт к ослаблению связи. Это особенно заметно на примере молекулы (117)2, где спин-орбитальное взаимодействие наиболее сильное. Эту особенность можно объяснить на основе подхода, предложенного в [25]. Рассмотрим скалярно-релятивистские орбитали. Для тяжёлых систем, с большой энергетической щелью между валентными в — р орбиталями, можно отбросить орбитали образованные валентными я атомными орбиталями, которые являются несвязанными. Оставшиеся десять валентных электронов от р5 оболочки дают конфигурацию п*Пд(Тд в которой только ад имеет чисто связывающий характер. В первом порядке спин-орбитальное взаимодействие расщепляет п орбитали uajz = ±|и ]2 = ±| компоненты, но для молекул с замкнутой оболочкой это не влияет на полную энергию молекулы. В более высоких порядках спин-орбитальное взаимодействие смешивает орбитали с одинаковым^ и одинаковой чётности. Поэтому занятая компонента = орбитали тгд смешивается с занятой орбиталыо ад, а занятая компонента jг = ±5 орбитали 7Гц смешивается с незанятой анти-связывающей сги орбиталыо. Такое смешение ведёт к повышению энергии, увеличивающемуся с увеличением длины связи, поскольку разность энергий 7ги и аи орбиталей становится меньше, а это ведёт к ослаблению связи.

В Главе 6 представлена эффективная реализация многоссылочного метода конфигурационного взаимодействия (М11С1) с вычислением матричных элементов на основе правил Слэтера. Было показано, что метод М11С1 позволяет явно коррелировать намного больше число электронов по сравнению с другими методами, основанными на использование теории унитарной группы, и рассчитывать состояния высокой мультиплетности. Программа позволяет коррелировать до 100 электронов (по сравнению с 20 электронами для большинства других программ) на компьютере с 1 СВ оперативной памяти, что является непревзойдённым показателем до настоящего времени. Её возможности были продемонстрированы в исследовании взаимодействия молекулы С#2 с поверхностью N¿(100), показывающим механизм каталитической активности металлов. Поверхность N¿(100) в расчёте моделировалась кластером N¿30 с тремя слоями атомов, расположенными в узлах кристаллической решётки Ni при экспериментальных расстояниях между атомами и с вложенным псевдопотенцналом. Геометрия этого кластера с молекулой СН2 на поверхности представлена на рис. 7. Расчёт включал 79 гауссовых

«о

■ч» "О

Рис. 7. Молекула СЯг на поверхности кластера N130-

базисных функций. Наинизшем спиновым состоянием кластера в принятой модели является квинтетное состояние с четырьмя открытыми оболочками. Наинизшими электронными состояниями С#2 являются триплетное и синглетное состояния. Триплетное состояние возникает из конфигурации с двумя открытыми оболочками, а синглетное состояние соответствует конфигурации с замкнутыми оболочками. Отсюда можно заключить, что наинизшими состояниями системы МзцСЯг являются септетное и квинтетное состояние с шестью открытыми оболочками и квинтетное состояние с четырьмя открытыми оболочками. Полные энергии этих состояний в зависимости от расстояния между СН-2 и поверхностью кластера УУгзо представлены в Табл. 6.

Полученные результаты показывают, ЧТО уровни энергий ле- Таблица 6. Полные энергии (в ат. ед.) №(100) -жат в узком энергетическом диа- си> комплекса в зависимости от расстояния Я ,т , (в ат. ед.) между поверхностью .№(100) и С Н?.

пазоне. Наблюдается пересечение

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

квинтетное состояние с четырьмя открытыми оболочками. Оценки дают, что доля корреляционной энергии в энергии адсорбции СН2 на поверхности N¿(100), равной примерно 4 эВ, составляет порядка 40%. Отсюда следует, что аккуратный учёт корреляционной энергии крайне важен для получения корректных оценок энергий взаимодействия молекул с поверхностью.

Хорошо известно, что метод МИС1 сводится в общем случае к нахождение собственных значений и соответствующих им собственных векторов обобщённой матричной задачи на собственные значения с симметричными матрицами. Поэтому эффективность метода КВ во многом определяется эффективностью метода, использованного для нахождения собственных значений матриц. В связи с этим в наших работах значительное внимание было уделено развитию итерационных методов нахождения собственных значения и соответствующих им собственных векторов обобщённой задачи на собственные значения с симметричными матрицами. Всего было предложено более 20 новых методов, среди которых выделим методы типа Ньютона и Ньютона-

Септет Квинтет (1) Квинтет (2)

Я (6 откр. обол.) (6 откр. обол.) (4 откр. обол.)

2.1 -199.030335 - -199.040820

2.3 -199.025548 - -199.041518

2.5 -199.005630 -199.003016 -199.038337

2.7 -199.004732 -199.003440 -199.032332

2.9 -199.004060 -199.003040 -199.024326

30.0 -198.893150 -198.890327 -198.864568

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

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

<2)

Метод Ньютона-Рафсона для нахождения всех нахождения экстремальных значений этого функционала заключается в определение корректирующего вектора ЗХ для некоторого исходного вектора X из уравнения Ньютона-Рафсона Н6Х = -С7, где Я - матраца гессиана вторых производных иб - вектор градиента. Вычисляя градиент и гессиан и подставляя их в уравнение Нью-тона-Рафсона, получим

{А - р{Х)В - ВХв+ - С(ВХ)+)6Х = (3)

уравнение Ньютона-Рафсона для корректирующего вектора 6Х, являющееся неоднородной системой линейных уравнений. Прямая проверка показывает, что X является решением однородной системы линейных уравнений, получаемой из (3). Это означает, что матрица гессиана данной задачи является вырожденной. В этом случае система линейных уравнений (3) имеет решение, если X ортогональна (7, что легко проверяется. Таким образом уравнение (3) даёт корректирующий вектор ортогональный текущему вектору X.

Уравнение (3) позволяет построить новый итерационный метод Нью-тона-Рэлея для вычисления выбранных собственных значений и собственных векторов обобщённой задачи на собственные значения. Он основан на решение этого уравнения методом Галёркина в подпространстве Крылова, базис в котором определяется системой векторов, получаемой из уравнения (3). Ясно, что в этом случае скорость сходимости метода Ньютона-Рэлея будет выше квадратичной вблизи решения, поскольку несколько корректирующих векторов Ньютоновского типа используются для нахождения нового вектора. При использовании других векторов скорость сходимости может быть ниже.

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

приближение матриц А и В, а. также диагональное приближение матрицы гессиана, приводит к простейшей формуле корректирующего вектора

___—__(4)

Ац — рВи — ВцХ{С{ — GiBiiX¡

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

Методы типа Ньютона-Лагранжа строятся исходя из эквивалентности обобщённой задачи на собственные значения и задачи нахождения всех минимумов функционала Лагранжа Ь(Х, А) = Х1АХ - А(Х*ВХ -1). Применение метода Ньютона-Рафсона к этому функционалу приводит к уравнению для

¿X и 5А ч

/ А - \В -ВХ \ (6Х\ _ (АХ - АВХ\

V -Х*В 0 ) \зх) и -Х1ВХ)~

Полагая, что X является В ортонормальным вектором, из второго уравнения этой системы получаем Х*В5Х = 0. Затем, из первого уравнения системы (5) следует

(А - АВ)6Х - 5ХВХ = -(А - АВ)Х . (6)

Умножение этого уравнения слева на (А — А В)"1 и затем на Х1В даёт следующее выражение для 6\

ХгВ(А - АВ)-1(Л - АВ)Х ( )

Х'В{А-\В)-1ВХ ' И

Корректирующий вектор Лагранжа получается после этого из решения системы линейных уравнений (6).

Использование выражений (6) и (7) для построения итерационных методов требует широких исследований. В нашей работе для оценки А использовались собственные значения Ритца, а в вычисление (А — АВ)-1 использовалось диагональное приближение для матриц Ар и Во- Уравнения (6) и (7) с такими приближениями позволяют построить метод Ньютона-Лагранжа и его блочное обобщение. Применение диагонального приближения в (6) позволяет построить диагональный метод Ньютона-Лагранжа и его блочное обобщение.

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

В Главе 7 представлены работы по оптимизации малых базисов из сгруппированных функций гауссова типа для атомов первого, второго и третьего периодов. Использование гауссовых функций в качестве базисных функций в неэмпирических квантовомеханических расчётах атомов и молекул было предложено Бойсом [26]. В реальных расчётах более удобным оказалось применение сгруппированных гауссовых функций, предложенных в [27, 28]. В последующие годы много различных сгруппированных наборов гауссовых функций было разработано для использования в неэмпирических расчётах.

В наших работах для атомов первого периода были построены наборы базисных функций образующиеся из (9я5р) примитивных гауссовых функций. При этом использовалась двухшаговая оптимизация. Вначале были оптимизированы экспоненты примитивных функций, а затем - коэффициенты группировки. Были построены базисы: (9з5р/5.чЗр), (9з5р/5к2р), (9з5р/4з2р), (9.ч5р/3я2р) и дополнительно базисы (10з5р/4й2р) и (10,ч5р/3з2р), в которых одна из примитивных в функций включалась в две различные сгруппированные функции. Для атома водорода базисы были также оптимизированы.

Для атомов первого и вто-

Таблнца 7. Полные Хартри-Фоковские энергии (в ат. ед.) с обратным знаком атомов от Н до Аг в базисах 3-21С, 3-213Р, 4-31С и 4-22БР.

poro периодов нами были построены малые валентно-расщеплённые базисы 3-21sp и 4-22sp, в которых для представления 2s и 2р атомных орбиталей используются гауссовы функции с одинаковыми значениями экспоненциальных параметров, что даёт значительное сокращение времени вычисления двух-элек-тронных интегралов. Построение этих базисов было обусловлено тем, что известные базисы 3-21G и 4-31G являются не полностью оптимизированными. В них функции описывающие валентные и остовные электроны оптимизированы раздельно. В базисах, предложенными нами, эти функции оптимизировались одновременно. Поэтому, полные энергии атомов в базисах 3-21sp и 4-22sp, представленные в Табл. 7, заметно ниже, чем в известных базисах. Тестирование построенных базисов в молекулярных расчётах также показало преимущество новых базисов по сравнению с известными 3-21G и 4-31G базисами.

Ат. Coc. 3-21G 3-21SP 4-31G 4-22SP

Н 0.49698 0.49698 0.49928 0.49928

Не 2.83568 2.83568 2.85516

Li 2S 7.38151 7.42675 7.43044

Be 1s 14.48682 14.55522 14.56844

В 24.38963 24.47394 24.49425 24.51853

С 3 p 37.48039 37.59611 37.63595 37.66980

N *s 54.10366 54.25499 54.32529 54.36910

О 3 p 74.39178 74.58856 74.70244 74.75842

F гР 98.84423 99.08933 99.26418 99.33243

Ne ls 127.80382 128.10069 128.35621 128.43679

Na Js 160.85404 161.77118 161.76301

Mg »S 198.46810 199.49071 199.48815

Al 2P 240.55101 241.69673 241.75546

Si 3p 287.34439 288.60217 288.73294 288.69592

P 4S 339.00003 340.38713 340.32016 340.51945

S 3 p 395.44743 396.97036 397.05120 397.26417

C1 2P 457.19662 458.89843 458.97140 459.19406

Ar '5 524.34296 526.17748 526.47000

В настоящее время, среди валентно-расщеплённых базисов, базис 6-31G наиболее часто используется в неэмпирических расчётах молекул, образованных атомами первого, второго и третьего периодов, включая атомы переходных металлов. Все эти базисы были построены по одной методологии путём минимизации функционала энергии соответствующих атомов. Однако, сравнение электронного строения атомов первого и второго периодов с атомами третьего периода, включая атомы переходных металлов, показывает, что последние имеют отличия от атомов первого и второго периодов, которые не учитываются при оптимизации всех базисов по одной методологии. Дело в том, что в атомах переходных металлов открытая d оболочка, согласно схеме группировки базисных функций, должна описываться как валентная. С другой стороны, она является внутренней, поскольку лежит ниже уровня Ферми (ниже s оболочки). Для других атомов третьего периода она также является внутренней. Поэтому, для учёта этой особенности, построение базисов для атомов третьего периода, включая атомы переходных металлов, должно отличаться от построения базисов для атомов первых двух периодов.

В связи с этим, для учёта особен- „,„„_. ,

Таблица 8. Энергии ионизации (в эи) ато-ностей атомов третьего периода, нами мов переходных металлов рассчитанные

было предложено использовать двух- методом ТФП с B3LYP функционалом, шаговую оптимизацию для построения базисов, в которой функции, описывающие d оболочку, оптимизировались отдельно от других sup функций. Таким образом были построены модифицированные базисы m6-31G для атомов третьего периода, включая атомы переходных металлов, а также поляризационные функции. Атомные расчёты показали их значительное преимущество по сравнению с базисами 6-31G. Например, среднее отклонение в энергиях возбуждения атомов переходных металлов, рассчитанных методом ТФП с B3LYP функционалом, составило 0.17 эВ для базиса m6-31G и 1.66 эВ для базиса 6-31G, а среднее отклонение в энергиях возбуждения для катионов - 0.06 и 1.21 эВ, соответственно. Энергии ионизации атомов переходных металлов, рассчитанные с новым базисом и представленные в Табл. 8, являются существенно более точными, чем энергии полученные с базисом 6-31G.

Новые базисы атомов третьего периода были детально тестированы в расчётах молекул, включавших все атомы третьего периода. Для примера, в Табл. 9 даны результаты расчётов метиленовых комплексов переходных

Атом Ионизация 6-31G m6-31G Эксп.

Sc d}sг -> dLsL 6.50 6.57 6.56

Ti cPs2 -> dV 6.73 6.61 6.83

V dV -+ d3s' 6.94 7.06 7.06

Cr dV -> d5 6.20 6.98 6.76

Mn dV -y dV 7.31 7.47 7.43

Fe dV dV 7.68 7.89 7.90

Co dV dV 8.01 8.26 8.28

Ni dV -v dV 8.34 8.62 8.67

Cu d'V -> d10 6.89 7.85 7.72

Среднее откл. -0.29 0.01

Ср. абс. откл. 0.29 0.08

МеСН? Сост. 6-ЗЮ** тб-ЗШ** Эксп. О.

г(Ме-Н) О. г(Ме-Н) ое

БсСМ Ч 1.859 87.5 1.866 89.3 97.9 ± 5.3

ТгСЩ 1.851 91.5 1.848 88.8 92.9 ± 2.1

УСЩ *в2 1.877 85.0 1.876 80.8 79.9 ± 1.4

СгСН} 4в, 1.863 77.2 1.854 63.0 51.7 ± 0.9

МпСН? 1.883 73.8 1.885 74.6 70.3 ±2.1

РеСНТ 4В! 1.844 84.2 1.850 88.0 83.6 ± 1.8

СоСНТ 'А, 1.796 91.2 1.827 61.8 77.9 ± 1.2

тент 2А, 1.773 98.5 1.805 89.4 75.1 ± 0.9

СиСЩ м, 1.746 115.0 1.852 72.9 63.3 ± 1.2

Ср. абс. откл. 15.0 8.2

металлов МеСЩ, показы- Таблица 9. ВЗЬУР длины связей Ме-С (в А) и энергии вающие ЧТО среднее аб- разрыва связи (на Ме+ + 3СН2 в ккал/моль) для катионов метиленовых комплексов МеСН?.

солютное отклонение энергий разрыва связи, полученных с базисом тб-ЗЮ*, от экспериментальных почти в два раза меньше, чем отклонение полученное с 6-ЗЮ* базисом. Сравнение геометрических параметров 38 молекул, содержащих атомы третьего периода, рассчитанных методом МР2, показывает, что среднее абсолютное относительное отклонение, полученное с базисом тб-ЗШ**, равно 3.87, а с базисом 6-ЗЮ** - 4.47. В тестовом расчёте из 9 молекул, образованных только атомами третьего периода, соответствующие отклонения составили 3.94 для базиса тб-ЗЮ** и 5.08 для 6-ЗШ**.

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

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

нического анализа аппроксимирующего полинома. Оптимальный аппроксимирующий полином ищется в виде разложения в ряд по классическим ортогональным полиномам Рк{Н)

П1111

п 6 [5,ЛГ]

пин И5>ад,)-

г 6 [0,ЛГ] *=°

и(Щ\

(8)

где - точки, в которых заданы значения потенциала 11(11$, [О, ЛГ] -

последовательность целых чисел от 0 до N, [5, ЛГ] - такая же последовательность чисел от 5 до N. Выбор минимальной степени аппроксимирующего полинома, равной 5, объясняется тем, что коэффициент Данхэма Уго, во втором порядке теории возмущений, выражается через производную четвёртой степени, а для правильной передачи её зависимости от Е необходимо увеличить степень аппроксимирующего полинома по крайней мере до 5.

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

№> 1Чь-г1 - к = 3,...,п . (9)

Соотношения (8) и (9) полностью определяют оптимальный аппроксимирующий полином для квази-эрмитовой аппроксимации функции заданной таблицей её значений.

Предложенный нами метод квазиэрмитовой аппроксимации был ИСПОЛЬ- Таблица 10. Коэф. разложения аппроксими-зован в разработке метода Данхэма РУ^го полинома в ряд по многочленам

, , тт Чебышёва второго рода 8, 9 и 10 степени,

вычисления коэффициентов Данхэма

и колебательно-вращательных уровней энергий. При этом оптимальный аппроксимирующий полином может искаться в виде ряда по многочленам Чебышева первого и второго рода и многочленов Лежандра. Работу алгоритма рассмотрим на примере аппрок-

Пол. 8 9 10

Ьо -0.78998+1 -0.78998+1 -0.78998+1

1>1 -0.13661-2 -0.13661-2 -0.13661-2

Ьг 0.24192-2 0.24192-2 0.24192-2

ь3 -0.36195-3 -0.36195-3 -0.36195-3

Ь4 0.37820-4 0.37820-4 0.37815-4

Ь5 -0.33386-5 -0.33388-5 -0.33388-5

¡>6 0.29571-6 0.29571-6 0.29570-6

Ьт -0.26372-7 -0.25604-7 -0.25604-7

Ьв 0.15138-8 0.15138-8 0.54131-8

Ьд 0.49092-9 0.49092-9

¡>10 0.17337-8

Ср. кв. ош. 0.48-10 0.45-10 0.35-10

состояния молекулы ЫН. Коэффициенты разложения полиномов 8, 9 и 10 степени в ряд по полиномам Чебышёва второго рода, полученные при приближении исходных точек методом наименьших квадратов, даны в Табл. 10. Среднеквадратичная ошибка для этих полиномов уменьшается с увеличением степени полинома, т. е. соглаено критерию (8) наилучшим является полином десятой степени. Однако из Табл.

10 видно, что для полиномов 8 и 9 степени условие (9) выполняется, а для полинома 10 степени - нет. Поэтому совместно условия (8) и (9) дают, что оптимальным аппроксимирующим полиномом является полином 9 степени.

Производные от оптимального аппроксимирующего полинома вычислялись по рекуррентным формулам. При этом для многочленов Лежандра были выведены новые рекуррентные формулы. Высокая степень оптимального аппроксимирующего полинома позволяет вычислять большое число коэффициентов Ук1. Например, в случае ЫН, оптимальный полином 9 степени позволил вычислить 25 коэффициентов Уы в 14 порядке теории возмущений. Знание большого числа коэффициентов Уы даёт возможность вычислять энергии колебательно-вращательных уровней энергий с высокой точностью. В Табл. 11 даны уровни энергий основного СОСТОЯНИЯ N2 ДЛЯ Таблица И. Колебательно-вращательные уровни энергии 3 = О И 3 = 10 вы- (в см_1) N1 для .1 = 0 и J — 10, вычисленные мето-

' гг дом Данхэма, ошибки (<5) с которыми они были вычислены численные методом Дан" и энергии рассчитанные конечно-разностным методом.

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

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

Часть II. Расчёты волновых функций малых и больших молекул и их исследование

В Главе 9 представлены расчёты и исследования волновых функций, которые были в основном нацелены на достижение предельной точности вычислений.

Так для #2 была поставлена задача получения численно точного решения уравнения Хартри-Фока в приближение Л К АО. Исследование показало, что МО уравнения Хартри-Фока и МО этого же уравнения в приближение ЛКАО отличаются, поскольку обычное приближение ЛКАО неодинаково описывает асимптотические свойства МО в пределах объединённого и разде-

Пор. ТВ 0 2 4 6-14 6 Числ.

J = 0

С(0) 1363.78 1361.37 1361.37 1361.37 0.00 1361.37

С(1) 4091.35 4067.09 4087.13 4067.13 0.00 4067.13

С (2) 6818.91 6750.99 6751.15 6751.16 0.00 6751.17

С(3) 9546.48 9413.04 9413.51 9413.53 0.01 9413.54

С(4) 12274.04 12053.26 12054.27 12054.31 0.04 12054.35

С(5) 15001.61 14671.64 14673.50 14673.57 0.08 14673.68

./ = 10

С(0) 1596.79 1593.56 1593.56 1593.56 1593.56

С(1) 4324.36 4297.77 4297.81 4297.81 4297.81

С(2) 7051.92 6980.14 6980.32 6980.33 6980.33

С(3) 9779.49 9640.68 9641.18 9641.19 9641.21

С(4) 12507.05 12279.38 12280.43 12280.47 12280.51

С(5) 15234.62 14896.25 14898.16 14898.23 14898.32

лённых атомов. В связи с этим нами было предложено модифицированное приближение ЛКАО, в котором асимптотические свойства МО описываются корректно. Использование этого приближения позволило вычислить Хартри-Фоковскую энергию Яг с 12 точными числами (Табл. 12), т. е. получить численно точное решение. Предыдущая точность была увеличена на 6 порядков.

Были проведены также неэмпирические расчёты состояний молекул 5гг и Ва,2, в которых получена наиболее точная оценка диссоци-онной энергии Ва,2 1090±30 см-1. Молекулы Л^Ог и А^Оз были рассчитаны с целью получения наиболее точного инфракрасного и Рама-новского спектров этих молекул, а также с целью исследования относительной энергетической стабильности основных конформеров этих молекул.

Молекулы Нв2 и Вег часто используются для тестирования наиболее совершенных методов учёта корреляционной энергии, поскольку вычислительная сложность таких методов быстро растёт с увеличением числа электронов. Расчёты основных состояний молекул Яег и Ве^, выполненных нами, преследовали две цели. Во первых, они были нацелены на достижение максимально высокой точности, и, во-вторых, на демонстрацию того, что высокоточные расчёты могут выполняться также методом МИС1. При этом, для расчёта В&2 применимость метода связанных кластеров является проблематичной. Дополнительный аргумент в пользу выбора метода М11С1 связан со структурой волновых функций слабо связанных состояний Не г и Ве 2- Хартри-Фоковкие однодетерминантные волновые функции Яег и Ве2 дают потенциальные кривые не имеющие связанных состояний. Поэтому возбуждённые КФС ответственны за образование химической связи в этих молекулах. В связи с этим возможность включения в список ссылочных конфигураций максимального числа конфигураций является важным фактором предопределяющим качество рассчитываемой волновой функции. По этому критерию М11С1 метод имеет заметное преимущество перед другими известными методами. Метод М11С1 не использовался ранее в расчётах Не2- Поэтому сравнение этого метода с другими высокоточными методами, на примере расчёта молекулы Яег, стало дополнительной целью работы.

Потенциальная кривая Х1^ Нег была рассчитана методом МР1С1 с ис-

Таблнца 12. Полная Хартри-Фоковская и орбитальная энергии (в ат. ед.) П2 при Я=1.4 ат. ед.

Молекулярный базис

Полная энер. Орб. энер.

Н

Не

218+203-1-205 218+205+203 228+215+215 235+225+225 235+228+225 245+235+238 258+245+245 268+253+255

258+245+245 268+255+258 268+258+258 268+255+255 275+268+265 278+268+268 278+26в+2б8 278+268+265

-1.133629571315 -1.133629571316 -1.133629571387 -1.133629571420 -1.133629571420 -1.133629571434 -1.133629571448 -1.133629571456

-0.5946585686 -0.5946585686 -0.5946585686 -0.5946585686 -0.5946585686 -0.5946585686 -0.5946585687 -0.5946585687

пользованием корреляционно-согласованных базисов сс-р'У^ и сс-р\^, дополненных одной, двумя и тремя присоединёнными функциями. Для улучшения сходимости ряда по КФС, волновые функции МЯС1 метода вычислялись с использованием псевдо-натуральных орбиталей, рассчитываемых в предварительных М11С1 вычислениях. Список ссылочных конфигураций включал примерно 1.5 * 103 конфигураций, описывающих возбуждения электронов на 2в, 2р, Зв и 3р орбитали Не. Все одно- и двукратно возбуждённые конфигурации, образующиеся из списка ссылочных конфигураций, включались в М11С1 волновую функцию. При этом размер конфигурационного пространства был равен примерно 4.6 * 106. Полученные МЯС1 полные энергии вначале были скорректированы, чтобы исключить базисную суперпозиционную ошибку, используя метод Бойза-Бернарди. После этого они были экстраполированы на базис бесконечного размера. Влияние релятивистских эффектов на потенциальную кривую было оценено в МЛС1 расчётах с использованием гамильтониана Дугласа-Крола.

Потенциальная кривая Не2 была рассчитана в интервале 2.0-80.0 ат. ед. в 62 точках. Значения полученного нерелятивистского и квази-релятивист-ского потенциала в некоторых точках даны в Таблице 13. Нерелятивистское

Таблица 13. Нерелятивистский и квази-релятивистский потенциал Нс2 (в К).

л (ат. ед.) МЯС1 н. р. МШЛ+рел. н. р. БАРТ (291 ЭАРТ+БМ |291 ССЗО(Т)+ЕС1 |30] МЯ-АСРЕ |31|. ЦМС [32]

3.0 3767.5035 3766.7816 3766.97 3767.681 - 3768.046 -

3.5 1110.6635 1110.4382 1110.38 1110.649 1112.430 1110.878 -

4.0 292.5836 292.5214 292.49 292.570 293.498 292.750 292.60

4.5 58.4157 58.4023 58.371 58.397 58.759 58.488 -

5.0 -0.4943 -0.4945 -0.4781 -0.4754 -0.395 -0.429 -

5.3 -9.1744 -9.1731 -9.1659 -9.1681 - - -

5.4 -10.2683 -10.2668 - - -10.255 -10.235 -10.998

5.6 -11.0023 -11.0004 -10.9996 -11.0037 -11.004 -10.980 -10.998

6.0 -9.6798 -9.6779 -9.6791 -9.6819 -9.686 -9.669 -

6.5 -6.8900 -6.8916 -6.8924 -6.89359 - -6.887 -

7.0 -4.6184 -4.6173 -4.6225 -4.62277 - -4.620 -

7.5 -3.0711 -3.0703 - - -3.068 -3.073 -3.073

8.0 -2.0647 -2.0642 -2.0669 -2.06669 - -2.066 -

9.0 -0.98886 -0.98906 -0.98984 -0.98977 -0.988 -0.990 -0.989

10.0 -0.51205 -0.51171 -0.51259 - -0.512 -0.513 -

12.0 -0.16657 -0.16679 -0.16592 - -0.165 -0.165 -0.166

15.0 -0.04334 -0.04272 _ -0.042 -0.043 -0.0423

Ае=5.607 ат. ед. было определено как точка минимума оптимального аппроксимирующего полинома восьмой степени. Коэффициенты этого полинома были получены по 14 точкам потенциальной кривой, используя методы разработанные нами и представленные ранее. Энергии диссоциации были оценены как разности полных энергий в точке минимума и при Я = 80.0 ат. ед. Значения Ое, полученные в нерелятивистском и квази-релятивистском расчётах, равны 11.0031 К и 11.0002 К, соответственно. Небольшое уменьшение Д, в

квази-релятивистском расчёте связано с тем, что волновая функция становится более компактной, что приводит к уменьшению поляризации атомов и, следовательно, £>е. Полученные результаты очень хорошо согласуются с литературными данными. Это позволяет заключить, что метод МГ1С1 можно успешно использовать для проведения высокоточных квантово-механических расчётов молекул.

Расчёты Вв2, выполненные аналогично Яег, привели к ряду важных результатов. Было показано, что особенность связи в Вв2 состоит в том, что на нижних четырёх колебательных уровнях в молекуле реализуется ковалент-ный тип связи, тогда как на более высоких уровнях связь образуется из-за поляризационного или дисперсионного взаимодействия. Впервые это было получено нами из анализа возбуждённых конфигураций, а экспериментально наличие двух типов взаимодействий на разных колебательных уровнях было получено позже (см. рис. 8). Полные электронные плотности Вв2 (рис. 9), рассчитанные при Я = Яе и Я, соответствующему пятому колебательно-

250

у+1/2

Рис. 9. ТФП/ВЗЬУР/сс-рУС^ электр. плот. Ве2: Рис. 8. Ве2. Эксп. ДС„+1/2 в (а) - Я=2.5 А, (б) - Я=4.0 А. Рис. (а) плот. - 0.025

зависимости от и [33]. е/А3; рис. (б) плот. - 0.025 е/А3 (а) и 0.005 е/А3 (б).

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

Исследование квази-релятивистских эффектов показало, что они нелинейным образом зависят от геометрии молекул (рис. 10 и 11), оказывая тем самым заметное влияние на значения колебательных уровней энергий. Кроме этого, проявление релятивистских эффектов зависит от типа химической связи (см. излом на кривой Е(МНС1)-Е(БК-МГ{,С1) на рис. 11). Расчёт потенциальной кривой Ве2, Ое, колебательных уровней энергий и сравнение их с экспериментальными значениями (см. Табл. 14) показывает, что расчёт Веч методом М11С1 был выполнен на пределе возможной точности неэмпирических

Рис. 10. Экстр. БЭвЕ потенц. Х1Т,+ Ве2 (см-1, А), Е(экстр)-Е(1;-а1Щ-сс-рУ57), Е(МКС1)-Е(БК-]УтС1).

-Е(экстр)

-Е(экстр) Е(|ШС|)

-Е(М!!С|) ЕрК-ШС!)

Рис. 11. Область 2.5-5.0 А рис. 10 в увеличенном масштабе.

расчётов, учитывая что в [39] использовался больший базис, а полученный потенциал является композиционным. Отклонение Ое от экспериментального составило 1.50 и 2.44 см~: для нерелятивистского и квази-релятивистско-го расчёта, а средние абсолютные ошибки в колебательных уровнях энергии находятся на уровне 4 см-1 для обоих расчётов. При этом наибольшие отклонения рассчитанных уровней энергий от экспериментальных наблюдается в области смены характера связи в Вв2: что вполне объяснимо. Расчёты показывают возможность существования также колебательного уровня с V = 11 в двух вращательных состояниях.

В Главе 10 представлены результаты полученные при решении задачи определения активных центров гем-цитохрома с (БТСДМКЗ) в реакции переноса электрона, возникающей в процессе биологической активности бактерии БкетапеНа one.idev.sis, и приводящей к окислению различных внешних акцепторов электронов, среди которых выделяются соединения тяжёлых металлов и актиноидов. На этом основана идея использования бактерий вкегиапеПа опег<1еп818 для биологической очистки вод от тяжёлых металлов, когда соединения в более высокой степени окисления являются менее растворимыми, например соединения и(У1) и и(1У). Данная работа стали возможными благодаря прогрессу, достигнутому в нами по разработке новых эффективных методов решения уравнения Хартри-Фока со свойствами линейной масштабируемости, которые открыли возможность прямого неэмпирического исследования больших молекулярных систем.

БТС образован 91 аминокислотой и четырьмя гем-группами и включает более 1500 атомов. Первые неэмпирические расчёты вТС и анализ волновых функций показал наличие неизвестного ранее физического явления и необыч-

Таблица 14. Колебательные уровни энергий Не? (в см -1), Я,. . в ат. ед., и в СМ

ММЛ БК- ммя БК-МИС1 МВ.С1 МЛ Ех1 Комп ЭАРТ СС+11 Эксп

ммл -Эксп -Эксп АСРР Сет мет. +Согг РС1

н. р. н. р. н. р. н. р. [34] |35] 136] [37] [38] [39| [33]

К 2.4411 2.4415 4.627 4.617 4.633 2.44 2.444 2.4536

А- 931.24 927.30 893 898 952.0 922.9 938.7 935 929.74

803.60 799.88 768.2 776.8 826.2 795.0 (812.4) 807.4

С(1)-С(0) 224.8 224.2 2.2 1.6 218.4 218.6 221.7 223.3 222.3 222.7 222.6

С(2)-С(0) 402.1 400.8 5.0 3.7 387.0 389.8 396.7 396.4 397.6 396.8 397.1

С(3)-С(0) 526.4 524.3 8.3 6.2 499.0 506.2 519.8 515.9 520.3 517.8 518.1

С(4)-С(0) 604.7 602.0 9.9 7.2 568.0 579.6 600.6 592.4 597.9 594.7 594.8

С(5)-С(0) 660.5 657.6 9.0 6.1 622.4 634.8 660.0 649.1 655.1 651.6 651.5

С(6)-С(0) 704.7 701.6 5.9 2.8 706.1 680.2 709.0 695.6 702.6 698.9 698.8

С(7)-С(0) 741.6 738.4 3.9 0.7 706.1 716.6 749.7 733.5 741.7 738.0 737.7

С(8)-С(0) 770.1 766.7 1.9 -1.5 734.8 744.3 782.0 762.9 772.4 768.6 768.2

С(9)-С(0) 789.6 786.1 -0.3 -3.8 754.2 763.4 805.3 783.1 794.3 790.4 789.9

С(10)-С(0) 800.2 796.5 -2.4 -6.1 764.6 773.6 793.7 807.1 803.1 802.6

С(11)-С(0) 803.4 799.7 768.0 776.7 811.9

Ср. аб. ош. 4.4 3.6 21.4 15.8 6.9 3.3

Ср. кв. ош. 5.5 4.3 25.1 18.2 8.9 4.2 3.4 0.3

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

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

Детально этот эффект и его влияние на относительную стабильность конформеров молекул исследовался на примере молекул 1УМ2 и 1СК\У. Было показано, что не учёт этого эффекта в методах силовых полей даёт неправильную относительную энергетическую стабильность различных конформеров. Его проявление в биомолекуле можно наблюдать сравнивая электростатический потенциал молекулы, рассчитанный неэмпирическим методом, учитывающим этот эффект естественным образом, с электростатическим потенциалом, полученным методом силовых полей, который этот эффект не учитывает (см. рис. 12). Можно заметить, что найденный эффект выделяет биологические молекулы из класса органических молекул, поскольку наблюдается в молекулах с внутренней структурой. Он указывает также на новые возможности создания материалов с необычными свойствами.

Исследование электронных свойств гем-группы - \MePIm2] привело к построению модели химической связи в этих комплексах, необходимой для по-

Рис. 12. 1VM2, 4-й конф. Электростатические потенциалы: (а) - с натуральными зарядами атомов, (6) - с зарядами атомов силового поля Amber 99 (б).

нимания, каким образом гем-группа интегрируется в тело цитохрома. Решение поставленной задачи осложнялось проблемой нарушения принципа заполнения орбиталей в методе ТФП, представленной ранее, которая не позволяет провести детальное исследование орбитальной структуры верхних заполненных орбиталей комплекса [МеР/тг] и его катионов. Однако, качественно решение задачи было получено из анализа заселённости волновых функций ТФП. Было найдено, что верхние занятые орбитали комплекса, его катиона и двойного катиона образованы атомными орбиталями норфиринового кольца. Отсюда следует, что атом переходного металла находится во втором окислительном состояние как в комплексе, так и в его катионах, а порфириновое кольцо может выполнять роль как акцептора, так и донора электронов. Это ясно следует из рис. 13, где представлены электростатические потенциалы

(а) (б) (в)

Рис. 13. Электростатические потенциалы комплекса \FePIm2] (а), его катиона (б) и дважды заряженного катиона (в), рассчитанные с натуральными зарядами на атомах.

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

Расчёт нейтрального БТС и его анионов в восстановленной и окисленной формах был выполнен методом Хартри-Фока в базисе 3-21ЭР и включал

Таблица 15. Полные суммы Малликеновских частичных зарядов положительно (<У+) и отрицательно заряженных аминокислот, а также заряды атомов Ре в гем-группах в восстановленной и окисленной формах ЭТС.

[ЙТС]" Восстановленная ЭТО [5ТС]3- [£.'ТС]5" [5ТС] Окисленная ЭТС [ЯТСр- [ЯТС]4'

0+ 1.5527 1.5422 1.3725 1.4316 1.4320 1.3810

ег -6.8673 -8.6365 -9.7966 -5.9609 -7.7340 -9.3671

Ре / Гем I 2.3374 2.3369 2.3364 2.3545 2.3537 2.3529

Ре / Гсм 111 2.1653 2.1725 2.0911 2.3603 2.3601 2.3599

Ре / Гем IV 2.3372 2.3370 2.3385 2.0357 2.0339 2.0322

Ре / Гем 11 2.1099 2.1186 2.0886 2.0562 2.0726 2.0665

около 8900 сгруппированных функций. Полученные волновые функции были детально проанализированы в терминах Малликеновских частичных зарядов аминокислот и гем-групп. Анализ показал, что аминокислоты в БТС взаимно поляризованы. Полные суммы частичных зарядов положительно и отрицательно заряженных аминокислот, а также заряды Ее в гем-группах в восстановленной и окисленной формах вТС даны в Табл. 15. При этом нейтральные аминокислоты не учитывались. поскольку они сохраняют свою нейтральность. Из представленных результатов следует, что верхние занятые ор-битали БТС, в восстановленной и окисленной формах, образованы орбита-лями атомов отрицательно заряженных аминокислот. Поэтому, они же являются активными центрами в реакции переноса электрона от вТС к внешнему акцептору электронов. Детальный анализ зарядов аминокислот указал также на наиболее вероятные аминокислоты и их атомы, которые могут быть активными центрами. Распределение зарядов в ЭТС можно понять, рассматривая её электростатический потенциал, пример которого представлен на рис. 14.

Рис. 14. Электростатический потенциал (а - спереди, б -сзади) нейтральной окисленной формы ЭТС рассчитанный с Малликеновскими частичными зарядами на атомах.

Основные результаты и выводы.

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

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

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

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

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

— эффективный метод параллелизации программ метода Хартри-Фока и ТФП.

— эффективную реализацию двухкомпонентной релятивистской программы ТФП с использованием ЯЕСР. Первый расчёт атома элемента (117) и его двухатомной молекулы.

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

— более 20 новых методов решения обобщённой задачи на собственные значения, включая методы типа Ньютона и Ньютона-Лагран-жа, их диагональные версии, а также их ленточные и блочные

обобщения.

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

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

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

3. Проведены высокоточные неэмпирические расчёты малых молекул, демонстрирующие расширенные возможности неэмпирических методов, в которых:

— получено численно точное решение уравнения Хартри-Фока для молекулы #2 с предложенным модифицированным приближением JIKAO корректно описывающим корреляционные свойства МО.

— методом КВ, используя оптимальную методологию, рассчитаны Х1Е+ потенциальные кривые Яег и Be-i на пределе точности неэмпирических расчётов: диссоционная энергия Вег отличается от экспериментальной менее чем 1.5 см-1, колебательные уровни энергии вычислены со среднеквадратичным отклонением порядка 4 см-1. Предсказано существование в состоянии Вс-2 колебательного уровня энергии с и = 11 в двух вращательных состояниях.

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

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

— рассчитаны наиболее точные ИК и Рамановский спектры молекул AI2O2 и AI2O3.

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

Основные публикации по теме диссертации

1. Митин А. В. Линейная экстраполяция в итерационном методе решения систем уравнений // Журнал вычислительной математики и математической физики. 1985. Т. 25, № 3. С. 325-334.

2. Mitin А. V. The Dynamical "Level Shift"Method for Improving the Convergence of the SCF Procedure // J. Comput. Chem. 1988. V. 9, N 2. P. 107 110.

3. Митин А. В. Метод псевдопотенциала для улучшения сходимости процедуры ССП // Журн. физ. химии. 1990. Т. 64, № 1. С. 219-222.

4. Митин А. В. Вычисление коэффициентов Yki двухатомной молекулы методом Данхэма по данным неэмпирических расчётов // Журн. физ. химии. 1990. Т. 64, № 7. С. 1837-1843.

5. Митин А. В. Вычисление колебательно-вращательных уровней энергий двухатомных молекул методом Данхэма по данным неэмпирических расчётов // Журн. физ. химии. 1990. Т. 64, № 8. С. 2041-2044.

6. Митин А. В. Оценка ошибок при вычислении коэффициентов Yki и колебательно-вращательных уровней энергий двухатомных молекул методом Данхема по данным неэмпирических расчётов // Журн. физ. химии. 1990. Т. 64, № 8. С. 2231-2233.

7. Митин А. В. Табулирование потенциальной кривой двухатомной молекулы в неэмпирическом расчёте // Журн. физ. химии. 1991. Т. 65, № 3. С. 716-721.

8. Митин А. В. Приближение ЛКАО в методе молекулярных орбиталей: влияние граничных условий // Журн. физ. химии. 1992. Т. 66, № 7. С. 1828-1835.

9. Mitin А. V., Hirsch G. Linear extrapolation in iterative methods // J. Math. Chem. 1994. V. 15, N 1. P. 109-113.

10. Mitin A. V. Iterative Methods for the Calculation of a Few of the Lowest Eigenvalues and Corresponding Eigenvectors of the AX = ABX Equation with Real Symmetric Matrices of Large Dimension // J. Comput. Chem. 1994. V. 15, N 7. P. 747-751.

11. Mitin A. V., Hirsch G., Buenker R. J. Accurate Atomic Gaussian Basis Functions for First-Row Atoms. Part I. Contracted Basis Sets Derived from 9s5p Primitives // J. Mol. Struct. (THEOCHEM). 1996. V. 362, N 3. P. 283-296.

12. Mitin A. V., Hirsch G., Buenker R. J. Accurate Small Split-Valence 3-21SP and 4-22SP Basis Sets for First-Row Atoms // Chem. Phys. Lett. 1996. V. 259, N 1-2. P. 151-158.

13. Mitin A. V., Hirsch G., Buenker R. J. Accurate Atomic Gaussian Basis Functions for Second-Row Atoms. I. Small Split-Valence 3-21SP and 4-22SP Basis Sets // J. Comput. Chem. 1997. V. 18, N 9. P. 1200-1210.

14. Mitin A. V. Calculation of the Rovibrational Energy Levels of Diatomic Molecules by the Dunham Method with a Potential Obtained from Ab Initio Calculations // J. Comput. Chem. 1998. V. 19, N 1. P. 94-101.

15. Mitin A. V. The Use of Symmetric Rank-One Hessian Update in Molecular Geometry Optimization // J. Comput. Chem. 1998. V. 19, N 16. P. 1877-1886.

16. Mitin A. V. New Methods for Calculations of the Lowest Eigenvalues of the Real Symmetric Generalized Eigenvalue Problem // J. Comput. Phys. 2000. V. 161, N 2. P. 653-667.

17. Mitin A. V. Exact Solution of the Hartree-Fock Equation for H2 Molecule in the Linear-Combination-of-Atomic-Orbital Approximation // Phys. Rev. A. 2000. V. 62, N 1. P. 010501 (R).

18. Mitin A. V. Large scale Hartree-Fock calculations. Compression of two-electron integrals and their indices // J. Molec. Struct. (Theochem). 2002. V. 592, N 1-3. P. 115-121.

19. Mitin A. V., Baker J., Wolinski K., Pulay P. Parallel Stored-Integral and Semidirect Hartree Fock and DFT Methods with Data Compression // J. Comput. Chem. 2003. V. 24, N 2. P. 154-160.

20. Mitin A. V., Baker J., Pulay P. An Improved 6-31G* Basis Set for First-Row Transition Metals //J. Chem. Phys. 2003. V. 118, N 17. P. 7775-7782.

21. Mitin A. V., Wüllen C. V. Two-Component Relativistic Density Functional Calculations of the Dimers of the Halogens from Bromine through Element 117 Using Effective Core Potential and All-Electron Methods // J. Chem. Phys. 2006. V. 124, N 6. P. 064305.

22. Mitin A. V., Merz K. M. An Improved 6-31G* Basis Set for Atoms Ga through Kr // Int. J. Quantum Chem. 2007. V. 107, N 15. P. 3028-3038.

23. Mitin A. V., Merz K. M. Erratum: An Improved 6-31G* Basis Set for Atoms Ga through Kr // Int. J. Quantum Chem. 2009. V. 109, N 5. P. 1158.

24. Mitin A. V. Performance of Multi-Reference CI Method in Calculations of Weakly Bonded Sr2, and Ba2 Molecules // Жури. физ. химии A. 2009. Т. 83, № 7. С. 1310-1314.

25. Mitin А. V., Kubicki J. D. Quantum Mechanical Investigations of Heme Structure and Vibrational Spectra: Effects of Conformation, Oxidation State, and Electric Field // Langmuir. 2009. V. 25, N 1. P. 548-554.

26. Mitin A. V. Large Scale Hartree-Fock Calculations with Conventional SCF Algorithm. Influence of Integral and Index Compression on Fock Matrix Construction // Жури. физ. химии. 2010. Т. 84, № 5. С. 912 919.

27. Mitin А. V. Accurate calculations of dissociation energies of weakly bonded He2 and Be2 molecules by MRCI method // Russ. J. Phys. Chem. A. 2010. V. 84, N 13. P. 2314-2319.

28. Митин А. В. Поляризация аминокислот и их взаимодействие в биомолекулах // Письма в ЖЭТФ. 2010. Т. 92, № 5. С. 398-402.

29. Mitin А. V. Ab initio calculations of weakly bonded He2 and Be2 molecules by MRCI method with pseudo-natural molecular orbitals // Int. J. Quantum Chem. 2011. V. Ill, N 11. P. 2560-2567.

30. Mitin A. V., Kubicki J. D., Merz К. M. Electronic Structure, Chemical Bonding, and Oxidation Numbers of Transition Metals in [MePIm2] Complexes and their Cations // Int. J. Quantum Chem. 2011. V. Ill, N 14. P. 3630-3642.

31. Mitin A. V. Lagrange Type Iterative Methods for Calculation of Eigenvalues of Generalized Eigenvalue Problem with Large Symmetric Matrices // Int. J. Quantum Chem. 2011. V. Ill, N 11. P. 2545-2554.

32. Mitin A. V. Effect of amino acid polarization in force field biomolecular calculations // Int. J. Quantum Chem. 2011. V. Ill, N 11. P. 2555-2559.

33. Dorofeeva О. V., Mitin A. V., Altova E. P. et al. Anomeric Effect in N-Azidomethylpyrrolidine: Gas-Phase Electron Diffraction and Theoretical Study // Phys. Chem. Chem. Phys. 2011. V. 13, N 4. P. 1490 1498.

34. Mitin A. V. Accurate theoretical IR and Raman spectrum of A1202 and A1203 molecules // Struct. Chem. 2011. V. 22, N 2. P. 411-418.

35. Paulus В., Mitin A. An Incremental Method for the Calculation of the Electron Correlation Energy in Metals // Recent Progress in Ccomputational

Sciences andD Engineering, Ed. by S. Т., M. G. Lecture Series on Computer and Computational Sciences, Vol. 7A-B. Leiden, The Netherlands: Brill Academic Publishers, 2006. P. 935-942. ISBN: 978-90-04-15542-8.

36. Mitin A. V. Accurate Dissociation Energies of He2 and Be2 Calculated by MR-CI Method // APS Conference Proceeding, Ed. by Т. E. Simos, G. Maroulis. Melville: APS, 2007. P. 231-234. ISBN: 10.1063/1.2836048.

37. Mitin A. V., Buenker R. J. Ab Initio Calculation of Adsorption Energy of the CH2 Molecule on the Ni Surface. Preliminary Results // in: Fourth European Workshop, Quantum Systems in Chemistry and Physics. Marly-Le-Roi, France: INJEP, 1999. P. B7.

38. Mitin A. V. Linear Scaling Orbital Based Super Direct SCF Method // in: Swedish Theoretical Chemistry 2006 (Svensk Teoretisk Kemi 2006), 4-5 May, 2006. Stockholm, Sweden: 2006. P. 28.

39. Mitin A. V. Linear Scaling Properties of One- and Two-Electron Integrals, Fock Operator, and Molecular Orbitals // in: 48th Sanibel Symposium, February 21-26, 2008. St. Simons Island, Georgia, USA: 2008. P. 2-21.

40. Mitin A. V. The Origin of Linear Scaling Phenomenon in Fock Matrix Calculations // in: Thirteen International Workshop on Quantum Systems in Chemistry and Physics, July 6-12, 2008. Michigan State University, East Lancing, MI 48824, USA: 2008. P. 30.

41. Митин А. В. Биологическая переработка соединений тяжёлых элементов. Теоретическое исследование реакционных центров СТЦ в реакции переноса электронов инициированной биологической активностью бактерии Shewanella oneidensis // Всероссийская конференция «Радиохимия - наука настоящего и будущего». Сборник тезисов. Москва, МГУ: 13-15 Апреля 2011. С. 110-111.

Цитированная литература

1. Clementi Е., Corongiu G., Detrich J. Н. et al. Parallelism in computational chemistry: Applications in quantum and statistical mechanics // Physica. 1985. V. 131B, N 1-3. P. 74-102.

2. Challacombe M., Schwegler E. Linear scaling computation of the Fock matrix // J. Chem. Phys. 1997. V. 106, N 13. P. 5526-5536.

3. Pepper M., Bürsten B. E. The electronic structure of actinide-containing molecules: a challenge to applied quantum chemistry // Chem. Rev. 1991. V. 91, N 5. P. 719-741.

4. Relativistic Electronic Structure Theory Part 1. Fundamentals, Ed. by P. Schwerdtfeger. Amsterdam, The Netherland: Elsevier, 2002. 946 p. ISBN: 978-0444512499.

5. Zaitsevskii A. V., van Wüllen C., Titov A. V. Relativistic pseudopotential model for superheavy elements: applications to chemistry of eka-Hg and eka-Pb // Russ. Chem. Rev. 2009. V. 78, N 12. P. 1173-1181.

6. Zaitsevskii A., van Wüllen C., Rykova E. A., Titov A. V. Two-component relativistic density functional theory modeling of the adsorption of element 114(eka-lead) on gold // Phys. Chem. Chem. Phys. 2010. V. 12, N 16. P. 4152-4156.

7. Stewart J. J. P. Optimization of Parameters for Semiempirical Methods I. Method // J. Comput. Chem. 1989. V. 10, N 2. P. 209-220.

8. Ren P., Ponder J. W. Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation // J. Phys. Chem. B. 2003. V. 107, N 24. P. 5933-5947.

9. White C. A., Head-Gordon M. Derivation and efficient implementation of the fast multipole method // J. Chem. Phys. 1994. V. 101, N 8. P. 6593-6605.

10. Strain M. C., Scuseria G. E., Frisch M. J. Achieving linear scaling for the electronic quantum Coulomb problem // Science. 1996. V. 271, N 5245. P. 51-53.

11. Challacombe M., Schwegler E., Almlöf J. Fast assembly of the Coulomb matrix: A quantum chemical tree code //J. Chem. Phys. 1996. V. 104, N 12. P. 4685-4698.

12. Strout D. L., Scuseria G. E. A quantitative study of the scaling properties of the Hartree-Fock method // J. Chem. Phys. 1995. V. 102, N 21. P. 8448-8452.

13. Dement'ev A. I., Stepanov N. F., Yarovoi S. S. Convergence Problems in the Solution of SCF Equations // Int. J. Quantum Chem. 1974. V. 8, N 1. P. 107-117.

14. Pongor G. On the application of Saunders' level shifting technique to CNDO/2 calculations // Chem. Phys. Lett. 1974. V. 24, N 4. P. 603-605.

15. Dennis J. E., More J. J. Quasi-Newton Methods, Motivation and Theory // SIAM Review. 1977. V. 19, N 1. P. 46 89.

16. Fiacco A. V., McCormick G. P. Nonlinear Programming. New York: John Willey and Sons, 1983. 210 p. ISBN: 0898712548, 9780898712544.

17. Baker J. Techniques for Geometry Optimization: A Comparison of Cartesian and Natural Internal Coordinates // J. Comput. Chem. 1993. V. 14, N 9. P. 1085-1100.

18. Pyykkö P. Relativistic Effects in Structural Chemistry // Chem. Rev. 1988. V. 88, N 3. P. 563-594.

19. Gropen O. Effective core potentials // in: Methods in Computational Chemistry. Vol 2. Relativistic Effects in Atoms and Molecules, Ed. by S. Wilson. New York: Plenum Press, 1988. P. 109-135. ISBN: 0306429462, 9780306429460.

20. Ermler W. C., Ross R. B., Christiansen P. A. Spin-Orbit Coupling and Other Relativistic Effects in Atoms and Molecules // Adv. Quantum Chem. 1988. V. 19. P. 139-182.

21. Hellmann H. A New Approximation Method in the Problem of Many Electrons // J. Chem. Phys. 1935. V. 3, N 1. P. 61.

22. Gornbas P. Über die metallische Bindung // Z. für Physik. 1935. V. 94, N 7-8. P. 473-488.

23. McMurchie L. E., Davidson E. R. Calculation of integrals over ab initio pseudopotentials // J. Comput. Phys. 1981. V. 44, N 2. P. 289-301.

24. Pitzer R. M., Winter N. W. Spin-orbit (core) and core potential integrals // Int. J. Quantum Chem. 1991. V. 40, N 6. P. 773 780.

25. van Lenthe E., Snijders J. G., Baerends E. J. The zero-order regular approximation for relativistic effects: The effect of spin-orbit coupling in closed shell molecules // J. Chem. Phys. 1996. V. 105, N 15. P. 6505-6516.

26. Boys S. F. Electronic Wave Functions. I. A General Method of Calculation for the Stationary States of Any Molecular System // Proc. Roy. Soc. London A. 1950. V. 200, N 1063. P. 542-554.

27. Whitten J. L. Gaussian lobe function expansions of Hartree-Fock solutions for the first-row atoms and ethylene // J. Chem. Phys. 1966. V. 44, N 1. P. 359-364.

28. Clementi E., Davis D. R. Barrier to Internal Rotation in Ethane // J. Chem. Phys. 1966. V. 45, N 7. P. 2593-2599.

29. Jeziorska M., Cencek W., Patkowski K. et al. Pair potential for helium from symmetry-adapted perturbation theory calculations and from supermolecular data //J. Chem. Phys. 2007. V. 127, N 12. P. 124303.

30. van Mourik Т., Dunning Т. H. A new ab initio potential energy curve for the helium dimer //J. Chem. Phys. 1999. V. Ill, N 20. P. 9248-9258.

31. Gdanitz R. J. Accurately solving the electronic Schrodinger equation of atoms and molecules using explicitly correlated (ri2-)MR-CI. VI. The helium dimer (He2) revisited // Mol. Phys. 2001. V. 99, N 11. P. 923-930.

32. Anderson J. B. Comment on "An exact quantum Monte Carlo calculation of the helium-helium intermolecular potential"[J. Chem. Phys. 115, 4546 (2001)] // J. Chem. Phys. 2004. V. 120, N 20. P. 9886-9887.

33. Merritt J. M., Bondybey V. E., Heaven M. C. Beryllium Dimer—Caught in the Act of Bonding // Science. 2009. V. 324, N 5934. P. 1548-1551.

34. Stark J., Meyer W. The ground state potential of the beryllium dimer // Chem. Phys. Lett. 1996. V. 258, N 3-4. P. 421-426.

35. Gdanitz R. J. Accurately solving the electronic Schrodinger equation of atoms and molecules using explicitly correlated (ri2-)MR-CI. The ground state of beryllium dimer (Be2) // Chem. Phys. Lett. 1999. V. 312, N 5 6. P. 578 584.

36. Roeggen I., Veseth L. Interatomic Potential for the XState of Be2, Revisited // Int. J. Quantum Chem. 2005. V. 101, N 2. P. 201-210.

37. Spirko V. Potential energy curve of Вег in its ground electronic state // J. Mol. Spectr. 2006. V. 235, N 2. P. 268-270.

38. Patkowski K., Spirko V., Szalewicz K. On the Elusive Twelfth Vibrational State of Beryllium Dimer // Science. 2009. V. 326, N 5958. P. 1382-1384.

39. Koput J. The ground-state potential energy function of a beryllium dimer determined using the single-reference coupled-cluster approach // Phys. Chem. Chem. Phys. 2011. V. 13, N 45. P. 20311-20317.

Подписано в печать:

01.07.2013

Заказ № 8612 Тираж -120 экз. Печать трафаретная. Типография «11-й ФОРМАТ» ИНН 7726330900 115230, Москва, Варшавское ш., 36 (499) 788-78-56 www.autoreferat.ru

Московский государственный университет имени М. В. Ломоносова

Химический факультет Федеральное государственное бюджетное учреждение науки Объединённый

институт высоких температур РАН

На правах рукописи

05201351715

Митин Александр Васильевич

Новые методы решения электронных уравнений квантовой химии

и их применение

02.00.17 - Математическая и квантовая химия

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

Москва - 2013

Содержание

Введение ................................... 6

Часть I. Новые методы решения уравнений теории электронной структуры молекул 13

Глава 1. Эффективные методы решения уравнения Хартри-Фока

со свойством линейной масштабируемости............14

1.1. Исследование асимптотических свойств числа ненулевых одно-и двух-электронных интегралов, наибольших матричных элементов матрицы Фока и наибольших ЛКАО коэффициентов молекулярных орбиталей......................16

1.2. Исследование масштабируемости по числу атомов сложности вычисления двух-электронной части матрицы Фока с хранимыми в файлах интегралами....................37

1.3. Метод ССП Хартри-Фока со сжатием данных и линейно масштабируемым вычислением матрицы Фока............46

1.4. Супер прямой линейно масштабируемый метод решения уравнения Хартри-Фока.........................57

Глава 2. Эффективные методы улучшения сходимости итераций метода самосогласованного поля...............63

2.1. Метод псевдопотенциала для улучшения сходимости ССП итераций .................................65

2.2. Метод динамического сдвига уровней для улучшения сходимости ССП итераций..........................67

2.3. Теория функционала плотности и нарушение принципа заполнения орбиталей...........................73

2.4. Лианеризованное уравнение для построения экстраполяцион-ных методов.............................77

2.5. Линейные методы экстраполяции.................78

2.6. Линейные проекционные экстраполяционные методы......80

2.7. Линейные методы экстраполяции, использующие специальные функции ...............................83

Глава 3. Эффективный метод оптимизации геометрии молекул 87

3.1. Алгоритм оптимизации геометрии молекул............89

3.2. Критерий остановки оптимизационного процесса........92

3.3. Результаты тестовых расчётов и обсуждение...........93

Глава 4. Эффективный метод параллелизации программ методов Хартри-Фока и ТФП.......................98

Глава 5. Двухкомпонентный релятивистский метод ТФП с RECP.

Расчёт потенциалов ионизации и энергий сродства к электрону атомов галогенидов от Вг до элемента (117) и их двухатомных молекул...............................102

5.1. Теоретические основания......................104

5.2. Детали расчёта атомов и молекул.................106

Глава 6. Развитие метода конфигурационного взаимодействия 115

6.1. Многоссылочный метод конфигурационного взаимодействия с явной корреляцией большого числа электронов.........117

6.2. Нахождение экстремальных собственных значений и соответствующих собственных векторов обобщёнными методами Кару ша и Ланцоша...........................121

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

6.4. Итерационные методы типа Ньютона-Лагранжа и другие блочные методы..............................137

Глава 7. Построение сгруппированных гауссовых базисных функций для неэмпирических расчётов.................150

7.1. Базисы сгруппированных функций гауссова типа для атомов первого периода построенные из наборов примитивных (9s5p) гауссовых функций.........................150

7.2. Малые валентно-расщеплённые базисы 3-21sp и 4-22sp сгруппированных гауссовых функций для атомов первого и второго периодов ............................... 155

7.3. Модифицированный базис сгруппированных гауссовых функций m6-31G* для атомов переходных металлов третьего периода167

7.4. Модифицированный базис сгруппированных гауссовых функций тб-ЗЮ* для атомов третьего периода............176

7.5. Поляризационные функции для базиса m6-31G атомов от Ga

до Кг.................................183

Глава 8. Численный алгоритм метода Данхэма решения колебательно-вращательного уравнения Шрёдингера для двухатомной молекулы...........................188

8.1. Вычисление коэффициентов Данхэма Уы и метод квази-эрми-товой аппроксимации потенциальной кривой заданной таблицей классическими ортогональными полиномами........190

8.2. Вычисление колебательно-вращательных уровней энергий методом Данхэма............................197

8.3. Оценка ошибок при вычисление коэффициентов Уы и колебательно-вращательных уровней энергий двухатомных молекул методом Данхэма..........................199

8.4. Табулирование потенциальной кривой двухатомной молекулы

в неэмпирическом расчёте .....................201

Часть И. Теоретические расчёты молекул и их исследования 205

Глава 9. Высокоточные неэмпирические расчёты малых моле

кул .....................................206

9.1. Приближение ЛКАО в методе молекулярных орбиталей .... 208

9.2. Точное решение уравнения Хартри-Фока для молекулы #2 в приближение ЛКАО.........................217

9.3. Расчёты потенциальных кривых основных состояний молекул

бУг и Ва,2...............................219

9.4. Расчёты потенциальных кривых основных состояний молекул Не2иВе2 ..............................224

9.5. Теоретический расчёт инфракрасного и Рамановского спектров молекул А/202 и А/гОз.......................242

Глава 10. Исследование многоэлектронных волновых функций

основных состояний биомолекул ..................251

10.1. Взаимная поляризация аминокислот и их взаимодействие в биомолекулах.............................253

10.2. Электронная структура комплексов [МеР/тг], их катионов и двухзарядных катионов.......................259

10.3. Реакционные центры СТЦ в реакции переноса электронов инициированной биологической активностью бактерии вкетапеИа oneidensis...............................272

Основные результаты и выводы ....................278

Литература..................................281

Введение

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

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

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

настоящей работе. Однако, из формы представления новых методов видно, что все они были реализованы в программах, что и позволило провести вычисления, представленные в работе.

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

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

матрицы Фока является эффективным методом решения уравнения Хартри-Фока, применимым для расчёта больших молекулярных систем. Был предложен также сверх-прямой метод ССП, все шаги которого линейно масштабируются по числу базисных функций. Исследовались проблемы сходимости ССП итераций, оптимизации геометрии молекул и параллелизации методов Хартри-Фока и теории функционала плотности (ТФП). Результатами этих исследований стал метод псевдопотенциала и динамического сдвига уровней для подавления осцилляций ССП итераций, метод оптимизации геометрии, основанный на симметричной формуле первого ранга для обновления матрицы гессиана, и эффективный метод параллелизации алгоритма ССП решения уравнений Хартри-Фока и ТФП. Была построена общая теория экстраполя-ционных методов в итерационных методах решения уравнений и предложено много новых экстраполяционных методов, включая проекционные экстрапо-ляционные методы, являющиеся наиболее мощными из них.

Хорошо известно, что применение нерелятивистских квантово-механи-ческих методов для молекулярных систем, образованных атомами элементов верхней части периодической таблицы Менделеева, вполне обосновано ввиду малости релятивистских эффектов в таких системах. Для молекулярных систем с атомами тяжёлых и сверх тяжёлых элементов ситуация прямо противоположная. Для таких систем использование релятивистских методов является предопределённым ввиду необходимости учёта релятивистских эффектов для правильного описания взаимодействий в таких системах. В связи с этим в нашей работе была представлена эффективная реализация релятивистского двухкомпонентного метода ТФП с использованием релятивистского эффективного остовного псевдопотенциала (ЛЕСР). Эта программа позволила провести расчёт потенциалов ионизации и энергий сродства к электрону атомов галогенидов от Вг до и их двухатомных молекул, а также соответствующие первые расчёты атома элемента (117) и его двухатомной молекулы.

Метод конфигурационного взаимодействия (КВ) является одним из основных методов учёта корреляционной энергии в неэмпирических расчётах.

Чтобы показать, что потенциал этого метода далеко не исчерпан, нами был разработан многоссылочный метод КВ, позволяющий явно коррелировать большое число электронов и рассчитывать состояния высокой мультиплетно-сти. Возможность коррелировать до 100 электронов на компьютере с 1 СВ оперативной памяти является непревзойдённой до настоящего времени. Этот метод был использован нами при исследование адсорбции молекулы С#2 на поверхности Ш.

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

Приближение линейной комбинации атомных орбиталей (ЛКАО) является одним из важнейших в теории теории электронной структуры молекул, поскольку подавляющее большинство неэмпирических расчётов проводятся с использованием этого приближения, и выбор системы базисных функций во многом определяет качество рассчитываемых параметров. В наших работах по разработке новых наборов базисных функций были оптимизированы малые сегментно-сгруппированные базисы для атомов первого периода и малые валентно-расщеплённые базисы для атомов первого и второго периодов, заметно превосходящие аналогичные базисы 3-2Ю и 4-31 С, а также базисы шб-ЗЮ* для атомов третьего периода, включая атомы переходных металлов, имеющих существенное преимущество перед аналогичными базисами б-ЗЮ*.

Учитывая это, метод 04 для термохимических расчётов соединений переходных металлов был построен с использованием модифицированных базисов тб-ЗЮ для атомов переходных металлов, разработанных нами.

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

Во второй части диссертационной работы представлены теоретические расчёты и исследования волновых функций, которые были связаны с достижением или использованием экстремальных результатов. Так в исследование молекулы Яг была поставлена задача получения численно точного решения уравнения Хартри-Фока в приближение Л К АО. Исследование показало, что молекулярные орбитали (МО) уравнения Хартри-Фока и орбитали этого же уравнения в приближение Л К АО отличаются, поскольку обычное приближение ЛКАО неодинаково описывает асимптотические свойства МО в пределе разъединённых атомов и пределе объединённого атома. В связи с этим

нами было предложено модифицированное приближение ЛКАО, в котором асимптотические свойства МО описываются корректно. Использование этого приближения позволило вычислить 12 правильных цифр в величине Хартри-Фоковской энергии Н2, т. е. получить численно точное решение. Предыдущая точность была увеличена на шесть порядков.

Неэмпирические расчёты основных состояний молекул Не2 и Ве2 преследовали две цели. Во первых, они были нацелены на расчёты максимально высокой точности, а, во-вторых, продемонстрировать, что высокоточные расчёты могут выполняться также методом КВ. Обе эти цели были успешно достигнуты. Результаты полученные нами для Не2 по крайней мере не уступают по точности результатам, полученными другими методами, а расчёты Ве2 выполнены на пределе возможной точности неэмпирических н�