CC BY
53
УДК 539.193:544.258
Вестник СПбГУ. Сер. 4. 2012. Вып. 1
В. С. Неверов, А. В. Комолкин
КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ЖИДКИХ КРИСТАЛЛОВ, СОДЕРЖАЩИХ БЕНЗОЛ С ЗАМЕСТИТЕЛЯМИ В пара-ПОЛОЖЕНИИ*
Введение. Изучение структуры и свойств жидких кристаллов можно производить экспериментальными методами (ЯМР [1, 2], рентгеноструктурный анализ [3], нейтронное рассеяние и др.) и посредством моделирования. Классический метод компьютерного моделирования молекулярной динамики (МД) [4] является одним из наиболее распространённых инструментов теоретического изучения структуры и динамики вещества. Заметную роль он играет и при установлении детальной структуры молекул и комплексов, определённой экспериментальными методами.
Одним из основных ограничений классического метода МД является недостаточная вычислительная мощность современных компьютеров. Количество вещества, которое можно исследовать с помощью молекулярной динамики, много меньше даже самых маленьких образцов, обычно изучаемых с помощью экспериментальных физических методов. Например, 1 мм3 жидкого кристалла содержит около 1018 молекул МББА (4,1-1019 атомов), в то время как типичные на сегодняшний день размеры исследуемых методом МД систем составляют 105 атомов. То же можно сказать и про интервалы времени, на протяжении которых изучается поведение модельных систем. Таким образом, исследовать с помощью МД многие физические процессы и заметное количество макромолекул оказывается чрезвычайно сложно, несмотря ни на эволюцию компьютерной техники, ни на разнообразные приёмы, применяемые для уменьшения количества необходимых вычислений. Например, в [5] проводится моделирование смеси жидких кристаллов на протяжении 100 нс. При этом исследуемая система (смесь четырёх жидких кристаллов) состоит из 256 молекул (около 10 тысяч атомов), чего явно недостаточно для получения достоверных физических характеристик при моделировании длинных молекул.
Многомасштабное моделирование молекулярной динамики. Однако существуют приёмы, позволяющие исследовать методом МД большое количество молекул на заметных временных интервалах. Во-первых, увеличивают шаг моделирования, который, вообще говоря, ограничен сверху характеристическим временем наиболее быстрых движений в системе (0,5-1 фс в случае модели с жёсткими химическими связями). Во-вторых, для ускорения моделирования уменьшают количество взаимодействующих центров в молекулах. Упрощённые, или «огрублённые», модели применяются в компьютерном моделировании давно. Они требуют меньше вычислительных ресурсов по сравнению с полноатомными моделями. Случаем крайнего «огрубления» молекул ЖК можно считать модель Гея—Берне [6], в которой молекула представляется единым центром взаимодействия, имеющим вытянутую форму. Примером менее «грубой» модели, применявшейся для жидких кристаллов, служит модель, использованная в работе Кросса и Фанга [7] при исследовании 4-н-пентил-4'-цианобифенила (5ЦБ). Другими примерами могут быть широко распространённая модель «объединённых атомов» [8], а также разнообразные модели «грубых зёрен» (coarse-grain) [9-11], где группы близко расположенных атомов (например, несколько метильных групп или все атомы
* Работа выполнена при финансовой поддержке РФФИ (грант № 09-03-01105-а) и внутреннего гранта СПбГУ.
© В.С.Неверов, А. В. Комолкин, 2012
бензольного кольца) заменяются одним центром взаимодействия, который обычно называют «суператомом». Кроме непосредственного уменьшения количества взаимодействующих центров в системе, «грубые» модели позволяют ещё больше увеличить шаг моделирования, так как минимальная масса «суператомов» в них, как минимум, на порядок больше, чем в полноатомной модели, соответственно больше и характеристическое время движений «суператомов».
Обратной стороной использования «грубых» потенциалов является недостоверность физических характеристик исследуемых систем, получающихся в результате расчётов, так как модель упрощается, из неё удаляются атомы водорода, расположенные в основном на периферии молекул и участвующие в межмолекулярных взаимодействиях. В частности, при моделировании с использованием «грубых» потенциалов невозможно исследовать свойства, зависящие от взаимного расположения атомов одной или разных молекул. Иными словами, для исследования молекулярных свойств вещества необходимо применять детальные (полноатомные) модели, а для увеличения времени эволюции системы — «грубые» модели.
В литературе предложено несколько способов преобразования полноатомных моделей в «грубые» и обратно [12-15]. В нашем исследовании использован метод MuScaT [16] (от multiscale transitions — переходы между масштабами) собственной разработки. Он основан на геометрических построениях и позволяет в рамках одного непрерывного моделирования осуществлять переходы между моделями разной детализации. Молекулы разных типов преобразуются каждая в свой «огрублённый» или, наоборот, детализированный аналог, но при этом процедура уменьшения/увеличения детализации (т. н. coarse-graining procedure и fine-graining procedure) — одинакова для всех типов молекул. Подробное описание метода и его сравнение с другими приведено в нашей статье [17].
Модель Big Ben. Бензол и его производные являются распространённым молекулярным фрагментом, встречающимся во многих классах молекул, в частности, в жидких кристаллах, макромолекулах и пр. Например, молекула термотропного жидкого кристалла 4-н-пентил-4'-цианобифенила (5ЦБ) включает в себя два бензольных кольца, атомы которых составляют более половины атомов всей молекулы (20 атомов из 38). Бензол состоит из 12 атомов в полноатомной модели, имеет практически сферическую симметрию и очевидным образом является хорошим кандидатом для представления в «грубой» модели.
Существует несколько моделей, разработанных специально для «грубого» представления бензола [7, 18, 20]. Кроме них можно использовать и общие подходы, например, потенциалы «объединённых атомов» [8] или MARTINI [14]. Но все эти подходы имеют свои плюсы и минусы. Бензол в модели «объединённых атомов» состоит из 6 «суператомов» (углерод + водород), что довольно много для эффективности расчётов. Кросс и Фанг в [7] используют для бензольного кольца одноточечную модель. Это очень эффективный подход, но он слишком «грубый» и приводит к плохому описанию свойств системы. В частности, теряется информация об ориентации кольца. То есть, во время процедуры «детализации» системы эта ориентация должна воссоздаваться случайным или каким-либо другим образом. Потенциалы MARTINI также представляют бензол в виде единственного центра взаимодействия, и недостатки тут те же самые. Четырёхточечная модель, описанная в [18], достаточно хорошо воспроизводит экспериментальные свойства бензола, но сравнительно дорога с точки зрения эффективности расчётов.
В предыдущей работе [19] мы представили новую «грубую» модель бензольного кольца с заместителями в пара-положении в молекуле («Big Ben»). Эта модель состоит
из двух тяжёлых центров взаимодействия («суператомов»), лежащих в плоскости кольца на линии, перпендикулярной к прямой, соединяющей заместители. Подобная двухточечная модель практически так же эф- _ фективна в расчётах, как и одноточечная, однако «за- - „ поминает» ориентацию кольца, что важно при «детализации» исследуемой системы. Модель бензольного кольца с заместителями в пара-положении показана
на рис. 1. Вес каждого «суператома» — 39 а. е. м., рис ^ Бензольное кольцо то есть половина общего веса всех атомов бензольно- с заместителями
го кольца. Суператомы соединены связями напрямую в пара-пож>жении
с заместителями. Четыре точки (две, принадлежащие в «крупнозернистой» модели бензольному кольцу, и две, принадлежащие замести- и модели «объединённых атомов» телям) образуют ромб. Расстояния между атомами и суператомами поддерживаются неизменными, параметры двугранных углов заданы таким образом, чтобы они могли немного изгибаться, имитируя тем самым движение заместителей друг относительно друга.
Для воспроизведения локальной структуры, формируемой бензольным кольцом, чтобы вычислить ван-дер-ваальсовы взаимодействия атомов двухточечной модели с остальными атомами и друг с другом используется модифицированный потенциал Леннард—Джонса «12—6»
Ег,
4е<,
"ю ДгЮ
12
ггЮ ДгЮ
(1)
Это — упрощённая форма известного потенциала WCA [21]. Смысл а^ и е^ здесь такой же, как в классическом выражении для энергии Леннард—Джонса: а обозначает расстояние между двумя атомами, на котором энергия их взаимодействия равна нулю, е — значение минимума этой энергии.
В нашем выражении для энергии (1) параметр Д сдвигает сферу Леннард—Джонса, радиус которой в классическом варианте определяется исключительно параметром а. Мы не обрезаем взаимодействия на расстоянии \/2а+Д, как это делается в \¥СА. Идея в данном случае заключается в том, чтобы задавать Д индивидуально для каждой пары типов взаимодействующих атомов (или даже для каждой пары атомов в молекуле), чтобы расположение атомов других молекул воспроизводило локальную структуру потенциалов «объединённых атомов».
При разработке двухточечной модели бензола мы считали потенциалы «объединённых атомов» эталоном, так как они широко используются и обладают хорошо известными достоинствами и недостатками. Параметры взаимодействия суператомов с прочими атомами молекулы подбирались таким образом, чтобы локальная структура вещества в «грубой» модели по возможности совпадала со структурой в модели «объединённых атомов». Кроме того, все атомы и группы атомов молекулы, не принадлежащие бензольным кольцам с заместителями в пара-положении, должны описываться с помощью этих потенциалов.
Ранее в [19] нами проведена проверка модели на трёх простых веществах, содержащих бензольное кольцо с заместителями в пара-положении: пара-ксилол, 1-этил-4-ме-тилбензол и 1-метил-4-пропилбензол. Рассмотрены два варианта соединения атомов внутри ромба двухточечной модели виртуальными связями и их влияние на кон-формации в алифатических цепях заместителей (1-этил-4-метилбензол). Кроме того,
6
а
а
рассмотрен случай длинной алифатической цепи. Было показано, что взаимодействия, возникающие между атомами такой цепи и суператомами двухточечной модели бензола, необходимо параметризовать особым образом. В качестве решения предложено разделять параметры, описывающие ван-дер-ваальсово взаимодействие между двумя атомами по модифицированному уравнению (1), на меж- и внутримолекулярные. В целом проверка двухточечной модели при расчёте простых веществ показала, что параметры систем, предсказываемые этой моделью, очень близки к тем, которые предсказываются моделью «объединённых атомов», и что замена 6 «объединённых атомов» на 2 суператома и модификация выражения для энергии Леннард—Джонса возможны и оправданы. Параметры взаимодействия всех типов пар атомов в трёх исследованных веществах приведены в табл. 1. Для пара-ксилола действительны только межмолекулярные параметры.
Таблица 1
Меж- и внутримолекулярные параметры леннард-джонсовского взаимодействия всех пар типов атомов в пара-ксилоле, 1-этил-4-метилбензоле и 1-метил-4-пропилбензоле
Пары Межмолекулярные Внутримолекулярные
атомов е, кДж/моль 0, Ä Д, Ä е, кДж/моль 0, Ä Д, Ä
CG-CG 0,227 3,75 -0,15 0,270 3,75 1,1
CG-CH2 0,163 3,827 -0,15 0,178 3,826 0,55
GG-СНз 0,199 3,827 -0,15 0,217 3,826 0,55
СНт—СНз 0,118 3,905 0,0 0,118 3,905 0,0
СНт—СНз 0,143 3,905 0,0 0,143 3,905 0,0
СН3-СН3 0,175 3,905 0,0 0,175 3,905 0,0
Задача и параметры. Задача настоящей работы состояла в том, чтобы продемонстрировать применимость двухточечной модели бензольного кольца к жидким кристаллам, исследуемым с помощью метода МД. Нами выбраны два нематических жидких кристалла: МББА и БОБТ. МББА очень хорошо изучен и может считаться эталонным веществом. БОБТ является изомером МББА и методом МД изучался впервые в работе [17]. Кроме того, сделана попытка представить с помощью двухточечной модели молекулу жидкого кристалла 5ЦБ, которая содержит бифенил, что усложняет задачу по «огрублению» бензольных колец.
Жидкие кристаллы являются хорошим тестом для модели из-за присущей им ори-ентационной упорядоченности. При моделировании анизотропного вещества молекулы должны располагаться друг относительно друга определённым образом и группироваться в определённых положениях. Воспроизведение жидкокристаллической фазы при моделировании молекул, описанных с помощью двухточечной модели бензольного кольца, будет означать принципиальную применимость этой модели для ускорения расчётов вязких жидкостей, состоящих из сложных молекул.
МББА и БОБТ. Вещества МББА и БОБТ принадлежат к классу оснований Шиф-фа или азометинам и являются изомерами. В молекулах обоих веществ к фрагменту бензилиденанилина с одной стороны алифатическая цепь присоединена напрямую, а с другой — через атом кислорода. Различие молекул в том, что в МББА к атому кислорода присоединена метильная группа, а к бензольному кольцу — цепочка С4Н9, а в БОБТ — наоборот. МББА и БОБТ различаются температурами фазовых переходов,
ограничивающих жидкокристаллическую фазу: у БОБТ эта фаза наблюдается при температурах от 339 до 341 К, у МББА — от 293 до 320 К.
В молекулах обоих веществ присутствуют два бензольных кольца с заместителями в пара-положении. Кольца соединены с основанием Шиффа с одной стороны и с алифатическими цепями — с другой. В полноатомной модели каждая молекула состоит из 41 атома. В модели «объединённых атомов» — из 20. Если представить молекулы этих веществ с помощью обсуждаемой двухточечной модели, то количество взаимодей-
ствующих центров сокращается до 12. На рис. 2 изображены молекулы обоих веществ в модели «объединённых атомов» и в двухточечной. Видно, что бензольное кольцо с заместителями представляется в виде ромба с одной связью внутри. Это виртуальная связь, которая призвана поддерживать правильную форму ромба. Стороны ромба не одинаковы, их длины зависят от заместителей. В табл. 2 приведены геометрические параметры модели в обеих молекулах (нумерация атомов соответствует приведённой на рисунке).
Таблица 2 Геометрические параметры двухточечной модели бензольного кольца в молекулах МББА и БОБТ
МББА БОБТ Длина, А
2-3 и 2-4 2-3 и 2-4 3,085
3-5 и 4-5 3-5 и 4-5 3,033
1-7 и 1-8 1-10 и 1-11 3,012
7-9 и 8-9 10-12 и 11-12 3,162
3-4 и 7-8 3-4 и 10-11 2,5
20 МББА
^ 19_18 14_12
17 11_N 4_6 10.
3 8 9/
12 15 13
10 7
О
17 15
20
18 16 10
БОБТ 12
14_N 4_6 10,_/
Л / Ч / 'и
2 3 8 О
9
13
Рис. 2. Структура молекул МББА и БОБТ в модели «объединённых атомов» (сверху)
и «грубой» (снизу)
4
4
В молекулярной динамике внутримолекулярная энергия рассчитывается как сумма компонент:
¿intra = ¿bond + ¿angle + ¿dih + ¿imp + ¿nb,
где ¿bond — энергия растяжения валентных связей; ¿angle — энергия искажения валентных углов; ¿dih — энергия вращения двугранных углов; ¿imp — энергия неправильных двугранных углов и ¿nb — энергия взаимодействия несвязанных атомов. Во всех наших расчётах для поддержания длин связей молекул применяется алгоритм SHAKE [22], что означает: ¿bond = 0. Энергия искажения валентных углов определяется по формуле
где К — коэффициент «упругости» валентного угла; ф0 — равновесное значение угла. Параметры валентных углов в случае исследуемых жидких кристаллов совпадают с аналогичными параметрами из набора потенциалов «объединённых атомов». Энергия неправильных двугранных углов равна нулю из-за отсутствия неправильных двугранных углов. Энергия несвязанных атомов вычисляется как сумма энергий ван-дер-ваальсового и электростатического взаимодействий. Энергия вращения двугранных углов вычисляется по формуле
Здесь Уп — значение энергетического барьера; 6 — двугранный угол; 6п — угол, при котором данный член суммы обращается в ноль.
Следует отметить, что выражение (2) имеет немного другой вид, нежели аналогичные выражения для вычисления энергии вращения двугранных углов, применяющиеся в других программах моделирования, например в GROMACS [23]. Поэтому параметры, приведённые из (2), отличаются от параметров из других источников, хотя и согласованы с ними. Причина — в способе реализации расчёта этой энергии в используемой нами программе моделирования [2].
Двугранные углы по связи между двумя суператомами двухточечного бензола (например, 1-10-11-12 в БОБТ) заданы таким образом, чтобы позволить молекуле немного сгибаться вокруг этой связи. Параметры двугранных углов обоих жидких кристаллов в двухточечной модели приведены в табл. 3. Величины У, 6 и п, указанные в таблице, использовались в уравнении (2).
Параметры расчётов. Модельные системы МББА и БОБТ были взяты нами из работы [17]. Обе системы состояли из 726 молекул каждая. Системы создавались в модели «объединённых атомов» и уравновешивались. В начальный момент времени молекулы располагались в узлах пространственной решётки 11 х 11 х 6 и были преимущественно ориентированы вдоль оси Z с отклонением не больше 35°. Плотность систем составляла 0,1 г/см3. Моделирование начиналось в ансамбле NVT (использовался термостат Берендсена [4]), в котором система сжималась до плотности, свойственной ей при нормальных условиях (МББА — 1,039 г/см3, БОБТ — 1,007 г/см3). Окончательное уравновешивание производилось в ансамбле NpT (по алгоритму Хувера [24]).
Первые 12 нс моделирование производилось в моделях «объединённых атомов» и полноатомной (см. [17]). Затем был произведён переход в «грубую» модель, моделирование в которой продолжалось до 19 нс. После был произведён переход в полноатомную модель, и в ней система уравновешивалась в течение 1 нс. Таким образом,
¿angle = K(ф - фс>)2,
(2)
n=1
Таблица 3
Параметры двугранных углов молекул МББА и БОБТ в двухточечной модели
МББА БОБТ
Угол V, кДж/моль е, ° п Угол V, кДж/моль е, ° n
9-7-8-1 87 180 1 12-10-11-1 87 180 1
9-8-7-1 87 180 1 12-11-10-1 87 180 1
2-3-4-5 87 180 1 2-3-4-5 87 180 1
2-4-3-5 87 180 1 2-4-3-5 87 180 1
3-2-5-6 0,48 90 2 3-2-5-6 0,48 90 2
4-2-5-6 0,48 90 2 4-2-5-6 0,48 90 2
7-1-9-10 0,48 90 2 2-5-6-7 3,5 180 1
8-1-9-10 0,48 90 2 2-5-6-7 2,9 180 3
7-9-1-2 10 45 4 5-6-7-8 2,9 180 3
8-9-1-2 10 45 4 5-6-7-8 1,0 180 1
9-1-2-5 80 180 2 6-7-8-9 2,9 180 3
3-5-2-1 80 0 2 10-12-1-2 10 45 4
4-5-2-1 80 0 2 11-12-1-2 10 45 4
1-9-10-11 3,5 180 1 12-1-2-5 80 180 2
2,9 180 3 3-5-2-1 80 0 2
9 10 11 12 2,9 180 3 4-5-2-1 80 0 2
1,0 180 1 - - - -
имеется возможность сравнить характеристики полноатомной системы до и после эволюции в «грубой» модели.
Кроме того, оба жидких кристалла были уравновешены в «грубой» модели при атмосферном давлении и различных температурах (МББА: 293, 300, 307, 314 и 320 K, БОБТ: 337, 340 и 343 K) в течение 2 нс (при каждой температуре), что позволило изучить температурные зависимости некоторых свойств.
Во всех расчётах использовалась кубическая ячейка с наложенными периодическими граничными условиями. Интегрирование уравнений движения производилось по методу Верле (leap-frog) [4]. Электростатические взаимодействия учитывались по методу Эвальда [25], радиус отсечения взаимодействий составлял 14,5 A, параметр а метода Эвальда равнялся 0, 27. Характеристическое время термостата Нозе—Хувера равнялось 0,5047 пс, время баростата Хувера — 5,047 пс. Для поддержания длины валентных связей применялся алгоритм SHAKE [22]. В модели «объединённых атомов» использовался список ближних соседей Верле [26], который обновлялся каждые 100 шагов. Шаг моделирования в «грубой» модели и в модели «объединённых атомов» составлял 2 фс, в полноатомной модели — 1 фс.
Результаты и обсуждение.
Параметр порядка жидких кристаллов. Жидкокристаллическая фаза обоих веществ воспроизводится в компьютерном эксперименте, когда бензольные кольца представлены в виде двухточечной модели. Параметры порядка и двуосности молекул в зависимости от температуры приведены в табл. 4. Там же представлены для сравнения эти же величины для полноатомной модели, полученные в результате исследования [17]. Эти данные имеются для меньшего количества точек на температурной шкале. Приведённые результаты рассчитаны как среднее за 2 нс для «грубой» и 1 нс для полноатомной системы.
Таблица 4
Параметры порядка МББА и БОБТ при разных температурах в «грубой» и полноатомной моделях
Вещество Т, К СО РА
$ху
293 0,643 0,056
300 0,645 0,078 0,649 0,064
МББА 307 0,628 0,064 0,666 0,056
314 0,625 0,061 0,620 0,058
320 0,628 0,080
337 0,605 0,053
БОБТ 340 343 0,618 0,606 0,059 0,054 0,697 0,056
Параметр порядка (Р2(еов 6), полином Лежандра второго порядка) рассчитывался следующим образом:
Параметр двуосности молекул рассчитывался по формуле
/ 3я1г2 '
<? — —
^уу
3 эпг 8 соэ 2ф ^
2
В приведённых формулах 6 обозначает угол между длинной осью молекулы и направлением директора жидкого кристалла; ф — угол между осью А" молекулярной системы координат и проекцией директора ЖК на плоскость ХУ молекулярной системы.
Приведённые в таблице значения позволяют сделать вывод, что жидкокристаллическая фаза исследованных жидких кристаллов воспроизводится моделью хорошо, вещество находится в анизотропной фазе. Однако температурная зависимость не наблюдается. На наш взгляд, это связано с тем, что даже 2 нс уравновешивания такой вязкой системы, какой является жидкий кристалл, недостаточно для того, чтобы параметр порядка пришёл к своему среднему значению. В данном случае мы приводим полученные значения, чтобы показать, что макроскопические характеристики жидкокристаллической фазы — параметр порядка и параметр двуосности — воспроизводятся «грубой» моделью так же хорошо, как и полноатомной. Однако эти величины нельзя считать характеристиками равновесной системы.
Коэффициенты самодиффузии. Коэффициенты трансляционной диффузии могут быть получены из формулы
Вп
(3)
где а = х, у или г определяет координатную ось, а На(Ь) — положение центра масс молекулы в момент времени Ь. В реальных условиях предел заменяют достаточно длинными конечными траекториями движения молекул.
На рис. 3 приведена зависимость параллельной (Бц) и перпендикулярной (Б±) директору ЖК компонент коэффициента самодиффузии в «грубой» модели МББА при температуре Т = 307 К от времени анализа траектории (Ь в выражении (3)). Нижняя кривая — зависимость перпендикулярной компоненты, верхняя — параллельной. На рисунке видно, что перпендикулярная директору ЖК компонента коэффициента самодиффузии приходит к равновесному значению после 3,5 нс, а параллельная — после 4,5 нс. Коэффициент самодиффузии рассчитывается как среднее по горизонтальной (равновесной) части кривой.
Кроме того, коэффициенты трансляционной диффузии могут быть рассчитаны из автокорреляционной функции линейной скорости центра масс молекул, однако этот метод даёт больший разброс значений при анализе результатов моделирования (см. [17]), и в данной работе мы его не используем.
Коэффициенты самодиффузии для двух исследованных жидких кристаллов приведены в табл. 5. Приведённые значения нельзя считать точными, это — верхние оценки.
5е-12
4,5е-12 4е-12 3,5е-12 3е-12 2,5е-12 2е-12 1,5е-12 1е-12 5е-13
1000 2000 3000 4000 5000 6000
3. Зависимость параллельной и перпендикулярной компонент коэффициента самодиффузии в МББА при Т = 307 K от времени: длина траекторий 6 нс
В «грубой» модели они рассчитывались по траекториям движения молекул, записанным в течение последних 1,6 нс из 2 нс моделирования при каждой температуре. В полноатомной модели — по траекториям, записанным в течение последних 0,8 нс из 1,6 нс при каждой температуре. На рис. 4, где приведена зависимость компонент коэффициента СД от времени анализа, видно, что такое время анализа не даёт возможности рассчитать истинное значение коэффициента для данной модели.
Из-за недостаточного времени анализа системы при каждой температуре температурная зависимость, характерная для коэффициентов самодиффузии, не наблюдается. Для получения достоверных данных время уравновешивания должно быть увеличено в два-три раза. Такой расчёт в полноатомной модели является весьма
0
Рис.
Таблица 5 Коэффициенты самодиффузии в МББА и БОБТ при разных температурах в «грубой» и полноатомной моделях, 10-12 м2/с
Вещество Г, К СО БА
Б* О! Б* Б| О!
293 2,7 3,9 2Д
300 2,7 4,3 2,0 3,2 3,3 2,7
МББА 307 3,1 4,7 2,3 4,5 6,7 3,5
314 4,1 6,1 3,1 10 15 7,8
320 4,8 7,8 3,3
337 6,9 10,4 5Д
338 7,0 10,5 5,3
339 6,4 10,7 4,3
БОБТ 340 341 342 343 9,5 8.7 8,3 7.8 15,0 12,4 12,7 11,7 6,8 6,9 6,0 5,8 13 21 8,3
8e-12
7e-12 6e-12 5e-12 4e-12 3e-12 2e-12 1e-12
0 200 400 600 800 1000 1200 1400 1600
Рис. 4- Зависимость параллельной и перпендикулярной компонент коэффициента самодиффузии в «грубой» модели МББА при T = 307 K от длины анализируемых траекторий: максимальная длина траекторий 1,6 нс
ресурсоёмким. Например, 1,5 нс уравновешивания требуют расчётов в течение 280 ч (12 сут) на 16 процессорах Intel Xeon E5520. В «грубой» модели провести такое уравновешивание проще. Такой расчёт мы сделали для МББА при температуре T = 307 K. На рис. 3 и 4 приведена одна и та же зависимость компонент коэффициента самодиффузии от времени анализа, рассчитанная по траекториям длиной 6 и 1,6 нс соответственно. В результате длительного анализа получились такие значения коэффициента: D = 1,5 • 10~12 м2/с, D\\ = 2,3 • 10~12 м2/с, D± = 1,1 • 10~12 м2/с. Если сравнить их со значениями из таблицы, полученными в результате анализа траекторий в течение 1,6 нс, то видно, что оценки в таблице оказываются завышенными в два раза.
Таким образом, чтобы правильно рассчитать коэффициенты самодиффузии в жидком кристалле, необходимо анализировать систему в течение, как минимум, 5 нс.
Экспериментальное значение коэффициента самодиффузии в МББА, полученное методом нейтронного рассеяния при 303 K [27]: D = 20 • 10~12 м2/с. Вычисленное значение при T = 307 K отличается от экспериментального для полноатомной модели — в 2 раза, для «грубой» — почти в семь раз.
Двумерные распределения. Чтобы оценить, как двухточечная модель воспроизводит локальную структуру вещества, мы построили цилиндрические функции распределения (ЦФР). ЦФР являются представлением двухчастичных функций распределения, как и широко применяемые функции радиального распределения (ФРР). ЦФР строятся в молекулярной системе координат и показывают расположение атомов молекул друг относительно друга. В начале координат находится центр масс молекулы.
Вертикальная ось совпадает с направлением оси Z тензора инерции молекулы, то есть, с длинной молекулярной осью. По горизонтальной оси откладывается абсолютная величина расстояния в плоскости ХУ молекулярной системы координат, т. е. изображения симметричны относительно оси Z. ЦФР показывают плотность вероятности расположения атома или группы атомов в определённом месте относительно центра масс молекулы. В отличие от сферически симметричных ФРР ЦФР дают информацию о вероятности найти атом в элементарном цилиндрическом слое данного радиуса и на данной высоте.
На рис. 5 и 6 изображены пространственные межмолекулярные распределения плотности центров масс молекул системы друг относительно друга в двух моделях: полноатомной и «грубой». На эти распределения на обоих рисунках наложены внутримолекулярные распределения атомов относительно центра масс молекулы. В случае полноатомной модели изображены распределения всех атомов, кроме атомов водорода. Ориентация бензилиденанилинового фрагмента обеих молекул идентична, атом кислорода находится в верхней части рисунков и обозначен буквой «О». Горизонтальные чёрные линии отделяют бензилиденанилиновый фрагмент от алифатических цепей. Рисунки выполнены в градациях серого, где белый цвет означает отсутствие атомов, серый — среднюю плотность центров масс или атомов, а чёрный — плотность, в два или более раз превышающую среднюю.
Белый цвет, на фоне которого построено внутримолекулярное распределение атомов, обозначает т. н. исключённый объём молекулы — объём, в котором не могут находиться атомы других молекул. Длина и ширина исключённых объёмов обеих молекул в двух сравниваемых моделях приведены в табл. 6.
На рисунках видно, что исключённые
Полноатомная модель
О
«Г Ч.
§ =N1
«Грубая» модель
О
I* | |
И Чм** г
Ы
Рис. 5. ЦФР в МББА в полноатомной и «грубой» моделях
Полноатомная модель
I -»о
«Грубая»
МОЧСЛЬ
¡О»
N V
К".
% Ф
Рис. 6. ЦФР в БОБТ в полноатомной и «грубой» моделях
Таблица 6
Размеры исключённых объёмов молекул МББА и БОБТ в полноатомной и «грубой» моделях
Параметры, А МББА БОБТ
ГА СО ГА СО
Длина 19,02 20,4 20,15 20,78
Ширина 4,64 5,9 5,02 6,15
объёмы молекул обоих веществ в «грубой» модели имеют цилиндрическую форму, схожую с той, которую они имеют в полноатомной модели. Из таблицы видно, что в «грубой» модели исключённый объём немного больше: длиннее и шире.
Избыточная ширина объясняется тем, что, хотя при подгонке параметров леннард-джонсовского взаимодействия суператомов двухточечной модели бензола с атомами других молекул используется весьма ограниченный набор простейших конфигураций, одновременно удовлетворительно воспроизвести зависимости энергии от расстояния во всех конфигурациях оказывается весьма непросто (подробнее о процедуре подгонки см. [19]). Провести подгонку по большему числу конфигураций практически невозможно. Однако, на наш взгляд, этого и не требуется, так как «грубая» модель призвана, в первую очередь, с удовлетворительной точностью воспроизводить только расположение молекул друг относительно друга, чтобы после перехода в полноатомную модель и недолгого моделирования получить хорошо уравновешенную систему.
Немного большая длина исключённого объёма молекул в «грубой» модели объясняется тем, что эта модель более жёсткая, в ней меньше связей и углов. Такая модель в среднем должна быть чуть вытянутее, чем полноатомная. Скорее всего, немного большая длина приводит лишь к небольшому завышению параметра порядка жидкого кристалла в «грубой» модели, однако в этой модели любые макроскопические параметры можно считать лишь приближёнными оценками.
Расположение молекул друг относительно друга модель воспроизводит удовлетворительно. На рисунках видно, что «грубая» модель приводит к более структурированному веществу, чем полноатомная, области концентрации центров масс молекул гораздо более выражены. Центры масс молекул сгруппированы в три слоя: возле атома азота, возле атома кислорода и возле первой группы в алифатической цепи, соединённой с бензольным кольцом напрямую. В полноатомной модели наблюдается подобная группировка, но значительно в меньшей степени. На наш взгляд, более сильная структурированность вещества в «грубой» модели не является недостатком, мешающим использовать её в многомасштабном моделировании, так как, во-первых, качественно модель воспроизводит структуру правильно, во-вторых, анализ структуры вещества необходимо делать в полноатомной модели.
Средний угол между длинными осями молекул в молекулярной системе координат составляет: в «грубой» модели — 38° у МББА, 35° у БОБТ, в полноатомной модели: 35,5° у МББА, 28,5° у БОБТ. Распределение плотности вероятности этого угла для двух веществ в полноатомной и «грубой» моделях показано на рис. 7.
Бифенил. 5ЦБ. Отдельной задачей является использование обсуждаемой двухточечной модели бензольного кольца для моделирования веществ, содержащих бифенил. Основной сложностью при замене бензольных колец в таких веществах на двухточечную модель является то, что одним из заместителей у обоих колец является другое кольцо, также подлежащее замене. Иными словами, модель в том виде, в котором она используется в молекулах МББА и БОБТ, не годится для представления бифенила, и её нужно адаптировать.
В качестве пробного вещества для этой задачи мы выбрали вещество 5ЦБ, первый жидкий кристалл, синтезированный в классе цианобифенилов в 1972 году. Он состоит из бифенила, к которому с одной стороны присоединена цианогруппа -С=^ а с другой — алифатическая цепь С5И11. Температура перехода 5ЦБ из кристаллического состояния в мезофазу — 291 К, а из мезофазы в изотропное состояние — 308 К.
В полноатомной модели молекула 5ЦБ состоит из 38 атомов. В модели «объединённых атомов» — из 19. В обсуждаемой двухточечной модели — из 11.
На рис. 8 изображено одно из возможных решений задачи адаптации модели для представления бифенила. Каждое кольцо представлено в виде двух суператомов, которые с одной стороны соединены виртуальными связями с заместителями. С другой
0,09 0,08 0,07 0,06 0,05 0,04 0,03 0,02 0,01
0 20 40 60 80 100 120 140 160 180
Рис. 7. Распределение плотности вероятности угла между длинными осями жидких кристаллов МББА и БОБТ в полноатомной и «грубой» моделях
стороны виртуальные связи соединяют суператомы не с заместителем — другим бензольным кольцом, а с ближайшим атомом за этим кольцом. Кроме того, добавлена ещё одна виртуальная связь между заместителями (связь 2-7), которая поддерживает правильную форму ромбов. Такая конфигурация позволяет каждому из колец вращаться независимо от другого. В табл. 7 приведены геометрические параметры модели, а в табл. 8 — параметры двугранных углов. Указанные в ней величины V, 8 и п, использовались в уравнении (2).
Расчёт. Наша модельная система 5ЦБ состоит из 2888 молекул. Она создавалась в модели «объединённых атомов». В начальный момент времени молекулы располагались в узлах пространственной решётки 19x19 х 8 и также были преимущественно ориентированы вдоль оси Z с отклонением не больше 35°. Плотность системы составляла 0,1 г/см3. Создание и первичное уравновешивание системы происходили так же, как в случае МББА и БОБТ, плотность после сжатия составляла 1,02 г/см3. Далее в течение 11,5 нс система эволюционировала в модели «объединённых атомов» в ансамбле NVT при атмосферном давлении и температуре 300 К. Затем был произведён переход в «грубую» модель, в которой на сегодняшний день продолжается расчёт.
Пока что мы не можем привести детальных характеристик системы 5ЦБ, представленной с помощью двухточечной модели. Это — тема отдельной публикации.
Таблица 7
Геометрические параметры
двухточечной модели бифенила в молекуле 5ЦБ
5ЦБ Длина, A
2- -3 и 2- -4 3,085
3- -7 и 4- -7 7,330
2- -6 и 2- -6 7,250
5- -7 и 6- -7 3,162
2-7 10,047
5 ЦБ
18 19
4_6 10_12 16
5 7 11 13
3 5
4
6
Рис. 8. Структура молекулы 5ЦБ в моделях «объединённых атомов» (сверху) и «грубой» (снизу)
Таблица 8 Параметры двугранных углов молекулы 5ЦБ в двухточечной модели
В настоящий момент можно сказать, что молекулу, содержащую бифенил, удалось представить с помощью этой модели описанным выше способом, что её возможно использовать в многомасштабном моделировании и что модель воспроизводит жидкокристаллическую фазу вещества. Всего промоделировано 13,5 нс эволюции системы, из них последние 2 нс — в «грубой» модели. Параметр порядка в «грубой» модели при T = 300 K равен 0,75.
Ускорение расчётов. Главной целью работы являлось создание модели, которая позволит ускорить вычисления. В проведённых моделированиях с использованием трёх вариантов предложенная нами «грубая» модель оказалась наиболее эффективной. В табл. 9 приведены значения времени, затраченного на расчёт одного шага моделирования, для всех использованных веществ в трёх моделях. При этом надо учитывать, что шаг интегрирования уравнений движения в «грубой» модели по сравнению с более детальной может быть увеличен. Значения в столбце «Ускорение» рассчитаны относительно полноатомной модели. Системы МББА и БОБТ состоят из 726 молекул, система 5ЦБ — из 2888. Измерение времени производилось на одном компьютере с 16 процессорами Intel Xeon E5520.
Выводы. Представлена новая модель бензольного кольца с заместителями в пара-положении для многомасштабного моделирования молекулярной динамики и показано её применение для моделирования ряда жидких кристаллов. Модель требует модификации уравнения Леннард—Джонса для вычисления энергии атом-атомного взаимодействия и введения в это уравнение дополнительного параметра сдвига для каждой пары типов взаимодействующих атомов отдельно для меж- и внутримолекулярных взаимодействий. Сама модель состоит из двух суператомов. Она весьма эффективна в расчётах и может сохранять ориентацию бензольного кольца в многомасштабном
Угол V, кДж/моль е, ° n
3-2-7-5 0,61 45 4
3-2-7-6 0,61 45 4
4-2-7-5 0,61 45 4
4-2-7-6 0,61 45 4
3-2-7-5 0,70 0 2
3-2-7-6 0,70 0 2
4-2-7-5 0,70 0 2
4-2-7-6 0,70 0 2
3-2-7-6 -1,20 0 0
5-2-7-8 0,48 90 2
6-2-7-8 0,48 90 2
2-7-8-9 2,9 180 3
3,5 180 1
2,9 180 3
7-8-9-10
1,0 180 1
2,9 180 3
8-9-10-11
1,5 180 1
3-2-7-4 8,7 180 4
5-7-2-6 8,7 180 4
Таблица 9
Длительность одного шага моделирования в разных моделях
Вещество Модель Количество атомов Длительность одного шага, мс Ускорение
FA 29766 1420 1,0
МББА UA 14520 550 2,6
CG 8712 164 8,6
FA 29766 1580 1,0
БОБТ UA 14520 383 4,1
CG 8712 160 9,9
FA 112632 15800 1,0
5ЦБ UA 57760 4130 3,8
CG 31768 2350 6,7
моделировании. Введение дополнительного параметра требует только однократной подготовки входных данных, но никак не замедляет сам расчёт.
Показано, что модель применима для ускорения расчётов жидких кристаллов, содержащих бензольные кольца с заместителями в пара-положении. Локальная структура вещества в ней близка к локальной структуре в полноатомной модели. Жидкокристаллическая фаза воспроизводится хорошо.
Кроме того, предложена модификация модели для представления молекул, содержащих бифенил. Показано, что такая модификация также применима для ускорения расчёта вещества в жидкокристаллической фазе.
Литература
1. Sandstrom D., KomolkinA. V., MaliniakA. Orientational order in a liquid crystalline mixture studied by molecular dynamics simulation and NMR //J. Chem. Phys. 1996. Vol. 104 P. 9620.
2. KomolkinA. V., LaaksonenA., MaliniakA. Molecular dynamics simulation of a nematic liquid crystal // J. Chem. Phys. 1994. Vol. 101. N 5. P. 4103-4116.
3. SocciE. P., FarmerB. L., Bunning T. J. et al. Molecular dynamics and x-ray scattering simulations of cyclic siloxane-based liquid crystal mesogens // Liquid Crystals. 1993. Vol. 13. N 6. P. 811-827.
4. Allen M. P., Tildesley D. J. Computer simulation of liquids. Oxford: Clarendon Press, 1987. 400 p.
5. PelaezJ., 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.
6. GayJ. G., BerneB. J. Modification of the overlap potential to mimic a linear site—site potential // J. Chem. Phys. 1981. Vol. 74. P. 3316.
7. Cross C. W., Fung B. W. Molecular dynamics simulations for cyanobiphenyl liquid crystals //J. Mol. Cryst. Liq. Cryst. 1995. Vol. 262. P. 507-524.
8. 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. 1988. Vol. 110. P. 1657-1666.
9. Izvekov S., VothG.A. A multiscale coarse-graining method for biomolecular systems // J. Phys. Chem. (B). 2005. Vol. 109. N 7. P. 2469-2473.
10. Zacharopoulos N., VergadouN., Theodorou D. N. Coarse graining using pretabulated potentials: Liquid benzene // J. Chem. Phys. 2005. Vol. 122. 244111.
11. WilsonM. R. Progress in computer simulations of liquid crystals // Int. Rev. Phys. Chem. 2005. Vol. 24. N 3-4. P. 421-455.
12. Milano G., Muller-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.
13. Santangelo G., Di Matteo A., Muller-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.
14. Marrink S. J., Fuhrmans M., RisseladaH. J., PerioleX. Coarse-Graining of Condensed Phase and Biomolecular Systems. The MARTINI force field. CRC press, 2008. 456 p.
15. Praprotnik M, Delle SiteL., Kremer K. Multiscale simulation of soft matter: From scale bridging to adaptive resolution // Ann. Rev. Phys. Chem. 2008. Vol. 59. P. 545-571.
16. КомолкинА. В. Органические и гибридные наноматериалы. Иваново. 2009. 300 с.
17. Неверов В. С., КомолкинА. В., Волкова Т. Г. Исследование влияния структурной изомерии на молекулярную подвижность жидких кристаллов методом молекулярной динамики // Вестн. С.-Петерб. ун-та. Сер. 4: Физика, химия. 2011. Вып. 1. С. 34-53.
18. DeVane R., KleinM. L., ChiuC. C. et al. Coarse-grained potential models for phenyl-based molecules: I. Parametrization using experimental data // J. Phys. Chem. (B). 2010. Vol. 114. P. 6386-6393.
19. Neverov V. S., KomolkinA. V. Coarse-grain model of the benzene ring with para-substituents in the molecule //J. Chem. Phys. 2012. (In print.)
20. Megariotis G., Vyrkou A., LeygueA., Theodorou D. N. Systematic coarse graining of 4-cy-ano-4-pentylbiphenyl // Ind. Eng. Chem. Res. 2011. Vol. 50. P. 546-556.
21. Andersen H. C., Weeks J. D., Chandler D. Relationship between the hard-sphere fluid and fluids with realistic repulsive forces // Phys. Rev. (A). 1971. Vol. 4. N 4. P. 1597-1607.
22. Ryckaert J. P., Ciccotti G., Berendsen H. J. C. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes //J. Comp. Phys. 1977. Vol. 23. N 3. P. 327-341.
23. HessB., Kutzner C., van der SpoelD., LindahlE. Gromacs 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation // J. Chem. Theory Comput. 2008. Vol. 4. N 3. P. 435-447.
24. MelchionnaS., Ciccotti G., HolianB. L. Hoover's style molecular dynamics for systems varying in shape and size // J. Mol. Phys. 1993. Vol. 78. P. 533.
25. EwaldP. Investigations of crystals by means of roentgen rays // Ann. Phys. 1921. Vol. 64. P. 253.
26. VerletL. Computer «experiments» on classical fluids. I. Thermodynamical properties of lennard-jones molecules // Phys. Rev. 1967. Vol. 159. P. 98-103.
27. CviklB., Dahlborg U., CepicM. 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.
Статья поступила в редакцию 4 октября 2011 г.
