Журнал физической химии, 2022, T. 96, № 7, стр. 1056-1061

Моделирование сжимаемости многокомпонентной жидкокристаллической фазы CHS1 методами молекулярной динамики

Г. И. Макаров a*, Р. В. Решетникова a, Е. В. Барташевич a

a Южно-Уральский государственный университет
454080 Челябинск, Россия

* E-mail: makarovgi@susu.ru

Поступила в редакцию 08.10.2021
После доработки 29.12.2021
Принята к публикации 10.01.2022

Полный текст (PDF)

Аннотация

Методами молекулярно-динамического моделирования исследована реакция хирального смектического жидкого кристалла CHS1 (смесь ахиральных молекул алкиловых эфиров 2-(4-гидроксифенил)-5-гидроксипиримидинов (LC1, LC2, LC3) и немезогенного хирального компонента: п-бис-[4-(гекс-2-ил-оксикарбонил)фенил]бензола (LC4)) на гидростатическое сжатие, рассчитаны объемный изотермический коэффициент сжимаемости и линейные изотермические коэффициенты сжимаемости этой фазы. Обнаружено, что при давлении ниже 50 МПа изучаемая фаза перераспределяет сжатие в направлении оси ориентации (Ox) в перпендикулярном направлении (Oz) за счёт отклонения мезогенов компонентов; при давлении выше 50 МПа фаза CHS1 сжимается равномерно по всем направлениям, что сопровождается увеличением упорядоченности системы.

Ключевые слова: молекулярная динамика, сегнетоэлектрические жидкие кристаллы, механические свойства

Жидкие кристаллы (ЖК) способны проявлять уникальные свойства в зависимости от их состава, в частности, они характеризуются выраженными анизотропными межмолекулярными взаимодействиями. ЖК широко применяются в бытовой технике: как дисплеи телефонов, телевизоров и компьютеров [1, 2], в лабораторном оборудовании, например, для электронных дисплеев для микролазеров [35] и как светомодулирующие среды [6]. Можно считать, что ЖК востребованы в широком диапазоне сфер нашей деятельности, начиная от биомедицины [7] и заканчивая различными разделами физики. Одна из важнейших областей применения ЖК – светомодулирующие среды, где они выполняют роль оптических переключателей. В основном, в данной области используются нематические ЖК, чье время отклика лежит в диапазоне нескольких миллисекунд. Однако для сегнетоэлектрических ЖК время отклика оказывается в микросекундном диапазоне, что более привлекательно для использования сегнетоэлектрических ЖК в устройствах, основанных на электрооптическом отклике.

Сегнетоэлектрические ЖК представляют собой смектические хиральные ЖК (Sm-C*), упорядоченные и ориентационно, и позиционно. Уменьшение времени отклика объясняется тем, что, благодаря такой упорядоченности они имеют симметрию, приводящую к спонтанной поляризации [8], которая крайне чувствительна к внешнему электрическому полю [9]. Чаще всего сегнетоэлектрические ЖК представляют собой многокомпонентные смеси, поскольку чистые соединения образуют SmC*-фазу в узком диапазоне температур, да и добиться проявления нужных физических свойств в таком случае может оказаться трудно. Обычно смесь состоит из четырех компонентов, варьирование соотношений которых позволяет задать требуемые физико-химические свойства. В большинстве случаев такая смесь состоит из ахиральных соединений и одного хирального допанта, который непосредственно влияет на проявление сегнетоэлектрических свойств. Такие смеси изучались по большей части экспериментально в отношении именно свойств хирального допанта и его влияния на поляризацию ЖК-фазы [1012]. Теоретических исследований, основанных на компьютерном моделировании структуры многокомпонентных смектических систем и прогнозировании их свойств, крайне мало. Имеющиеся исследования в основном полагаются на изучение влияния свойств хиральной добавки путем моделирования методом Монте-Карло [1315]. Часть работ посвящена исследованию физических свойств ЖК и фазовых переходов методами классической молекулярной динамики [16, 17], но, к сожалению, такие исследования крайне мало распространены, и совершенствование подходов, достоверно воспроизводящих особенности структуры и свойств многокомпонентных ЖК-фаз с хиральными допантами, очень актуально. Также следует отметить, что исследований многокомпонентных смектиков не проводилось вплоть до 2018 г.

Основная особенность жидкокристаллических тел – их структурная анизотропия, приводящая к векторному характеру свойств – зависимости ряда свойств от направления в фазе. Электрооптические свойства жидких кристаллов определяются их структурной организацией, которая, в свою очередь, может включать дефекты. Поэтому существенный интерес представляют дефекты ЖК, обусловленные присутствующими в них примесями или воздействием деформаций. Исследование дефектов ЖК, возникающих при его деформации, требует внимательного учёта связанных с этим процессом характеристик фазы, таких, как изотермическая сжимаемость. Именно изучению изотермической сжимаемости многокомпонентной SmC*-фазы CHS1 – свойства, необходимого для корректного моделирования этой системы – посвящена настоящая работа, которая является первым шагом к исследованию структурных дефектов ЖК-систем различного происхождения, на данный момент недостаточно изученных.

Мы воспользовались разработанной ранее [8] моделью многокомпонентной хиральной смектической жидкокристаллической фазы (CHS1), которая была успешно смоделирована методами классической молекулярной динамики и уравновешенной метадинамики с обменом потенциалами. Эта структурная модель верно воспроизводит экспериментально измеренную спонтанную поляризацию. Сегнетоэлектрическая жидкокристаллическая смесь CHS1 (рис. 1) состоит из ахиральных молекул алкиловых эфиров 2-(4-гидроксифенил)-5-гидроксипиримидинов (LC1, LC2, LC3) и немезогенного хирального компонента: п-бис-[4-(гекс-2-ил-оксикарбонил)фенил]бензола (LC4).

Рис. 1.

Структура и состав сегнетоэлектрической жидкокристаллической смеси CHS1.

МЕТОДЫ

Мы используем метод моделирования молекулярной динамики (МД) для расчета многокомпонентной SmC*-фазы CHS1, пребывающей в основном низкополяризованном состоянии. Расчет сжимаемости основывался на явном изучении зависимости объема от давления, изучаемая система моделировалась МД при каждом выбранном давлении, для каждой полученной траектории рассчитывался средний объем симуляционной ячейки, причём для усреднения использовался конечный участок траектории, на котором объем переставал направленно изменяться и флуктуировал вблизи некоторого постоянного значения. Физико-математические принципы подхода излагаются ниже.

Моделирование МД и анализ полученных траекторий выполнялись с использованием программного обеспечения GROMACS [18] версии 5.1.4. Шаг интегрирования составлял 2 фс во всех симуляциях, координаты записывались в файл траектории каждые 4 пс, а физические величины, характеризующие систему, сохранялись каждые 0.02 пс. Протяжённость каждой полученной траектории составляла 2 нс, анализировался участок со 100 пс по 2 нс. Уравнения движения интегрировались с помощью алгоритма стохастической динамики. Алгоритм LINCS [19] использовался для ограничения длины связей с атомами водорода. Моделирование проводилось при постоянной температуре (298.15 К), поддерживаемой термостатом масштабирования скоростей с добавочным стохастическим членом [20] с временем привязки 2 пс. Использовались периодические граничные условия при постоянном анизотропном давлении, контролируемом баростатом Берендсена [21] с временем привязки 2 пс; при этом для нормальных напряжений устанавливалась сжимаемость в 1 ГПа–1, а для касательных напряжений сжимаемость устанавливалась равной нулю. Алгоритм сети частиц Эвальда с шагом сетки 0.125 нм и интерполяцией шестого порядка использовался для обработки дальнодействующих электростатических [22] и ван-дер-ваальсовых [23] взаимодействий. Молекулярно-механические модели компонентов жидкокристаллической фазы строились с использованием силового поля GAFF [24] и модели точечных зарядов RESP [25] на основе оптимизированных трехмерных структур и молекулярных электростатических потенциалов, рассчитанных методом HF/6-31G*.

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

SmC*-фаза CHS1 – анизотропная структура, т.е. следует ожидать разных реакций на сжатие по разным направлениям. Поэтому при ее МД-моделировании следует использовать так называемый анизотропный контроль давления, заключающийся в независимом расчёте давления, действующего в разных направлениях моделируемой ячейки, и его независимом приведении к заданным величинам через масштабирование соответствующих параметров ячейки. Этот процесс управляется алгоритмом баростатирования, необходимым параметром которого являются изотермические сжимаемости, определенные для каждого из трех направлений симулируемой ромбической ячейки.

Для вычисления этих величин методами МД-моделирования рассмотрим зависимость объема моделируемой системы от давления. Пусть объем системы V – функция давления p и температуры T, тогда приращение объема выражается дифференциальной формой уравнения состояния:

(1)
$dV = {{\left( {\frac{{dV}}{{dp}}} \right)}_{T}}dp + {{\left( {\frac{{dV}}{{dT}}} \right)}_{p}}dT,$
(2)
$d\ln V = \beta dT - {{\kappa }_{T}}dp.$

В формуле (2) κT – изотермическая сжимаемость, которую мы в дальнейшем будем также называть объемной изотермической сжимаемостью, а β – коэффициент температурного расширения. Предположим, что наша система подвергается сжатию при постоянной температуре. Тогда член уравнения (2), связанный с приращением температуры, обратится в ноль:

(3)
$d\ln V = - {{\kappa }_{T}}dp.$

Предположим, давлению p0 соответствует объем системы V0, а давлению p – объем V. Проинтегрируем уравнение (3) с этими граничными условиями:

(4)
$\mathop \smallint \limits_{{{V}_{0}}}^V \,d\ln V = - {{\kappa }_{T}}\mathop \smallint \limits_{{{p}_{0}}}^p \,dp,$
(5)
$\ln \frac{V}{{{{V}_{0}}}} = - {{\kappa }_{Т}}(p - {{p}_{0}}).$

Используя уравнение (5), можно измерить или рассчитать методами молекулярной динамики объемный коэффициент изотермической сжимаемости κT: для этого требуется установить зависимость объема от давления при постоянной температуре, преобразовать координаты и затем линейно аппроксимировать зависимость логарифма нормированного объема от разности давлений.

Аналогичным образом мы можем рассмотреть параметры симуляционной ячейки как функции температуры и давления тем же образом, что и объем системы. Например, для параметра a справедливо:

(6)
$d\ln a = \frac{1}{a}{{\left( {\frac{{da}}{{dp}}} \right)}_{T}}dp.$
Обозначим μa величину, аналогичную объемному коэффициенту изотермической сжимаемости κT:
(7)
${{\mu }_{a}} = \frac{{ - 1}}{a}{{\left( {\frac{{da}}{{dp}}} \right)}_{T}}.$
Тогда, обозначив a0 величину параметра а при давлении p0, проинтегрируем уравнение (6) аналогично уравнению (4):
(8)
$\ln \left( {\frac{a}{{{{a}_{0}}}}} \right) = - {{\mu }_{a}}(p - {{p}_{0}}).$
Для других параметров ячейки получаются аналогичные уравнения.

Коэффициенты μa, μb и μc назовем частными изотермическими сжимаемостями или линейными изотермическими сжимаемостями. Они характеризуют реакцию анизотропной системы на сжатие в соответствующем направлении при постоянной температуре. Для ромбической ячейки, используемой нами при моделировании SmC*-фазы CHS1, верно следующее:

(9)
$V = abc.$
Разделим обе части на объем:
(10)
$d\ln V = \frac{{dc}}{c} + \frac{{db}}{b} + \frac{{da}}{a} = d\ln c + d\ln b + d\ln a.$
Суммируя все уравнения для линейных изотермических сжимаемостей (10) и учитывая уравнение объема ромбической ячейки, получаем:
(11)
$\begin{gathered} \ln \frac{a}{{{{a}_{0}}}} + \ln \frac{b}{{{{b}_{0}}}} + \ln \frac{c}{{{{c}_{0}}}} = - ({{\mu }_{a}} + {{\mu }_{b}} + {{\mu }_{c}})(p - {{p}_{0}}) = \\ = ln\frac{{abc}}{{{{a}_{0}}{{b}_{0}}{{c}_{0}}}} = \ln \frac{V}{{{{V}_{0}}}} = - {{\kappa }_{T}}(p - {{p}_{0}}). \\ \end{gathered} $
Следовательно, сумма линейных изотермических сжимаемостей равна коэффициенту объемной изотермической сжимаемости:

(12)
${{\kappa }_{T}} = {{\mu }_{a}} + {{\mu }_{b}} + {{\mu }_{c}}.$

Исходя из изложенных выше соображений, мы использовали следующий алгоритм расчета коэффициентов изотермической сжимаемости. Полученное ранее низкополяризованное состояние смеси CHS1, уравновешенное при анизотропном контроле давления, было дополнительно уравновешено на протяжении 2 нс при давлении 0.1 МПа. При анализе этой траектории были получены V0, a0, b0, c0 как средние значения соответствующих величин. Исходя из полученного состояния выполнялся следующий расчет равновесной молекулярной динамики при более высоком давлении, полученное состояние служило исходным для следующего расчета при еще более высоком давлении и так далее. Всего мы выполнили три серии расчетов, исходящих из одного и того же состояния, отвечающего давлению 0.1 МПа: от 0.1 до 4.5 МПа включительно с шагом в 0.5 МПа, от 5 до 45 МПа включительно с шагом в 5 МПа и от 50 до 500 МПа включительно с шагом в 50 МПа. Расчеты проводили при анизотропном контроле давления, которое было одинаковым со всех сторон. Для каждого давления регистрировали средний объем, начиная со 100 пс. По полученным траекториям для каждого из накладываемых давлений также рассчитывались значения угла отклонения мезогена (центральной конформационно жесткой части молекулы) от оси ориентации Θ и ориентационной энтропии мезогена S, определяемой как

(13)
$S = - \sum {P({{\varphi }},\theta )\ln ~P({{\varphi }},\theta )} ,$
где φ и θ – углы, описывающие направление мезогена в полярных координатах.

На рис. 2 можно видеть, что логарифм нормированного объема зависит от разницы давлений близко к ожидаемой линейной зависимости в диапазоне давлений от 50 до 500 МПа включительно, тогда как в диапазоне давлений от 0.5 до 45 МПа включительно система реагирует на избыточное давление сильно нелинейно, причем в наибольшей степени такое поведение характерно не для объема, а для параметров ячейки. Так, в этом диапазоне давлений ребро ячейки a, соответствующее направлению Ox, в котором ориентированы молекулы фазы, сокращается сильнее, чем можно было бы ожидать, исходя из уравнения (8). Напротив, ребро ячейки c, соответствующее направлению Oz, увеличивается при обжатии ячейки. Мы полагаем, что это явление следует интерпретировать так: при равномерном сжатии ячейки со всех сторон давлением менее 50 МПа увеличивается отклонение молекул фазы от оси ориентации в плоскости xOz, отчего их проекции на направление Ox сокращаются, а на направление Oz – увеличиваются, что и приводит к выраженному сокращению ребра a при одновременном росте ребра c. В пользу такого истолкования свидетельствует увеличение среднего угла отклонения мезогена Θ (рис. 2,3): в этом диапазоне давлений он увеличивается на 3 градуса против обычного. Возникающее отклонение пространственной укладки фазы от ее обычного состояния регистрируется также по увеличению ориентационной энтропии.

Рис. 2.

Зависимости логарифмов нормированных объема и параметров ячейки от изменения давления для жидкокристаллической смеси CHS1.

Рис. 3.

Зависимости среднего угла отклонения мезогена Θ и ориентационной энтропии мезогена S от давления.

При давлении больше 50 МПа эта внутренняя упругость системы перестает противостоять сжатию в направлениях Oy и Oz, и система начинает сжиматься равномерно, без увеличения отклонения компонентов от оси ориентации, что регистрируется по уменьшению среднего угла отклонения мезогена Θ к прежнему значению. С ростом давления уменьшается ориентационная энтропия, что закономерно – чем сильнее сжата фаза, тем меньше колеблются уложенные в ней молекулы. В целом, можно заключить, что с ростом давления восстанавливается пространственная укладка фазы. Из-за выявленного необычного поведения системы при давлении менее 50 МПа мы использовали для расчета коэффициентов изотермической сжимаемости точки в диапазоне от 50 до 500 МПа включительно, аппроксимируя их линейными зависимостями в координатах $\ln \frac{V}{{{{V}_{0}}}}$ и $(p - {{p}_{0}})$, основываясь на уравнениях (5) и (8). Полученные значения коэффициентов изотермической сжимаемости (ТПа–1) равны: μa = 52 ± 4, μb = 83 ± 6, μc = 56 ± 4, μa + μb + μc = 190 ± 10, κT = 193 ± 5.

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

Таким образом, используя метод моделирования молекулярной динамики, мы вычислили объемную и линейные изотермические сжимаемости SmC* жидкокристаллической смеси CHS1. Полученные значения найдут применение при дальнейших исследованиях реакции этой фазы на деформации. Также было выявлено, что при давлении ниже 50 МПа фаза CHS1 реагирует на сжатие нелинейно, пытаясь релаксировать приложенное давление посредством изменения своей структуры.

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

  1. Bahadur B. Liquid Crystals – Applications and Uses. Volume 3. Singapore: World Scientific, 1992. 424 p.

  2. Demus D., Goodby J., Gray G.W., Spiess H., Vill V. Handbook of Liquid Crystals. N.Y.: Wiley-VCH, 1998. 897 pages.

  3. Humar M., Ravnik M., Pajk S. et al. // Nat. Photonics. 2009. V. 3. № 10. P. 595.

  4. Peddireddy K., Jampani V.S.R., Thutupalli S. et al. // Opt. Express. 2013. V. 21. № 25. P. 30233.

  5. Humar M., Muševič I. // Opt. Express. 2010. V. 18. № 26. P. 26995.

  6. Chigrinov V.G., Srivastava A. // Mol. Cryst. Liq. Cryst. 2014. V. 595. № 1. P. 39.

  7. Woltman S.J., Jay G.D., Crawford G.P. // Nat. Mater. 2007. V. 6. № 12. P. 929.

  8. Makarov G.I., Bartashevich E.V., Khnykina K.A. et al. // J. Mol. Liq. 2019. V. 283. P. 630.

  9. Miller J.S., Drillon M. Magnetism: Molecules to Materials. N.Y.: Wiley-VCH, 2003. 430 pages.

  10. Lemieux R.P. // Acc. Chem. Res. 2001. V. 34. № 11. P. 845.

  11. Eelkema R., Feringa B.L. // Org. Biomol. Chem. 2006. V. 4. № 20. P. 3729.

  12. Proni G., Spada G.P. // Enantiomer. 2001. V. 6. № 2–3. P. 171.

  13. Cook M.J., Wilson M.R. // J. Chem. Phys. 2000. V. 112. № 3. P. 1560.

  14. Earl D.J., Wilson M.R. // J. Chem. Phys. 2003. V. 119. № 19. P. 10280.

  15. Wilson M.R., Earl D.J. // J. Mater. Chem. 2001. V. 11. № 11. P. 2672.

  16. Quevillon M.J., Whitmer J.K. // Materials (Basel). 2018. V. 11. № 1. P. 64.

  17. Frezzato D., Saielli G. // J. Phys. Chem. B. 2016. V. 120. № 9. P. 2578.

  18. Abraham M.J., Murtola T., Schulz R. et al. // SoftwareX. 2015. V. 1–2. P. 19.

  19. Hess B., Bekker H., Berendsen H.J.C. et al. // J. Comput. Chem. 1997. V. 18. № 12. P. 1463.

  20. Bussi G., Donadio D., Parrinello M. // J. Chem. Phys. 2007. V. 126. № 1. P. 014101.

  21. Berendsen H.J.C., Postma J.P.M., Van Gunsteren W.F. et al. // J. Chem. Phys.1984. V. 81. № 8. P. 3684.

  22. Darden T., York D., Pedersen L. // J. Chem. Phys. 1993. V. 98. № 12. P. 10089.

  23. Wennberg C.L., Murtola T., Hess B. et al. // J. Chem. Theory Comput. 2013. V. 9. № 8. P. 3527.

  24. Wang J., Wolf R.M., Caldwell J.W. et al. // J. Comput. Chem. 2004. V. 25. № 9. P. 1157.

  25. Bayly C.I., Cieplak P., Cornell W.  et al. // J. Phys. Chem. 2002. V. 97. № 40. P. 10269.

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