ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ МНОГОСТАДИЙНЫХ БИМОЛЕКУЛЯРНЫХ ФОТОРЕАКЦИЙ В ЖИДКОСТИ: АЛГОРИТМЫ И ПРОГРАММНЫЕ СРЕДСТВА Текст научной статьи по специальности «Физика»

www.volsu.ru

DOI: https://doi.org/10.15688/mpcm.jvolsu.2021.1.4

УДК 519.688

Дата поступления статьи: 15.01.2021 Дата принятия статьи: 26.02.2021

ББК 22.193

см о

CS

03

<J

m о К ■я о со

е

©

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ МНОГОСТАДИЙНЫХ БИМОЛЕКУЛЯРНЫХ ФОТОРЕАКЦИЙ В ЖИДКОСТИ: АЛГОРИТМЫ И ПРОГРАММНЫЕ СРЕДСТВА1

Доктор физико-математических наук,

профессор кафедры информационных систем и компьютерного моделирования,

Волгоградский государственный университет

serguei.feskov@volsu.ru

https://orcid.org/0000-0001-5020-8211

просп. Университетский, 100, 400062 г. Волгоград, Российская Федерация

Аннотация. В работе рассмотрены инструменты численного моделирования многостадийных бимолекулярных фотореакций, обусловленных диффузионной подвижностью молекул-реагентов в вязких растворах. Для описания этих процессов использован аппарат дифференциальной теории встреч (differential encounter theory, DET), расширенной для учета когерентной динамики отдельных степеней свободы (например, электронных спинов реагентов и интермедиатов). Модель включает в себя диффузию реагентов в среде и многостадийные / многоканальные физико-химические процессы, протекающие как в отдельных молекулах, так и столкновительных комплексах. Предложены методы численного решения модельных уравнений, связанных с расчетом эволюционных операторов. Подробно представлен алгоритм вычисления квантового пропагатора для матрицы плотности на основе расщепления Троттера. Разработан программный пакет, ориентированный на численное моделирование многостадийных фотореакций в рамках предложенных подходов. Представлена структура основных компонентов пакета, приведены примеры их использования для построения программной модели реакции. Проведено тестирование, показавшее соответствие численных результатов и точных решений модельных уравнений в частных случаях. На примере фотореакции с участием 9,10-диметилантрацена и 1,3-дицианобензола в растворе ацетонит-рила показаны основные процедуры конфигурирования и численного исследования многостадийных бимолекулярных фотопроцессов. Показана важность учета когерентной эволюции электронных спинов реагентов на радикальной стадии реакции.

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

Сергей Владимирович Феськов

Введение

Многостадийные химические реакции с участием мобильных реагентов в жидкости лежат в основе некоторых важных биологических процессов, таких, например, как дыхание и фотосинтез. Изучение этих реакций важно для понимания их физических механизмов, а также для некоторых практических целей, например, разработки методов диагностики и лечения некоторых видов заболеваний, поиска способов управления химическими процессами, создания эффективных устройств преобразования и хранения солнечной энергии. Среди множества бимолекулярных процессов особый интерес представляют реакции фотоиндуцированного межмолекулярного переноса электрона (ФММ-ПЭ), что обусловлено возможностями их экспериментального исследования с высоким временным разрешением [12; 14; 17]. Объем экспериментов в этой области растет, но возможности интерпретации спектроскопических данных остаются ограниченными. В большинстве случаев нестационарные спектры дают только косвенную информацию о кинетике отдельных стадий. Эффективным в этом случае может быть использование подходов, основанных на сопоставлении результатов натурного и численного экспериментов [1;8]. Как показывают примеры, такой подход помогает определить неизвестные энергетические и динамические характеристики фотосистемы и восстановить картину реакции [6; 7; 9; 10; 13].

Одной из особенностей ФММПЭ является образование радикал-ионных пар на некоторых этапах реакции, к примеру, при фотохимическом разделении зарядов Мх + + М2 ^ М* + М2 ^ Мг+ + М2-, где молекула Мх выступает в качестве фотосенсибилизатора и донора электрона, а М2 — в качестве акцептора. Наличие у М[+ и М2-неспаренных электронов и сильное обменное взаимодействие между радикалами ведет к спиновой селективности рекомбинации зарядов, то есть зависимости скорости рекомбинации от спинового состояния реагентов. Синглетная ионная пара (с суммарным спином ^ = 0) в таких процессах, как правило, является реакционноспособной, тогда как рекомбинация триплетных пар (в = 1) оказывается фактически запрещена [18]. Вместе с тем спиновое состояние пары может изменяться на временах, сопоставимых с временем диффузионной встречи реагентов в жидкости. Среди физических механизмов, обеспечивающих «переворот» электронного спина молекулы, часто выделяют Дд-механизм во внешнем магнитном поле, сверхтонкое взаимодействие с ядерными спинами, релаксационную спиновую динамику [18]. Учет спиновых и магнитных эффектов необходим для построения корректной модели многостадийного переноса электрона.

Интерес к когерентным эффектам в фотохимии обусловлен также новыми экспериментальными данными о реакциях с участием 9,10-диметилантрацена и 1,3-дициано-бензола (фталонитрила) в полярных растворителях [13]. В этих реакциях, в частности, исследовалось влияние среды и постоянного внешнего магнитного поля на кинетику и квантовый выход триплетных продуктов рекомбинации зарядов. Оказалось, что некоторые детали кинетических кривых, зафиксированные в растворителе с высокой полярностью (ацетонитрил), могут быть объяснены только при учете когерентной гамильтоновой спиновой эволюции ион-радикальных пар. С другой стороны, кинетическое описание синглет-триплетных переходов не дало удовлетворительного результата [13].

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

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

1. Модель

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

F + Q — F* + Q ^ F*...Q — F•+...Q~. (1)

Здесь и далее мы обозначим символом F*...Q столкновительный комплекс, образованный при контактном взаимодействии частиц F* и Q в жидкости, а символами F'+ и Q'- — ион-радикалы с неспаренными электронами на внешних орбитах. Будем считать, что при t = 0 молекулы флуорофора и тушителя равномерно распределены по всему объему жидкости и их концентрации известны: [F*] = с/, [Q\ = cq. Электронное состояние пары реагентов F* + Q, находящихся на расстоянии г друг от друга, обозначим №о>.

Учтем далее, что при разделении зарядов по схеме (1) ионные пары могут образовываться в различных электронных и спиновых состояниях |"фп> (п = 1,N). Выделение этих подсистем необходимо для корректного описания процессов заселения энергетических уровней Еп, так как механизмы электронных переходов вида |^0> ^ |^га> для разных п могут сильно различаться. Используем математический аппарат дифференциальной теории встреч [4; 5] и запишем уравнение для кинетики N(t) возбужденного состояния в виде

N(t) = -cqkCS(t)N(t) - N(t)/tfi, N(0) = 1, (2)

где Tfi — среднее время жизни молекулы F* относительно излучательного перехода в основное состояние (флуоресценции), а kcs(t) — скорость тушения флуоресценции по механизму переноса электрона на молекулу Q. Величина kcs(t) в уравнении (2) может быть записана в виде суммы по нескольким каналам ФРЗ

kcs(t) = £ Ш(t) = W W^(r)v(r,t) d3r. (3)

n n

Функция v(r,t) в этом соотношении описывает распределение возбужденных пар по расстояниям между молекулами; wj(^)(r) — парциальная скорость реакции для пары реагентов, находящихся на фиксированном расстоянии г друг от друга. Интегрирование в уравнении (3) формально выполняется по всему пространству, однако в реальных условиях оно ограничено достаточно узкой областью вблизи контактного радиуса, где W^(r) заметно отлична от нуля (ширина этой области — порядка одного ангстрема).

Динамика распределения пар в растворе v(r,t) может быть найдена из уравнения

= - Е wcD (r)v(r, t) + Lv(r, t), v(r, 0) = 1, (4)

где Ь — оператор Смолуховского, описывающий взаимную радиальную диффузию реагентов в жидкости

£ = " Р Г' (5)

г2иг иг

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

= Е Ц)Рп - [н, ц(г,г)] + 2(ж(г), ц(г,г)} + щг,г). (6)

п

Здесь 1-е слагаемое в правой части описывает генерацию ионных пар в ходе разделения зарядов, 2-е и 3-е слагаемые представляют соответственно когерентную и некогерентную динамику квантовых переходов, последнее слагаемое учитывает диффузионную мобильность реагентов в жидкости. Квадратными и фигурными скобками в уравнении (6) обозначены соответственно коммутатор и антикоммутатор двух операторов; Н — эффективный спиновый гамильтониан ионной пары, Рп — оператор проектирования в п-е состояние. В уравнении (6) учтем кулоновское взаимодействие между заряженными частицами, в результате чего оператор Смолуховского Ь преобразуется к виду

Ь = #г2 е-и/квТ^еи/к*Т, (7)

г2 иг иг

где и (г) — кулоновский потенциал; кв — постоянная Больцмана; Т — температура среды и И — коэффициент взаимной диффузии частиц.

Уравнения для у(г,Ь) и ц(г, ¿) удовлетворяют отражающим граничным условиям при г = а (а — контактный радиус двух реагентов)

Jvv

Anr2D ^ ог

0, г, t) =4 ш 2 De-и/к*Т еи/квТ£

г=а ОГ

0. (8)

Здесь Jv и ^ — диффузионные потоки нейтральных реагентов и ион-радикалов через сферу радиуса .

Оператор Ш удобно представить в виде суммы

™(Г) = Е пп' (г)Рп + Щп'п(г)Рп', (9)

п,п'

где Шпп/ (г) — скорость распада состояния |-фп) в результате необратимого перехода в состояние |*фп/), а Шп>п(г) — скорость заселения конечного состояния. Сохранение общего числа частиц в системе обеспечивается за счет выполнения условия №пп! (г) = = ^п'п(г). Населенности состояний |-фп) рассчитываются по формуле

Рп(Ъ = Тг(Рп,ц(г, 1)}£г. (10)

Система уравнений (2)—(10) полностью определяет кинетику многостадийного ФММПЭ с учетом когерентной эволюции электронных спинов в рамках дифференциальной теории встреч. В качестве основных «входных» параметров используемой модели выделим следующие:

г=п

г=а

г=а

1) глобальные параметры системы, включая температуру растворителя, его вязкость и полярность, среднее время флуоресценции Р* и др.;

2) число квантовых состояний интермедиатов |*фп), принимающих участие в реакции;

3) характеристики пространственного движения каждого типа реагентов, включая коэффициент взаимной диффузии И(г), контактный радиус а, потенциал межчастичного взаимодействия и (г);

4) набор квантовых состояний, выступающих в качестве первичных продуктов ФРЗ, а также скорости переходов в эти состояния в виде функций

5) характеристики кинетических (некогерентных) переходов между состояниями |"фп) в форме функций Шпп/(г) и с указанием типа перехода (обратимый/необратимый, одно-/двунаправленный и др.);

6) матрица спинового гамильтониана Н в выбранном базисе;

7) параметры численной схемы, включая границы расчетной области, параметры пространственной и временной сетки и др.

Основными результатами численных расчетов являются населенности квантовых состояний системы N(1) и Рп(1), а также нестационарные спектры люминесценции и поглощения, вычисляемые на основе этих населенностей и функций Vг, Ь) и ц(г, ¿). Полученные в ходе расчетов данные в дальнейшем могут использоваться для анализа результатов спектроскопических экспериментов.

Компьютерное моделирование в рамках указанной модели опирается на численное решение дифференциальных уравнений (4) и (6) для плотностей пар реагентов и интермедиатов в различных квантовых состояниях. С вычислительной точки зрения наибольшую сложность представляет численное решение уравнения для ц(г, ¿), так как оно включает в себя разнородные компоненты: с одной стороны, диффузию в заданном потенциале, с другой — гамильтонову динамику квантовых степеней свободы и некогерентные электронные переходы. Далее в этой статье основное внимание уделяется численному решению уравнения (6). Отметим, что расчет собственно населенностей N(1) и Рп(1) на основе (3) и (10) не представляет сложности, так как в итоге сводится к вычислению пространственных интегралов от рассчитанных ранее функций Vг, Ь) и

Мг, Ь).

2. Численный метод и программная реализация

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

Т ~Тн •Тщ •Ть, (11)

где Тн(г, ¿), Тщ(г, Ь) и Ть(г, Ь) отвечают соответственно за когерентную спиновую динамику, некогерентные электронные переходы и диффузию вдоль координаты г. В общем случае соотношение (11) является приближенным и справедливым только для достаточно малых ¿. Условием применимости разложения следует считать корректный

выбор временного шага АЪ. Его величину можно оценить из соотношения Л^ АЪ ^ 1, где Л^ — верхняя граница спектра эволюционного оператора Т.

По определению, оператор эволюции Т(г, АЪ) позволяет рассчитать матрицу плотности квантовой системы в момент времени Ь + АЬ по ее известному значению в момент Ь

ц(г, г + АЬ) = Т(Г, АЬ)ц(г, г) ~ТН(г, АЬ)Тф(г, М)ТЬ(г, Аг)ц(г, Ь). (12)

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

В данной работе мы подробно рассмотрим метод численного расчета пропагаторов ТН и Т^, принимая во внимание, что вычисление Т^ рассматривалось ранее другими авторами (см., например, [3; 15]).

2.1. Пропагатор чистых состояний

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

* = $ £). (13)

Волновая функция чистого состояния в общем случае представима в форме

ф(г) = а(г)ф! + Ь(1)ф2 = , (14)

где а(Ь) и Ь(Ь) — комплекснозначные коэффициенты, удовлетворяющие условию нормировки |а({)12 + |Ь({)12 = 1. Согласно квантовой теории, эволюция ф(£) задана уравнением Шредингера

Ш ^ = (15)

Для гамильтониана вида (13) известно точное решение уравнения (15), однако в данном случае нас интересует метод расчета оператора эволюции Т({), который был бы удобен для численной реализации и мог бы быть легко обобщен на более сложные квантовые объекты. Запишем формальное решение уравнения Шредингера (15) в виде

ф(1) = ехр (-ф(0) = Т(1)ф(0), (16)

где Т(1) = ехр(—гНЬ/К), а ф(0) — волновая функция системы в начальный момент времени. Таким образом, Т(Ь) — искомый пропагатор двухуровневой системы, представленный в форме матричной экспоненты. В общем случае расчет Т(1) может быть выполнен с помощью разложения экспоненты в ряд Маклорена. Однако этот способ неэффективен с вычислительной точки зрения, так как требует большого объема матричных операций и его производительность быстро падает с ростом размерности матрицы Н. Кроме того, при ограничении числа суммируемых членов ряда оператор Т(1) теряет свойство унитарности, что в свою очередь приводит к таким побочным эффектам, как потеря нормировки волновой функции в ходе компьютерного моделирования.

Рассмотрим альтернативный вариант вычисления оператора T(t), основанный на известной формуле Троттера

е A)t = lim (еМ/mеЫ'тГ , (17)

m^-x V )

где А и В — произвольные (в общем случае некоммутирующие) операторы. Разложим матрицу гамильтониана на диагональный и внедиагональный компоненты

»=$ =(о- D+(¿т)=*>(18)

Согласно формуле Троттера (17)

T(mAt) ~ (Ti(At)T2(At)) ...Т^Т2, (19)

где TT— = exp(—гHlAt/fi) и Т2 = exp(—iH2At/К). Смысл выполненного разложения состоит в том, что матрицы операторов rTl и Т2 могут быть легко вычислены аналитически. Прямое суммирование бесконечных рядов в данном случае дает следующий точный результат

H (At) _ V (-AtУ Нп _ [ex?— EiAt0 ^

Tim _ ¿Д Hl Л 0 exp (-E At ß))

H (Af\_^ Í-iAt V Hrn _( cos(QAt) -i sin (üAt)\

l2(At)-^n{ n ) H _ [-г sin(QAt) cos(QAt) ) '

выраж

моделирования At

(20) (21)

\ П I \— ■/, Siül^üLl CUNIH^Ll I

п=0

Полученные выражения позволяют сформулировать критерий выбора шага численного

At < max(E1/n,E1/%, tt). (22)

Легко видеть, что операторы Т1 и Т2 являются унитарными, а значит предложенная схема консервативна по построению. Метод может быть напрямую обобщен на квантовые системы большей размерности: для этого достаточно использовать в разложении (18) недиагональной части гамильтониана Н такие компоненты Нкк', в которых отличны от нуля только два элемента Нкк' = Нк'к = 0, расположенные симметрично относительно главной диагонали.

2.2. Пропагатор смешанных состояний

Применим теперь технику расщепления к описанию квантовых систем, находящихся в смешанном состоянии и содержащих неэрмитовы компоненты в гамильтониане. Для компактности снова используем двухуровневую модель, в которой когерентная динамика задана оператором (13), а необратимый распад — матрицей IV с отличными от нуля элементами Шц = -Я\ и W22 = -К2. Введем оператор плотности V(t)

v(t) = е р-тьх^I = (Р2$) Р12Ш, (23)

где р!! и р22 — населенности состояний ф! и ф2, удовлетворяющие условию р!! + р22 = 1, а р!2 и р2! — комплексные числа, указывающие на связь между базисными состояниями (когерентности). Для матрицы р имеем следующее уравнение эволюции

f=- Т [й • й+2 •

(24)

Перепишем это уравнение в представлении Лиувилля. Для этого представим р в виде вектора-столбца, например, с помощью преобразования

(Р11 Р12^ \Р21 Р22)

(?1l\ Р12 Р21

\Р22 J

(25)

Уравнение (24) в этом представлении можно записать в форме

^=— *£№)•

(26)

где £ — супероператор (лиувиллиан системы), представляющий собой матрицу размером 4 х 4

£ = h

(-% R1 -Q

+Q 0

-Q

+Q

0

АЕ/Т - iЕR/2 0 +Q

0 -АЕ/Т - iE R/2 -Q

+Q -Q -iR2J

(27)

Здесь мы использовали обозначения АЕ = Е! — Е2 и ЕД = Д! + Д2.

Уравнение Лиувилля (26), так же как и матричное уравнение Шредингера (15), — это система линейных ДУ 1-го порядка, поэтому техника расщепления может быть также применена к расчету квантового пропагатора для вектора р(1). Получаем

Z(At) = e-liAt/h ~ Д е

- i£kA t/h

(28)

где £ — элементы разложения матрицы £, допускающие прямое суммирование бесконечного ряда (аналогично формулам (20), (21)).

Отметим некоторые особенности структуры лиувиллиана £ по выражению (27). К примеру, недиагональная часть матрицы £р содержит только недиагональные элементы Н, то есть параметры, определяющие когерентную динамику переходов в системе. При этом матрица £р сохраняет симметричность, что дает возможность использовать формулу (21) для вычисления элементарного пропагатора И^АЬ). Другая особенность связана с появлением мнимых слагаемых на главной диагонали, содержащих параметры Д!, Д2 и ЕД. Как легко видеть из формулы (20), такие слагаемые нарушают унитарность пропагатора, а следовательно, ведут к потере нормировки матрицы плотности. Однако это свойство является прямым следствием используемой модели, а именно наличием в системе некогерентных процессов, описываемых матрицей IV.

3. Программный код FLUT

В этом разделе мы представляем программную разработку, предназначенную для численного решения уравнений математической модели (2)-(10). Основным результатом расчетов являются кинетики населенностей флуоресцентного и ионных состояний F*, F+, Q-, а также парные функции распределения реагентов и интермедиатов в растворителе v(г, t) и \inn(г, t). Разработанный код предполагается использовать для анализа экспериментальных данных о многостадийных бимолекулярных фотореакциях, включая процессы в биологических системах. Методика фитирования нестационарных спектров люминесценции/поглощения с фемто- и пикосекундным временным разрешением ранее уже применялась для исследования подобных систем и определения неизвестных значений параметров реагентов и среды [1; 6—10].

Особенностью пакета FLUT является учет когерентной динамики отдельных степеней свободы молекулярной системы в рамках теории матрицы плотности и уравнения фон Неймана. Эта особенность позволяет в частности моделировать квантовые эффекты, возникающие при спин-спиновом взаимодействии флуорофора и тушителя в состоянии ион-радикалов (F + и Q-). Также программная модель включает в себя каналы равновесного (термического) и неравновесного (горячего) переноса заряда между реагирующими частицами, характеристики импульса фотовозбуждения, внутреннюю конверсию и колебательную релаксацию. Перенос электрона «через растворитель» и влияние диэлектрической релаксации среды (динамический эффект растворителя) моделируются в рамках теории Зусмана [19] и теории возмущений по величине донорно-акцепторной связи.

Код FLUT организован в виде системы классов С++, представляющих базовые понятия теории межмолекулярного переноса заряда (см. рис. 1). В качестве компонента верхнего уровня выступает объект класса TModel, который также служит контейнером для объектов более низких уровней. Такая архитектура позволяет построить программную модель ФММПЭ, как композицию объектов нескольких существующих классов внутри агрегата TModel. С точки зрения пользователя, конфигурирование программы заключается в добавлении необходимых компонентов в TModel и установлении связей между ними. Рассмотрим эти компоненты более подробно.

Одним из основных понятий предметной области является понятие диффузионного пространства (diffusion space [11]), которое в данном случае используется для описания распределения пар молекул-реагентов в жидкости. Диффузионных пространств в модели может быть несколько, так как каждое из них привязано к определенному состоянию реагентов и интермедиатов. В общем случае состояние представлено набором электронных, колебательных и спиновых квантовых чисел. Понятие диффузионного пространства реализовано классом TDS, показанным на рисунке 1. Элементы данных (поля класса) дают пользователю возможность задать граничные условия для оператора диффузии, начальное распределение цпп(г, 0), потенциал межчастичного взаимодействия U(г), околоконтактную зависимость коэффициента D(r) и некоторые другие. Параметры растворителя входят в модель несколькими способами, в частности, вязкость влияет на коэффициент пространственной диффузии реагентов, а полярность отвечает за экранирование молекулярного взаимодействия и энергетику переноса заряда (через энергию реорганизации и свободную энергию образования продукта).

Для описания переходов между квантовыми состояниями в уравнении (6) используются функции связи Wnn/(г). В коде FLUT понятие о связи реализовано с помощью

класса TCoupling. Выделяется несколько типов связей по отношению к типам реагентов и продуктов, среди которых отметим (а) связи между диффузионными пространствами (например, описывающие высвечивание фотона F*...Q ^ F...Q и другие мономолекулярные процессы в растворе) и (б) связи между диффузионным пространством и локализованным состоянием пары — продуктом реакции. К связям 2-го типа можно отнести, например, ассоциацию двух молекул при их столкновении в растворе, или наоборот, распад молекулы на две части. Для моделирования процессов различных типов в разработанном пакете реализована система родственных классов, часть которой показана на рисунке 2. Элементы этой структуры являются наследниками базового класса TCoupling, и реализуют различные специфические связи в системе. Краткое описание предназначения некоторых из этих элементов дано на рисунке 2.

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

Для тестирования алгоритма (28) нами проведено компьютерное моделирование динамики квантовой системы, включающей три связанных состояния с гамильтонианом

( 0 V \2

H= I V12 0 У2з\ . (29)

\г>1з V23 0 )

Расчеты проводились для матрицы плотности с начальным условием p11(t) = 1. Сравнение с точными решениями уравнения фон Неймана в частных случаях и независимыми кодами показало совпадение в пределах вычислительной погрешности в широком диапазоне параметров vnn>. Пример тестового сравнения показан на рисунке 3, где изображена динамика диагональных элементов (панель А) и недиагональных элементов (панель Б) матрицы плотности для случая v12 = v23 = 10 нс-1, v13 = 0. Результаты численного моделирования показаны точками, точное решение — сплошными линиями. Расчеты выполнены с шагом At = 1 пс.

4. Фотохимия слабых растворов 9,10-диметилантрацена и фталонитрила

в полярной среде

Рассмотрим процедуру конфигурирования программного пакета на примере реакции тушения флуоресценции фотовозбужденных молекул 9,10-диметилантрацена (DMeA) по механизму переноса электрона на фталонитрил (PN) в полярном растворителе (ацетонитрил, ACN). Экспериментальные исследования этой фотореакции методом накачки-зондирования проводились недавно в Женевском университете (см. [13]). Одним из результатов, полученных в эксперименте, стало сильное влияние внешнего магнитного поля на интенсивность спектральной линии, соответствующей триплетному состоянию 3DMeA*. Такое влияние может быть объяснено, с одной стороны, сильной зависимостью скорости рекомбинации зарядов от спинового состояния ион-радикальной

пары 2ОМеА'+...2Р№ , а с другой — эволюцией электронных спинов молекул 2ОМеА'+ и 2Р№- в ходе диффузионных блужданий в жидкости.

Программная модель фотореакции

TModel

Class

0 Fields

0 couplings

^ description

* dt

modeljiame

M N_norm_dev

# num_couplings

^ num_states

^ states

^ time

3 Methods

© ~TModel

© add_coupling

© add_state

© advance_timestep

© calc_populations

© calc_propagators

© ¡nit

© load_configuration_from

© save_configuration_to

© TModel

Оператор эволюции для диффузионного пространства

TDSPropagator

Class

0 Fields

a

^ a_balance В

# ds_nu

Ф numjntervals ** num_polynomi... ^ Omega

# P & size

# X

0 Methods © -TDSPropagator © calc © get_balance © in it © map

© TDSPropagator

TGIobals Л

Class

0 Fields Ф С

* dT 0 kT

в mapjen 0 num_intervals

* Rmax

* Rmin

* skipped_points

* T

* tau_D

* tau_L 0 Methods

© load_from © save_to © TGIobals

TMapping

Class

0 Fields

interval_c enter interval_length mapjen

« M

0 Rmax

* Rmin 0 Methods

© --TMapping

© get_mapped_R

© get_mapped1_R

© init_with

© TMapping

Диффузионное пространство

TDS «

Class

0 Fields

# bcsjnner & bcs_outer

# D_f unction

# DO

M description

# ics

# index M K_ass ** K_dis

# K_esc

# K_rec ^ name

# num_intervals

# population

# propagator 0 QO

** Q1

0 Rc

^ V_f unction 0 Methods

© ~TDS

© calc_propagator

© get_ds_nu

© get_ds_populations_all

© get_escaped

© get_Q

© get_recombined

© get_Rho

© get_survived

© get_transition_rate

© get_V3D

© inii

© init_profile

© load_configuration_from

© propagate

© save_configuration_to

© set_Q

© TDS

public

TQuantum3x3DS «

Class ->TDS

0 Fields

Ф A

M omega_A

^ omega_D

0 rho_0

& rho_1

* U 0 Methods

© ~TQuantum3x3DS

© calc_propagator

© calc_quantum_propagator

© in it

© init_profile

© propagate

© TQuantu m3x3 DS

Квантовое (когерентное) диффузионное пространство

Рис. 1. Основные классы программной модели

Рис. 2. Иерархия классов, представляющая различные типы связей между состояниями

в системе

Рис. 3. Динамика в трехуровневой квантовой системы (29). Точки — результаты численного моделирования с помощью схемы (28), сплошные линии — точный результат. На панели А показаны населенности базисных состояний р„„(£), на панели Б — когерентности

6„„' = г/\/2 (р„„' - р„'„)

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

1т-* , 1

F "+1 Q диффузия [ 1р * 1Q]

h v

(t=0)

1F+1Q

разделение зарядов

спиновая конверсия

Д

[1F...1Q]

1 [2F+...2Q-] Л,

3 [2F+... 2Q"] ^ ]

^7<s) / рекомбинация зарядов (синглетная и триплетная)

столкновительныи комплекс

диффузия 2р+ + 2 Q"

д свободные с ионы

t

Рис. 4. Схема фотохимических реакций с участием 9,10-диметилантрацена (флуорофор Р) и фталонитрила (тушитель Q) в вязком полярном растворителе. Синглет-триплетная конверсия ион-радикальных пар обусловлена Дд-механизмом и сверхтонким взаимодействием

с ядерными спинами

Основной этап начинается с определения квантовых состояний частиц, принимающих участие в реакции. В нашем случае начальным состоянием будем считать син-глетное возбужденное состояние пары [1Р*...1Q], где Р и Q — молекулы флуорофора и тушителя. Согласно схеме на рисунке 4, промежуточными состояниями являются также синглетное и триплетное состояния ион-радикальной пары, образующиеся при фотоиндуцированном разделении зарядов. На схеме они обозначены как 1[2Р+ ...2Q-] и 3[2Р+...^-]. Триплетное состояние с нейтральными реагентами [3Р...^] и основное состояние [1Р...1Q] также принимают участие в этих процессах, и заселяются в ходе

рекомбинации зарядов в ионных парах.

Во всех указанных состояниях партнеры находятся в несвязанных парах и совершают диффузионное движение в жидкости. Для описания динамики распределения пар реагентов введем соответствующие диффузионные пространства. В частности, для син-глетного возбужденного состояния программный код C++ может выглядеть следующим образом (здесь statel — новый объект класса TDS)

TDS *state1 = new TDS("S1", "excited"), state1->D0 = D;

state1->bcs_inner = bc_reflecting; state1->bcs_outer = bc_constant_c; state1->ics = ic_uniform; state1->population = 1.0; model->add_state(state1);

// создаем новый объект TDS

// коэффициент диффузии

// внутренняя граница

// внешняя граница

// начальное распределение

// начальная населенность

// добавляем объект в модель

В соответствии с уравнением (4), в качестве начального условия для Vг, Ь) задаем равномерное распределение; внутреннюю границу (на контактном радиусе а) делаем отражающей, а внешнюю границу задаем прозрачной для частиц. Отметим также, что реагенты в возбужденном состоянии остаются физически нейтральными (заряд каждой частицы равен нулю), поэтому для данного состояния не вводим дополнительного потенциала и (г) взаимодействия между частицами.

Аналогичным образом может быть сконфигурировано квантовое диффузионное пространство для ионных пар. В общем случае размерность этого пространства связана с числом связанных квантовых состояний, и в случае радикальной пары (2 спина 5 = 1/2) и одного ядерного спина ( 5 = 1/2) это число равняется 8. Таким образом, матрица плотности для этого случая имеет размеры 8 х 8, а лиувиллиан системы — 64 х 64. Вместе с тем, используя свойства симметрии гамильтониана Н, можно существенно уменьшить размерность решаемой задачи [13].

state2 = new TQuantumDS("CS", "radical-ions"); state2->D0 = D;

state2->bcs_inner = bc_reflecting; state2->bcs_outer = bc_absorbing; state2->ics = ic_zero; state2->V_function = V_func; model->add_state(state2);

Отметим здесь изменение начальных и граничных условий (на внешней границе), а также поле V_function, устанавливающее потенциал взаимодействия заряженных частиц — ионов F + и Q-. Используемый компонент V_func — это объект класса для расчета экранированного кулоновского взаимодействия, который, однако, может настраиваться дополнительно.

Коснемся теперь конфигурирования объектов TCoupling для установки связей в программной модели. Следующий код формирует компонент для описания переноса электрона из возбужденного состояния в состояние ион-радикальных пар

cs_coupling = new TC_UTCoupling("charge separation", "S_1->CS"); cs_coupling->bind(state1, state2, W_func); cs_coupling->set_fluorescence_destination(state5); model->add_coupling(cs_coupling);

Функция bind устанавливает связь между диффузионными пространствами statel (возбужденные пары) и state2 (ионы) через компонент W_func, который представляет соответствующую функцию связи Wcs(О- Выбор типа объекта W_func и спецификация его полей позволяет использовать различные модельные зависимости Wcs (г), включая такие, которые учитывают динамический эффект растворителя, адиабатическое насыщение скорости переноса электрона и другие эффекты. Настройка объекта cs_coupling требует также указания направления излучательных переходов из возбужденного состояния (функция set_fluorescence_destination).

С помощью дополнительных объектов-наследников класса TCoupling аналогичным образом настраиваются программные компоненты для процесса внутренней конверсии, а также для рекомбинации зарядов в синглетных и триплетных ионных парах (см. общую схему реакции). Глобальные параметры системы заносятся в поля класса TGlobal, характеристики неравномерной пространственной сетки вдоль координаты г задаются в компоненте TMapping.

Ключевыми элементами всех классов диффузионных пространств и связей является функция calc_propagator(), выполняющая расчет матрицы эволюционного оператора в соответствии с численными схемами, описанными выше. Вызов этих функций производится неявно при инициализации программной модели model->init(). Собственно численное моделирование выполняется путем последовательного вызова функции model->advance_timestep(), которая производит расчет нового состояния молекулярной системы на заданном временном шаге At. Выгрузка результатов моделирования — на-селенностей квантовых состояний и профилей распределения частиц в пространстве — производится через вызов соответствующих методов класса TModel.

Рис. 5. Сопоставление результатов компьютерного моделирования (сплошные линии) и экспериментальной кинетики (точки) бимолекулярной реакции с участием 9,10-диметилантрацена и 1,3-дицианобензола согласно схеме 4. Светло-серым цветом (красным

в электронной версии статьи) обозначены результаты в сильном (насыщающем) магнитном поле, темно-серым (синим в электронной версии статьи) — в отсутствие поля [13]. Панель А — кинетика ион-радикальных пар, панель Б — кинетика нейтрального триплетного состояния

Рассмотрим некоторые из результатов моделирования, изображенные на рисунке 5. На левой панели (А) показана кинетика Р2(£) накопления и распада ионных пар, обусловленная конкуренцией между разделением зарядов и их рекомбинацией. Величина Р2(¿) в данном случае рассчитывалась как сумма населенностей синглетного и три-

плетного состояний !\2Р+ ...2(^-} и 3[2F+...2Q-]. Светло-серым и темно-серым цветом (красным и синим в электронной версии статьи соответственно) выделены результаты в насыщающем магнитном поле и в отсутствие поля. Для сравнения на рисунке 5 показаны данные экспериментов, в которых кинетика Р2(1) фиксировалась через изменение интенсивности спектральной линии, характерной для катион-радикала ЭМеА+. Аналогичные результаты для кинетики Р3(Ь) нейтральных триплетных пар [3^..!(] показаны на правой панели (Б).

Отметим хорошее соответствие численных и опытных данных на панели А, а также слабую зависимость Р2(£) от магнитного поля. Как показывают расчеты, основным каналом распада ионного состояния является синглетная рекомбинация, тогда как три-плетный канал и спиновые переходы оказывают на него незначительное влияние. Однако магнитное поле сильно влияет на скорость заселения триплетных продуктов рекомбинации, что хорошо заметно на рисунке 5Б. Действительно, образование пар [3^..!(] возможно только при условии достаточно быстрых синглет-триплетных переходов, протекающих на временах, сравнимых с длительностью диффузионной встречи ионов в жидкости. Удовлетворительное совпадение расчетных и экспериментальных результатов для Р3(1) наблюдается только на начальных этапах реакции (длительностью примерно до 10 нс). Отметим, тем не менее, что воспроизведение этого начального этапа стало возможно только при использовании когерентной модели спиновой конверсии, в которой были учтены одновременно А -механизм и сверхтонкие взаимодействия.

Заключение

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

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

Описана структура пакета, приведены результаты ряда тестовых расчетов. Проведено численное исследование фотохимических процессов с участием 9,10-диметилантрацена и фталонитрила в полярном растворителе, результаты расчетов соотнесены с известными экспериментальными данными.

ПРИМЕЧАНИЕ

1 Работа выполнена при поддержке РФФИ и Администрации Волгоградской области (проект 19-47-340020).

СПИСОК ЛИТЕРАТУРЫ

1. Вотэ, Э. Объединение теории и эксперимента для понимания процесса сверхбыстрого фотоиндуцированного переноса заряда [На англ.] / Э. Вотэ // Математическая физика и компьютерное моделирование. — 2020. — Т. 23, № 2. — C. 91-99. — DOI: 10.15688/mpcm.jvolsu.2020.2.8.

2. Феськов, С. В. Алгоритмы численного моделирования сверхбыстрого фотоиндуцированного межмолекулярного переноса заряда в жидкостях / С. В. Феськов, С. С. Хохлова // Вычислительные методы и программирование. — 2020. — Т. 21. — C. 27-40. — DOI: 10.26089/NumMet.v21r103.

3. Allen, M. P. Computer Simulation of Liquids / M. P. Allen, D. J. Tildesley. — New York : Clarendon Press, 1989. — 408 p.

4. Burshtein, A. I. Non-Markovian theories of transfer reactions in luminescence and chemiluminescence and photo- and electrochemistry / A. I. Burshtein // Adv. Chem. Phys. — 2004. — Vol. 129. — P. 105-418. — DOI: 10.1002/047168077X.ch3.

5. Burshtein, A. I. Unified theory of photochemical charge separation / A. I. Burshtein // Adv. Chem. Phys. — 2000. — Vol. 114. — P. 419-587. — DOI: 10.1002/9780470141731.ch6.

6. Feskov, S. V. Double-channel photoionization followed by geminate charge recombination/Separation / S. V. Feskov, A. I. Burshtein // J. Phys. Chem. A. — 2009. — Vol. 113. — P. 13528-13540. — DOI: 10.1021/jp901863t.

7. Feskov, S. V. Kinetics of nonequilibrium electron transfer in photoexcited ruthenium(II)-cobalt(III) complexes / S. V. Feskov, A. O. Kichigina, A. I. Ivanov // J. Phys. Chem. A. — 2011. — Vol. 115. — P. 1462-1471. — DOI: 10.1021/jp108607t.

8. Feskov, S. V. Non-equilibrium effects in ultrafast photoinduced charge transfer kinetics / S. V. Feskov, V. A. Mikhailova, A. I. Ivanov // J. Photochem. Photobiol. C. — 2016. — Vol. 29. — P. 48-72. — DOI: 10.1016/j.jphotochemrev.2016.11.001.

9. Hot recombination of photogenerated ion pairs / V. Gladkikh, A. I. Burshtein, S. V. Feskov, A. I. Ivanov, E. Vauthey // J. Chem. Phys. — 2005. — Vol. 123. — P. 244510. — DOI: 10.1063/1.2140279.

10. Ionkin, V. N. Numerical simulations of ultrafast charge separation dynamics from second excited state of directly linked zinc-porphyrin-imide dyads and ensuing hot charge recombination into the first excited state / V. N. Ionkin, A. I. Ivanov // J. Phys. Chem. A. — 2009. — Vol. 113. — P. 103-107. — DOI: 10.1021/jp806589m.

11. Krissinel', E. B. Spherical symmetric diffusion problem / E. B. Krissinel', N. Agmon // J. Comput. Chem. — 1996. — Vol. 17. — P. 1085-1098. — DOI: 10.1002/(SICI)1096-987X(19960715)17:9<1085::AID-JCC1>3.0.CO;2-O.

12. Kumpulainen, T. Ultrafast elementary photochemical processes of organic molecules in liquid solution / T. Kumpulainen, B. Lang, E. Vauthey // Chem. Rev. — 2017. — Vol. 117. — P. 10826-10939. — DOI: 10.1021/acs.chemrev.6b00491.

13. Magnetic field effect on ion pair dynamics upon bimolecular photoinduced electron transfer in solution / S. V. Feskov, M. V. Rogozina, A. I. Ivanov, A. Aster, M. Koch, E. Vauthey

// J. Chem. Phys. - 2019. - Vol. 150, № 2. - P. 024501. - DOI: 10.1063/1.5064802.

14. Martin, C. Ultrafast and fast charge separation processes in real dye-sensitized solar cells / C. Martin, M. Ziolek, A. Douhal // J. Photochem. Photobiol. C. - 2016. - Vol. 26. -P. 1-30. - DOI: 10.1016/j.jphotochemrev.2015.12.001.

15. Pines, E. Geminate recombination in excited state proton transfer reactions: Numerical solution of the Debye - Smoluchowski equation with backreaction and comparison with experimental results / E. Pines, D. Huppert, E. Agmon // J. Chem. Phys. - 1988. -Vol. 88. - P. 5620-5630. - DOI: 10.1063/1.454572.

16. Potovoi, V. V. Theory of non-thermal electron transfer / V. V. Potovoi, A. I. Ivanov // Chem. Phys. - 1999. - Vol. 247. - P. 245-259. - DOI: 10.1016/S0301-0104(99)00197-4.

17. Rosspeintner, A. Ultrafast photochemistry in liquids / A. Rosspeintner, B. Lang, E. Vauthey // Annu. Rev. Phys. Chem. - 2013. - Vol. 64. - P. 247-271. - DOI: 10.1146/annurev-physchem-040412-110146.

18. Spin Polarization and Magnetic Effect in Radical Reactions / K. M. Salikhov, Yu. N. Molin, R. Z. Sagdeev, A. L. Buchachenko. - Amsterdam : Elsevier, 1984. - 419 p.

19. Zusman, L. D. Outer-sphere electron transfer in polar solvents / L. D. Zusman // Chem. Phys. - 1980. - Vol. 49, № 2. - P. 295-304. - DOI: 10.1016/0301-0104(80)85267-0.

REFERENCES

1. Vote E. Obyedinenie teorii i eksperimenta dlya ponimaniya protsessa sverkhbystrogo fotoindutsirovannogo perenosa zaryada [Na angl.] [Combining Theory and Experiment for Understanding of Ultrafast Photoinduced Charge-Transfer Processes]. Matematicheskaya fizika i kompyuternoe modelirovanie [Mathematical Physics and Computer Simulation], 2020, vol. 23, no. 2, pp. 91-99. DOI: 10.15688/mpcm.jvolsu.2020.2.8.

2. Feskov S.V., Khokhlova S.S. Algoritmy chislennogo modelirovaniya sverkhbystrogo fotoindutsirovannogo mezhmolekulyarnogo perenosa zaryada v zhidkostyakh [Algorithms for the Computer Simulation of Ultrafast Photoinduced Intermolecular Charge Transfer in Liquids]. Vychislitelnye metody i programmirovanie [Numerical Methods and Programming], 2020, vol. 21, pp. 27-40. DOI: 10.26089/NumMet.v21r103.

3. Allen M.P., Tildesley D.J. Computer Simulation of Liquids. New York, Clarendon Press, 1989. 408 p.

4. Burshtein A.I. Non-Markovian Theories of Transfer Reactions in Luminescence and Chemiluminescence and Photo- and Electrochemistry. Adv. Chem. Phys., 2004, vol. 129, pp. 105-418. DOI: 10.1002/047168077X.ch3.

5. Burshtein A.I. Unified Theory of Photochemical Charge Separation. Adv. Chem. Phys., 2000, vol. 114, pp. 419-587. DOI: 10.1002/9780470141731.ch6.

6. Feskov S.V., Burshtein A.I. Double-Channel Photoionization Followed by Geminate Charge Recombination/Separation. J. Phys. Chem. A, 2009, vol. 113, pp. 13528-13540. DOI: 10.1021/jp901863t.

7. Feskov S.V., Kichigina A.O., Ivanov A.I. Kinetics of Nonequilibrium Electron Transfer in Photoexcited Ruthenium(II)-Cobalt(III) Complexes. J. Phys. Chem. A, 2011, vol. 115, pp. 1462-1471. DOI: 10.1021/jp108607t.

8. Feskov S.V., Mikhailova V.A., Ivanov A.I. Non-Equilibrium Effects in Ultrafast Photoinduced Charge Transfer Kinetics. J. Photochem. Photobiol. C, 2016, vol. 29, pp. 48-72. DOI: 10.1016/j.jphotochemrev.2016.11.001.

9. Gladkikh V., Burshtein A.I., Feskov S.V., Ivanov A.I., Vauthey E. Hot Recombination of Photogenerated Ion Pairs. J. Chem. Phys., 2005, vol. 123, pp. 244510. DOI: 10.1063/1.2140279.

10. Ionkin V.N., Ivanov A.I. Numerical Simulations of Ultrafast Charge Separation Dynamics From Second Excited State of Directly Linked Zinc-Porphyrin-Imide Dyads and Ensuing Hot Charge Recombination Into the First Excited State. J. Phys. Chem. A, 2009, vol. 113, pp. 103-107. DOI: 10.1021/jp806589m.

11. Krissinel' E.B., Agmon N. Spherical Symmetric Diffusion Problem. J. Comput. Chem., 1996, vol. 17, pp. 1085-1098. DOI: 10.1002/(SICI)1096-987X(19960715)17:9<1085::AID-JCC1>3.0.CO;2-O.

12. Kumpulainen T., Lang B., Vauthey E. Ultrafast Elementary Photochemical Processes of Organic Molecules in Liquid Solution. Chem. Rev., 2017, vol. 117, pp. 10826-10939. DOI: 10.1021/acs.chemrev.6b00491.

13. Feskov S.V., Rogozina M.V., Ivanov A.I., Aster A., Koch M., Vauthey E. Magnetic Field Effect on Ion Pair Dynamics Upon Bimolecular Photoinduced Electron Transfer in Solution. J. Chem. Phys., 2019, vol. 150, no. 2, pp. 024501. DOI: 10.1063/1.5064802.

14. Martin C., Ziolek M., Douhal A. Ultrafast and Fast Charge Separation Processes in Real Dye-Sensitized Solar Cells. J. Photochem. Photobiol. C, 2016, vol. 26, pp. 1-30. DOI: 10.1016/j.jphotochemrev.2015.12.001.

15. Pines E., Huppert D., Agmon E. Geminate Recombination in Excited State Proton Transfer Reactions: Numerical Solution of the Debye — Smoluchowski Equation with Backreaction and Comparison with Experimental Results. J. Chem. Phys., 1988, vol. 88, pp. 5620-5630. DOI: 10.1063/1.454572.

16. Potovoi V.V., Ivanov A.I. Theory of Non-Thermal Electron Transfer. Chem. Phys., 1999, vol. 247, pp. 245-259. DOI: 10.1016/S0301-0104(99)00197-4.

17. Rosspeintner A., Lang B., Vauthey E. Ultrafast Photochemistry in Liquids. Annu. Rev. Phys. Chem., 2013, vol. 64, pp. 247-271. DOI: 10.1146/annurev-physchem-040412-110146.

18. Salikhov K.M., Molin Yu.N., Sagdeev R.Z., Buchachenko A.L. Spin Polarization and Magnetic Effect in Radical Reactions. Amsterdam, Elsevier, 1984. 419 p.

19. Zusman L.D. Outer-Sphere Electron Transfer in Polar Solvents. Chem. Phys., 1980, vol. 49, no. 2, pp. 295-304. DOI: 10.1016/0301-0104(80)85267-0.

NUMERICAL SIMULATION OF MULTISTAGE BIMOLECULAR PHOTOREACTIONS IN LIQUIDS: ALGORITHMS AND SOFTWARE

Sergey V. Feskov

Doctor of Physical and Mathematical Sciences, Professor, Department of Information Systems and Computer Simulation, Volgograd State University serguei.feskov@volsu.ru https://orcid.org/0000-0001-5020-8211

Prosp. Universitetsky, 100, 400062 Volgograd, Russian Federation

Abstract. The paper considers tools for numerical simulation of multistage bimolecular photoreactions assisted by diffusive mobility of reactant molecules in viscous solutions. To describe these processes, the differential encounter theory (DET) is extended to account for the coherent dynamics of certain degrees of freedom (for example, electronic spins of reactants and intermediates). The model involves diffusion of reactants in solution and multistage / multichannel physicochemical processes proceeding both at the level of individual molecules and encounter complexes. Algorithms for numerical solution of model equations are proposed, which are related to the evaluation of evolution operators. An algorithm for computing the quantum propagator for the density matrix based on the Trotter splitting is presented. A software package for simulation of multistage photoreactions has been developed using the suggested numerical approaches. The author presents a structure of key software components and examples of the program model construction. A software testing has been carried out, showing

good correspondence between the numerical results and exact solutions of the model equations in certain particular cases. As an example, a photoreaction with participation of 9.10-dimethylanthracene and 1.3-dicyanobenzene in the acetonitrile solution has been considered, and basic procedures for configuring and simulating multistage bimolecular photoprocesses are shown. The importance of coherent description of the electronic spin evolution at the radical stage is shown.

Key words: computational chemistry, numerical methods, simulation software, Trotter splitting, bimolecular reactions, photochemistry, coherent dynamics.