Первые результаты
В связи с новизной поставленной задачи физического моделирования химической кинетики мы отступаем от принятого в научной литературе порядка изложения, когда сначала показывают исходные посылки, методы и используемые инструменты, а уж затем предъявляют первые полученные результаты. Мы попытаемся заинтересовать читателя на первых же шагах изложения этой части нашей работы. Поэтому мы начинаем изложение предъявлением первых практических результатов. Это облегчит труд читателю, поскольку при решении задачи нами было задействовано несколько научных дисциплин – химия, физика, моделирование и вычислительная математика. Начинать изложение с подробностей примененных и изобретенных математических приемов было бы лучшим способом потерять любого читателя.
Мы предъявляем для ознакомления результаты трёх компьютерных экспериментов, поставленных с целью прогноза химической кинетики.
Пример 1. Кинетика реакций изотопного обмена в плотном газе водорода
Пример 2. Кинетика накопления комплексов ван дер Ваальса в плотном газе пропилена
Пример 3. Моделирование экстракции с помощью клеточного автомата
Все подробности выбора объекта каждого эксперимента и средств моделирования находятся в нижеследующих разделах этой части работы.
Пример 1. Кинетика реакций изотопного обмена в плотном газе водорода
Для компьютерного эксперимента была выбрана модель сверхплотного газа водорода. Начальное состояние системы организовано в компьютере с помощью клеточного автомата. Задан кубический объем, состоящий из 8 тысяч кубических ячеек, каждая с ребром в 2 Å. Клеточный автомат в одну ячейку может поместить лишь одну молекулу водорода. Программа случайным образом размещает по ячейкам 400 молекул H2 и 400 молекул HD. Ожидается появление молекул D2 в результате изотопных обменов в ходе столкновений молекул. Популяция молекул D2 вначале равна 0. Температура T=300 К. Программа распределяет молекулы по энергиям в соответствии с законом Больцмана и с заданной температурой. Эти условия в среднем соответствуют такой ситуации, когда клетка, занятая молекулой, окружена объемным слоем пустых клеток. Это соответствует такой плотности газа, которая остается неизменной при переходе охлажденного газа в жидкость без дальнейшего сжатия газа. В жидкости любая пара соседних молекул будет совершать колебательные движения в потенциальной яме ван дер Ваальса,.
Далее модель поступает на вход нашей оригинальной программы молекулярной динамики, которая не навязывает парам молекул колебательный характер движения, а просто следит за свободными поступательными движениями и соударениями молекул. При встрече двух молекул HD программа с вероятностью 0.25 заменяет эти две молекулы молекулами H2 и D2 независимо от скоростей молекул. При встрече молекул H2 и D2 программа заменяет их двумя молекулами HD. Скорости и направления движения новых молекул определяются законами сохранения энергии и импульса.
Результаты работы программы показаны на рисунке 1.
Рис. 1. Кинетика реакций изотопного обмена в сверхплотном газе водорода. По оси абсцисс отложено время процесса в секундах. По оси ординат отложены численности популяций в системе. 1 – популяция HD, 2 – популяция H2, 3 – популяция D2.
Выводы по данным результатам
Показана возможность прогнозирования кинетики реакций изотопного обмена на основании чисто физических представлений о динамике хаотического движения молекул и о квантовых законах перестройки молекулярных структур в процессе соударения молекул. При этом не использовались никакие данные о константах скоростей прямых и обратных реакций. Также не составлялись и не решались никакие кинетические уравнения. Вместо таких химических предположений и экспериментальных данных использовалась только техника составления и решения уравнений свободного движения для каждой молекулы ансамбля совместно с законами обмена скоростями в парных столкновениях. Вероятности элементарных актов изотопного обмена определялись на основе квантовых представлений о прохождении частицы под потенциальным барьером.
Следует обратить внимание на прихотливый характер полученных кривых. Начало процесса имеет явно экспоненциальный характер, что вполне привычно для экспериментатора и легко моделируется решением систем кинетических уравнений. Но затем наступает более длительный период линейного роста продуктов изотопного обмена. Далее следует вполне ожидаемое динамическое равновесие.
Заметны колебания популяции молекул D2. Колебания начинаются после достижения динамического равновесия в системе. Можно предположить, что это не флуктуации, хотя бы потому, что амплитуда колебаний больше средней статистической флуктуации, равной 0.9. Предположение связано и с характером всей кривой, не согласующимся с картиной флуктуаций.
Также хорошо видно, что колебания популяций молекул HD и H2 строго согласованы с колебаниями популяции молекул D2. Этого, впрочем, следовало ожидать, поскольку в замкнутой системе полное число частиц должно оставаться неизменным.
Напрашивается содержательная аналогия с поведением популяций в экологической системе Хищники-Жертвы-Ресурсы. Такая система описывается дифференциальными уравнениями Лоттки-Вольтерра. В простейшем случае, когда в системе имеется один тип хищников и один тип жертв, решения уравнений носят либо стационарный, либо регулярный колебательный характер. У нас же колебания получились не совсем регулярными, что вполне согласуется со случайностями каждого парного столкновения частиц в нашей системе и с вероятностным характером акта изотопного обмена в столкновении.
Ясно, что невозможно составить и решить такую систему уравнений кинетики реакций второго порядка, которая предсказывала бы нерегулярный колебательный характер популяций после достижения динамического равновесия в системе. Дело в том, что дифференциальные уравнения составляются для макроскопических приращений популяций, а наша программа молекулярной динамики, дополненная условиями прохождения обменных реакций, оказывается ближе к описанию реальных событий в химической системе.
В следующей части данной работы, посвященной уже биологической эволюции, мы вернёмся к дифференциальным уравнениям Лоттки-Вольтерра. Мы составим систему уравнений для экологической системы, где присутствуют несколько типов хищников и жертв, дополним систему уравнений возможностью описывать случайные события типа эпидемий и других катастроф, и только тогда сумеем предсказать адекватную картину развития системы. Эта картина описывает нерегулярные колебания популяций в системе, что соответствует наблюдениям за реальными экологическими системами. Тогда станет более ясно, что в примере кинетики изотопного обмена картина нерегулярных колебаний не может быть списана на неточности компьютерной имитации этого процесса.
Пример 2. Кинетика накопления комплексов ван дер Ваальса в плотном газе пропилена
Мы полагаем, что химическая кинетика определяется на микроскопическом уровне следующими факторами.
- Вступающие в реакцию молекулы должны найти друг друга в пространстве, перебирая различные сочетания типов молекул в случайных парных столкновениях.
- В ходе столкновения две подходящие молекулы должны расположиться друг относительно друга удачным для реакции образом. Благоприятны лишь такие расположения молекул в предварительно возникающем комплексе, когда реакция может осуществиться путем минимальных перемещений минимального количества атомов. Тогда возникает реакционный центр во временном комплексе.
- Колебательные состояния обеих молекул должны обеспечить высокую вероятность структурной перестройки ради получения желаемого продукта.
- Продукты реакции, если это не единственная молекула, должны иметь возможность разойтись в пространстве, чтобы тут же не произошла обратная перестройка структур в исходные, иначе факт элементарного акта реакции не будет заметен в макроскопических измерениях.
Была поставлена конкретная задача отделить первые два фактора и выявить их роль в общей кинетической картине.
Предполагалось смоделировать газ, состоящий из одинаковых молекул, которые при столкновениях способны образовывать реакционный комплекс. Но решено для физической ясности при благоприятных столкновениях, когда две молекулы находятся в подходящих взаимных положениях, образовывать не готовый реакционный комплекс, а некий пред-комплекс, когда молекулы просто держатся друг друга в потенциальной яме ван дер Ваальса.
Выбрана модель комплекса Пропилен-Пропилен. Такие промежуточные комплексы ведут к цепочке последовательных реакций, на выходе которой получаются олигомеры, заготовки полипропилена. Механизмы всех промежуточных реакций описаны в монографии [9].
В данном упрощенном компьютерном эксперименте выполнена имитация самого первого этапа процесса полимеризации. Упрощение картины состоит в том, что комплексы ван дер Ваальса до конца эксперимента не распадаются и не переходят к акту химической реакции. Такое упрощение оправдано тем, что комплексу ван дер Ваальса требуется некоторое время, в течение которого две частицы совершают сложные колебания друг относительно друга, приспосабливаясь наилучшим способом к конфигурации реакционного промежуточного комплекса. Характер таких колебаний будет нами проанализирован в соответствующем разделе этой части работы.
В кубическом объеме с идеально отражающими стенками случайно разместили 720 молекул пропилена, придали им случайные скорости, распределённые по закону Больцмана при температуре 400 К, и поручили их программе молекулярной динамики. В результате была получена картина накопления во времени комплексов Пропилен-Пропилен, показанная на рисунке 2.
Рис. 2. Убывание популяции мономеров пропилена 1 и накопление временных комплексов типа ван дер Ваальса из двух молекул пропилена 2. По оси абсцисс отложено время, выраженное в количестве вычислительных шагов программы. Шагом считается такой промежуток времени, на котором состоялось одно продуктивное столкновение двух молекул пропилена, столкновение, приведшее к образованию комплекса. Программа анализирует условия столкновения и решает, состоится комплекс или нет.
Картина вполне ожидаемая, известная. Это экспоненциальный рост продукта и такое же убывание исходного вещества. В данном имитационном эксперименте мы убедились, что для формирования такой кинетической зависимости было достаточно учета двух первых физических факторов. Учет третьего фактора привел бы только к изменению параметров экспонент.
Выводы по данным частным результатам
В простой физической ситуации, когда при столкновении молекул образуются комплексы ван дер Ваасльса, но не случаются химические реакции, принятая методика моделирования дала качественно верные результаты. При этом мы не привлекали в имитационный эксперимент никакой химической информации, кроме сведений о возможности образования подобного комплекса в процессе олигомеризации пропилена. Мы также не составляли никаких кинетических уравнений.
Краткие сведения о методике проведения эксперимента, общей для всех описываемых здесь имитаций
Для вычисления характеристики конфигурации благоприятного столкновения построена модель комплекса, показанная на рисунке 3. При регистрации факта парного столкновения молекул пропилена программа молекулярной динамики проверяет соответствие конфигураций возникающего случайно и заданного стандартного комплекса.
Рис. 3. Предполагается, в соответствии с [9], что при сближении показанных молекул к общему центру масс винильные группировки сблизятся до нужного расстояния, и возникнет четырехчленный цикл, который дает начало дальнейшему процессу олигомеризации пропилена.
Возможное использование полученных теоретических кинетических кривых
Можно вписать в эти кривые аналитические зависимости, пользуясь методом наименьших квадратов. Продифференцировать аналитические кривые и получить скорость реакции Пропилен+ Пропилен –> комплекс. Тогда можно проверить, укладывается ли полученная константа «реакции» в представления химиков о кинетических уравнениях и их природе.
Пример 3. Моделирование экстракции с помощью клеточного автомата
Это короткое исследование потребовало самой простой вычислительной техники – программы клеточного автомата. В данном исследовании проверялась возможность моделирования процесса экстракции чисто физическими средствами. Поэтому здесь использована весьма абстрактная модель. Тем не менее результаты этой работы могут представлять определенный методический интерес.
Задаем в программе клеточного автомата фазовое пространство из nCells ячеек. Задаем число частиц первого типа nA. Задаем число частиц второго типа nB. Организуем бесконечный цикл. В цикле распределяем случайно все частицы по ячейкам. Накладываем ограничение: в ячейку может попасть не больше 2 частиц.
На каждом шаге цикла проверяем факты попадания в одну и ту же фазовую ячейку двух частиц разного сорта. Такие две частицы имеют химическую возможность вступить в реакцию, образовать новое соединение такого типа, что сама химическая система удалит новое соединение вовне. Это и будет моделировать экстракцию.
Построим график зависимости популяций частиц от числа шагов. Это и будет конечным результатом работы программы.
Выясняется, что при малом числе частиц кинетические кривые получаются ступенчатыми. При увеличении числа частиц ступеньки исчезают, но проявляются флуктуации. При достаточно большом числе частиц кривые получаются совершенно гладкими. И при этом расчет занимает мало машинного времени.
Приведем результаты для системы со сравнительно большим числом частиц. nCells = 100000; nA = 5000; nB = 1000.
Рис. 4. Популяции молекул экстрагируемого вещества 1 и экстрагента 2. В имитации экстракции задано очень большое число ячеек по сравнению с числом отслеживаемых частиц. Это соответствует возможному присутствию в жидкости каких-то нейтральных частиц, не имеющих отношения к экстракции.
Получается почти гладкая картина, показанная на рисунке 4. И ясно, что огибающая графика уменьшения популяции каждого компонента системы это падающая экспонента. Но экспонента в качестве аппроксимации не очень хороша, поскольку она отвечает теоретическому предположению, что объекты nA и nB имеют дифференциальную природу. То есть, могут изменяться на бесконечно малую величину. Теоретическая экспонента nB(t) тянется в бесконечность и никогда не падает до нуля.
Клеточный автомат преодолевает это сильное предположение о природе объектов. И находим в имитационном эксперименте, что при очень малом числе частиц несколько десятков шагов теплового движения приводят к полной экстракции частиц В. А когда частиц стало больше, то понадобилось 230 шагов для полного удаления частиц В из системы.
Мы видим, что клеточные автоматы не только успешно решают конкретные прикладные задачи. Они еще могут плавно заполнять непреодолимый для математики пробел между дискретной и непрерывной моделью подобного процесса, где сама Природа считает на штуки. С помощью клеточного автомата мы можем посчитать на штуки элементарные акты реакций, а потом так увеличить статистику, чтобы кривые стали гладкими. Для таких гладких кривых можно подобрать эмпирические формулы, а затем получить из них дифференциальные уравнения. И, в отличие от нормы, когда коэффициенты уравнений приходится брать из опыта, клеточный автомат позволяет прямо связать коэффициенты дифференциальных уравнений с характеристиками элементарных процессов.
