Математическое моделирование и численное исследование динамики внутримолекулярных таутомерных превращений без учета процессов релаксации протона. Часть I. описание математической модели Текст научной статьи по специальности «Математика»

МАТЕМАТИКА

94

УДК 519.6+539.194

Г. В. Квитко, Э. Л. Кузин, К. С. Латышев А. И. Иванов, Д. В. Шоть

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ И ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ДИНАМИКИ

ВНУТРИМОЛЕКУЛЯРНЫХ ТАУТОМЕРНЫХ ПРЕВРАЩЕНИЙ БЕЗ УЧЕТА ПРОЦЕССОВ РЕЛАКСАЦИИ ПРОТОНА

Часть I. Описание математической модели

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

There is a quantum statistical model of prototropic tautomeic transformation caused by vertical transition between different adiabatic potentials with two minima divided by potential barrier offered. Formula for the probability function defining time of tautomer transformation without regard to vibronic relaxation is received.

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

Key words: mathematical model, hydrogen binding, proton, Schrodinger equation, adiabatic model potentials, statistical operator, probability of tautomerny transformation.

Введение

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

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

© Квитко Г. В., Кузин Э. Л., Латышев К. С., Иванов А. И., Шоть Д. В., 2013

Вестник Балтийского федерального университета им. И. Канта. 2013. Вып. 4. С. 94-101.

Постановка задачи

Каждая из таутомерных форм (назовем их I и II) относится к невырожденному состоянию одного и того же точного электронно-ядерного гамильтониана молекулы Н . Поэтому возможные наборы {Т^} и {Т«1)} стационарных состояний каждой из форм (I) и (II) являются решениями одного и того же стационарного уравнения Шредингера. Не исключено, что среди этих наборов состояний могут оказаться вырожденные пары. Однако такое вырождение реализуется только случайно и в общем случае отсутствует. В противном случае различить таутомеры спектроскопическими методами было бы невозможно.

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

В простейшей модели движение мостикового атома водорода формы I (или II) рассматривается (в рамках приближения Борна — Оппенгеймера) в одномерном эффективном адиабатическом потенциале двух донорно-акцепторных центров В(В') и А(А'). Здесь элементарный акт прототроп-ного таутомерного превращения I ^ II моделируется реакцией

в,—н а ^ в; н—а;.

Точки означают внутримолекулярную водородную связь, знак «;» — движение протона в истинном адиабатическом потенциале заменяется его движением в эффективном потенциале, создаваемом центрами В (В ) и А (А ), а также фрагментами молекулы, не указанными на моделях I и II.

Очевидно, что форма адиабатических потенциалов моделей

В;---Н....А; (I) и В;...Н----А, (II) неодинакова. Различие в длинах

связей и валентных состояний атомов форм I и II приводит к тому, что в адиабатическом приближении электронная энергия молекул I и II по-разному зависит от координаты протона мостиковой связи и, следовательно, в этом приближении формам I и II соответствуют разные приближенные эффективные гамильтонианы H< и Н« протона. Поэтому при описании таутомерного перехода приходится рассматривать не один, как это делалось, например, в работе [4], а два адиабатических потенциала V<(х) и У«(х). Здесь под координатой, вообще говоря, понимается так называемая нормальная координата протона. Как правило, данные потенциалы имеют два минимума, разделенных потенциальным барьером. Каждый из двух минимумов этих потенциалов соот-

96

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

Положим, что под влиянием внешнего воздействия форма II образуется из I с вероятностью М(Ь), зависящей от времени, и пусть Е — некоторая наперед заданная вероятность обнаружения новой формы. Тогда время таутомерного превращения т можно определить из уравнения

М (т) = Е. (1)

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

Определение вероятности обнаружения таутомерной формы

Пусть в начальный момент времени t = 0 молекула (таутомер) формы I находится в макроскопическом объеме О замкнутой равновесной системы, содержащей ансамбль из N таких молекул. В этом случае движение протона в молекуле формы I можно приближенно описать с помощью статистического оператора (матрицы плотности) р(0) [2]:

(3(0) = ехр^И / кТ] / Бриг(ехр[-И1 / кТ]) где И — эффективный гамильтониан протона формы I:

Н = - 2- т-г+х), (2)

2— а х

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

вида У1(х) = ^ Vк!) хк, где V— постоянные коэффициенты, опре-к=0

деляющие конкретную геометрию этого модельного потенциала.

Отметим, что задачи для стационарного и нестационарного уравнений Шредингера с такими потенциалами при п = 4 и п = 6 были численно решены и исследованы в работах [5; 6].

Существует другая, дираковская, форма представления статистиче-

N

ского оператора р(0) = ^ | Фт > рт < Фт |, где | Фт > — собственные

т=0

вектора эффективного гамильтониана И1, то есть решения спектральной задачи И | Ф(;т) > = Е{(т) | Ф^т) >, удовлетворяющие условиям полноты и ортонормированности, а рт — так называемые статистические ве-

са (или уровни заселенности) рт = ехр(-Е[т / кТ)/ ^ ехр(-Е^^ | кТ),

к=0

N

удовлетворяющие естественному требованию ^ рт = 1.

т=0

Положим, что в момент і = 0 исследуемая модель формы I испытывает некоторое внешнее воздействие, например фотовозбуждение пико- или фемтосекундного диапазона.

Смоделируем это воздействие некоторым оператором возбуждения РУ(х, і), полагая, что именно оно переводит эффективный гамильтониан протона Н,, соответствующий таутомерной форме I, в гамильтониан Н,, таутомерной формы II, то есть меняет адиабатический потенциал V, (х) на V,, (х). Тогда для любого момента і > 0 Нп = Н, + РУ(х, і).

В работе [2] рассматривается и обсуждается достаточно сложная ситуация, когда оператор РУ(х, і) имеет явную зависимость от времени, и отмечены проблемы математического и расчетного характера, неизбежно возникающие в такой ситуации.

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

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

В этом случае оператор возбуждения моделируется следующим об-

п

разом: РУ(х, і) = 9(і) • №(х), где №(х) = ^ ^кхк, а 9(і) — функция Хе-

к=0

Г0, если і є (-да, 0), висайта (скачка): 9(і) = <

[1, если і є [0, да).

Постоянные коэффициенты цк полинома №(х) определяются разностью соответствующих коэффициентов в эффективных потенциалах

уп(х) и у, (х): ^ = ^П) - ^П).

В любой момент времени і > 0 состояние протона в возбужденной молекуле в гейзенберговском представлении можно описать статистическим оператором р(і) = и(і) р(0) и + (і) [6], где и(і) — оператор эволюции (или пропагатор), позволяющий учесть динамику процесса перехода I ^ II. Знак «+« над пропагатором означает эрмитовское сопряжение. и (ґ) — это унитарный оператор, который в данной постановке задачи имеет следующий вид: іи(і) = ехр(-(і / Й)НПі). Здесь і — мнимая

98

единица, а h — постоянная Планка. С учетом соотношений (2) и эрми-товости оператора Нп статистический оператор p(t ) принимает вид

p(t) = exp(x, x + Ax) • [exp(-Hn / kT) / Sp(exp(-iH:l / kT))]exp((z' / h)Hnt).

Мгновенно изменение гамильтониана Hl ^ Нп не меняет исходного состояния системы [7]. Поэтому матрица плотности р(0) принята в качестве начального состояния в последующей эволюции системы.

Для дальнейшего рассмотрения существенно, что процесс перехода протона от возбужденного к новому равновесному состоянию связан с изменением вероятности его пребывания в любой из произвольно выбранных областей [x, x + Ax], связанных с движением протона.

Определим проекционный оператор Р как оператор, проектирующий любой вектор гильбертова пространства H на подпространство собственных векторов гамильтониана Нп, соответствующего таутомер-ной форме ll, имеющего размерность N, следующим образом:

N

р = ^|ф<;1 ) хф^і. (3)

п=0

Известно, что в случае дискретного спектра оператор (3) совпадает с единичным: Р = I. Поэтому левую часть последнего соотношения называют разложением единицы по отношению к собственным векторам в данном случае оператора Нп . В этом случае имеем следующее соот-

N

ношение Р I Ф,™) > = Фц) і Фіт) > і Фі") >, где выражения < Ф5") | Ф5т) >

n=0

a

представляют собой интегралы вида < Ф^) | Ф[m) > = J Фп)(x) Ф5т)(x)dx.

- a

В общем случае эти интегралы перекрывания могут быть несобственными: a ^ +со. Итак, основываясь на вышесказанном, в операторном виде вероятностную функцию M(t ) можно представить формулой

M(t) = Sp jexp(--Hut) —exp(-Hi/ kT )— exp(— Hnt )P | w F[ h ll ' Sp[exp(-h,/kT)] VKh 11 ’ J.

Ясно, что численная реализация решения (1) для такой сложной функции возможна только при переходе к матричному представлению.

Выбор базиса

При выборе базиса, в котором будет проводиться численное решение поставленной задачи, будем полностью следовать идеям и технологиям, предложенным и численно реализованным в работах [5; 6]. Спектральная задача для уравнений Шредингера с гамильтонианами Н, и Нп решается вариационным методом Ритца. Приближенное решение ФЩ задачи на собственные значения представляется в виде разложения по некоторому конечному набору полной системы функций, обра-

зующих базис в рассматриваемом пространстве. В работе [6] было использовано следующее представление:

N

Ф<">(х, ю) =£ спкик(х, ю), 7 = 1,11. (4)

к=0

где и^ (х, ю) — собственные функции некоторого гармонического осциллятора с собственной частотой (параметром) ю:

ик(х, ю) = (2к к!Х(ю)>/л)-1/2 ехр[-(х /Х(ю))2 /2]И^ (х /Х(ю)).

Здесь Нк(^) — полиномы Эрмита; ^ = х /Х(ю), а Х(ю) = (ют/Й) 1/2. Свойства базисных функций ^ (х, ю) хорошо изучены. Они ортонор-

мированы и удовлетворяют условию полноты.

Используя известный операторный метод, разработанный Дираком, в явном виде были найдены матричные элементы для каждого из гамильтонианов И1 и Нп . Для полинома степени т = 6 они имеют следующий вид [6]:

Нпп =4 +(Ч +1/2)( п+1/2)+3/4-Х4(2п2 +2п+1) + 5/8-Х6(4п3 +6п2 +8п+3), Нп ,п+1 = Нп+і,п = [X 1 + 3/2 X з(п +1) + 5/4 X 5(2 п2 + 4 п + 3)][(п +1)/2]1/2, Нп,п+2 = Ни+2,и = [(X1 -1/2)+Х4(2п+3) + 5/4X5(п2 + 3п+3)] [(п+1)(п+1)/4]1/2, Нп ,и+3 = Нп+Зп = [X 3 + 5/2 X 3(п + 2)][(п + 1)(п + 2)(п + 3)/8]1/2,

Нп ,п+4 = Нп+4,п = [X 4 + 3/2 -X 6(2п + 5)][(п + 1)(п + 2)(п + 3)(п + 4)/16]1/2, Нп,п+5 = Нп+5,п = X 5 [(п + 1)(п + 2)(п + 3)(п + 4)(п + 5)/ 32]1/2 ,

Нп,п+6 = Нп+6,п = X 6 [(п +1)(п + 2)(п + 3)(п + 4)(п + 5)(п + 5)/ 64] 1/2 .

Все остальные матричные элементы Н^ равны нулю. Видно, что матрица 11Н^ 11 имеет ленточный вид и является симметрической: Н^ = Н^к . Если положить, что X5 = X6 = 0, то получаем матричные элементы для полинома степени т = 4 . Постоянные коэффициенты Xk (к = 0,1,...,6), которые входят в матричные элементы Нкп, это обез-

размеренные коэффициенты в полиномах V (х), Уп (х), определяющие их геометрию. Они связаны со старыми коэффициентами и(и1; ^п) следующим соотношением: Xk = 4ик(тю / Й)-к/2 / Йю, где ю — частота (оптимального) гармонического осциллятора.

Вариационный метод Ритца обеспечивает получение приемлемых решений Ф7п) стационарного уравнения Шредингера за счет соответствующего оптимального выбора коэффициентов спк в разложении (4) искомых решений по базису собственных функций и^ (х, ю) гармонического осциллятора. Однако сами базисные функции зависят от парамет-

99

100

ра ю — частоты колебаний осциллятора. Понятно, что при ограниченном (обрезанном) базисе N величина этого параметра может существенно повлиять на решения Ф^п) стационарного уравнения Шредингера.

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

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

Перепишем вероятностную функцию M (t ) в матричной форме, позволяющей в конечном итоге провести ее численный расчет. Для этого перейдем к описанию в базисе собственных векторов оператора Н(, то есть | Ф^ >, которые (так же как и собственные векторы оператора Нц ) удовлетворяют условию ортогональности, нормированы на единицу и удовлетворяют условию полноты S | Ф<т) > < Ф<т) |= I, где I —

m=0

единичный оператор. Тогда после серии несложных, но громоздких преобразований получаем:

M(t) = £< Ф(іт) I Р(t)Р| Ф(,т) > = SМп, expH^f),

m n ,s

Mn,s = S < Ф(т) | Ф(і) > < Ф(і) | Ф?) > < Ф(,) | Ф(п) > < Ф(п) | Ф(т) >, (5)

m ,q

Pq = exp(-E((q) /kT)/Sexp(—E((r) /kT), Q™ = (Ej>) -Е(П))/h.

r

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

Ims(®j, ю2) =< Ф((т) | Ф^ > = | Фi(m) (x,) • Ф(і) (xю2)dx.

—да

В части II настоящей работы будут представлены, проанализированы и обсуждены результаты численных расчетов, полученных на основе предложенной здесь квантово-статистической модели.

Работа выполнена при финансовой поддержке РФФИ по проектам № 12-01-00477а, № 11-01-00098а и № 11-01-00558а.

Список литературы

1. Квитко Г. В., Кузин Э. Л., Шоть Д. В. Математическая модель внутримолекулярного таутомерного превращения и процессы релаксации протона / / Вестник Российского государственного университета им. И. Канта. 2009. Вып. 10. С. 104 — 111.

2. Квитко Г. В., Кузин Э. Л., Новиков В. И. Квантовая статистическая модель внутримолекулярного таутомерного превращения // Теоретическая и экспериментальная химия. 1975. Т. 11, № 6. С. 754—761.

3. Грибов Л. А., Дементьев В. А., Завалий М. В. и др. Компьютерное моделирование изомеризации сложных молекул с использованием суперкомпьютера типа МВС 1000 / / Журнал структурной химии. 2005. Т. 46, № 2. С. 303—310.

4. Brickmann J., Zimmermann H. Lingering time of proton in well of double-minimum potential of hydrogen bonds // The Journal of Chemical Physics. 1966. Vol. 50, № 4. P. 1608 — 1618.

5. Квитко Г. В., Кузин Э. Л., Шоть Д. В. Численное решение уравнения Шре-дингера с полиномиальными потенциалами (Ч. I) // Вестник Балтийского федерального университета им. И. Канта. 2011. Вып. 5. С. 115 — 119.

6. Квитко Г. В., Кузин Э. Л., Шоть Д. В. Численное решение уравнения Шре-дингера с полиномиальными потенциалами (Ч. II) // Вестник Балтийского федерального университета им. И. Канта. Вып. 10. С. 79 — 93.

7. Ландау Л. Д., Лифшиц Е. М. Квантовая механика. М., 1963.

Об авторах

Геннадий Васильевич Квитко — канд. физ.-мат. наук, доц., Балтийский федеральный университет им. И. Канта, Калининград.

E-mail: [email protected]

Эдуард Леонидович Кузин — канд. хим. наук, доц., Балтийский федеральный университет им. И. Канта, Калининград.

E-mail: [email protected]

Константин Сергеевич Латышев — д-р физ.-мат. наук, проф., Балтийский федеральный университет им. И. Канта, Калининград.

E-mail: [email protected]

Иванов Алексей Иванович — д-р физ.-мат. наук, проф., Балтийский федеральный университет им. И. Канта, Калининград.

E-mail: [email protected]

Дмитрий Владимирович Шоть — ассист., Балтийский федеральный университет им. И. Канта, Калининград.

E-mail: d. [email protected]

About authors

Dr Gennady Kvitko — ass. prof., I. Kant Baltic Federal University, Kaliningrad.

E-mail: gkvitko.univ@gmail. com

Dr Edward Kuzin — ass. prof., I. Kant Baltic Federal University, Kaliningrad.

E-mail: [email protected]

Prof. Konstantin Latyshev — I. Kant Baltic Federal University, Kaliningrad.

E-mail: [email protected]

Prof. Aleksey Ivanov — I. Kant Baltic Federal University, Kaliningrad.

E-mail: [email protected]

101

Dmitry Shott — ass., I. Kant Baltic Federal University, Kaliningrad. E-mail: [email protected]