CC BY
124
В. С. Неверов, А. В. Комолкин, Т. Г. Волкова
ИССЛЕДОВАНИЕ ВЛИЯНИЯ СТРУКТУРНОЙ ИЗОМЕРИИ НА МОЛЕКУЛЯРНУЮ ПОДВИЖНОСТЬ ЖИДКИХ КРИСТАЛЛОВ МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ*
Введение. Свойства жидких кристаллов (ЖК), т. е. состояния вещества, характеризующегося, с одной стороны, текучестью, а с другой — анизотропией свойств, определяются его молекулярной структурой [1, 2]. При этом структурные изомеры, отличающиеся порядком заместителей в центральном жёстком фрагменте, могут иметь
сильно различающиеся свойства. В качестве примера таких веществ, образующих ЖК-фазу, можно привести МББА (4-метоксибен-зилиден-4/-н-бутиланилин) [3] и БОБТ (4--н-бутилоксибензилиден-4/-толуидин) [4]. На рис. 1 приведены схематические изображения обеих молекул, а в табл. 1 — температуры фазовых переходов.
Таблица 1 Температуры фазовых переходов
Вещество Тк^м, К Тм^і, К АТ, К
МББА 293 320 27
БОБТ 339 341 2
40
37 35,
(16
'38 36 21
БОБТ
30 31
XX Ь 3 /4 6 XX
2^4 У^Г28
57
32
33
39
23 25
Рис. 1. Схематические изображения молекул МББА и БОБТ
* Работа выполнена при финансовой поддержке проекта Министерства образования и науки РФ «Развитие механизмов интеграции Ивановского государственного университета и Института проблем химической физики РАН» (РНП.2.2.1.1.2820), а также частично поддержана грантом РФФИ № 09-03-01105-а
© В. С. Неверов, А. В. Комолкин, Т. Г. Волкова, 2011
Оба вещества принадлежат к классу оснований Шиффа. Основания Шиффа (азоме-тины, имины) — органические соединения общей формулы RR1C=NR2, где ИиИ1 — водород, алкил или арил, R2 — алкил или арил. Многие азометины обладают свойствами жидких кристаллов [6] и находят широкое применение в качестве температурных датчиков, индикаторов, модуляторов света, для создания ЖК-дисплеев [5-11].
В обеих молекулах к фрагменту бензилиденанилина с одной стороны алифатическая цепь присоединена напрямую, а с другой — через атом кислорода. Различие заключается в том, что в МББА к атому кислорода присоединена метильная группа, а к бензольному кольцу — цепочка С4Н9, в БОБТ — наоборот.
Целью работы является выяснение физических основ столь различного поведения этих веществ путём сравнительного анализа конформационной структуры и динамики молекул, локальной структуры веществ в жидкокристаллической фазе. Исследование выполнялось методом компьютерного моделирования молекулярной динамики (МД).
Ограничения классического метода МД. Классический метод компьютерного моделирования молекулярной динамики [12] является одним из наиболее важных инструментов теоретического изучения структуры и динамики ЖК, биологических макромолекул, жидкостей, газов и других молекулярных систем. Важную роль он играет и при установлении детальной структуры молекул и комплексов, определённой методами ЯМР [13, 14] и рентгеноструктурного анализа [15]. Метод МД имеет очень большую разрешающую способность, которая позволяет в каждый момент времени знать координаты и скорости всех атомов в исследуемой модельной системе. Таким образом, компьютерное моделирование обеспечивает получение всесторонней информации практически о любых микроскопических характеристиках исследуемой системы.
В основе метода МД лежит расчёт классических траекторий движения атомов в фазовом пространстве координат и скоростей. При вычислении сил, действующих на атомы, учитываются межмолекулярные и внутримолекулярные взаимодействия [12].
Известно, что одним из основных ограничений классического метода МД является недостаточная вычислительная мощность современных компьютеров. Количество вещества, которое можно исследовать с помощью молекулярной динамики, много меньше даже самых маленьких образцов, обычно изучаемых с помощью физических экспериментов. Например, 1 мм3 жидкого кристалла содержит около 1018 молекул МББА (4,1 • 1019 атомов), в то время как типичные на сегодняшний день размеры исследуемых методом МД систем составляют 105 атомов. То же можно сказать и про интервалы времени, на протяжении которых изучается поведение модельных систем. Кроме того, при моделировании вязких жидкостей, в частности нематических жидких кристаллов, нельзя выполнять исследования на малых интервалах времени, так как есть вероятность попадания системы в локальный минимум энергии, в котором она будет находиться достаточно долго. В этом случае результаты короткого моделирования могут получиться статистически недостоверными. Для ЖК, которые являются вязкими системами, желательное время моделирования составляет десятки наносекунд.
Сегодня исследовать с помощью МД многие физические процессы и заметное количество макромолекул оказывается чрезвычайно сложно, несмотря на развитие компьютерной техники и разнообразные приёмы, применяемые для уменьшения количества необходимых вычислений. Например, в [16] проводится моделирование смеси жидких кристаллов на протяжении 100 нс. При этом исследуемая система (смесь четырёх жидких кристаллов) состоит из 256 молекул (около 10 тысяч атомов), чего явно недостаточно для получения достоверных физических характеристик при моделировании длинных молекул.
Многомасштабное моделирование молекулярной динамики. Существуют различные приёмы, позволяющие исследовать методом МД большое количество молекул на заметных временных интервалах. Во-первых, часто увеличивают шаг моделирования, который, вообще говоря, ограничен сверху характеристическим временем наиболее быстрых движений в системе. Это время в случае моделей с жёсткими химическими связями составляет 0,5-1 фс (10_ 15 с), но шаг может быть увеличен до 2-3 фс при использовании модели «объединённых атомов» [17], в которых атомы водорода объединены с атомами углерода в единые центры взаимодействия с массой 13-15 а. е. м. Иногда ради ускорения вычислений используют шаг моделирования 2 фс в полноатомной модели, в которой присутствуют атомы водорода с массой 1 а. е. м., (см., например [16]) однако это приводит к существенно менее точному предсказанию свойств системы.
Во-вторых, для ускорения моделирования часто уменьшают количество взаимодействующих центров в молекулах. Упрощённые, или «огрублённые», модели применяются в компьютерном моделировании давно. Они требуют меньше вычислительных ресурсов по сравнению с полноатомными моделями, поэтому можно сказать, что компьютерное моделирование началось с таких моделей и затем постепенно перешло к полноатомным. Случаем крайнего «огрубления» молекул ЖК можно считать модель Гея—Берне [18], в которой молекула представляется единым центром взаимодействия, имеющим вытянутую форму. Примером менее «грубой» модели, применявшейся для жидких кристаллов, можно считать модель, использованную в работе Кросса и Фанга [19] при исследовании 4-н-пентил-4'-цианобифенила. Другими примерами могут служить широко распространённая модель «объединённых атомов» [17], а также разнообразные модели «грубых зёрен» (соагее^гаш) [20-22], в которых группы близко расположенных атомов (например, несколько метильных групп или все атомы бензольного кольца) заменяются одним центром взаимодействия, который обычно называют «суператомом». Использование модели «объединённых атомов» позволяет уменьшить количество взаимодействующих центров примерно в два раза, что приводит к четырёхкратному выигрышу во времени вычислений. Модели «грубых зёрен» ещё больше сокращают это время.
Кроме непосредственного уменьшения количества взаимодействующих центров в системе, «грубые» модели позволяют ещё больше увеличить шаг моделирования, так как минимальная масса «суператомов» в них, как минимум, на порядок больше, чем в полноатомной модели, соответственно больше и характеристическое время движений «суператомов».
В качестве примера того, как «грубая» модель сокращает время расчёта, рассмотрим молекулу жидкого кристалла МББА. В полноатомной модели она состоит из 41 атома. В приближении «объединённых атомов» — из 20. Если же заменить два бензольных кольца на два «суператома» в приближении «грубых зёрен», то получится молекула, состоящая всего из 10 «атомов». Количество вычислений на каждом шаге моделирования пропорционально квадрату числа атомов, соответственно расчёт с использованием приближения «грубых зёрен» будет идти примерно в 17 раз быстрее, чем в полноатомной модели.
Обратной стороной использования «грубых» потенциалов является недостоверность физических характеристик исследуемых систем, получающихся в результате расчётов, так как модель упрощается, из неё удаляются атомы водорода, расположенные в основном на периферии молекул и участвующие в межмолекулярных взаимодействиях. В частности, при моделировании с использованием «грубых» потенциалов невозможно исследовать свойства, зависящие от взаимного расположения атомов одной или разных молекул. Например, реакционная способность молекулы связана с её конформационной
структурой, однако за счёт объединения атомов само понятие конформации может оказаться неопределённым. В отсутствие атомов водорода не будут наблюдаться эффекты, возникающие из-за водородных связей, и т. д. Полноатомные потенциалы не обладают этими недостатками, но требуют слишком больших вычислительных ресурсов. Таким образом, для исследования молекулярных свойств вещества необходимо применять детальные (полноатомные) модели, а для увеличения времени эволюции системы — «грубые» модели.
В литературе предложено несколько способов преобразования полноатомных моделей в «грубые» и обратно. Например, в работах Милано и Мюллера—Платэ [23, 24] при исследовании виниловых полимеров потенциалы взаимодействия «суператомов» разрабатываются и оптимизируются по данным моделирования фрагментов полноатомных молекул. Можно, напротив, использовать унифицированные и огрублённые модели, например потенциал MARTINI [25], предлагающий пять типов частиц для представления любых больших групп атомов. При исследовании растворов больших молекул применяют также адаптивные модели (для растворителя), в которых детализация модели зависит от расстояния до большой молекулы исследуемого вещества, так что в определённом диапазоне расстояний молекула растворителя представляется двумя моделями: более и менее грубой [26].
В данной работе мы используем метод MuScaT [27] (от multiscale transitions — «переходы между масштабами») собственной разработки. Он позволяет в рамках одного непрерывного моделирования осуществлять переходы между моделями разной детализации. Молекулы разных типов преобразуются каждая в свой «огрублённый» или, наоборот, детализированный аналог, но при этом процедура уменьшения/увеличения детализации (т. н. coarse/fine-graining procedure) одинакова для всех типов молекул.
Переходы осуществляются таким образом, что каждая молекула в системе, промоделированной с использованием одного типа потенциалов, заменяется на молекулу, описанную с помощью другого типа потенциалов. При этом координаты и скорости атомов, совпадающих в обоих типах потенциалов, сохраняются, а координаты недостающих атомов рассчитываются для каждой молекулы относительно совпадающих атомов. Скорости недостающих атомов задаются в каждом случае особым образом. Переход от одной модели к другой происходит с изменением числа атомов (силовых центров) в молекуле.
В случае перехода от более к менее детализированной модели «лишние» атомы убираются, а остальные замещаются на соответствующие силовые центры с новыми параметрами взаимодействия. Геометрия каждой молекулы при этом сохраняется такой, какая получилось в результате моделирования, то есть переход затрагивает имеющиеся в системе конформации только в том случае, если из молекулы удаляются атомы, непосредственно формирующие ту или иную конформацию. Скорости замещающих атомов совпадают со скоростями замещаемых. Кроме того, в «грубых» моделях могут присутствовать «суператомы», координаты которых не соответствуют никакому атому из исходной системы. Например, атомы бензола могут замещаться одним «суператомом», расположенным в геометрическом центре кольца. Скорость такого нового «суператома» будет равна средней скорости шести атомов углерода.
Более сложным случаем является обратный переход от «грубой» модели к детализированной, так как здесь новые атомы появляются не взамен старых, а добавляются к уже имеющимся. Пространственное расположение таких атомов определяется с помощью специальной z-матрицы, параметры в которой задаются относительно положений атомов с уже определёнными координатами таким образом, чтобы длины связей нового
атома и валентные и двугранные углы, составленные с его участием, принимали свои равновесные значения. Это позволяет минимизировать неизбежный скачок энергии на первых шагах моделирования после перехода. Скорость каждого добавленного атома обычно задаётся равной скорости того замещённого атома, с которым у нового атома есть связь. За счёт этого минимизируется вероятность разрыва связей на первых шагах после преобразования системы.
Сразу после такого перехода свойства системы должны быть как можно ближе к равновесным. Происходит неизбежный скачок термодинамических параметров, но время приведения системы к равновесию после перехода много меньше времени, необходимого для уравновешивания системы из начальной конфигурации. В частности, сразу после перехода распределение энергии по степеням свободы будет сильно отличаться от равновесного, так как начальные скорости новых атомов внесут вклад исключительно в трансляционную кинетическую энергию.
В данном исследовании жидких кристаллов мы использовали только два вида потенциалов: OPLS-AA [28] и OPLS [17], то есть полноатомную модель и модель «объединённых атомов», хотя метод MuScaT не накладывает никаких ограничений на то, между какими моделями происходит преобразование. Использование более «грубых» моделей для моделирования эволюции системы и последовательная их «детализация» для изучения свойств являются темой отдельной работы.
Детали моделирования. Для уточнения параметров моделей исследуемых молекул мы провели их квантово-химическое моделирование. Это было сделано, во-первых, для того, чтобы определить угол между плоскостями бензольных колец (двугранный угол a вдоль связи азот-углерод бензольного кольца), во-вторых, чтобы уточнить заряды атомов в полноатомной модели обеих молекул. Квантово-химический расчёт производился в программе GAMESS [29]. Оптимизация геометрии молекул выполнялась в базисе 6-31G**++. Равновесные значения углов a и b приведены в табл. 2 вместе с другими па-
раметрами двугранных углов [14]. Фрагмент -N=CH- лежит в плоскости анилинового бензольного кольца. Угол между плоскостями бензольных колец равен 45°. В оптимальной геометрии вычислены частичные заряды на атомах по методу Мил-ликена в базисе STO-3G, которые ближе к значениям, рекомендованным разработчиками потенциала OPLS-AA [28]. Заряды атомов, имеющиеся в потенциале OPLS-AA,
взяты из него, а заряды отсутствующих в OPLS-AA атомов — из результатов квантово-
химического расчёта. В табл. 3 приведены значения зарядов атомов обеих молекул, использованные нами при моделировании.
Исследуемые системы МББА и БОБТ состояли из 726 молекул каждая. Использовались периодические граничные условия. Интегрирование уравнений движения производилось по методу Верле (leap-frog) [12]. Электростатические взаимодействия
Таблица 2 Параметры двугранных углов, использованные при моделировании
БОБТ МББА
Угол 9, ° V, ккал/моль п 9, ° V, ккал/моль п
а 45 10 4 45 10 4
b 0 80 2 0 80 2
с 180 1 2 180 1 2
d 180 1Д 3 -
180 1,7 1 180 1,3 1
е 180 1,7 1 180 -0,05 2
180 1,7 1 180 0,2 3
180 1,3 1 180 1,3 1
/ 180 1,3 1 180 -0,05 2
180 1,3 1 180 0,2 3
Значения зарядов атомов, использованные при моделировании
№ атома БОБТ МББА № атома БОБТ МББА
атом заряд, е атом заряд, е атом заряд, е атом заряд, е
1 N -0,26 N -0,26 21 Н 0,06 Н 0,06
2 С 0,13 С 0,13 22 Н 0,115 Н 0,115
3 С -0,02 С -0,02 23 н 0,115 н 0,115
4 С 0,0 С 0,0 24 н 0,115 н 0,115
5 с 0,0 с 0,0 25 н 0,115 н 0,115
6 с 0,0 с 0,0 26 н 0,06 н 0,06
7 с 0,0 с 0,0 27 н 0,06 н 0,06
8 с 0,15 с 0,15 28 н 0,06 н 0,06
9 о -0,26 о -0,26 29 н 0,06 н 0,115
10 с 0,12 с 0,0 30 н 0,06 н 0,115
11 с 0,0 с ОД 31 н 0,06 н 0,115
12 с 0,0 с 0,0 32 н 0,06 н 0,115
13 с 0,0 с 0,0 33 н 0,06 н 0,06
14 с 0,10 с 0,0 34 н 0,06 н 0,06
15 с 0,0 с 0,0 35 н 0,115 н 0,06
16 с 0,0 с 0,04 36 н 0,115 н 0,06
17 с 0,0 с 0,12 37 н 0,115 н 0,06
18 с 0,0 с 0,0 38 н 0,115 н 0,06
19 с 0,04 с 0,0 39 н 0,06 н 0,06
20 с 0,0 с 0,0 40 н 0,06 н 0,06
- 41 н 0,06 н 0,06
учитывались по методу Эвальда [30], радиус отсечения взаимодействий составлял 14,5 А, параметр а метода Эвальда равнялся 0,27.
В модели «объединённых атомов» использовался список ближних соседей Вер-ле [31], который обновлялся каждые 100 шагов. Шаг моделирования в модели «объединённых атомов» составлял 2 фс, в полноатомной модели — 1 фс.
Система создавалась в модели «объединённых атомов» и уравновешивалась. В начальный момент времени молекулы располагались в узлах пространственной решётки 11 х 11 х 6 и были преимущественно ориентированы вдоль оси Z с отклонением не больше 35°. Плотность системы составляла 0,1 г/см3. Моделирование начиналось в ансамбле NVT (использовался термостат Берендсена [12]), в котором система сжималась до плотности, свойственной ей при нормальных условиях (МББА — 1,039 г/см3, БОБТ — 1,007 г/см3). Окончательное уравновешивание производилось в ансамбле NpT (по алгоритму Хувера [32]).
Затем производилось длительное моделирование (5 нс или 2,5 млн шагов) в ансамбле NVT с использованием термостата Нозе-Хувера [33, 34]. Молекулы жидкого кристалла перемешивались, теряя первоначальную «неестественную» трансляционную упорядоченность. Через 5 нс осуществлялся переход в полноатомную модель, в которой система уравновешивалась в течение 0,2 нс, после чего производился анализ траекторий атомов в течение 1,4 нс. Общее время моделирования в полноатомной модели составило
1,6 нс. После этого был произведён переход обратно в модель «объединённых атомов»,
в которой система эволюционировала в течение ещё 5 нс. Затем снова происходил переход в полноатомное приближение, в котором снова после уравновешивания анализировались траектории движения атомов. Общее время моделирования обеих систем составило 13,2 нс.
Расчёты осуществлялись на вычислительных кластерах физического факультета СПбГУ. Реальное время моделирования двух обсуждаемых систем — около двух месяцев. Расчёт траекторий длительностью 1 нс в полноатомной модели происходил приблизительно в 10 раз дольше, чем в модели «объединённых атомов».
Отдельно следует упомянуть температуры, при которых производилось моделирование. Так как диапазон существования жидкокристаллической фазы БОБТ крайне мал (2 градуса), расчёт производился при температуре из середины этого диапазона (340 К). Поэтому для моделировая МББА была взята температура 307 К, которая находится в середине диапазона существования жидкокристаллической фазы этого вещества. Кроме того, МББА был промоделирован дополнительно при температурах 300 и 314 К для оценки зависимости от температуры некоторых параметров.
Результаты и обсуждение.
Усреднённая форма конформационно подвижных молекул. Для получения информации о конформациях молекул, их форме и взаимном расположении мы использовали распределения плотности вероятности по двугранным углам, статистику нахождения элементов алифатических цепей в той или иной конформации и цилиндрические функции распределения. Цилиндрические функции распределения (ЦФР) [35] являются представлением двухчастичных функций распределения, как и широко применяемые функции радиального распределения (ФРР). ЦФР строятся в молекулярной системе координат и показывают расположение атомов молекул друг относительно друга. В начале координат находится центр масс молекулы. Вертикальная ось на всех рисунках совпадает с направлением оси Z тензора инерции молекулы, то есть с длинной молекулярной осью. По горизонтальной оси отложена абсолютная величина расстояния в плоскости ХУ молекулярной системы координат, и рисунки симметричны относительно оси Z. ЦФР показывают плотность вероятности расположения атома или группы атомов в определённом месте относительно центра масс молекулы. В отличие
от сферически симметричных ФРР, ЦФР дают информацию о вероятности найти атом в элементарном цилиндрическом слое данного радиуса и на данной высоте.
На рис. 2 изображено пространственное межмолекулярное распределение плотности центров масс молекул системы друг относительно друга. На это распределение наложено внутримолекулярное распределение атомов относительно центра масс молекулы. Ориентация бензилиденанилинового фрагмента обеих молекул идентична, атом кислорода находится в верхней части рисунка (ср. с рис. 1). Горизонтальные чёрные линии отделяют на рисунке бензилиденанилиновый фрагмент от алифатических цепей. Рисунок выполнен в градациях серого, где белый цвет означает отсут-Рис. 2. ЦФР центров масс молекул ствие атомов или центров масс молекул, серый
Рис. 3. Распределение плотности вероятности угла между длинной осью молекулы и директором ЖК
цвет — среднюю плотность центров масс или атомов, а чёрный — плотность, в два или более раз превышающую среднюю.
Белый цвет, на фоне которого построено внутримолекулярное распределение атомов, обозначает т. н. исключённый объём молекулы — объём, в котором не могут находиться атомы других молекул. Видно, что формы исключённого объёма БОБТ и МББА несколько различаются: у БОБТ он более вытянутый (19 А против 17,5 А у МББА), у МББА — более широкий (6 А против 5 А у БОБТ).
Рассмотрим сначала, как молекулы располагаются друг относительно друга в исследуемых веществах. На рисунке у обеих молекул видно, как структурируется ближайшее окружение. Заметна высокая плотность (центров масс) около фрагмента —М=СН—. Кроме того, видно, что плотность низка около длинных алифатических цепей молекул, обладающих высокой подвижностъю и как бы «расталкивающих» другие молекулы. Такое поведение молекул было найдено и в 4-н-пентил-4/-цианобифениле [36]. Около коротких цепей, напротив, наблюдается увеличение плотности. Высокая плотность центров масс молекул БОБТ на оси Z в «торце» молекулы (на рис. 2 — непосредственно под короткой алифатической цепью) означает, что там сосредоточены молекулы, не сонаправленные с молекулярной осью Z (которая преимущественно направлена вдоль директора ЖК).
На рис. 3 изображено распределение молекул по углу между длинной осью молекулы и директором. Это распределение позволяет судить о преимущественной ориентации молекул в системе. Видно, что молекулы ориентированы либо по, либо против
директора, что и должно наблюдаться в нематической фазе жидкого кристалла. Молекул, ориентированных перпендикулярно директору, в системе не наблюдается.
Учитывая сказанное, можно сделать вывод, что молекулы БОБТ, расположенные непосредственно под короткой алифатической цепью, отклонены от директора на сравнительно небольшие углы. Это справедливо и для МББА, однако распределение молекул в ближайшем окружении в данном случае более равномерное практически по всему контуру молекулы.
Рассмотрим короткие алифатические цепи молекул: группу -СН3 в БОБТ и -О-СНз в МББА. Атомы водорода метильной группы БОБТ и метокси-фрагмента МББА по-разному локализуются в пространстве (см. рис. 2). Атомы водорода в БОБТ равномерно распределяются по кругу, то есть группа СН3 вращается вокруг связи угле-род-ароматический углерод. Распределение плотности вероятности двугранных углов по этой связи показывает, что группа СН3 вращается свободно.
Для молекулы МББА наблюдается иная картина. Распределение плотности вероятности двугранного угла вдоль связи кислород-ароматический углерод показывает, что наиболее вероятны углы ± 90о, то есть группа СН3 может находиться с обеих сторон по отношению к плоскости бензольного кольца. Это не совпадает с экспериментальными результатами, приведёнными в [37], где из анализа спектров ЯМР 13С получено значение этого угла 9о. Интерпретация в [37], однако, производилась в предположении, что данный угол фиксирован, распределение плотности вероятности не учитывалось.
Кроме того, видны две отдельные области наиболее вероятной локализации атомов водорода группы СН3, которые мы наблюдаем из-за того, что группа СН3 вращается вокруг оси, находящейся под углом к плоскости бензольного кольца (и плоскости рисунка). Иными словами, на рисунке мы наблюдаем отдельно атом водорода, двугранный угол для которого по связи «кислород-углерод СН3» составляет 180о, и отдельно — группу, состоящую из двух оставшихся атомов водорода.
Перейдём теперь к длинным алифатическим цепям: -С4Н9 в МББА и -О-С4Н9 в БОБТ. На рис. 4 изображены ЦФР отдельных групп СН2 и СН3 бутильной цепи. В качестве фона серым тоном изображены ЦФР атомов водорода всей длинной алифатической цепи. Чёрным тоном (высокая плотность вероятности) выделен вклад в общую ЦФР от атомов водорода той или иной группы атомов. Номера атомов в молекулах подписаны для каждого случая (ср. с рис. 1). Сравниваемые группы БОБТ и МББА размещены на одном уровне. Штрихами обозначены положения центров масс молекулы. Алифатические цепи изображены повёрнутыми одинаково (т. е. молекула БОБТ развёрнута относительно положения, показанного на рис. 2). В сочетании с распределениями плотности вероятности
БОБТ МББА
32, 33, 34 39, 40, 41
Рис. 4. Пространственное рас-
пределение отдельных групп атомов водорода алифатических цепей
1,6 1,4 1,2 1
0,8 0,6
0,4
0,2 0
-200 -150 -100 -50 0 50 100 150 200
Рис. 5. Плотность вероятности распределения двугранных углов е в алифатических цепях
двуграных углов в алифатической цепи такие ЦФР позволяют разделить вклады, вносимые различными группами атомов в конформационные переходы.
На рис. 5 изображены плотности вероятности распределения двугранных углов, расположенных на одном расстоянии от бензольного кольца, в алифатических цепях (/ для МББА и е для БОБТ). Минимумы функции плотности вероятности между значениями угла 180 и ± 60° позволяют произвести анализ поведения цепей в рамках модели вращательных изомеров [38]. Видно, что у молекул БОБТ шире диапазон углов, соответствующих состояниям «гош+» и «гош—», в то время как у молекул МББА шире диапазон углов, соответствующих состоянию «транс». Площади под графиками для обеих молекул совпадают. Вероятности реализации состояний «гош+» и «гош—» первого двугранного угла у БОБТ примерно одинаковые и составляют в сумме около 45 %. У МББА сумма вероятностей этих же состояний равна приблизительно 7,5 %.
Для получения общей картины распределения состояний алифатических цепей в исследуемых молекулах подсчитали вероятности комбинаций различных конформаций в этих цепях. Были взяты два двугранных угла е и / (см. обозначения углов на рис. 1).
Вероятности реализации той или иной комбинации приведены в табл. 4. Первое обозначение в столбце «Конформации» относится к углу е, второе — к углу / Сумма вероятностей реализации состояний «транс—транс», «гош—транс» и «гош+транс» у БОБТ равна 88 %, у МББА — 86 %. Полученные
Таблица 4 Вероятность реализации конформаций в алифатических цепях, %
Конформации БОБТ МББА
g- g- 1,8 ОД
g+, g+ 1,8 0,3
g- * 22,7 3,6
g+, * 24,3 3,5
g- g+ 1,4 ОД
g+^ g- 1,5 ОД
g- 2,7 6,8
g+ 2,5 6,3
1 41,4 79,2
данные не противоречат результатам ранних исследований бутильной цепи гомолога МББА — 4-этоксибензилиден-4'-н-бутиланилина [39].
Вероятность «вытянутого» состояния «транс-транс» для молекул МББА примерно в два раза выше, чем для молекул БОБТ. Это можно объяснить тем, что в молекулах БОБТ алкильная цепочка соединяется с бензольным кольцом через кислород, а не напрямую.
Рассмотрим подробнее пространственные распределения отдельных групп длинных алифатических цепей. Атомы водорода 33 и 34 молекулы МББА (номера атомов указаны на рис. 1) рассматривать не имеет смысла, так как в молекуле БОБТ нет атомов, находящихся в аналогичном положении. Атомы водорода 35 и 36 молекулы МББА можно сравнить с атомами водорода 26 и 27 молекулы БОБТ, так как эти группы СН2 алифатических цепей находятся в схожих положениях и отстоят от бензольного кольца на один атом. Если взглянуть на отдельные пространственные распределения этих атомов водорода, то видно, что в молекуле БОБТ вращение вокруг связи «кислород-углерод» (угол 1) не происходит — практически всегда атомы водорода первой группы СН2 алифатической цепи локализованы немного ближе к центру молекулы, чем атом углерода этой группы, то есть двугранный угол вдоль этой связи находится в основном в положении «транс». В молекуле МББА можно говорить о небольшом вращении вокруг связи «углерод-углерод» (угол е), что подтверждается видом функций распределения этих двугранных углов. В молекуле МББА более широкое распределение угла вокруг положения «транс» и максимумы функции распределения, соответствующие конформациям «гош±», смещены к значениям ± 75°. Можно говорить о немного большей подвижности двугранного угла вокруг этой связи в МББА.
Если сравнить атомы 37 и 38 в МББА с атомами 28 и 29 в БОБТ, то видно, что форма области пространства, занятого этими атомами, схожа у обеих молекул. Однако в случае МББА атомы занимают весь объём пространства равномерно, тогда как в распределении для БОБТ заметны области более и менее вероятной локализации.
Наконец, сравнивая распределение атомов водорода метильных групп в МББА и БОБТ, можно видеть, что форма области пространства, занятого этими атомами у МББА, шире, чем у БОБТ. Несмотря на то, что в модели вращательных изомеров молекула МББА имеет явно предпочтительную вытянутую транс-конформацию (см. табл. 4), распределения двугранных углов совместно с ЦФР показывают, что бутильная цепь молекулы МББА имеет большую свободу отклонения от длинной молекулярной оси, что приводит к увеличению диаметра исключённого объёма молекулы.
Внутримолекулярные функции радиального распределения. Для определения эффективной длины молекул МББА и БОБТ были проанализированы внутримолекулярные ФРР между различными атомами и группами атомов.
На рис. 6 изображены ФРР между атомами водорода крайних метильных групп, что соответствует полной длине молекулы. Длина обеих молекул варьируется в одинаковом диапазоне от 11,5 до 18,1 А, однако профили ФРР различны. У МББА наблюдается два относительно нерезких максимума на 15,2 и 16,2 А, в то время как у БОБТ наблюдается один резкий максимум на 16,95 А.
0,1 0,06 0,02
0 5 10 15 20 25 30
Рис. 6. Внутримолекулярные ФРР полной длины молекул
1 Л 1 1
/Д :
: МББА
\ БОБТ -
■ ■ ■ ■ '
Рис. 7. Внутримолекулярные ФРР длинных алифатических цепей МББА и БОБТ
Бензилиденанилиновый фрагмент представляют собой жёсткую конструкцию длиной 11,9 ± 0,1 А в обеих молекулах.
ФРР атомов водорода метильной группы относительно атома кислорода БОБТ имеет два резких максимума (рис. 7): первый на расстоянии 6,2, второй — 6,9 А. ФРР атомов водорода метильной группы относительно первого атома углерода бутильной цепи в МББА похожа: два максимума находятся на расстоянии 5,15 и 5,9 А (нужно принять во внимание расстояние до атома кислорода в БОБТ). Более интенсивные и размытые первые максимумы соответствуют гош-конформациям цепи, а вторые максимумы — транс-конформации. Видно, что у МББА транс-конформация более вероятна, чем у БОБТ (см. табл. 4), в то время как длинная алифатическая цепь БОБТ находится в основном в гош-конформациях. Однако полная длина молекулы МББА получается меньше, чем длина БОБТ, из-за подвижности её метильной группы.
Таким образом, эффективная длина молекулы БОБТ больше, чем МББА, следовательно, если рассматривать эти молекулы как жёсткие стержни, то температуры фазовых переходов БОБТ должны быть выше, чем соответствующие температуры МББА. Подобная закономерность наблюдается для жёстких цилиндрических молекул ряда замещённых п-три-, п-тетра-, п-пента- и п-гексафенилов [40]. Это может объяснять факт наличия различных температурных диапазонов существования ЖК-фазы исследуемых молекул.
Угол между бензольными кольцами. Распределение количества молекул по углу между бензольными кольцами позволяет проанализировать конформационные переходы в центральном жёстком фрагменте молекул. На рис. 8 изображены четыре распределения по двугранным углам а и Ь для МББА и БОБТ, построенные по траекториям на конечном этапе моделирования.
В начальный момент времени все молекулы находились в одинаковой конформации. Угол а = +135°, а Ь = 0°. На рисунке видно, что для угла а существуют четыре вероятных положения: ± 45° и ± 135°. За время моделирования небольшая часть бензольных колец повернулась из положения +135° в положения +45° и —135° (поворот на 90°), и ещё меньшая часть — в положение —45° (два поворота на 90°), несмотря на то, что все эти положения равновероятны. За 13 нс моделирования произошло всего несколько переходов, поэтому мы не можем оценить характеристическое время таких «перескоков», а можем лишь констатировать, что бензольные кольца вращаются вокруг своих связей а прыжками на 90 градусов.
Вращение по углу Ь может происходить прыжками по 180°. На момент окончания моделирования зафиксировано небольшое число молекул БОБТ, у которых произошёл такой конформационный переход. В молекулах МББА такого перехода не обнаружено.
Ориентационный порядок жидкого кристалла. Для обоих жидких кристаллов был рассчитан параметр порядка Бгг = Р2(сов 0), полином Лежандра 2-го порядка:
0,35
0,25
0,15
0,05
1 MBBA 1-2-BBA 2-1-11 BOBT 1-2- 1 3-4 ! .
M -12 3-4 ! ) 1
B OBT 2-1-14- 15 I 1 ) 1 1 ) II
1 1 1 1 1 f 1 1 ) 1
$ 1 ? ) 1
! 1 1 I ;Т i 1 J II
V 1 1 1 1 i 1
1 1 1 1 1 !
I 1 1 1 1 I !
i -н ь- 1 1 1 : 1 / А \ - i ! 1 4j
-200
-150
-100
50
0
50
100
150
200
Рис. 8. Распределение угла между бензольными кольцами
3 1
Szz = P2(cos0) = -cos20 — (1)
а также Pj(cos 0), полином Лежандра 4-го порядка:
Р4(cos 0) = — (35 cos4 0 — 30 cos2 0 + 3^. (2)
Кроме того, был рассчитан параметр двуосности молекул:
с, с, _3sin20cos2c|> ^
&ХХ дуу — 2 ,
а также построены функции ориентационного распределения, которые можно вычислить на основании траекторий для различных векторов и анализировать, используя формулу
f (cos0) = c • exp(a2P2(cos 0) + a4P4(cos0)+a6p5(cos 0)), (4)
где c — нормировочная постоянная, а коэффициенты щ — соответствующие параметры порядка жидкого кристалла (Szz, P4 и т. д.)
Рис. 9. Ориентационные функции распределения для МББА и БОБТ
Таблица 5
Параметры порядка и двуосности ЖК
Вещество Т, К Szz Pi
МББА 300 К 0,647 0,065 0,267
МББА 307 К 0,665 0,056 0,289
МББА 314 К 0,667 0,056 0.288
БОБТ 340 К 0,697 0,029 0,321
Параметры порядка Szz и P4, а также параметр двуосности Sxx -Syy для исследуемых жидких кристаллов приведён в табл. 5. Из неё видно, что все вещества находятся в ЖК-фазе. Однако температурная зависимость параметра порядка МББА (увеличение Szz с температурой) показывает, что эта величина не успевает прийти к равновесию за 1,5 нс, в течение
которых уравновешивались системы в данной работе. Вычисленные параметры для МББА соответствуют экспериментально исследованным в [41, 42].
На рис. 9 приведены ориентационные функции распределения f (cos 0), где 0 — угол между осью Z молекулы и направлением директора жидкого кристалла. Все они дают острый максимум около 0 = 0, что совпадает с тем, как ведут себя эти функции в нематической фазе жидкого кристалла [14].
Трансляционная и вращательная диффузия молекул. Коэффициенты трансляционной диффузии могут быть получены из формулы
N
Da
Иш(1/2*)(^ R(t) - Да(о)]'
j = 1
(5)
где а = х, у или г определяет координатную ось, а йа(і) — положение центра масс молекулы в момент времени і. Кроме того, коэффициенты трансляционной диффузии могут быть рассчитаны из автокорреляционной функции линейной скорости центра масс молекул.
В жидких кристаллах коэффициенты трансляционной диффузии анизотропны. В одноосной фазе можно выделить две компоненты тензора диффузии: и Вхх = Вуу = П±, ориентированные соответственно параллельно и перпендикулярно директору ЖК. Движение молекул вдоль направления директора в жидких кристаллах происходит быстрее, чем перпендикулярно ему, поэтому следует ожидать, что будет больше .
Рассчитанные анизотропные коэффициенты самодиффузии для БОБТ и МББА, уравновешенных при температурах, соответствующих середине их жидкокристаллического температурного диапазона, приведены в табл. 6 вместе с изотропным значением. Кроме того, для МББА приведены также коэффициенты самодиффузии при температурах 300, 307 и 314 К, что позволяет оценить температурную зависимость.
Таблица 6
Коэффициенты трансляционной диффузии, м2/с
Вещество Т, К О х 10~12 £>ц х 10~12 х 10~12
МББА 300 3,2 3,3 2,7
МББА 307 4,5 6,7 3,5
МББА 314 10 15 7,8
БОБТ 340 13 21 8,3
Значения для МББА хорошо совпадают с экспериментальным изотропным значением О = 20 • 10-12 м2/с, полученным при 303 К [43] методом нейтронного рассеяния.
Для обоих веществ коэффициент самодиффузии вдоль направления директора больше, чем в перпендикулярных направлениях: для МББА в 2 раза при 307 К и в 1,25 раза при 300 К, для БОБТ — в 2,5 раза.
Рассчитанные коэффициенты показывают, что подвижность БОБТ по сравнению с МББА выше, особенно вдоль направления директора ЖК. Если сравнивать коэффициенты БОБТ с коэффициентами МББА из середины температурного диапазона (307 К), то видно, что они различаются в 2,5-3,5 раза, если сравнивать с коэффициентами МББА при температуре 314 К, различия меньше — на 20-25
Другой способ определения коэффициентов самодиффузии связан с вычислением автокорреляционной функции линейной скорости цетра масс молекул:
СЮ
Аш = \ J (г;«(°) ‘ Та, (6)
0
где угловые скобки означают среднее по ансамблю; V — линейная скорость молекулы; М — молекулярная масса; кв — постоянная Больцмана, а та — время корреляции, равное интегралу от корреляционной функции.
На рис. 10 приведена нормированная автокорреляционная функция линейной скорости центра масс молекул. Коэффициенты трансляционной диффузии, рассчитанные из автокорреляционной функции линейной скорости центра масс, приведены в табл. 7. Из-за неточности интегрирования автокорреляционной функции этим способом нами были вычислены только изотропные коэффициенты самодиффузии.
Рис. 10. Автокорреляционная функция линейной скорости центра масс
Таблица 7
Коэффициенты диффузии, вычисленные из автокорреляционных функций
Вещество Г, К Агапв X 10~12, М2/С Ас* X 10~12, рад2/с
МББА 300 25 67
МББА 307 40 100
МББА 314 - -
БОБТ 340 70 120
Аналогично автокорреляционная функция скорости вращения длинных осей молекул, приведённая на рис. 11, позволяет оценить коэффициенты вращательной диффузии. Вычисленные таким образом коэффициенты и приведённые в таблице хорошо совпадают с экспериментальными: БГ01 = 230 х 10~12 рад2/с [43].
На рис. 12 представлена автокорреляционная функция ориентаций длинных осей молекул, позволяющая оценить, насколько молекулы жидких кристаллов в равновесном состоянии сохраняют первоначальную ориентацию. Видно, что молекулы БОБТ сохраняют её немного лучше, чем молекулы МББА, и что в обоих случаях можно говорить о дальнем ориентационном порядке в системе.
Выводы. Моделирование молекулярной динамики методом МиБеаТ позволило наблюдать медленные процессы в жидких кристаллах, такие как скачкообразные переориентации бензольных колец и медленная диффузия молекул. Моделирование двух изомеров позволило установить причины их различных свойств в жидкокристаллической фазе.
Рис. 11. Автокорреляционная функция скорости вращения длинных осей молекул
1. Оба исследуемых вещества — МББА и БОБТ — являются жидкими кристаллами, и их ЖК-фаза воспроизводима в компьютерном эксперименте, несмотря на малый температурный диапазон существования этой фазы у БОБТ. Это видно из параметров порядка (см. табл. 5), полученных в результате расчётов, а также из автокорреляционных функций ориентации длинных осей молекул. Такая воспроизводимость ЖК-фазы в компьютерном эксперименте может быть связана с несоответствием «компьютерной» температуры температуре реальной. Эта известная проблема компьютерного моделирования обсуждается, в частности, в [44].
2. Заметно различаются конформации бутильной и бутокси-цепей исследованных молекул. В бутильной цепи МББА в два раза вероятнее полностью вытянутое транссостояние, чем в бутокси-цепи БОБТ.
3. Из межмолекулярных цилиндрических функций распределения видно, что свободный объём молекул БОБТ длиннее и уже, чем свободный объём молекул МББА, что может объяснять более высокие температуры фазовых переходов стержнеобразных молекул ЖК.
4. Ближние окружения молекул ЖК различаются. Молекулы БОБТ имеют три области более вероятной локализации центров масс соседей: около атомов кислорода и азота, а также на оси Z около группы -СНз, в то время как соседние молекулы МББА распределяются более равномерно вдоль «стенки цилиндра» и практически не локализуются около группы —О—СН3.
5. В обеих молекулах наблюдаются скачкообразные вращения бензольных колец на 90 и 180°. Кроме начального состояния, в котором угол между кольцами равен 135°, в конце моделирования наблюдались ещё состояния с углами 45, —45 и —135°.
6. Конформационная структура алифатических цепей оказывает влияние на свойства жидких кристаллов.
Литература
1. Чандрасекар С. Жидкие кристаллы. М., 1980.
2. Де Жен П. Физика жидких кристаллов. М., 1977.
3. Haller I. Elastic constants of the nematic liquid crystalline phase of p-methoxybenzylidene-p-n-butylaniline (MBBA) // J. Chem. Phys. 1972. Vol. 57. N 4. P. 1400-1405.
4. Marinov Y., Simova P. Hydrodynamic flow and surface tension temperature dependence during the nematic-isotropic phase transition // J. Phys. (D). 1992. Vol. 25. P. 1495-1499.
5. Evans O. R., Lin W. Crystal engineering of nlo materials based on metal-organic coordination networks // Accounts of Chemical Research. 2002. Vol. 35. N 7. P. 511-522.
6. Magnetism: A Supramolecular Function // NATO Science Series C. Proc. NATO Advanced Research Workshop / Ed. by O. Kahn. Carcans-Maubuisson. 1995.
7. Jezierska A., Panek J. J. First-principle molecular dynamics study of selected schiff and mannich bases: application of two-dimensional potential of mean force to systems with strong intramolecular hydrogen bonds // J. Chem. Theory Comput. 2008. Vol. 4 N 3. P. 375-384.
8. Ohshima A., Momotake A, Arai T. Photochromism, thermochromism, and solvatochromism of naphthalene-based analogues of salicylideneaniline in solution // J. Photochem. Photobiol. 2004. Vol. 162. N 2-3. P. 473-479.
9. Filarowski A., Koll A., Gowiak T. Proton transfer equilibrium in the intramolecular hydrogen bridge in sterically hindered schiff bases // J. Mol. Struct. 2002. Vol. 615. N 1-3. P. 97-108.
10. Hadjoudis E., Dziembowska T., Rozwadowski Z. Photoactivation of the thermochromic solid di-anil of 2-hydroxy-5-methyl-isophthalaldehyde in beta-cyclodextrin // J. Photochem. Photobiol. (A). 1999. Vol. 128. N 1-3. P. 97-99.
11. Ogawa K., Harada J. Aggregation controlled proton tautomerization in salicylideneani-
lines // J. Mol. Struct. 2003. Vol. 647. N 1-3. P. 211-216.
12. Allen M. P., Tildesley D. J. Computer simulation of liquids. Oxford, 1987.
13. Sandstrom D., Komolkin A. V., Maliniak A. Orientational order in a liquid crystalline mixture studied by molecular dynamics simulation and NMR // J. Chem. Phys. 1996. Vol. 104. P. 9620.
14. Komolkin A. V., Laaksonen A., Maliniak A. Molecular dynamics simulation of a nematic liquid crystal // J. Chem. Phys. 1994. Vol. 101. N 5. P. 4103-4116.
15. Socci E. P., Farmer B. L., Bunning T. J., Pachter R., Adams W. W. Molecular dynamics
and x-ray scattering simulations of cyclic siloxane-based liquid crystal mesogens // Liq. Crystals. 1993. Vol. 13. N 6. P. 811-827.
16. Pelaez J., Wilson M. Molecular orientational and dipolar correlation in the liquid crystal mixture e7: a molecular dynamics simulation study at a fully atomistic level // Phys. Chem. Chem. Phys. 2007. Vol. 9. N 23. P. 2968-2975.
17. Jorgensen W. L., Tirado-Rives J. The opls potential functions for proteins. energy minimizations for crystals of cyclic peptides and crambin // J. Am. Chem. Soc. Protein & Peptide. 1988. Vol. 110. P. 1657-1666.
18. Gay J. G., Berne B. J. Modification of the overlap potential to mimic a linear site-site potential // J. Chem. Phys. 1981. Vol. 74. P. 3316.
19. Cross C. W., Fung B. M. Molecular dynamics simulations for cyanobiphenyl liquid crystals // J. Mol. Cryst. Liq. Cryst. 1995. Vol. 262. P. 507-524.
20. Izvekov S., Voth G. A. A multiscale coarse-graining method for biomolecular systems // J. Phys. Chem. (B). 2005. Vol. 109. N 7. P. 2469-2473.
21. Zacharopoulos N., Vergadou N., Theodorou D. N. Coarse graining using pretabulated potentials: Liquid benzene // J. Chem. Phys. 2005. Vol. 122. P. 244111.
22. Wilson M. R. Progress in computer simulations of liquid crystals // Int. Rev. Phys. Chem. 2005. Vol. 24. N 3-4. P. 421-455.
23. Milano G., Mueller-Plathe F. Mapping atomistic simulations to mesoscopic models: a systematic coarse-graining procedure for vinyl polymer chains // J. Phys. Chem. (B). 2005. Vol. 109. P. 18609-18619.
24. Santangelo G., Di Matteo A., Mueller-Plathe F., Milano G. From mesoscale back to atomistic models: A fast reverse-mapping procedure for vinyl polymer chains // J. Phys. Chem. (B). 2007. Vol. 111. P. 2765-2773.
25. Marrink S. J., Fuhrmans M., Risselada H. J., Periole X. The MARTINI force field, chapter 2. CRC press, 2008.
26. Praprotnik M., Delle Site L., Kremer K. Multiscale simulation of soft matter: From scale bridging to adaptive resolution // Annual Rev. Phys. Chem. 2008. Vol. 59. P. 545-571.
27. Комолкин А. В. Органические и гибридные наноматериалы. Иваново, 2009.
28. Jorgensen W. L., Maxwell D. S., Tirado-Rives J. Development and testing of the opls allatom force field on conformational energetics and properties of organic liquids // J. Am. Chem. Soc. 1996. Vol. 117. P. 11225-11236.
29. Gordon M. S., Schmidt M. W. Theory and Applications of Computational Chemistry: the first forty years. Amsterdam, 2005.
30. Ewald P. Investigations of crystals by means of roentgen rays // Ann. Phys. (Leipzig). 1921. Vol. 64. P. 253.
31. Verlet L. Computer “experiments” on classical fluids. I. Thermodynamical properties of lennard-jones molecules // Phys. Rev. 1967. Vol. 159. P. 98-103.
32. Melchionna S., Ciccotti G., Holian B. L. Hoover’s style molecular dynamics for systems varying in shape and size // J. Mol. Phys. 1993. Vol. 78. P. 533.
33. Melchionna S. Constrained systems and statistical distribution // Phys. Rev. (E). 2000. Vol. 61. P. 6165-6170.
34. Hoover W. G. Canonical dynamics: Equilibrium phase-space distributions // Phys. Rev. (A). 1985. Vol. 31. P. 1695-1697.
35. De Vries A. Some comments on cylindrical distribution functions, with special reference to X-ray studies of liquid crystals // Acta Crystallographica Section A. 1972. Vol. 28. N 6. P. 659-660.
36. Komolkin A. V., Maliniak A. Local structure of anisotropic systems determined by molecular dynamics simulation. application to a nematic liquid crystal // Mol. Phys. 1995. Vol. 84. N 6. P. 1227-1238.
37. Pines A., Chang J. J. Study of the isotropic-nematic-solid transitions in a liquid crystal by carbon-13-proton double resonance // Phys. Rev. (A). 1974. Vol. 10. N 3. P. 946-949.
38. Helfer C. A., Mattice W. L. The Rotational Isomeric State Model // Physical Properties of Polymers Handbook. N.-Y., 2007. P. 43-57.
39. Komolkin A. V., Molchanov Y. V., Yakutseni P. P. Computer simualtion of a real liquid
crystal // Liq. Crystals. 1989. Vol. 6. N 1. P. 39-45.
40. Demus D., Demus H., Zaschke H. Fluessige Kristalle in Tabellen. Leipzig., 1974.
41. Chang R. Orientational order in MBBA from optical anisotropy measurements // Mol. Cryst. Liq. Cryst. 1975. Vol. 30. P. 155-165.
42. Jen S., Clark N. A., Perchan P. S., Priestley E. B. Polarized Raman scattering studies
of orientational order in uniaxial liquid crystalline phases // J. Chem. Phys. 1977. Vol. 66. N 10.
P. 4635-4661.
43. Cvikl B., Dahlborg U., Cepic M. et al. Investigation of restricted molecular reorientation in MBBA in the nematic phase by quasielastic incoherent neutron scattering // Physica Scripta. 1991. Vol. 44. P. 63-68.
44. Неверов В. С., Комолкин А. В. Исследование методом молекулярной динамики структурных и термодинамических свойств воды // Хим. физика. 2010. Т. 29. Вып. 3. С. 33-42.
Статья поступила в редакцию 28 сентября 2010 г.
