Журнал вычислительной математики и математической физики, 2019, T. 59, № 1, стр. 135-142

1. ВВЕДЕНИЕ

Эта статья посвящена теоретическому и численному исследованию математической модели сравнения двух молекул в задаче, возникающей в структурной химии. Эта модель сводится к сравнению двух объектов из $N$ упорядоченных точек с жесткой геометрией, ведущих себя как твердые тела в ${{R}^{3}}$. Принцип сравнения в этой модели основан на оптимизации совмещения этих двух объектов путем сдвигов и поворотов. Оптимальное совмещение этих объектов предполагает минимизацию по сдвигам и углам поворота некоторой функции сравнения геометрии объектов, представляющую собой сумму расстояний между точками двух объектов с одинаковыми индексами. Для минимизации [1], [2] функции сравнения используется метод нулевого порядка Розенброка [3], [4]. Полученные результаты использованы для исследования геометрии реальных молекул.

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

2. ОПИСАНИЕ МАТЕМАТИЧЕСКОЙ МОДЕЛИ

В двух следующих разделах приведем некоторые определения и теоремы.

Определение 1. Геометрической структурой, или короче – структурой, будем называть жесткую геометрическую конструкцию из $N$ упорядоченных точек в ${{R}^{3}}$ с координатами $({{x}_{i}},{{y}_{i}},{{z}_{i}})$, $i = 1,\; \ldots ,\;N$, движение которой в ${{R}^{3}}$ осуществляется как движение твердого тела.

Будем считать, что каждой $i$-й точке структуры сопоставлен весовой коэффициент ${{w}_{i}} \geqslant 0$, причем $\sum\nolimits_{i = 1}^N \,{{w}_{i}} = W > 0$. Пусть ${{i}_{k}}$, $k = 1,\; \ldots ,\;K$ – все индексы ненулевых весовых коэффициентов ${{w}_{i}}$. Тогда

Пусть даны две геометрических структуры, каждая из которых состоит из $N$ точек с координатами соответственно $({{x}_{{1,i}}},{{y}_{{1,i}}},{{z}_{{1,i}}})$ и $({{x}_{{2,i}}},{{y}_{{2,i}}},{{z}_{{2,i}}})$, $i = 1,\; \ldots ,\;N$. При этом при каждом фиксированном значении индекса $i$ точка с координатами $({{x}_{{1,i}}},{{y}_{{1,i}}},{{z}_{{1,i}}})$ первой структуры соответствует точке $({{x}_{{2,i}}},{{y}_{{2,i}}},{{z}_{{2,i}}})$ второй структуры.

Пусть геометрия двух структур одинакова. Возьмем две любые точки, задаваемые векторами ${{{\mathbf{r}}}_{{1,0}}}$ и ${{{\mathbf{r}}}_{{2,0}}}$, которые одинаково расположены относительно соответствующей геометрической структуры (в частности, в качестве ${{{\mathbf{r}}}_{{1,0}}}$ и ${{{\mathbf{r}}}_{{2,0}}}$ можно взять точки соответствующих структур с одинаковыми номерами). Тогда путем вращений можно совместить одинаковые точки до полного совпадения двух структур. Если рассматриваемые структуры имеют различную геометрию, то естественно поставить задачу оптимального “совмещения” структур путем всевозможных сдвигов и вращений одной структуры относительно другой как твердых тел.

Критерий оптимальности совмещения двух геометрических структур будем связывать с минимумом по сдвигам и углам вращения (углам Эйлера [6]) функции сравнения вида

где векторы ${{{\mathbf{r}}}_{{1,i}}}$ и ${{{\mathbf{r}}}_{{2,i}}}$ определяют положение точек в первой и второй структуре, вектора ${{{\mathbf{r}}}_{{1,0}}}$ и ${{{\mathbf{r}}}_{{2,0}}}$ определяют смещения первой и второй структуры в соответствующие точки, $Q = Q(\varphi ,\theta ,\psi )$ – матрица вращения на углы Эйлера ($\psi $ – угол прецессии, $\theta $ – угол нутации, $\varphi $ – угол собственного вращения):

Таким образом, функция $U$ вида (2) представляет собой сумму расстояний между “одноименными” точками двух геометрических структур с весами ${{w}_{i}}$ после совмещения точек, определяемых векторами ${{{\mathbf{r}}}_{{1,0}}}$ и ${{{\mathbf{r}}}_{{2,0}}}$ и поворота второй структуры относительно первой. Исследуем некоторые проблемы, возникающие при минимизации функции сравнения двух структур вида (2).

3. МИНИМИЗАЦИЯ ФУНКЦИИ СРАВНЕНИЯ ГЕОМЕТРИЧЕСКИХ СТРУКТУР

Покажем, что, совместив некоторые характерные точки двух геометрических структур, можно свести задачу минимизации функции сравнения по полному набору переменных к поиску минимума функции $U$ по углам поворота $\varphi $, $\theta $, $\psi $.

Здесь мы даем более полные и строгие формулировку и доказательство идеи, высказанной в работе [7].

Теорема 1. Минимум функции $U$ реализуется в точке, соответствующей “центру масс” двух геометрических структур и определяемой векторами

Доказательство. В функции (2) зафиксируем углы $\varphi $, $\theta $, $\psi $, а значит, и матрицу $Q$. Таким образом, будем рассматривать функцию (2) относительно переменных ${{{\mathbf{r}}}_{{1,0}}}$ и ${{{\mathbf{r}}}_{{2,0}}}$. Относительно этих переменных функция является выпуклой и квадратичной, поэтому ее минимум по ${{{\mathbf{r}}}_{{1,0}}}$, ${{{\mathbf{r}}}_{{2,0}}}$ реализуется в тех точках, в которых производная по ним равна нулю. Таким образом, получаем уравнения

Осуществляя дифференцирование, имеем В силу невырожденности матрицы ${{Q}^{{\text{т }}}}$ уравнение (7) эквивалентно (6). Решая (6) относительно ${{{\mathbf{r}}}_{{1,0}}}$, получаем Для того, чтобы скобка в (8) была равна нулю, выберем ${{{\mathbf{r}}}_{{2,0}}}$ из условия Тогда из (8) вытекает

Как отмечалось выше, уравнения (6) и (7) эквивалентны, поэтому для полученных точек (9) и (10) справедливо (5) и, значит, построенные ${{{\mathbf{r}}}_{{1,0}}}$ и ${{{\mathbf{r}}}_{{2,0}}}$ дают минимум функции (2) при фиксированных углах $\varphi $, $\theta $, $\psi $. Теорема доказана.

Таким образом, минимум функции (2) соответствует сдвигу “центров масс” геометрических структур, определяемыми формулами (10) и (9) в начало координат. Функцию (2) можно теперь рассматривать как функцию углов поворота

где ${{{\mathbf{r}}}_{{1,i}}}$, ${{{\mathbf{r}}}_{{2,i}}}$ – заданные координаты геометрических структур, а “центры масс” ${{{\mathbf{r}}}_{{1,0}}}$ и ${{{\mathbf{r}}}_{{2,0}}}$ вычисляются по формулам (10) и (9). Минимум функции (11) определяется по углам поворота $\varphi $, $\theta $, $\psi $.

Пусть минимум функции $U = U(\varphi ,\theta ,\psi )$ вида (11) достигается в точке $(\mathop {{{\varphi }_{0}},{{\theta }_{0}},\psi }\nolimits_0 )$. Мерой близости двух геометрических структур назовем величину

Эту величину можно рассматривать как количественную характеристику близости геометрических структур, так как она представляет собой усредненное расстояние между точками с одинаковыми индексами в двух структурах после “совмещения”.

Две структуры будем называть приблизительно равными, если

где ${{s}_{0}}$ – заданная величина (в приложениях определяется конкретной практической ситуацией). Неравенство $s \leqslant {{s}_{0}}$ будем называть критерием близости структур. При проведении численных расчетов величина $s$ может содержать вычислительную погрешность.

Отметим, что если точка минимума не единственна и существует точка $(\mathop {{{\varphi }_{1}},{{\theta }_{1}},\psi }\nolimits_1 )$ такая, что $U(\varphi \mathop {_{1},{{\theta }_{1}},\psi }\nolimits_1 ) = U(\mathop {{{\varphi }_{0}},{{\theta }_{0}},\psi }\nolimits_0 )$, то величина $s$ при этом не меняется и на критерий близости структур неединственность точки минимума влияния не оказывает.

Тест на сравнение геометрии двух геометрических структур теперь можно разбить на три этапа: 1) смещение центра масс каждой геометрической структуры в центр координат; 2) минимизация функции (11) по углам; 3) вычисление величины $s$ вида (13) и вывод о близости структур.

Две геометрические структуры назовем равными, если их можно совместить точками с одинаковыми индексами путем различных перемещений этих структур, как твердых тел. Как известно [8], такое перемещение можно осуществить линейным сдвигом и ортогональным поворотом.

Приведенный здесь алгоритм сравнения двух геометрических структур 1)–3) можно считать корректным, если две равные структуры с несовпадающими координатами ${{{\mathbf{r}}}_{{1,i}}}$ и ${{{\mathbf{r}}}_{{2,i}}}$ совмещаются при реализации минимума функции (11), то есть мера близости двух равных структур $s = 0$.

Теорема 2. Неотрицательная функция (11) достигает минимума, равного нулю, тогда и только тогда, когда две геометрические структуры равны.

Доказательство. Пусть имеется две равные структуры. Это значит, что они могут быть совмещены точно точками с одинаковыми индексами и их центры масс могут быть совмещены с началом координат. В этом случае точки с одинаковыми индексами этих структур имеют одинаковые координаты, определяемые векторами ${{{\mathbf{R}}}_{i}}$, $i = 1,\; \ldots ,\;N$. Одинаковые структуры с разными координатами точек с одинаковыми индексами можно получить, “разведя” эти структуры, повернув вторую структуру относительно первой и сместив каждую структуру, а именно. Повернем вторую структуру относительно второй на углы $\bar {\varphi }$, $\bar {\theta }$, $\bar {\psi }$. Сместим центр масс первой структуры в точку, определяемую вектором ${{{\mathbf{R}}}_{{1,0}}}$, а центр масс второй структуры сместим в точку, определяемую вектором ${{{\mathbf{R}}}_{{2,0}}}$. Координаты первой и второй структур будут иметь значения

где ${{Q}_{0}}(\bar {\varphi },\bar {\theta },\bar {\psi })$– матрица вида (3). Полученные в (14) координаты ${{{\mathbf{r}}}_{{1,i}}}$, ${{{\mathbf{r}}}_{{2,i}}}$ очевидно, реализуют любое расположение двух одинаковых структур в ${{R}^{3}}$.

Будем считать теперь, что в качестве исходных данных нам даны координаты двух одинаковых структур с, вообще говоря, различными координатами (14) одноименных точек, то есть $\sum\nolimits_{i = 1}^N \,{{w}_{i}}{{\left| {{{{\mathbf{r}}}_{{1,i}}} - {{{\mathbf{r}}}_{{2,i}}}} \right|}^{2}} \geqslant 0$. Имея в виду, что для координат двух одинаковых структур справедливо представление (14), применим к этим двум геометрическим структурам алгоритм сравнения 1)–3).

1) Вычислим координаты центров масс двух одинаковых структур. Согласно (10) и (14) для центра масс первой структуры получаем

так как ${{{\mathbf{R}}}_{i}}$ – координаты точек структуры относительно “центра масс” и поэтому $\sum\nolimits_{i = 1}^N \,{{w}_{i}}{{{\mathbf{R}}}_{i}} = 0$. Аналогично для центра масс второй структуры из (9) и (14) получаем

2) Функция сравнения (11) с учетом (14), (15) и (16) приобретает вид

Равный нулю минимум функции $U,$ в соответствии с (1) достигается, очевидно, при Отсюда следует, что $Q{{Q}_{0}} = E$ и $Q = Q_{0}^{{ - 1}}$. Так как матрица вращений ${{Q}_{0}}$ ортогональна, то обратная к ней матрица совпадает с транспонированной $Q_{0}^{{ - 1}} = Q_{0}^{{\text{т }}}$ и, таким образом, получаем матрицу поворота $Q$, реализующую равный нулю минимум функции $U$ вида $Q(\varphi ,\theta ,\psi ) = Q_{0}^{{\text{т }}}(\bar {\varphi },\bar {\theta },\bar {\psi })$. Так как для матриц $Q$ и ${{Q}_{0}}$ справедливо представление (3), то нетрудно видеть, что $Q_{0}^{{\text{т }}}(\bar {\varphi },\bar {\theta },\bar {\psi }) = Q(\pi - \bar {\psi },\bar {\theta },\pi - \bar {\varphi })$ и, тогда Отсюда следует, что для двух равных структур равный нулю минимум функции $U$ достигается в точке $(\varphi ,\theta ,\psi ) = (\pi - \bar {\psi },\bar {\theta },\pi - \bar {\varphi })$. При этом структуры совмещаются сдвигами (15), (16) и поворотом (18).

Обратное утверждение теоремы очевидно, так как из $U = 0$ для каждой точки $i$ с весом ${{w}_{i}} \ne 0$ в (11) следует, что координаты точек первой структуры совпадают с координатами точек второй структуры, то есть точки и структуры в целом совмещаются сдвигом на ${{{\mathbf{r}}}_{{1,0}}}$ первой и сдвигом на ${{{\mathbf{r}}}_{{2,0}}}$ второй структуры и поворотом, определяемым матрицей $Q$, второй структуры на углы $\varphi $, $\theta $, $\psi $. Теорема доказана.

Таким образом, согласно теореме 2 алгоритм сравнения 1)–3), включающий минимизацию функции сравнения, идентифицирует одинаковые структуры и, таким образом, может давать объективную оценку разницы в геометрии “близких” структур.

Осуществляя сравнение двух структур в соответствии с алгоритмом 1)–3), определение минимума функции сравнения $U(\varphi ,\theta ,\psi )$ по углам поворота после совмещения центров масс структур будем проводить численно.

4. ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ ОПТИМИЗАЦИИ

В этом разделе приводятся результаты вычислительного эксперимента по сравнению геометрических структур путем минимизации функции сравнения (11). Поскольку функция сравнения (11) не является выпуклой, для ее численной минимизации применялся метод нулевого порядка Розенброка [3], [4], хорошо зарекомендовавший себя при решении задач структурной химии [7]. В программе, реализующей алгоритм сравнения 1)–3) для минимизации функции сравнения (11) методом Розенброка, используется программа библиотеки оптимизации [9].

Эффективность предложенного и обоснованного выше алгоритма сравнения 1)–3) покажем на примерах приложений в структурной химии. В качестве геометрической структуры в этом случае будем рассматривать молекулы с упорядоченной структурой атомов-точек. При сравнении молекул нас будут интересовать различия в пространственном строении молекул при одной структурной формуле.

В рассмотренных ниже примерах приведены результаты сравнения молекул структуры лактида C₆H₈O₄, исследованной в работе [10], состоящих из $N = 10$ основных атомов (4 атома кислорода, 6 атомов углерода). Координаты атомов водорода не участвовали в расчетах, так как точность их определения ниже, чем у других атомов и для поставленной задачи совмещения молекул их присутствие несущественно. Геометрия этих молекул изображена на фиг. 1. Координаты молекул получены на основе данных рентгеноструктурного анализа. В рамках проведенного здесь эксперимента совмещаемым атомам в молекулах присваивается вес ${{w}_{i}} = 1$.

Пример 1. Для проверки состоятельности алгоритма сравнения 1)–3) и эффективности работы программы рассмотрим пример совмещения двух одинаковых молекул лактида. Начальные координаты ${{{\mathbf{r}}}_{{1,i}}}$, ${{{\mathbf{r}}}_{{2,i}}}$, $i = 1,\; \ldots ,\;10$, атомов двух молекул получены на основе формул (14), (3) смещением исходных координат одной и той же молекулы для формирования координат “первой” молекулы и смещением и поворотом исходных координат молекулы на углы $(\bar {\varphi },\bar {\theta },\bar {\psi }) = (\pi {\text{/}}3,\pi {\text{/}}6,\pi {\text{/}}2)$ для формирования координат “второй” молекулы. В табл. 1 приведены значения координат в ангстремах первой и второй молекулы. Точка минимума функции $U$ определена как $({{\varphi }_{0}},{{\theta }_{0}},{{\psi }_{0}}) = (90.00^\circ ,30.00^\circ ,120.00^\circ )$, что соответствует (18). Значение величины меры близости (которая в химии называется характеристикой близости), рассчитанное в соответствии с формулой (12) составило $s = 8.19 \times {{10}^{{ - 8}}}$ Å. По отношению к минимальному расстоянию ${{R}_{{{\text{min}}}}}$ между атомами молекулы

это составляет $s{\text{/}}{{R}_{{{\text{min}}}}} = 6.83 \times {{10}^{{ - 8}}}$, что можно считать практическим нулем. Таким образом, алгоритм 1)–3) одинаковые молекулы совмещает полностью.

Пример 2. Сравнение трех реальных молекул лактида. В табл. 2 приведены координаты в ангстремах атомов кислорода ${{O}_{j}}$, $j = 1,\; \ldots ,\;4$, и углерода ${{C}_{j}}$, $j = 1,\; \ldots ,\;6$, трех симметрически независимых (присутствуют в одном кристалле) молекул лактида C₆H₈O₄, одна из которых изображена на фиг. 1.

Применение алгоритма сравнения 1)–3) для сравнения молекул 1 и 2, 1 и 3, 2 и 3 дает результаты, приведенные в табл. 3. Полученные для этих расчетов точки минимума по углам также приведены в табл. 3. Имеющийся практический опыт работы по изучению конформации молекул, основанный на результатах сравнения значительного числа структур [11], привел в следующей условной классификации: $s \leqslant {{s}_{0}} = 0.1$ Å – молекулы приблизительно, или практически, равны, $0.1\;{{\AA}} < s \leqslant 0.2$ Å – молекулы близки, $s > 0.2$ Å – молекулы различны.

Анализ величин невязок $\Delta {{r}_{i}} = \left| {{{{\mathbf{r}}}_{{1,i}}} - {{{\mathbf{r}}}_{{2,i}}}} \right|,\;i = 1, \ldots ,10$, то есть расстояний между одноименными атомами в двух молекулах после совмещения (в точке минимума функции сравнения $U$), и величины $s$ при ${{s}_{0}} = 0.1$ Å позволяет сделать вывод, что первая и третья молекулы, а также вторая и третья практически равны по геометрии, так как в том и другом случае $s < {{s}_{0}}$. Наибольшие отличия имеются между молекулами первой и второй $(s = 0.11\;{{\AA}} > {{s}_{0}})$; эти молекулы можно считать близкими. Как видно из табл. 3, максимальные невязки имеют атомы заместителей (атомы вне цикла). Для более выпуклого выявления этих различий проведен дополнительный расчет: при совмещении молекул атомами циклов (расчет с весовыми множителями ${{w}_{i}} = 0$ для заместителей). В этом случае получили $s < 0.04$ Å. Это означает, что циклы молекул практически одинаковы. Выявленные различия в положении атомов в первой и второй молекуле при совмещении циклами видны на фиг. 2, который выполнен по координатам молекул в минимуме функции сравнения $U$.

Пример 3. Установление собственной симметрии молекулы.

На основе алгоритма сравнения 1)–3) проведена проверка собственной симметрии молекулы лактида (фиг. 1). Предполагаемая ось симметрии второго порядка проходит через центр цикла вертикально. Для применения алгоритма сравнения 1)–3) сформирована “вторая” молекула, в которой по сравнению с исходной молекулой изменен порядок нумерации атомов. Как видно на фиг. 1, в силу предполагаемой симметрии атом № 1 второй молекулы должен соответствовать атому № 2 первой, атом № 5 – атому № 7 и т.д. Общая перенумерация атомов во “второй” молекуле указана в табл. 4. В этой же таблице приведены результаты сравнения “двух” молекул. Поскольку характеристика сравнения в этом случае $s = 0.009\;{{\AA}} < {{s}_{0}} = 0.1$ Å, молекулы считаем равными, а следовательно, исходная молекула обладает собственной симметрией второго порядка с высокой точностью.

Приведенные выше расчеты проведены авторами с использованием программы COMPARISON на языке $C\sharp $, написанной Н.В. Бахтадзе, которой авторы выражают искреннюю признательность.