CC BY
36
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
Вестн. Ом. ун-та. 2012. № 2. С. 200-204.
УДК 519.876.5, 535.379, 544.431.22, 546.21, 519.85
А.А. Юнусов, М.Ю. Овчиннков, С.Л. Хурсан, И.М. Губайдуллин
МОДЕЛИРОВАНИЕ ГЕНЕРАЦИИ
СИНГЛЕТНОГО КИСЛОРОДА
ПРИ РАЗЛОЖЕНИИ ДИМЕТИЛДИОКСИРАНА
Синглетный кислород играет существенную роль в органическом синтезе и биологических процессах. Распад диоксиранов, катализированных рядом анионов, сопровождается высокоэффективной генерацией синглетного кислорода. Сложностью верификации возможных моделей является то, что единственной величиной, получаемой из эксперимента, является интенсивность хемилюминесценции. Для поиска параметров модели было предложено использовать генетический алгоритм, распараллеленный с использованием ОрепМР. В ходе исследования были получены кинетические параметры, адекватно описывающие процесс и находящиеся в соответствии с теоретическими прогнозами.
Ключевые слова: глобальная минимизация, химическая кинетика, обратная кинетическая задача, кинетические параметры, генетический алгоритм, синглетный кислород, ОрепМР.
Введение
Синглетный кислород Ю2 на протяжении многих лет привлекает внимание исследователей ввиду его существенной роли в органическом синтезе и биологических процессах [1]. Немаловажное значение имеет участие 102 и в хемилюминесцентных реакциях [2]. Ранее была обнаружена реакция, сопровождающаяся высокоэффективной генерацией 102, а именно распад диоксиранов [3], катализированный рядом анионом (С1-, Вг-, I-, ОН- и 02-) [4]. Выход синглетного кислорода, измеренный в этой пероксидной системе методом ИК-ХЛ, достаточно высок.
Данные о закономерностях протекания этой реакции будут полезными для более глубокого понимания процессов генерации возбужденных состояний в реакциях пероксидов, а также факторов влияющих на стабильность диоксиранов в растворе. Высокий выход 102 позволяет рассматривать систему диоксиран-нуклеофильный ион как перспективную в химических лазерах и органическом синтезе.
1. Реакция индуцированного анионом хлора распада диметил-диоксирана (ДМД)
На основе результатов исследования каталитического процесса генерации 102 [5], а также «темнового» процесса окисления хлорид-иона, приводящего к выводу катализатора из системы, была предложена схема взаимодействия ДМД с хлорид-ионом, изображенная на рис. 1.
о—С1
(I)
Н3С о Н3С,
о—с
(И)
HjC о Н3С
X
о'
\ Н3С О- Н,С О--------------------------------------С1
Н3С о Н3С о Рис. 1. Исследуемая схема
© А.А. Юнусов, М.Ю. Овчиннков, С.Л. Хурсан, И.М. Губайдуллин, 2012
Реакции (I) и (II) относятся к каталитическому направлению процесса, а реакция переноса электрона (III), протекающая параллельно с (II), инициирует цепочку превращений, ведущих к необратимому расходованию катализатора С1-. Преимущественно реализуется экзотермическая реакция (II) с образованием Ю2[5], однако параллельно происходит эндотермический процесс переноса электрона (III), который дает небольшой вклад в суммарную реакцию. Один из продуктов этого взаимодействия - анион-радикал ДМД - расходуется, по-видимому, в реакции с растворителем. Другой продукт окислительно-восстановительной реакции (III) - 2-хлороксиизопропильный радикал -приводит к переводу катализатора в неактивную форму. Данный процесс не был исследован.
2. Прямая и обратная задача
2.1. Прямая задача
Прямая задача химической кинетики состоит в нахождении состава реагирующей смеси по заданной кинетической модели. Процессы в химической кинетики в общем случае описываются следующей системой обыкновенных нелинейных дифференциальных уравнений [6]:
&х ____ N
= р, 1 = 1, М; р. = Ё *№■ (1)
&
]=1
(2)
где XI - концентрации участвующих веществ (моль/л), & - скорость изменения концентрации (моль/(сек*литр)), М - число веществ, N - число стадий, щ - матрица стехиометрических коэффициентов, Ш] - скорость ]-й стадии, к/ - константа скорости стадии. Таким образом, прямая задача состоит в решении СОНДУ.
В нашем случае СОНДУ имеет следующий вид:
&х
—1 = -к1х1 х2 - (к2 + к3)х1 х3,
&
ёх,
&
2 = -к1 х1 х2 + к2 х1 х3,
&х
—3 = к1 х1 х2 - (к2 + к3) х1 х3 &
&х1
&
= к2 х1 х3,
где XI($ - концентрация диметилдиоксирана (моль/л), Х2(^ - концентрация катализатора С1 - (моль/л), хз(Ц - концентрация 2-хлороксиизопропилат-иона (моль/л), Х4@) -концентрация синглетного кислорода (моль/л), к.1($ - константа скорости первой стадии (л/(моль*с)), к2@) - константа второй стадии (л/(моль*с)), кз(р - константа «темно-вого» канала реакции (л/(моль*с)).
2.2. Обратная задача
Обратная задача химической кинетики
- это задача восстановления по экспериментальным данным вида кинетической модели и её параметров. Для этого минимизируется критерий отклонения экспериментальных и расчетных данных, который в традиционном случае выглядит следующим образом:
К М
гК 'у'Е
р\
(3)
р=хц ХК
1=1 j=l
где Хцк - расчётные значения концентраций наблюдаемых веществ (моль/л); Хф - экспериментальные данные по наблюдаемым веществам (моль/л); К - количество точек эксперимента; М - количество наблюдаемых веществ, участвующих в реакции.
Особенность данной работы заключается в том, что нет информации о концентрациях, минимизировался критерий отклонения интенсивностей хемилюминесценции:
N
р= 11 I
к=1
- г\
с1 1с1Ь
(4)
где 1оР - экспериментальные данные интенсивности хемилюминесценции, 1о1° - расчётные значения, N - число точек эксперимента.
Интенсивность ХЛ зависит от квантового выхода фо1, константы скорости каталитического процесса к и концентраций анионов хлора и диметилдиоксирана:
!с1 = <Рс№ = (Рск [С1—] [ ВМВ ] (5)
Величины фо1 и к являются искомыми параметрами модели, которые не меняются во время процесса, но неизвестны до опыта, значения концентраций известны перед экспериментом.
Для того чтобы упростить поиск параметров модели, мы избавились от фо1, заменив критерий (5) на следующий:
1' = Ь*. = п = к[С1—] [ ВМВ] (6)
Рс1
Дополнительной особенностью задачи являлось наличие помех, вызванных недостатками оборудования и тем, что смешивание исходных веществ происходит не за бесконечно малое время. Из-за этих помех невозможно одновременно с одинаковой точностью определить все константы одновременно.
3. Методы решения
3.1. Метод решения обратной задачи Для решения СОНДУ первоначально использовался классический метод Рунге-Кутты 4-го порядка с фиксированным шагом. Но при дальнейших исследованиях обнаружилось, что задача слишком жёсткая, ошибки накапливались слишком быстро, и был использован явный вложенный метод Дормана-Принса 9(8).
Число процессоров
Рис. 2. Ускорение и эффективность параллельной программы
3.2. Метод решения обратной задачи
Решение обратной задачи химической кинетики, являющейся, по сути, задачей оптимизации, целевой функцией которой является (4), где вместо Icl взято (6).
Существует большое количество методов оптимизации: покоординатного спуска, Монте-Карло, имитации отжига, эволюционные алгоритмы, генетические алгоритмы, дифференциальной эволюции, муравьиный алгоритм, метод роя частиц и т. д. До недавнего времени в основном использовался метод покоординатного спуска, с развитием параллельных технологий стали активно применяться методы Монте-Карло, генетические алгоритмы [7], методы роя частиц, а также индексный метод глобальной оптимизации [8]. В данной работе использовался генетический алгоритм, изображенный на рис. 2.
4. Программная реализация и распараллеливание
Явный вложенный метод Дормана-Принса 9(8) был взят из библиотеки GNU Scientific Library, распространяющейся по лицензии GPL.
Этапы генетического алгоритма были реализованы следующим образом:
1. Начальная популяция представляла собой N случайным образом выбранных наборов констант изучаемой системы.
2. Скрещивание производилось следующим образом: — раз выполнялась следующая процедура.
3. Генерируется 2 случайных целых числа I и т на промежутке 1, N. Они генерируются по закону равномерного распределения.
4. Генерируется новый набор констант: к^ = °к1 + (1 - а)кт, где а е (0;1) - случайное число, также генерирующееся по закону равномерного распределения.
5. Новый набор констант добавляется к популяции.
• Во время мутации ко всем членам популяции применялась следующая операция:
к1п = К& +r, г е (-а, а X а > °
где г - случайное число, генерируемое по закону нормального распределения, а - заранее известное число, параметр генетического алгоритма.
• Вычисление целевой функции представляет собой решения для каждого экземпляра популяции прямой задачи - решения СОНДУ и вычисления выражения (4).
• При селекции популяция сортировалась по значениям целевой функции, и из неё удалялась половина наборов с наихудшими значениями.
Профилирование первоначального последовательного варианта показало, что наиболее ресурсоёмким этапом является этап вычисления целевой функции, более 95 % времени тратилось на данный этап. Соответственно в первую очередь решалась задача ускорения данного этапа.
Первоначально функция вычисления целевой функции выглядела следующим образом:
void computeFitnessFunction(vector<vector<double> > population,
InitialData initial_data, vector<double> &fitnessFunction, const vector<double> &x0)
{
unsigned int i;
fitnessFunction.resize(population.size()) ; for (i = 0; i < population.size(); i++) {
double jdiff;
jdiff = jDifference(population[i], initial_data, x0); fitnessFunction[i]= jdiff;
}
}
Ввиду очевидного параллелизма по распараллеливание с использованием техноданным (data parallelism), было предложено логии OpenMP:
void computeFitnessFunction(vector<vector<double>> population,
InitialData initial data, vector<double> &fitnessFunction, const vector<double> &x0)
{
unsigned int i;
fitnessFunction.resize(population.size());
#pragma omp parallel for shared(fitnessFunction, population, \ initial data, x0) private (i)
for (i = 0; i < population.size(); i++) {
double jdiff;
jdiff = jDifference(population[i], initial data, x0) fitnessFunction[i]= jdiff;
}
}
С использованием функции OpenMP omp_set_num_threads, определяемой как void omp_set_num_threads(int num_threads), были получены следующие графики ускорения и эффективности распараллеливания программы (рис. 2). Все расчёты проводились на персональном компьютере с двумя процессорами AMD Phenom II X4 940, че-тырмя гигабайтами памяти.
5. Результаты и заключение
С примением вышеописанного программного комплекса были найдены кинетические константы, адекватно описывающие экспериментальные данные для двух температур - 9,5 °С (283 K) и 14,5 °С (288 K), в целом совпадающие с результатами, полученными методом квазистационарных приближений.
Согласно теоретическому анализу [5], при температуре 288 K константа
[Cl ] •104, моль/л к, kJ к2 -102
1,0 125,58 1,4
1,0 125,7 1,12
1,2 119,34 1,44
2,0 108,98 1,65
констант
к = 1,3 ± 0,3 10 л • моль • с , отношение
• 10-2, начальная диметилдиоксирана равна растворителем являлся аце-
k3 / k2 = 1,1 ± 0,2 •
концентрация 6,0 -10-3 моль/л,
тон. В ходе данного исследования были получены следующие значения данных параметров.
При малой концентрации С1- константа каталитического канала определяется плохо, ввиду того, что он слабо влияет на процесс, и поэтому константы к2 и кз определялись при обработке таких экспериментов. Для экспериментов, в которых концентрация С1- была наибольшей, к1 определялась достаточно конкретно, к2 и кз могли иметь различные значения. На рис. 3 показано сравнение экспериментальных и расчётных значений интенсивностей.
Главным результатом данной работы следует считать не то, что в целом были подтверждены значения, полученные методом квазистационарных приближений, а следующее - применение параллельных вычислений позволяет:
1) разделить константы к1, к2, кз, что не было возможным при старом подходе, при котором была получена зависимость отношения кз/к2;
2) с более высокой точностью определять константы, чем при методе квазиста-ционарных приближений, который упрощает систему;
0 _I________I_________I________I_________I________I_________I________I_________I________
0 2 4 6 8 10 12 14 16 18 20
t, сек.
10-4 л
Рис. 3. T = 283K, [Cl—] = 4,8--------------
(моль • с)
3) масштабировать данный подход на другие системы, поведение которых описывается более сложными системами дифференциальных уравнений, которые не имеют решения в аналитическом виде.
В ходе дальнейших исследований предполагается также обработать массив данных для 316 K и получить активационные параметры реакции.
ЛИТЕРАТУРА
[1] Wasserman H. H., Murray R. W. Organic Chemistry: Singlet Oxygen. N.Y.: Acad. Press, 1979. V. 40.
[2] Kazakov D. V, Kazakov V. P., Maistrenko G. Ya., Mal’zev D. V, Schmidt R. // J. Phys. Chem. A. 2007. V. 111. P. 4267.
[3] Adam W, Saha-Muller C.R., Zhao C.-G. // Acc. Chem. Res. 1989. V. 22. № 6. P. 205.
[4] Adam W, Kiefer W, Schlucker S., Saha-Muller C., Kazakov D. V, Kazakov V. P., Laty-
pova R. R. Bioluminescence and Chemilumines-cence: Progress and Current Applications / eds. P. E. Stanley, L. J. Kricka. Singapore : World Scientific Publishing, 2002. P. 129
[5] Ovchinnikov M. Yu., Khursan S. L., Kazakov D. V., Adam W. // J. Photochem. Photobiol. A. Chem. 2010. V. 210. P. 100.
[6] Слинько М. Г. Основы и принципы математического моделирования каталитических процессов. Новосибирск : Ин-т катализа им. Боре-скова СО РАН, 2004. 488 с.
[7] Линд Ю. Б., Губайдуллин И. М., Мулюков Р. А. Методология параллельных вычислений для решения задач химической кинетики и буровой технологии // Системы управления и информационные технологии. 2009. № 2/36. С. 44-49.
[8] Губайдуллин И. М., Тихонова М. В., Рябова В. В. Применение индексного метода глобальной оптимизации при решении обратных задач химической кинетики // Вычислительные методы и программирование. 2011. Т. 12. № 1. С. 127-135.
