Содержание к диссертации
Введение
1 Введение 4
2 Теоретический анализ фотодиссоциации озона в состояниях, отвечающих полосе Хартли 11
2.1 Цели и задачи работы 11
2.2 Квантовая динамика 14
2.3 Расчет спектральных характеристик 21
2.4 Понижение размерности в задачах квантовой динамики 25
2.5 Эффективный одномерный гамильтониан в задаче фотодиссоциации 27
2.6 Многоуровневое обобщение и проблема диабатических состояний 35
2.7 Электронные расчеты 39
2.8 Вычисление эффективного гамильтониана: аппроксимация диабатических кривых и гармонических частот вдоль реакционного пути 46
2.9 Детали динамического расчета: программы, вычисления и сходимость результатов 52
2.10 Обсуждение результатов 56
2.11 Вращательное уширение 62
3 Классическая динамика системы O^D) + Н2О —> 20Н 67
3.1 Цели и задачи работы 67
3.2 Молекулярная динамика. Общие вопросы 69
3.3 Неэшшрическая молекулярная динамика 75
3.4 Формулировка Борна-Оппенгеймера 79
3.5 Формулировка Кара-Паринелло 81
3.6 Информация о системе (Н)2(0)2 84
3.7 Использованные модельные подходы 87
3.8 Анализ поверхности потенциальной энергии 88
3.9 Динамика Борна-Оппенгеймера 92
3.10 Динамика Кара-Паринелло 103
Основные результаты и выводы 106
- Квантовая динамика
- Многоуровневое обобщение и проблема диабатических состояний
- Молекулярная динамика. Общие вопросы
- Динамика Борна-Оппенгеймера
Введение к работе
Объект исследования и актуальность темы
Ряд химических систем на сегодняшний день представляют особый интерес для теоретического рассмотрения. Одна из них — система, состоящая из фотовозбужденной молекулы озона и воды. Она интересна не только с теоретической, но и с прикладной точки зрения, поскольку ее компоненты составляют одну из наиболее важных частей стратосферы — озоновый экран, динамика которого с 80-х годов XX века привлекает пристальное внимание.
В диссертационной работе рассмотрены два основных процесса, протекающие в этой системе: непосредственно фотодиссоциация электронно-возбужденной молекулы озона и реакция продуктов диссоциации с молекулой воды.
Фотодиссоциация молекулы озона была предметом многочисленных теоретических и экспериментальных исследований, но тем не менее и по сей день часть вопросов, ее касающихся, остаются открытыми Так, например, не до конца решена проблема неадиабатического взаимодействия различных электронных состояний, которое обуславливает существование двух экспериментально наблюдаемых каналов диссоциации (с образованием синглетных и триплетных продуктов). Остается открытым вопрос о необходимости учета этого взаимодействия в теоретическом моделировании спектра поглощения молекулярного озона в ультрафиолетовом диапазоне (полоса Хартли).
Высокая реакционная способность синглетного атома кислорода — одного из продуктов фотодиссоциации — приводит к тому, что он может реагировать практически с любой молекулой, находящейся в атмосфере, за исключением, пожалуй, только кислорода, азота, углекислого газа и благородных газов В частности, реакция синглетного атома кислорода с молекулами воды приводит к появлению в атмосфере гидроксильных радикалов, играющих ключевую роль в естественном цикле разрушения озона. Знание динамических характеристик реакции образования гидроксил-радикалов позволит лучше понять механизмы, регулирующие количество озона в стратосфере.
I РОС НАЦИОНАЛЬНАЯ I
| БИБЛИОТЕКА I
Цель работы
Целью работы было теоретическое исследование процесса фотодиссоциации молекулы озона с учетом неадиабатического взаимодействия электронных состояний и изучение влияния этого взаимодействия на профиль полосы поглощения, моделирование реакционной динамики синглетного атома кислорода с молекулой воды.
В соответствии с целью диссертационной работы были сформулированы следующие задачи
На основании неэмпирических расчетов высокого уровня получить информацию о характере и величине неадиабатического взаимодействия электронных состояний.
Предложить варианты построения и определить границы применимости методов, использующих понижение размерности в задачах квантовой реакционной динамики.
Используя полученные динамические уравнения, а также матричные элементы неадиабатических взаимодействий, методом волновых пакетов рассчитать и сопоставить с экспериментом относительный выход синглетных и триплетных продуктов в процессе диссоциации фотовозбужденной молекулы озона.
Оценить влияние взаимодействия электронных состояний на профиль полосы Хартли и определить, насколько необходим его учет при теоретическом моделировании ультрафиолетового спектра молекулярного озона.
Методами неэмпирической классической реакционной динамики изучить превращения в системе НгО + 0(XD). На основании полученной картины обсудить применимость метода и определить его место среди других вариантов изучения реакционной динамики.
Оценить характеристические времена превращений в системе НгО + 0(*D) и проанализировать распределение избыточной энергии по степеням свободы продуктов.
Научная новизна
В диссертационной работе впервые:
Классифицированы различные подходы понижения размерности в приложении к задачам квантовой динамики.
Найдены неадиабатические матричные элементы взаимодействия возбужденных синглетных электронных состояний молекулы озона в одномерных сечениях на основании результатов неэмпирических квантово-химических расчетов высокого уровня.
Теоретически рассчитан выход синглетных продуктов в процессе фотодиссоциации молекулы озона в ультрафиолетовом диапазоне.
Оценено влияние неадиабатического взаимодействия состояний на профиль полосы поглощения озона.
Проанализирована применимость методов неэмпирической классической динамики к системе НгО + 0(XD).
Определены времена жизни промежуточных метастабильных состояний в системе H20 + 0(1D) и проанализированы распределения трансляционной, колебательной и вращательной энергии продуктов этой реакции
Научная и практическая ценность
В настоящей работе проведено систематическое квантово-динамическое моделирование фотодиссоциации молекулы озона с использованием подходов с понижением размерности и квазидиабатических потенциальных кривых Применение указанных приближений впервые позволило произвести теоретический расчет выхода синглетных и триплетиых продуктов в реакции фотодиссоциации, а также оценить влияние неадиабатических взаимодействий на профиль полосы поглощения Результаты работы могут быть использованы как для дальнейшего развития методов моделирования фотодиссоциации озона, так и для получения новых данных по выходам различных продуктов в этом процессе Кроме того, ряд подходов и выводов работы может быть применен в дальнейшем развитии систематических методов моделирования
динамики химических реакций с учетом квантовых эффектов Анализ реакционной динамики системы 1120 + 0() выявляет особенности колебательно-вращательного распределения энергии по степеням свободы продуктов реакции, которые могут быть обнаружены экспериментально
Публикации и апробация работы
Квантовая динамика
Согласно постулатам квантовой механики, эволюция квантовой системы описывается временным уравнением Шредингера с начальным условием Ф(Х,). _ = $0(Х). Волновая функция Ф(Х,), зависящая от совокупности пространственных переменных X и времени, полностью определяет состояние системы, т.е., зная ее, можно вычислить любое интересующее свойство А по формуле Символ А обозначает эрмитов оператор, отвечающий физической величине А. В уравнении (2.3) оператор Н — это оператор Гамильтона системы, или гамильтониан. Оператор Гамильтона тесно связан с полной энергией системы, в частности, если время не входит в него явно как переменная, среднее значение (ФЯФ), вычисленное на волновой функции Ф(Х, t) какой-либо системы, есть в точности полная энергия этой системы. Таким образом, наиболее естественное и подробное описание динамики атомно-молекулярной системы может быть получено при решении уравнения (2.3). Однако возможны ситуации, когда исследователя интересует не детальное описание самого процесса, а лишь состояние системы после его окончания. Более того, не всегда точно известно начальное квантовое состояние или их заведомо несколько. В этом случае используется стационарный подход, оперирующий понятием S-матрицы, которая связывает начальное и конечное асимптотические состояния ФІП причем от времени они явно не зависят и являются суперпозицией решений стационарного уравнения Шредингера Вероятность перехода системы из состояния г в состояние / {Pi- /) дается квадратом модуля матричного элемента S Такое описание процессов, происходящих в квантовых системах носит название теории рассеяния. Подробности можно найти в книгах [15, 16]. Связь решений уравнении (2.3) и (2.6) для случая не зависящего от времени гамильтониана может быть получена методом разделения переменных. Оба подхода, и временной, и стационарный, могут быть применены к задачам квантовой реакционной динамики. Выбор определяется особенностями рассматриваемого процесса, а также вычислительными и программ- ными ресурсами. Хорошее сравнение преимуществ и недостатков обоих методов можно найти в [17]. Кратко перечислим их. Временные методы особенно эффективны в случае, когда точно известно начальное состояние реагирующей системы. При этом получаемые вероятности реакций оказываются суммированы по всем конечным состояниям, однако возможен расчет вероятности каждого перехода (stateo-state), что дает одну строку матрицы S. Вычислительная сложность асимптотически оценивается в ЛГ2, где N — условный параметр, характеризующий размер задачи.
Стационарные подходы позволяют получить информацию о каждом из возможных переходов, однако это достигается ценой возросших вычислительных затрат, связанных с диагонализациеи матриц и асимптотически равных 0(N3). В ряде случаев возможно их уменьшение применением итерационных схем. Как уже отмечалось, временные подходы дают более ясную физическую картину процесса и подобны методу классических траекторий (см. главу 3), поскольку исследователь выбирает начальные условия, а далее изучает эволюцию системы. Задача фотодиссоциации — типичный пример такого процесса, поскольку начальное состояние квантовой системы известно. Кроме того, реализация временных методов часто проще стационарных. По причинам, приведенным выше, в нашей работе был использован временнбй формализм, основанный на уравнении (2.3). Мы приведем краткий обзор способов его решения. Представим волновую функцию системы в момент времени t как результат действия некоторого линейного оператора — оператора эволюции 17 — на волновую функцию в начальный момент времени Поскольку уравнение (2.3) является линейным дифференциальным уравнением в частных производных, такое представление возможно. Тогда (2.3) можно переписать как уравнение на оператор эволюции ih -U = HU (2.9) с начальным условием U(0) = , где — единичный оператор. Формальное решение уравнения (2.9) может быть записано в виде операторной экспоненты где экспонента понимается в смысле сходящегося по операторной норме ряда: Множество операторов U(t) при различных t образует абелеву группу: U(ti)U(t2) = U{t\ + іг) чем мы воспользуемся ниже. Сложив уравнение (2.8) с эрмитово-сопряженным к нему, используем унитарность оператора эволюции: где At — шаг интегрирования, a tn±i = tn ± At. Разложив экспоненту в ряд, приходим к простейшей трехточечной схеме решения уравнения (2.3): Метод, использующий выражение (2.13), подвержен влиянию численных неустойчивостей, накоплению ошибок и требует коротких шагов. Оценку максимальной длины шага можно найти в [18]. Формула (2.13) применима и в случае явно зависящего от времени гамильтониана. Кроме того, метод очень прост в реализации. К похожей форме — формуле Франка-Николсона [19] —- можно придти, используя групповое свойство оператора U(t) (2.14) Вновь разлагая экспоненту в ряд, получаем искомую формулу: Метод с ее использованием сильно устойчив, однако требует обращения матриц. Обращение матрицы — вычислительно трудная операция, хотя оно может быть выполнено только один раз в самом начале вычислений в случае гамильтониана, не зависящего от времени.
Теоретически возможно использование метода и для зависящего от времени гамильтониана, однако обращение тогда необходимо делать на каждом шаге. Матричную экспоненту (2.10) можно представить не только разложением в ряд Тейлора, но и в виде ряда по другим ортогональным системам, например, полиномам Чебышева [18]: Етіп и Етах — минимальная и максимальная энергия системы в выбранном методе вычислений, CQ — 1; Сп = 2, п 1, a Jn — функция Бесселя первого рода (мнимого аргумента). Основное преимущество (2.16) состоит в том, что функция Бесселя Jn(x) убывает экспоненциально при х п. Кроме того, связь многочленов Чебышева с многочленами наилучших равномерных приближений дает указание на то, что ошибка аппроксимации равномерно распределена на интервале (—1,1), т.е. в терминах исходного гамильтониана распределена равномерно по всему энергетическому спек-тру. Результат действия Тп( iHnorm)ty(t) = фп может быть вычислен по рекуррентному соотношению с начальными условиями ф0 = Ф(), ф\ = — iHnorm {t), Высокая точность аппроксимации (2.16) и возможность ее априорной оценки делает этот метод исключительным. Однако он неприменим к системам с явной зависимостью гамильтониана от времени [18]. Несколько иной способ приближения оператора эволюции представляет метод расщепления оператора. U(t) записывают в виде произведения ехр (-ІІА) = ехр (-±0) ехр (-Ltv) exp (-±tf) +0{Є). (2.21) Существуют регулярные способы выписывания таких произведений с точностью до любого заданного порядка, однако их сложность быстро растет, и на практике ограничиваются формулами, точными до порядка 2 и 3 включительно. Общую теорию и многочисленные формулы можно найти в [20, 21]. Основное достоинство (2.21) заключается в возможности выбирать такие представления, в которых операторы Г или V диагопальны, хотя и не одновременно (в силу их некоммутативности). Вычисление же любой функции диагонального оператора — тривиальная задача. В координатном представлении (значения волновой функции на пространственной сетке) диагоналей оператор потенциальной энергии. Оператор кинетической энергии диагоналей в импульсном представлении, которое связано с координатным преобразованием Фурье. Основное достоинство этого подхода в том, что переход от координатного к импульсному представлению и обратно, может быть выполнен при помощи быстрого преобразования Фурье [22], сложность которого оценивается как 0(N In iV), где N - число точек пространственной сетки. Наш обзор не исчерпывает всех способов приближенного решения уравнения (2.3). Подбор примеров преследовал цели иллюстрации использования различных свойств точного решения для построения той или иной формы его приближения. Обзор других подходов можно найти в [18,17,19]. В нашей работе мы выбрали метод разбиения оператора эволюции (2.21) по нескольким причинам:
Многоуровневое обобщение и проблема диабатических состояний
Задача фотодиссоциации включает себя динамику системы в нескольких электронных состояниях. Для этого нам необходимо многоуровневое обобщение гамильтониана (2.71) и учет взаимодействия электронных и ядерных степеней свободы. Многоуровневое обобщение гамильтониана реакционного пути (2.71) может быть получено, если считать его компоненты матрицами. Таким образом, можно написать Определенные таким образом диабатические состояния еще называют строго диабатическими. Диабатические состояния в общем случае не являются собственными для электронного гамильтониана, и матрица U в (2.83) не диагональна. Условие (2.85) приводит к следующей матричной системе дифференциальных уравнений, определяющих матрицу Т: где F — блок пхп матрицы F, отвечающий п интересующим состояниям. Решение уравнения (2.86) существует тогда и только тогда, когда для любой пары ядерных координат Rj, #& выполнено равенство где нижний индекс у величин F указывает на дифференцирование по соответствующей переменной R в определении их матричных элементов (2.81). Следует отметить, что условие (2.87) в общем случае не выполняется для многоатомных молекул. Однако его выполнение может быть обеспече- но на какой-то кривой в (3/V — 6)-мерном пространстве внутренних координат системы. Таким образом, невозможность в общем случае построения строго диабатических состояний многоатомной системы не исключает построение таких состояний, матричные элементы взаимодействия которых Fap(R) равны нулю вдоль некоторой кривой в пространстве внутренних переменных. Построение подобных состояний хорошо согласовано с концепцией реакционного пути. Систематическое рассмотрение диабатических состояний было проведено в работе [36], где проанализированы основные идеи, касающиеся диабатических состояний двухатомных молекул. В той же работе обсуждались и иные определения диабатических состояний как таких состояний, которые гладко пересекают области (квази)пересечений. Это и другие определения приобрели некоторое распространение в литературе. Так, например, (квази)диабатические состояния конструировали на основе орбитальной картины [37, 38], подобия по конфигурациям (configu rational uniformity) [39], и по критерию минимального изменения дипольных моментов [40, 41]. Все эти методы фактически опираются на идею о том, что вдали от областей квазипересечения волновые функции диабатичесих и адиабатических состояний совпадают, а в областях квазипересечений параметры, описывающие волновые функции диабатических состояний, ведут себя гладким образом.
Единственным примером построения квазидиабатических состояний в молекуле озона можно назвать работу [42]. В ней были вычислены матричные элементы неадиабатического взаимодействия FQp(R) при различных геометрических конфигурациях для первых двух низколежащих элек- тронных состояний, на основе градиентного кода в методе конфигурационного взаимодействия с набором исходных конфигураций (MR-СІ, см. также ниже раздел 2.7). В этой работе были получены и проанализированы функциональные зависимости матричных элементов неадиабатического взаимодействия вблизи точек квазипересечения. Мы также придерживались определения квазидиабатических состояний как таких, которые гладко проходят области квазипересечений, а в областях, удаленных от этих областей, совпадают с адиабатическими состояниями. Таким образом, процедура построения диабатических состояний состояла в идентификации частей диабатических кривых по обе стороны от области квазипересечения и гладкой сшивки этих частей в единую кривую (см. раздел 2.8). Многоуровневое обобщение эффективного гамильтониана (2.73) использует общий для всех выбранных электронных состояний вид реакционного пути /(т/) (см. стр. 34). 2.7. Электронные расчеты Квантово-химические расчеты, представленные в нашей работе, были выполнены для 8 низших сингл етных состояний. Среди низкол ежащих возбужденных электроЕшых состояний молекулы озона находится также ряд триплетных состояний, однако они рассмотрены не были ввиду малой интенсивности синглет-триплетных переходов в молекулах, содержащих только легкие атомы. Предположительно, триплетные состояния играют роль в формировании ближнего инфракрасного спектра озона [43]. Молекула озона — популярный объект для проверки методов расчета электронной структуры в силу многокофигурационного характера своего основного состояния [44].
Поскольку нам был необходим расчет сечений ППЭ нескольких электронных состояний, рассмотрение результатов одно-конфигурационных подходов можно в основном опустить. Отметим лишь несколько работ, содержащих результаты применения одноконфигураци-онных подходов к молекуле озона [44, 45, 46, 47]. Расчет вертикальных энергий возбуждения, а также участков ППЭ или даже глобальных поверхностей выбранных электронных состояний проводился неоднократно разными авторами. Вертикальные энергии возбуждения оценивали методом уравнений движения в приближении связанных кластеров (ЕОМ-СС) [48], а потенциалы ионизации — многокофигу-рационным спин-тензорным пропагаторным методом [49, 50]. Глобальная ППЭ основного состояния была построена в 2001-2003 гг с использованием метода конфигурационного взаимодействия с множеством исходных конфигураций (MR-СІ) с поправкой на черырехкратные несвязные возбуждения (поправка Давидсона) [51, 52]. Отдельные участки вблизи различных локальных минимумов на поверхности основного состояния изучали модифицированным методом функционала ПЛОТЕЮСТИ [53], обобщенным методом валентных связей [54] и методом конфигурационного взаимодействия с поправками по теории связанных кластеров (MR-AQCC) [33]. Ряд низколежащих электронных состояний молекулы озона, а также фрагменты их потенциальных поверхностей в окрестностях минимумов, отвечающих геометрическим конфигурациям симметрии C2v, были рассчитаны в работах [55, 56]. Наши расчеты во многом были основаны на опыте этих двух работ, поэтому осветим их основные результаты подробнее. В обеих работах использовали многоконфигурационные обобщения теории возмущений второго порядка. Хотя формулировки многоконфигурационной теории возмущений различны, обе они базируются на исходной волновой функции многоконфигурационного метода самосогласованного поля с полным активным пространством (CASSCF). Подход, использованный в [55], сформулирован в работах [57, 58] и получил название CASPT2. Он использует идеологию "сначала диагонализация, потом возмущение", являясь, таким образом, теорией возмущений для выбранного электронного состояния [59]. Энергии возбужденных состояний в работе [56] рассчитывались на основе метода MRMP2 [60]. Этот подход также сформулирован для одного электронного состояния в идеологии "сначала диагонализация, потом возмущение", но отличается от CASPT2 выбором гамильтониана нулевого приближения и другими техническими деталями.
Молекулярная динамика. Общие вопросы
В этой части под термином молекулярная динамика мы будем подразумевать методы, позволяющие предсказывать реакционную способность и направление химических реакций. Основным подходом, используемым для этих целей, служит метод классических траекторий. Этот метод основан на изучении множества классических траекторий движения ядер молекулы в потенциале, создаваемом другими ядрами и распределением электронной плотности. Примером программного кода, решающего эту задачу, служит Venus99 [74]. Получение количественных результатов по реакционному рассеянию требует изучения эволюции целого пучка траекторий с различными начальными условиями. Причин, обуславливающих необходимость использования набора начальных условий, несколько: с одной стороны, в химическую реакцию вступает множество молекул, каждая из которых может находиться в своем состоянии — имеется равновесие между разрешенными по энергии состояниями с различными значениями координат и импульсов частиц. Таким образом, даже в рамках классической механики следует учитывать статистическую природу начальных условий. Существуют несколько методов, позволяющих построить начальные распределения, находящиеся в статистическом равновесии. Перечислим их. Выборка по ортаптам [75]. Метод может быть использован для многоатомных молекул в случае, когда полная размерность фазового пространства велика и число траекторий в выборке много меньше числа ор-тантов — частей фазового пространства с определенным знаком у каждой переменной 1. На гиперповерхности в фазовом пространстве где Н — функция Гамильтона системы, a VQ — минимум потенциала, определяются максимальные и минимальные допустимые значения координат и импульсов: x ,xf и р .,р ., соответственно. Искомые начальные значения импульсов и координат получаются как где vi — случайный вектор в фазовом пространстве системы с компонентами О или 1, pSCaie — множитель, определяющий отношение средней потенци-альной энергии к средней кинетической (для гармонического осциллятора — 1), х\ — координаты частиц системы в равновесии.
Процедура также включает вычитание скорости центра масс системы и придания нужного углового момента и повторяется до тех пор, пока "исправленные" значения Нетрудно проверить, что всего таких частей 2п, где п — размерность пространства. Так, в двумерном пространстве их 4 и называются они квадрантами, в трехмерном пространстве — 8 и называются они октантами, соответственно. скоростей и импульсов не будут отвечать заданному значению углового момента и покоящемуся центру масс. Сходимость обычно достигается за несколько шагов. Микроканоническое распределение по нормальным модам [76]. Энергии п классических гармонических осцилляторов в микроканоническом ансамбле определены выражением где х, рх — вектор-столбцы с компонентами декартовых координат и импульсов соответственно, а М — матрица масс (метрический тензор) в декартовых координатах. Данный метод описывает точное микроканоническое распределение в системе гармонических осцилляторов. Как и метод выборки по ортантам, он лишь приближенно описывает распределение энергии в системе ангармонических осцилляторов. Термическое распределение. Энергия системы осцилляторов описывается распределением Больцмана при заданной температуре ТУІЬ. Плотность вероятности найти осциллятор в состоянии с квантовым числом щ для ъ-той моды дается выражением = (- )1 (-)] (314) Энергия, приходящаяся на каждую колебательную моду, есть ТИОІ{Щ + \). Статистический вклад каждого начального условия полностью определен заданным набором щ и распределением Больцмана, а значения координат и импульсов системы выбирают так же, как в микроканоническом распределении. Помимо теплового равновесия, другой причиной, требующей процедуры усреднения по траекториям, служит сама квантовая природа начального состояния. В частности, данеє когда (жесткая) молекула при средних и низких температурах находится в своем основном состоянии, ее ядра колеблются {нулевые колебания), что требует построения такого набора классических состояний, которые бы учитывали этот эффект. Кроме того, в ряде нетермически активируемых реакций, в частности фотоинду-цированных процессов, реагенты находятся в возбужденных электронно-колебательных состояниях, и имеют существенно отличное от равновесного (больцмаиовского) распределение энергии по степеням свободы.
Выбор начальных условий, в частности, может производится при помощи функции Вигнера [77]. Функция Вигнера — функция в фазовом пространстве систе- мы, играющая роль (псевдо)плотности вероятности распределения классических координат и импульсов для квантовой системы, описываемой волновой функцией Ф(хі,..., хп) х exp (f (рії/і + РпУпїїіїуі...dyn. Функция Вигнера не является плотностью вероятности в строгом смысле слова, поскольку она может приобретать отрицательные значения. Выра-ЖЄЕІИЄ (3.15), будучи проинтегрировано по всем доступным значениям импульса, дает всюду неотрицательно определенную плотность вероятности в конфигурационном пространстве: а интеграл по всему конфигурационному пространству дает всюду неотрицательную плотность вероятности для импульсов в системе являясь, таким образом, источником корректных плотностей вероятности координат и импульсов по отдельности. Использование функции Вигнера позволяет получить такой набор начальных условий (и их весов), который отвечает, скажем, нулевым колебаниям в системе или описывает возбужденное состояние. В качестве примера применения подходов, основанных на использовании функции Вигнера, приведем работу [78]. Один из основных выводов работы можно сформулировать так: использование начальных условий, полученных из квантового рассмотрения системы, значительно улучшает согласие результатов метода классических траекторий с экспериментом. Несколько слов следует сказать об алгоритмах вычисления классических траекторий. Уравнения движения в форме Ньютона (иногда Ла-гранжа или Гамильтона) решаются обычными методами, пригодными для решения обыкновенных дифференциальных уравнений. Простейшие алгоритмы интегрирования, алгоритм Верлета, алгоритм "прыжок лягушки" и алгоритм Верлета по скоростям используют для расчета классических траекторий (равно как и поиска оптимальных структур) в программе CHARMM [79]. Более сложные, но и более точные и численно устойчивые методы — многошаговые методы Адамса с начальными шагами по методу Рунге-Кутты [80] — используются в уже упоминавшейся программе классической динамики [74]. Основные требования к алгоритмам: точность, вычислительная простота и устойчивость. Многошаговые методы типа Адамса (Адамса-Моултона) удовлетворяют требованиям по точности и устойчивости, требуя однако несколько больших вычислительных затрат, нежели алгоритмы Верлета и "прыжок лягушки".
Динамика Борна-Оппенгеймера
Борн-оппенгеймеровская динамика системы Н О + O D) была исследована нами в двух приближениях: методом CASSCF и обобщенным методом валентных связей (GVB — Generalized Valence Bond). Причиной выбора служило то, что оба метода способны корректно описывать процессы гемолитической диссоциации (образование радикальных продуктов из синглетных реагентов) и перестройку волновой функции молекулы при внутренних движениях. Оба метода отвечают также вычислительным требованиям, предъявляемым к способам решения динамической задачи, поскольку возможно вычисление как значений энергии системы при заданной геометрии, так и сил, действующих на ядра. Каждому из методов присущи свои преимущества и недостатки. Так, выбор модельного пространства в методе CASSCF гораздо менее трудоемок, что дает возможность использовать его в самых различных случаях без тщательного и долгого анализа особенностей волновой функции при различных геометриях ядерного остова. В то же время неспецифичный выбор модельного пространства обычно ведет к необоснованному увеличению его размера и, соответственно, вычислительных затрат. Метод GVB, напротив, требует тщательного выбора орбитального базиса и геминалей, однако вычислительные затраты существенно ниже, чем в методе CASSCF. Мы провели тесты для ряда модельных пространств в обоих методах. В случае CASSCF в модельное пространство были включены все двукратно занятые орбитали, кроме остовных ls-орбиталей атомов кислорода. Размер модельного пространства увеличивали путем включения в него незанятых (виртуальных) молекулярных орбиталей. Последовательно расширяемые модельные пространства 14е/8 (1 виртуальная орбиталь), 14е/10 (3 виртуальных орбитали) и 14е/12 (5 виртуальных орбиталей) дальше будут обозначены lv, 3v и 5v, соответственно. Последнее пространство, наилучшим образом описывающее ППЭ системы, слишком велико для использования в динамических расчетах. Для его уменьшения мы включили в остов и обе 2є-орбитали атомов О, что при незначительном снижении точности приводит к существенному уменьшению размера модельного пространства метода CASSCF. Это дополнительное пространство —14е/10 (5 виртуальных орбиталей, 4 остовных) обозначено как 5v . В методе GVB также был рассмотрен ряд "модельных" пространств, включающих 1, 3 и 5 геминалей (пар коррелирующих орбиталей). Выбранные геминальные пары оказались хорошо согласованными с парами натуральных орбиталей метода CASSCF, что с химической точки зрения говорит о "локализованных (двухцентровых) связях" в этой системе.
В дальнейшем различные варианты расчета в методе GVB будут различаться по количеству геминалей и обозначаться 1р, Зр, 5р для одной, трех и пяти геминалей, соответственно. Несколько слов следует сказать о выборе базисного набора. Он производился на основании оценок энергии диссоциации молекулы перекиси водорода на два гидроксильных радикала в методе CASSCF для модельного пространства lv. При переходе от базиса 6-31-b+G(p,d) к корреляционно-согласованному трехэкспонентному базису aug-cc-pVTZ оценки улучшаются лишь незначительно (на 3 % от абсолютной величины энергии диссоциации). Оба использованных метода имеют тот недостаток, что, правильно описывая диссоциационный предел, завышают энергии систем с замкнутой оболочкой. Поэтому мы оценили энергию диссоциации перекиси водорода (см. выше) для всех использованных модельных пространств как в методе CASSCF, так и в методе GVB. Результаты ( в процентах от экспериментального значения, составляющего 0.0807 ат. ед. (2.2 эВ)) [97] приведены в табл. 3.1. Из таблицы видно, что расширение пространства CASSCF приводит к систематическому улучшению описания перекиси водорода в равновесной конфигурации, в то время как качество описания методом GVB остается примерно на одном уровне. Однако метод GVB требует намного меньших вычислительных затрат. Поэтому в динамическом моделировании мы проанализируем результаты обоих методов. На основании сравнения всего массива данных, а также с учетом предварительных расчетов для динамического моделирования были выбраны: в CASSCF — модельное пространство 5v (14е/10 с замороженным остовом из Is и 2s орбиталей обоих атомов кислорода), а в методе GVB — расчет с использованием всех 5 геминалей. Результаты вместе с экспериментальными данными сведены в табл. 3.2. Как следует из сравнения результатов расчетов и эксперимента, метод GVB сильно недооценивает прочность связи 00 в молекуле перекиси водорода, что может привести к уменьшению времени жизни промежуточного высоко возбужденного реакционного комплекса с геометрией, близкой к геометрии молекулы перекиси водорода. Метод CASSCF также несколько недооценивает частоту колебаний ОО связи, однако в количественном отношении отклонение значительно меньше. Кроме того, в методе GVB нарушен порядок следования по энергии частот валентных колебаний ОН связей. К сожалению, решение полной динамической задачи превращений ЩО + О и) наталкивается на серьезные трудности.
Нам не удалось подобрать схему спаривания в методе GVB для реагентов, поэтому вся траектория была рассчитана только методом CASSCF. Начальная геометрия был а получена путем выделения фрагмента структуры Н2О -f О из кластера Оз(ЩО)5 [72]. Начальная величина и направление вектора скорости были выбраны на основе результатов моделирования фотодиссоциации индивидуальной молекулы озона. Согласно этим данным, около 50% от энергии превращения Од —» 02(1Ag) -4-0() уносит атом кислорода. Направление его вылета примерно соответствует направлению слабой водородной связи, образованной исходной молекулой озона с кластером воды. Мы оценили величину кинетической энергии атома в 0.018 ат. ед. ( 0.49 эВ). Шаг интегрирования уравнений ядерного движения был выбран в 0.1 фс ( 4.1 ат. ед.). Показано, что длина шага близка к оптимальной — при ее увеличении до 0.15 фс, наблюдаются значительные вариации полной энергии системы (которая должна сохраняться). Три стадии процесса — формирование молекулы 0(0), изомеризация в Н2О2 и распад перекиси водорода — могут быть качественно прослежены на кривой зависимости полной кинетической энергии частиц системы от времени. Соответствующая зависимость приведена на рис. 3.3. Из рисунка можно заключить, что молекула Н20(0) существует в промежутке от 50 до 300 фс, молекула Н2О2 существует вплоть до 450 фс, и в дальнейшем распадется с образованием гидроксильных радикалов. Поскольку известны координаты и скорости частиц системы, можно разложить полную кинетическую энергию каждого из фрагментов ОН на поступательную, колебательную и вращательную компоненты. Нами была написана программа, вычисляющая эти компоненты на основании данных динамического расчета. На рис. 3.5 приведены кинетическая часть колебательной энергии и вращательная энергия фрагментов ОН. Отметим, что продукты реакции "наследовали" исходную асимметрию, свойственную реагентам. Так, верх- няя половина рисунка отвечает тому гидроксильному радикалу, который образовался из налетающего атома О, в то время как нижняя половина — тому радикалу ОН, фрагмент которого наличествовал в исходной молекуле воды. Колебательные распределения характеризуются относительно узким спектром доступных колебательных волновых чисел продуктов: первый гидроксил имеет среднее значение кинетической части колебательной энергии в ( )0.65 колебательного кванта ( возникла из того, что для гармонического осциллятора Т = V = у), т.е. следует ожидать распределения с примерно равной долей ОН радикалов в состоянии си = 0ии = 1, а второй гидроксил имеет очень низкое значение колебательного возбуждения, что говорит о преимущественном заселении колебательного уровня с v = 0.
