CC BY
36НАНОСИСТЕМЫ
DOI: 10.17725/rensit.2020.12.069
Описание вращательных движений молекул в компьютерных моделях воды с помощью кватернионов Волошин В.П., Наберухин Ю.И.
Институт химической кинетики и горения им. В.В. Воеводского СО РАН, http://www.kinetics.nsc.ru/ Новосибирск 630090, Российская Федерация Е-mail: voloshin@kinetics.nsc.ru, naber@ngs.ru
Поступила 05.01.2020,рецензирована 10.01.2020, принята 15.01.2020
Аннотация. Движение жёсткой молекулы в компьютерной модели может быть представлено как перемещение её центра масс и вращение вокруг оси, проходящей через этот центр. Для описания вращательного движения в данной статье используется алгебра кватернионов. Она позволяет изучать вращательное движение почти так же легко, как поступательное, так как кватернион содержит ось вращения и угол поворота вокруг неё практически в явном виде. Этот метод использован для изучения вращения молекул в молекулярно-динамической модели воды. Исследованы распределения углов между векторами пространственного перемещения и осями вращения молекул, а также углов обоих этих векторов с внутренними векторами молекул (дипольным моментом, вектором нормали, HH-вектором). Тем самым установлены корреляции между пространственными и вращательными перемещениями молекулы воды и определён общий характер её движения. В приложении даны конкретные формулы, позволяющие удобно работать с кватернионами в нашей задаче. Ключевые слова: молекулярно-динамическое моделирование, вода, пространственное смещение, угловое смещение, корреляция перемещения и вращения
УДК 532.74
Благодарности: Работа поддержана грантом РФФИ № 18-03-00045.
Для цитирования: Волошин В.П., Наберухин Ю.И. Описание вращательных движений молекул в компьютерных моделях воды с помощью кватернионов. РЭНСИТ, 2020, 12(1):69-80; DOI: 10.17725/rensit.2020.12.069._
Description of molecule rotation in computer models of water using quaternions
Vladimir P. Voloshin, Yuri I. Naberukhin
Voevodsky Institute of Chemical Kinetics and Combustion SB RAS, http://www.kinetics.nsc.ru/ Novosibirsk 630090, Russian Federation Е-mail: voloshin@kinetics.nsc.ru, naber@ngs.ru.
Received January 05, 2020; peer reviewed January 10, 2020; accepted January 15, 2020
Abstract. The motion of a rigid molecule in a computer model can be considered as the movement of its center of mass and rotation around an axis passing through this center. This article uses quaternion algebra to describe rotational motion. It provides an opportunity to study rotational motion almost as easily as translational, since the quaternion contains the axis and the angle of rotation in almost explicit form. We use this method to study the rotation of molecules in the molecular dynamics model of water. The distributions of angles between the translational displacement vectors and the axes of rotation of the molecules, as well as between these vectors and the internal vectors of the molecules (dipole moment, normal vector, and HH vector), are investigated. Correlations between the translational and rotational movements of the water molecule are established and the general character of its movement is determined. Appendix contains exact formulas and detailed descriptions of the procedures for using quaternions to describe the rotation of molecules.
70
' w ВОЛОШИН В.П., НАБЕРУХИН Ю.И.
НАНОСИСТЕМЫ
Keywords: molecular dynamics simulation, water, translational displacement, angular displacement, correlation of displacement and rotation axis PACS 45.20.dc, 47.11.-j
Acknowledgments: This work was supported by RFBR grant No. 18-03-00045.
For citation: Vladimir P. Voloshin, Yuri I. Naberukhin. Description of molecule rotation in computer models of water using quaternions. RENSIT, 2020, 12(1):69-80; DOI: 10.17725/rensit.2020.12.069.
Содержание
1. Введение (70)
2. Модели (72)
3. поступательные и угловые перемещения и их корреляции (72)
3.1. расчёт поступательных и угловых перемещений (72)
3.2. Углы между вектором перемещения и осью вращения (73)
4. ориентация внутренних векторов воды (731)
4.1. ориентация вектора дипольного момента (73)
4.2. ориентация HH-ВEктоPA (74)
4.3. ориентация вектора нормали молекулы воды (75)
5. заключение (75)
6. Приложение (77)
6.1. общая информация о кватернионах
(77)
6.2. Процедура 1. Сложение вращений (77)
6.3. Процедура 2. Вращение вектора (78)
6.4. Процедура 3. определение кватерниона текущей ориентации (78)
6.5. Процедура 4. Вычисление кватерниона углового смещения (79)
Литература (79)
1. ВВЕДЕНИЕ
Движение любой молекулы можно представить как поступательное движение её центра, обычно центра масс, и вращение вокруг этого центра. Поступательное движение активно изучается как экспериментально, так и с помощью компьютерного моделирования. Совершенно другая ситуация с изучением вращения. Так, в реальном эксперименте обычно рассматривается не полное вращение молекулы, а переориентация лишь одного
её вектора, например, вектора дипольного момента. В компьютерной модели мы знаем координаты каждого атома молекулы в любой момент времени и потому можем описать вращение молекулы в целом. Если молекула жёсткая, то есть изменением расстояний между отдельными её атомами можно пренебречь, то её вращение можно описать как вращение твёрдого тела в механике [1], то есть как поворот молекулы как целого вокруг оси вращения, проходящей через её центр. Именно этот случай мы будем рассматривать в этой статье.
Текущее положение молекулы в момент времени / описывает радиус-вектор, проведённый из начала координат в точку с текущими координатами её центра. Перемещение этого центра за некоторый интервал времени описывается вектором, равным разности радиус-векторов конечного и начального положений
г = г — г . (1)
2-12 1 V/
Определив координаты центра масс и вычтя их из координат каждого атома молекулы, мы получим координаты этих атомов относительно центра молекулы. Изменения относительных координат атомов во времени можно описать как вращение молекулы вокруг её неподвижного центра. Наиболее часто для этого используются углы Эйлера, матрицы вращения или кватернионы [2]. Параметры выбранного способа описания вращения образуют систему координатных осей некоторого пространства ориентаций, которое можно использовать для описания ориентации и вращения молекул подобно тому, как для описания положения и перемещения этих молекул мы используем декартову систему координат "обычного" пространства положений. Для углов Эйлера пространство ориентаций
имеет три координатных оси, для кватернионов четыре, для матриц вращений девять.
Опыт молекулярно-динамических расчётов показывает, что язык кватернионов очень удобен для расчёта вращений [2-4]. Нужные нам свойства кватернионов и расчётные процедуры мы даём в Приложении, а здесь наметим только общую идеологию их применения для изучения вращательных движений молекул. Основное удобство использования кватернионов заключается в том, что они допускают представление вращения через самые наглядные параметры — вектор оси вращения u и угол поворота ф вокруг этой оси, так как кватернион можно записать в виде
q = cos(ф/2) + u•sin(ф/2). (2)
Для однозначного описания ориентаций необходимо задать нулевую (базовую) ориентацию молекулы, играющую ту же роль, что и начало координат в пространстве положений. Кватернион вращения q переводящего молекулу из этой нулевой ориентации в текущую, используется в качестве координат этой текущей ориентации в пространстве ориентаций, то есть ql — аналог радиус-вектора г1 положения молекулы. Сравнивая координаты двух текущих ориентаций молекулы, мы определяем кватернион вращения q2 переводящий молекулу из одной ориентации в другую — аналог вектора поступательного перемещения г21. Кватернион вращательного перемещения q21 вычисляется из кватернионов начальной и конечной ориентации, q1 и q2, по особой формуле
42-1=ч 2 • я;> (3)
отличной от формулы (1) для подсчёта пространственных перемещений. Алгоритмы определения кватернионов текущей ориентации и кватернионов вращения описаны в Приложении.
Описание жёсткой молекулы с помощью координат центра масс и кватерниона ориентации абсолютно полно и однозначно. Преобразование координат атомов нулевой (базовой) конфигурации с помощью
кватерниона ориентации молекулы позволяет получить текущие относительные координаты её атомов, а прибавление к ним координат центра масс этой молекулы восстанавливает исходные молекулярно-динамические
координаты её атомов. Таким образом, вместо 9 координат для 3 атомов молекулы воды нам достаточно знать 3 координаты центра масс и 4 компоненты кватерниона. Для более крупных молекул экономия будет ещё больше: этих 7-ми параметров хватит, например, для замены 36 координат атомов молекулы бензола.
Кватернионы были впервые предложены Уильямом Гамильтоном в 1843 году [5], с тех пор активно изучались и к настоящему времени подробно описаны во многих математических учебниках [6]. Ныне они активно используются для описания вращения в компьютерных играх, и в интернете существует множество публикаций на эту тему [7-9]. Используются кватернионы и в механике [10]. Применение кватернионов позволило создать альтернативный подход к моделированию динамики жёстких молекул [3], [11], который, впрочем, в настоящее время используется очень редко. Сравнение методов анализа вращательной подвижности с помощью кватернионов, углов Эйлера и матриц вращения [2] показало, что именно использование кватернионов даёт наиболее устойчивые решения. Крайне полезными кватернионы оказались также при определении ориентаций молекул, не являющихся жёсткими [12-14]. Однако ни в одной из упомянутых работ нет простой и ясной инструкции, показывающей по шагам, как с помощью кватернионов описывать вращение жёстких молекул в молекулярно-динамических моделях. Мы надеемся, что такой инструкцией могут служить процедуры, представленные нами в Приложении. Они легко преобразуются в компьютерные программы. Использование этих инструкций позволило очень просто описать вращения молекул воды в ходе молекулярно-динамического моделирования. Результаты этих расчетов изложены в разделах 3 и 4 данной статьи.
' ^ ВОЛОШ11Н В.П., НАБЕРУХ11Н Ю.11.
НАНОСИСТЕМЫ
2. МОДЕЛИ
В работе использована молекулярно-дпнампческая модель воды, приготовленная с помощью пакета молекулярно-дпнампческого моделирования ЬАММРБ [15]. Модель содержала 8000 молекул с потенциалом взаимодействия Т1Р4Р/2005 [16] в кубическом боксе с периодическими граничными условпямп, шаг моделирования 2 фс, давление 1 бар, температура 300 К. После предварительной релаксации в течение 1 не рассчитывалась траектория длительностью 4 не, для которой через каждые 200 фс записывались мгновенные конфигурации, то есть всего 20 001 конфигурация.
Для каждой молекулы каждой записанной конфигурации были рассчитаны координаты центров масс молекул, а также кватернионы пх ориентации Процедуры вычисления кватернионов орпентацпй описаны в Приложении. Расчёт пространственных п угловых перемещений производился для интервалов времени 0.2, 0.4, 0.6, 1, 2, 4 п 6 пс. Для каждого интервала в качестве начальной п конечной конфигураций использовались все возможные пары, разделённые таким интервалом; таким образом, для наиболее короткого интервала было использовано по 20 000 наборов смещений для каждой пз 8000 молекул, а для наиболее длинного по 19 970 наборов. Мы полагаем, что тщательное усреднение на траектории, длительность
которой существенно превышает время жпзнп даже самых долговременных корреляций в моделях использованного размера, позволило получить очень надёжно усреднённые хорошо воспроизводимые результаты.
3. ПОСТУПАТЕЛЬНЫЕ И УГЛОВЫЕ ПЕРЕМЕЩЕНИЯ И ИХ КОРРЕЛЯЦИИ 3.1. Расчёт поступательных и угловых
перемещений
Для расчёта поступательных перемещений на каждом интервале времени определялся вектор перемещения центра масс каждой молекулы, а также его длина. Эта информация сохранялась для последующего использования. Расчёт угловых перемещений в свою очередь заключался в определении кватернионов вращения на основе кватернионов начальной п конечной орпентацпй по формуле (3), пз которых при помощи формулы (2) извлекались п сохранялись для дальнейшего использования углы поворотов п направления осей вращения. Описание процедур определения орпентацпй п вычисления кватернионов вращений приведено в Приложении. Вероятностные распределения величин поступательных п угловых перемещений за разные интервалы времени от 0.2 до 6 пс показаны на Рис. 1, слева поступательные перемещения, справа — угловые.
Формы всех распределений очень похожи. Отношение величины перемещения
Рис. 1. Распределения поступательных (а) и угловых (Ь) перемещений молекул воды при Р = / бар, Т = 300 К за разные промежутки времени (длительность интервапов указана около кривых). Срый пунктир — результат аппроксимации
формулой у — ах:-ехр(-Ьх:).
0.2 ПС
Поступательное перемещение
НАНОСИСТЕМЫ
к длительности интервала времени есть средняя скорость молекул на этом интервале. Поэтому распределения пространственных перемещений хорошо описываются распределением Максвелла (Рпс. 1а). Но распределения угловых перемещений существенно отклоняются от этой формы (Рпс. 1Ь).
Поскольку невозможно различить вращение на угол ф п на дополнительный к нему угол 2тт — ф, мы всегда выбирали меньший угол поворота, не превышающий тт. Для больших интервалов времени это приводит к обрывам на распределениях углов. Значение распределения вблизи обрыва формируется не только вращениями, имеющими такой угол поворота, но также вращениями с углом, дополняющим данный до 2тт. Впрочем, заметное количество молекул с такими углами было только при самом большом интервале, п поэтому существенного влияния на основные результаты это не оказало.
3.2. Углы между вектором перемещения и осью вращения.
Помимо угла вращения, кватернион вращения содержит в явном виде ось данного вращения, а значит мы можем рассмотреть углы между осью вращения п другими векторами в процессе движения молекулы. На Рис. 2 показаны распределения косинусов углов между вектором перемещения п осью
1.34
>
о т 1.2-
>3
т 1.1 -
Т
>-.
Ц
о 1.0-
ф
X 0.9-
(II
э 0.8-
С)
х
н о 0.7-
0.6-
а
..... .,»* 1* / / / / * / * / / М = 0.6 пс \ *» ..... \ Ч •« \ ч \ ч \ ч \ N -быстрые \ ----средние \
-1.00 -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 Косинус угла между вектором перемещения и осью вращения
1.00
вращения в модели воды. Распределение нормировано на распределение для случайно ориентированных векторов, то есть оно показывает, насколько чаще/реже данный угол реализовался в модели по сравнению со случайной ориентацией. Мы разделили молекулы по величине пх пространственного перемещения, поместив в каждую группу по 1 /3 от полного числа молекул, п изобразили отдельное распределение для каждой группы.
На Рпс. 2а представлены распределения для интервала 0.6 пс, а на Рпс. 2Ь для интервала, длительность которого больше в 10 раз. Они показывают, что наиболее часто вектор перемещения п ось вращения взаимно перпендикулярны друг другу. Но преимущество этой ориентации не слишком большое: даже для быстрых молекул перпендикулярная ориентация превышает случайную на 20%, а параллельная ниже случайной на 30%. Увеличение интервала времени уменьшает отклонение от случайного для всех групп.
4. ОРИЕНТАЦИИ ВНУТРЕННИХ ВЕКТОРОВ МОЛЕКУЛ ВОДЫ 4.1. Ориентация вектора дипольного момента
Кватернион текущей ориентации содержит всю информацию об ориентации молекулы, а значит содержит п информацию о направлениях её внутренних векторов.
Рис. 2. Распределение косинусов ума между векторами перемещения и осями вращения. Жирная линш относится к наиболее быспрым молекулам, пунктир — к средним, точки — к наиболее медленным. Труппы наиболее быстрых, средних и наиболее медленных молекул содержат по 1/3 от полного числа молекул. Распределения нормированы
на случайное.
74
' ^ ВОЛОШ1IH В.П., НАБЕРУХ1IH Ю.11.
НАНОСИСТЕМЫ
Например, текущее направление вектора дппольного момента, который мы в нулевой ориентации воды направили вдоль осп X, можно получить, преобразовав вектор {1,0,0} кватернионом ориентации согласно Приложению 6.3. Вычислив углы, во-первых, между направлением диполя п вектором перемещения, а, во-вторых, между ним же п осью вращения, можно найти преимущественные ориентации дппольного момента относительно векторов обоих этих движений. Однако движения происходят за некоторый интервал времени, в течение которого направление диполя может измениться. Мы полагаем, что наиболее разумно использовать внутренние вектора молекулы пз середины этого интервала. Определив угол поворота п ось вращения на всём интервале, мы вычисляли срединную ориентацию молекулы, повернув её начальную ориентацию на половину угла этого поворота вокруг той же самой осп. На Рис. 3 показано распределение косинусов углов между диполем пз этой срединной ориентации с вектором перемещения (по горизонтали) п с осью вращения (по вертикали), на этот раз без деления молекул на группы. Как оказалось,
Рис. 3. Двухмерное распределение косинусов углов между диполем из средней ориентации и вектором перемещения (по горизонт сит) и тем же диполем и осью вращения (по вертикали). Распределение нормировано на случайное: серые заштрихованные полоски — области, заполненные с (светло-серые) от тех,
(тёмно-серые и
наиболее вероятным является направление дппольного момента, перпендикулярное как осп вращения, так п вектору перемещения, причём перпендикулярность дппольного момента п вектора перемещения более жёсткая. Поскольку вектор перемещения п ось вращения наиболее часто также перпендикулярны друг другу, то, видимо, все три этих вектора наиболее часто взаимно перпендикулярны.
Интересно, что помимо главного максимума присутствуют два заметных максимума выше п ниже него, для которых дппольный момент по-прежнему перпендикулярен вектору перемещения, но параллелен осп вращения. Смысл этих максимумов рассмотрим в Заключении.
4.2. Ориентация НН-вектора
Рассмотрим ориентацию вектора НН, соединяющего центры атомов водорода молекулы воды. В нулевой ориентации он направлен вдоль осп У. Данный вектор перпендикулярен дппольному моменту, п вместе с ним лежит в плоскости молекулы. Вновь построим двухмерное распределение косинусов углов между данным вектором в средней ориентации интервала времени, п вектором перемещения центра масс молекулы, а также между ним п осью её вращения на этом интервале. Как следует пз Рис. 4, вектор НН
меньше случайной случайной
Рис. 4. Двухмерное распределение косинусов углов между НН-вектором и вектором перемещения (по горизонт сит) и НН-вектором и осью вращения (по вертиксиш). Нормировка, раскраска более и менее вероятных областей, а также границы между ними как на Vue. 3.
наиболее часто параллелен осп вращения. От угла между НН-вектором п вектором перемещения на данном интервале времени распределение почти не зависит. 4.3. Ориентация вектора нормали молекулы воды
На Рис. 5 показаны распределения косинусов углов между вектором нормали к плоскости молекулы воды п вектором перемещения (по горизонтали) п тем же вектором нормали п осью вращения (по вертикали). Вектор нормали определялся как перпендикуляр к дппольному моменту п вектору, проведённому от первого водорода молекулы воды ко второму, то есть совпадал с направлением осп Z нулевой ориентации. Наиболее вероятны движения, при которых вектор перемещения параллелен нормали, а ось вращения перпендикулярна ей, причём эта вторая закономерность более сильная.
5. ЗАКЛЮЧЕНИЕ
Итак, сравнивая направления вектора перемещения, осп вращения, а также собственных векторов молекулы пз срединной ориентации интервала, мы обнаружили: 1. Вектор перемещения п ось вращения преимущественно перпендикулярны друг другу.
Косинус угла между нормалью At - 1 пс и вектором перемещения
Рис. 5. Двухмерное распределение косинусов умов между вектором нормапи к плоскости молекулы воды и векторам перемещения (по горизонта/ш) и осью вращения (по вфтикали). Нормировка, раскраска более и менее вероятных областей, а также границы между ними как на Vue. 3.
2. Вектор дппольного момента в большинстве случаев перпендикулярен вектору перемещения. Наиболее часто при этом он перпендикулярен п осп вращения, однако присутствует заметная доля молекул, для которых они параллельны.
3. Вектор НН у большинства молекул параллелен осп вращения, п почти не зависит от ориентации вектора перемещения.
4. Вектор перемещения параллелен вектору нормали, а ось вращения перпендикулярна ему, причём эта вторая закономерность более сильная.
Эти закономерности наблюдаются на любых рассмотренных интервалах времени, однако наиболее сильно они проявляются для наиболее подвижных молекул на интервале движения порядка 1 пс.
На Рис. 6 показана вероятная картина совместного перемещения центра
масс молекул воды п её вращательного движения, соответствующая обнаруженным закономерностям. Здесь схематически изображена ориентация молекулы воды в середине интервала перемещения, показаны векторы движения п собственные векторы молекулы, в том числе направление диполя этой срединной ориентации. На данном рисунке молекула "качается" на двух донорных водородных связях. Остальных связей нет, плп они существенно слабее.
Впрочем, все перечисленные выше закономерности будут наблюдаться п в том случае, если обе связи акцепторные, но тогда ось вращения будет расположена со стороны неподелённых пар электронов. Чтобы проверить, правильно ли мы представили схему наиболее вероятного движения, построим ещё раз распределение углов между диполем п вектором смещения, а также с осью вращения, аналогичное Рпс. 3, используя на этот раз диполь пз начальной ориентации молекулы на интервале. Если наше предположение правильное, то преимущественный угол между начальным диполем п вектором смещения должен быть тупым, если же ось вращения
ВОЛОШИН В.П.,НАБЕРУХИНЮ.И.
НАНОСИСТЕМЫ
Рис. 6. 1предполагаемый виО движения молекулы удовлетворяющего всем выявленным корреляциям величин и умов между векторами. Изображена молекула воды в срединной ориентации интервала движения: О — атом кислорода этой молекулы, Н — её атомы водорода, О' и О" — кислороды молекул, с которыми данная молекула образует донорные водородные связи. Тонкая спгрелка вверх указывает направление вектора дппольного момента в срединной ориентации, пунктирная стре>жа влево вверх соответствует диполю в начале интервала. Вектор перемещения ценпгра масс (серая спгрелка с чёрным контуром) параые^чен вектору нормапи (ч'ёрная жгфная стре^жа), НН-вектор (ч'ёрная стре>жа, проходящая через ценпгры водородов) параыелен оси вращения (пунктир с круговой спгрелкой,
проходит со стороны неподеленных пар, то угол будет острым. Представленное на Рис. 7 двухмерное распределение полностью подтверждает наше первоначальное предположение: главный максимум этого распределения соответствует углу 140° между вектором перемещения п начальным дппольным моментом. Раскачивание молекулы на двух сильных донорных Н-связях происходит значительно чаще, чем на двух акцепторных связях.
Вспомним про дополнительные максимумы на Рпс. 3, соответствующие орпентацпям дппольного момента, перпендикулярным
Рис. 7. Двухмерное распределение косинусов углов между диполем в начаче интервала перемещения и вектором перемещения (по горизонтами) и этим же диполем и осью вращения (по вертикачи). Нормировка, раскраска более и менее вероятных областей, а также границы между ними как на Vue. 3. Основной максимум соответствует углам 90° между начальным диполем и осью вращения и 140° между этим диполем и вектором перемещения.
вектору перемещения п параллельным осп вращения. Как оказалось, на Рпс. 7 они остались на тех же местах. Их присутствие можно объяснить тем, что иногда наиболее сильными связями, на которых раскачивается молекула воды, являются одна донорная п одна акцепторная. Тот факт, что на обоих рисунках дппольный момент для этих максимумов почти параллелен осп вращения п перпендикулярен вектору перемещения означает, что со стороны неподелённой пары акцепторной связи ось вращения расположена значительно ближе к центру масс, чем со стороны водорода донорной связи. Если бы эти расстояния были одинаковы, дополнительные максимумы располагались бы при значениях косинусов +0.7 по каждой осп. Если Н-связей больше, или они одинаково сильные, движение становится более скованным, величины перемещений уменьшаются, п молекула выходит за пределы основных максимумов распределений.
Попробуем определить, насколько далеко от центра масс располагается ось вращения. Считая, что перемещение центра масс, перпендикулярное осп вращения,
НАНОСИСТЕМЫ
происходит исключительно в ходе вращения молекулы согласно схеме на Рис. 6, можно оценить радиус вращения по формуле R = L|2sin(ф|2), где Ь—длина перпендикулярного к оси перемещения, а ф — угол поворота. Для интервала 0.2 пс получается среднее значение К02пс = 1.28 А. По мере увеличения интервала растёт и средний радиус: К1 пс = 1.59 А, а Я = 1.65 А. Любое из этих
2 пс
значений соответствует расположению оси заметно дальше центров водородов, почти около центров кислородов, с которыми рассматриваемая молекула имеет водородные связи. Увеличение радиуса с ростом интервала времени означает, что движение центра масс молекулы на более длинных интервалах нельзя свести исключительно к её вращению. Каждая молекула одновременно участвует и в других движениях, в том числе в коллективных движениях разных масштабов. Однако полученные значения в основном подтверждают схему, приведённую на Рис. 6.
6. ПРИЛОЖЕНИЕ
6.1. общая информация о кватернионах
Будем рассматривать только те свойства кватернионов, которые необходимы для описания ориентации и вращения жёстких молекул в компьютерных моделях. Кватернион представляет собой гиперкомплексное число с одной действительной частью и тремя мнимыми с разными мнимыми единицами: q = w + ix + )у + kz, (4)
где ь ), к — мнимые единицы. Квадрат любой из них равен —1, а их умножение друг на друга производится аналогично векторному умножению единичных векторов вдоль осей трёхмерной правой декартовой системы координат:
¡¡ = _у = кк = -1;
у = -^ = к; 'к = -к'= ¡; к = -¡к = '
Четыре компоненты кватерниона позволяют представить его в виде точки в четырёхмерном
(5)
пространстве: q = (ш, х, у, z}. В данной работе мы будем использовать кватернионы, нормированные на единицу, то есть такие, для которых ш2 + х2 + у2 + z2 = 1.
Кватернион также можно представить в виде пары, состоящей из скаляра № и вектора V = (х, у, z}, проекции которого на мнимые оси {, ), к считаются направленными вдоль осей X, У и Z, соответственно: q = № + V. Нормированный кватернион, записанный в такой форме, наиболее наглядно описывает вращение. Если молекула повернулась на угол ф по часовой стрелке вокруг оси, направленной вдоль единичного вектора и, то это вращение описывается нормированным кватернионом вращения
q = cos(ф/2) + и^т(ф/2). (6)
Комплексно сопряжённый ему кватернион имеет вид
q, = ш — ix — )у — kz = cos(ф/2) — и^т(ф/2).
Он описывает обратное вращение, то есть вращение на тот же угол вокруг оси, ориентированной в противоположном направлении.
6.2. процедура 1. Сложение вращений
Пусть молекула последовательно совершила два вращения, кватернион первого q1 = (ш1, х1, у1, z1}, второго q2 = (ш2, х2, у2, z2}. Кватернион суммарного вращения равен произведению второго кватерниона на первый: q1+2 = q2•q1. Умножение кватернионов не коммуникативно — менять порядок сомножителей нельзя. Для вычисления произведения запишем кватернионы в виде сумм (4): q1+2 = (ш2 + ix2 + )у2 + кг2)^(ш1 + ix1 + )у1 + kzl). Раскрыв скобки с учётом правил перемножения мнимых единиц (5), получим компоненты итогового кватерниона q1+2:
ш = ш •ш — х •х — у •у — z•z:
1 + 2 2 1 2 1 У 2 -М 2 1'
х = ш •х + х •ш + у •z — z •у ;
1+2 2 1 2 1 2 1 2 1
у = ш •у — х •z + у •ш + z •х ;
1+2 2 1 2 1 2 1 2 1
1+2
= + - у2^х1 +
После вычислений нужно произвести перенормировку итогового кватерниона, чтобы устранить возможную неточность записи компонент исходных кватернионов и
78 ВОЛОШИН В.П., НАБЕРУХИН Ю.И.
НАНОСИСТЕМЫ
ошибки численного вычисления.
6.3. процедура 2. Вращение вектора
Если q = ш + ix + )у + kz — кватернион вращения, в результате которого вектор ^ = ix0 + )у0 + kz0 обратился в вектор ^ = ix1 + )у1 + kz1, то для них выполняется равенство ^ = q•v0•q,. Вычислим это произведение в два этапа: вначале рассчитаем вспомогательный кватернион q2 = = ^ + )у0 + kz0)•(w -ix — )у — kz). Раскрыв скобки и сгруппировав однотипные компоненты, получим:
= Vх + УоТ +
x2 = хо^—Уo•z +
у2 = V21 + —Vх; ^ = —хоТ + утх + V™
Компоненты итогового вектора ^ = q•q2 = (ш
+ ix + )у + kz)•(w2 + ix2 + )у2 + kz2) равны: x1 = + х^2 + у^2 — z•y2; у1 = ш-у2 — х^2 + у-ш2 + z•x2;
Z1 = w•z2 + хТ2 — Ух2 + ^2.
Действительная компонента этого
произведения тождественно равна нулю. Неточность представления сомножителей и ошибки численного вычисления могут слегка изменить длину вектора, поэтому, если данная процедура используется многократно, итоговый вектор следует нормировать заново.
6.4. процедура 3. определение кватерниона текущей ориентации
Подобно тому, как текущее положение молекулы описывается радиус-вектором, проведённым от начала координат к точке расположения её центра масс, текущую ориентацию молекулы будем описывать кватернионом, представляющим вращение от нулевой (базовой) ориентации к текущей.
Ориентацию жёсткой молекулы можно описать при помощи двух векторов, проведённых из центра масс молекулы к центрам двух её атомов, если эти векторы не лежат на одной прямой. В наших работах нулевую ориентацию любой молекулы мы определяли по одной и той же методике: первый выбранный атом (атом типа 1) располагался относительно центра масс в направлении оси X, а второй (атом типа 2) располагался на
плоскости ХУ при положительных значениях У Все прочие атомы для определения ориентации не использовались. В молекуле воды атомом типа 1 мы считали кислород, а атомом типа 2 — первый из водородов. Таким образом, в нулевой ориентации дипольный момент молекулы воды оказался направлен вдоль оси X, а линия, соединяющая центр второго и первого водородов, вдоль оси У. Ось Z при этом представляла собой нормаль к плоскости молекулы.
Для описания текущей ориентации молекулы мы использовали кватернион вращения, переводящего её из нулевой ориентации в текущую. Таким образом, компоненты этого кватерниона играли роль "текущих координат молекулы в пространстве ориентаций". Нулевая ориентация
описывается кватернионом нулевого вращения q0 = (1,0,0,0}. Для любых других ориентаций такой поворот и соответствующий ему кватернион нужно вычислять. Оказалось, что проще рассчитать кватернион, описывающий обратное вращение, а затем взять комплексно сопряжённый от него.
Составим это обратное вращение из двух простых последовательных вращений. Пусть в текущей ориентации направление от центра молекулы к центру её атома типа 1 описывается вектором единичной длины г10 = (х10,у10^10}. Первое вращение должно превратить его в вектор (1,0,0}. Это можно получить, вращая молекулу по часовой стрелке вокруг вектора = (0, z, —у10}, представляющего собой векторное произведение начального и конечного вектора, а потому перпендикулярного им обоим. Угол поворота определим из скалярного произведения тех же векторов: ф1 = агссоэ^^). Определив длину вектора вдоль оси вращения
= у]^о + Уо , запишем единичный вектор вдоль этой оси ^ = (0, z1, —у1}, где z1 = z10/d1 и у1 = y10/d1. Таким образом, кватернион первого вращения будет иметь вид ql = (^(ф^^О^т^^-у^т^^)}.
Пусть исходное положение атома типа 2 относительно центра масс описывается
НАНОСИСТЕМЫ
единичным вектором г20 = {х20,У20^20}. Во время первого вращения этот вектор изменится на г2 = {х2,у2^2}. Проекции этого вектора можно рассчитать по алгоритму, описанному в Процедуре 2 (6.3). Второе вращение будем производить вокруг оси X (чтобы не сместить атом типа 1) таким образом, чтобы центр атома типа 2 расположился на плоскости ХУ при положительном значении У Для этого достаточно повернуть молекулу на угол (р2 = агссо8(у2 / ) вокруг оси {1,0,0}.
Здесь + ¿2 — длина проекции вектора
г2 на плоскость УХ, перпендикулярную оси вращения X, а у2, в свою очередь, её проекция на ось У. Поворот именно на этот угол уложит вектор на плоскость ХУ. Таким образом, кватернион второго поворота должен быть равен q2 = {cos(ф2/2), sin(ф2/2), 0, 0}. Итоговый кватернион поворота данной молекулы из текущей ориентации в нулевую равен q2'q1. В качестве "координат ориентации" молекулы будем использовать кватернион обратного вращения, из нулевой в текущую, то есть кватернион текущей ориентации будет равен q
= ^^У = q,l•q,2.
6.5. процедура 4. Вычисление кватерниона углового смещения
Пусть q1 = х1, у1, z1} — кватернион
начальной ориентации молекулы, а q2 = х2, У2, z2} — кватернион её конечной ориентации. Кватернион вращения, переводящий молекулу из начальной ориентации в конечную, представляет собой их "разницу": q21 = q2•q1,
= К + Ц + 1У2 + kz2)'(wl - Ц - 1У1 - КУ
Раскрывая скобки и суммируя однотипные слагаемые, получим компоненты кватерниона
q2-l:
^-1 = W2'W1 + Х2'Х1 + У2'У1 + ^
х = —w 'х + х ^ — у ^ + z 'у ;
2-1 2 1 2 1 У 2 1 2 /1'
У2-1 = -^Т1 + Х2 'Z1 + У2^1 - ^^ Z2-1 = - Х2'У1 + У2'Х1 + Z2■W1.
Если данная процедура используется многократно, итоговый кватернион следует заново нормировать на единицу.
Отметим, что разницу ориентаций тождественно описывают два кватерниона:
q — cos(w /2) + v и а — —cos(w /2) —
*trot \rrof' rot >-2rot Уг rot' '
v . Первый описывает поворот на угол ф
rot ^ J- J ' rot
вокруг оси v а второй на угол 2п — вокруг оси —■v . В результате вычислений может получиться любой из них, и без привлечения промежуточных ориентаций мы не можем определить, какой угол является правильным. Однозначное определение угла поворота возможно только в том случае, если он заведомо меньше п. Для такого угла cos((prJ2) > 0. Потому, если после вычисления кватерниона углового смещения действительная часть оказывалась отрицательной, мы меняли знаки всех компонент кватерниона на обратные.
ЛИТЕРАТУРА
1. Савельев ИВ. Курс общей физики, Том I, Механика, колебания и волны, молекулярная физика. М., Наука, 1970.
2. Zhao F, van Wachem BGM. A novel Quaternion integration approach for describing the behaviour of non-spherical particles. Acta Mechanica, 2013, 224:3091-3109.
3. Rapaport DC. Molecular dynamics simulation using quaternions. J. Computational Physics, 1985, 60:306-314.
4. Rapaport DC. The art of molecular dynamics simulation. Cambridge Univ. Press, 2004.
5. Hamilton WR. On a new species of imaginary quantities connected with a theory of quaternions. Proc. Royal Irish Acad, 1844, 2:424434. URL: http: //www.maths.tcd.ie/pub/ HistMath/People/Hamilton/Quatern1.
6. Кантор ИЛ, Солодовников АС. Гипркомплексные числа. М., Наука, 1973, 144 с.
7. Mukundan R. Quaternions: From Classical Mechanics to Computer Graphics, and Beyond.
Proceedings of the 7th Asian Technology Conference in Mathematics, 2002.
8. Vince J. Quaternions for Computer Grcphics. Springer-Verlag, London, 2011. DOI 10.1007/978-0-85729-760-0.
9. Serrano H. Quaternions in Computer Graphics. URL: https://www.haroldserrano.com/blog/ quaternions-in-computer-graphics, 2015.
10. Голубев ЮФ. Алгебра кватернионов
80 ВОЛОШИН В.П., НАБЕРУХИН Ю.И.
НАНОСИСТЕМЫ
в кинематике твердого тела. Препринт ИПМ им. М.В.Келдыша, № 39, 2013, 23 с., URL: http://library.keldysh.ru/preprint. asp?id=2013-39.
11. Karney CFF. Quaternions in molecular modeling, URL: http://arxiv.org/abs/physics/0506177v, 2006.
12. Kneller GR. Superposition of molecular structures using quaternions. Molecular Simulation, 1991, 7:113-119. DOI: 10.1080/08927029108022453.
13. Chevrot G, Calligari P, Hinsen K, Kneller GR. Least constraint approach to the extraction of internal motions from molecular dynamics trajectories of flexible macromolecules. J. Chem. Phys., 2011, 135:084110. DOI: 10.1063/1.3626275.
14. Chevrot G, Hinsen K, Kneller GR. Model-free simulation approach to molecular diffusion tensors. J. Chem. Phys, 2013, 139:154110. DOI: 10.1063/1.4823996.
15. Plimpton S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comp. Phys., 1995, 117:1-19. URL: http://lammps.sandia. gov.
16. Abascal JLF, Vega C. A general purpose model for the condensed phases of water: TIP4P/2005. J. Chem. Phys., 2005, 123:234505. DOI: 10.1063/1.2121687.
Волошин Владимир Петрович
к.ф.-м.н, с.н.с.
Институт химической кинетики и горения им.
В.В. Воеводского СО РАН
3, ул. Институтская, Новосибирск630090, Россия
voloshin@kinetics.nsc.ru.
Наберухин Юрий Исаевич
д.х.н., г.н.с.
Институт химической кинетики и горения им. В.В. Воеводского СО РАН, 3, ул. Институтская, Новосибирск630090, Россия naber@ngs.ru.
