CC BY
5УДК 547.583.66
Д. О. Москвин1, В. В. Соцкий1, А. С. Кашицын12, Н. В. Усольцева1
МОДЕЛИРОВАНИЕ НАДМОЛЕКУЛЯРНОЙ ОРГАНИЗАЦИИ И ФАЗОВОГО
СОСТОЯНИЯ ГЕКСАГЕПТИЛОКСИТРИФЕНИЛЕНА МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ
MOLECULAR DYNAMICS SIMULATION OF SUPRAMOLECULAR ORGANIZATION AND PHASE STATE OF HEXAHEPTYLOXYTRIPHENYLENE
!НИИ наноматериалов ФГБОУ ВПО «Ивановский государственный университет», 153025 Иваново, ул. Ермака, д. 37. E-mail: Moskvin1987@mail.ru.
2Ивановский государственный университет, Шуйский филиал, 155908, Ивановская область, Шуя, ул. Кооперативная, д. 24.
В работе приводятся результаты расчета и интерпретации функций радиального распределения, полученных методом молекулярной динамики для мезогенного 2,3,6,7,10,11-гексагептилокситрифенилена. Описываются особенности аппаратно-программных средств, применяемых для проведения вычислительного эксперимента. Показано, что в рамках используемой модели полученные результаты вычислительного эксперимента согласуются с экспериментальными данными поляризационной оптической микроскопии и рентгено-структурного анализа.
Ключевые слова: молекулярная динамика, графические устройства, жидкие кристаллы.
The results on calculation and interpretation of radial distribution functions for 2,3,6,7,10,11-hexaheptyloxytriphenylene obtained by the molecular dynamics simulation are presented. The features of hardware and software resources applied in the calculation experiments are described. Within the framework of the used model the obtained results of the calculation experiments are in a good agreement with the experimental data of polarizing microscopy and X-ray diffraction.
Key words: molecular dynamics, graphics processing units, liquid crystals.
Введение
Жидкие кристаллы (ЖК) в настоящее время находят широкое применение во многих областях деятельности человека (электроника, энергетика, медицина и т. д.). В связи с этим актуальной является проблема получения ЖК, обладающих необходимым для конкретной области набором свойств. В последнее время для выбора стратегии получения ЖК-соединений используются прогностические методы, основанные на построении различных математических моделей [1-5].
Так как соединения, обладающие мезоморфными свойствами, как правило, имеют достаточно сложную структуру, для их описания необходима модель, позволяющая на атомном уровне отслеживать явления, происходящие в системе за наносекундные интервалы времени. Одним из методов, удовлетворяющих данным требованиям, является метод молекулярной динамики (МД), идея которого заключается в численном интегрировании уравнений движения для каждого атома моделируемой
© Москвин Д. О., Соцкий В. В., Кашицын А. С., Усольцева Н. В., 2013
системы [5-12]. При этом для ЖК размер моделируемой системы, требуемый для получения достоверных результатов, составляет десятки тысяч атомов, поэтому МД-моделирование даже в наносекундных интервалах наблюдения за системой обдает большой вычислительной сложностью (в худшем случае вычислительная сложность растет пропорционально N, где N - количество атомов моделируемой системы [6, 8, 12, 13]). В то же время самые затратные вычисления метода МД (расчет сил несвязных взаимодействий) могут быть выполнены независимо (параллельно) для каждого атома моделируемой системы. В связи с этим целесообразно внедрение в МД-моделирование технологий параллельных вычислений, способных значительно ускорить вычислительный эксперимент [13-18].
Графические устройства (GPU), обладая массивно-параллельной архитектурой, на настоящий момент являются наиболее мощными ускорителями вычислительных экспериментов [4, 13-22].
Материалы и методы
В качестве объекта исследования был выбран 2,3,6,7,10,11-гексагептил-окситрифенилен (I) (рис. 1).
OC7H15
C7H15° 1 I) 1
C7H15O T 11 ^Y^OC7Hi5 OC7H15
Рис. 1. Графическая модель молекулы соединения I
По данным ранее выполненных исследований [23], известно, что при охлаждении в температурном интервале от 68,6 до 93,0 °С соединение I находится в мезофазе, тип которой колончатая гексагональная. Кроме того, в [23] приведены результаты рентгеноструктурного анализа соединения I.
Для выполнения вычислительного эксперимента по МД-моделированию исследуемого соединения была использована программа «МББтОпё» [24, 25], выполняющая вычисления на графических устройствах. Потенциалы атом-атомного взаимодействия в программе представляются в следующем виде:
v v v V
U = Zkb (rb - rb0 )2 + ZК (« - «0 f + Zy i1 + cos(^ + %)) +
be B
+
i j * i
12 6 rr
VV У V У
+ k
4i4j
r
У у
где U - потенциальная энергия системы, ^, Vf - константы упругости валентных связей, валентных углов и торсионных углов соответственно, гь - длина связи в текущий момент времени, гЬ0 - равновесная длина связи, B - множество всех валентных связей системы, а - величина угла в текущий момент времени, а0 - равновесный угол, ^ - значение торсионного угла в текущий момент времени, п - кратность торсионного барьера, ф0 - начальная фаза, Ai, Bj - параметры потенциала Ван-дер-ваальсовых взаимодействий, г у - расстояние между ^ым и ^ым атомами, qi, q - их заряды, k - коэффициент пропорциональности.
Интегрирование уравнений движения осуществлялось методом Эйлера. Для поддержания температуры и давления системы применялись термостат и баростат Берендсена [8]. На систему были наложены периодические граничные условия.
Вычислительный эксперимент по МД-моделированию пространственной структуры и фазового состояния соединения I был выполнен по следующей схеме: моделируемая система была составлена из 360 молекул (58320 атомов), внешнее давление составило 1 атм., интервала наблюдения - 10 не. Эксперимент был выполнен при трех различных значения температуры: 97 °С (температура изотропной жидкости), 72 °С (температура мезофазы) и 62 °С (температура кристаллической фазы). По завершении моделирования были получены функции радиального распределения (ФРР), рассчитанные по центрам масс молекул, проанализировав которые мы сможем сделать вывод о пространственной структуре и фазовом состоянии исследуемой системы.
Результаты и их обсуждение
Так как моделируемая система достаточно сложная, то при анализе ФРР, полученных в результате МД-моделирования, возникает ряд трудностей.
Для анализа ФРР нами была разработана модель, позволяющая интерпретировать её для анизотропных систем (получить зависимость количества соседей от расстояния для идеального кристалла с гексагональной упаковкой молекул).
В рамках модели были сделаны следующие допущения: каждая молекула системы была аппроксимирована цилиндром, центр которого совпадает с центром масс молекулы, диаметр - с шириной молекулы, высота - с толщиной молекулы, центр системы координат перенесен в центр молекулы, для которой рассчитывается координационное число, центры масс молекул лежат в сфере некоторого радиуса (центр сферы совпадает с началом координат) (рис. 2).
Очевидно, что центры ближайших цилиндров, лежащих в той же плоскости, что и начало координат (центр цилиндра для которого рассчитывается координационное число), будут располагаться на окружности радиуса г, равного 2^ (где, Яс - радиус цилиндра), от него (т. е. будут лежать на окружности, отсекаемой от сферы данной плоскостью); в противном случае произойдет перекрывание молекул (при г меньше 2Яс) пли разуплотнение гексагональной упаковки (при г больше 2Ис). Центры соседей, лежащих в плоскостях, параллельных плоскости, содержащей начало координат, будут лежать на окружностях с центрами, смещенными относительно начала координат на расстояние, кратное И, в направлении вектора a (рис. 2), а также в центрах этих окружностей.
Рис. 2. Проекция сферы, содержащей центры цилиндров (обозначены точками), на плоскость
При вычислении координационного числа мы будем двигаться от самого нижнего соседа (выколотая точка на рис. 2) до начала координат (где располагается центр цилиндра, для которого мы рассчитываем координационное число) в направлении вектора а с шагом, равным И. В силу симметрии вторую половину сферы можно не обходить (окружности первой половины, центры которых не совпадают с началом координат, учесть дважды, окружности с центром в начале координат - один раз). Математически данную модель можно описать следующим образом:
(1о(Ю$(Ю ¡г(г)&) -1
(2)
где
1о(Ю:
[1, если Я е О I 0, иначе
д(Я) =
11, если Я = 0
12, иначе
О = {Я: - Я + {Я8/Н} < Я < 0 а (Я - (Я8/Н})/Н е 2},
, . т , . , . . Ц если А(р) е Г _ ,
У(г) = £!г(А(Р)), 1г(А(р)) = Г , Г = {А(р): р = г},
[0, иначе
где Яц - радиус сферы, И - высота цилиндра А(р) - центр цилиндра, р - расстояние от центра рассматриваемой окружности до центра цилиндра (обобщенная координата точки А), Я/И} - дробная часть от Я/И.
В формуле (2) единица вычитается для того, чтобы не учитывать цилиндр, расположенный в начале координат (сам у себя он соседом быть не может), д(Я) - для того, чтобы точки, лежащие на окружностях с центром в начале координат, были учтены один раз, точки, лежащие на всех остальных окружностях - два раза.
Анализируя ФРР, полученные в результате МД-моделирования (рис. 3), можно увидеть, что структуры кристаллической фазы и мезофазы не имеют значимых различий в области до 20 А, соответствующей параметрам колонки. На основании этого можно сделать вывод о сохранении структуры колонки при фазовых превращениях.
Рис. 3. ФРР для соединения (I) (здесь для наглядности, соответствующая температуре мезофазы, перемещена на 2 единицы по оси ординат, кривая, соответствующая температуре кристаллической фазы, - на 4)
Пики, соответствующие первой координационной сфере, находятся на расстоянии 3,5 А, что на 0,1 А отличается от значения (3,6 А), полученного с помощью рентгеноструктурных исследований [23] (относительная погрешность составляет 2,9 %). Также на этих ФРР имеются выраженные пики на расстоянии 7 А и 10,5 А, соответствующие второй и третьей координационным сферам. Кроме того, на исследуемых кривых наблюдается пик, отвечающий ширине колонки, на расстоянии 21,9 А, что на 0,3 А меньше величины, приведенной в [23] (22,2 А), однако величина относительной погрешности результата составляет 1,4 %.
Что касается ФРР, отвечающей температуре изотропной жидкости, то пики, соответствующие первой, второй и третьей координационным сферам, а также пик, отвечающий ширине колонки (3,5 А, 7 А, 10,5 А и 21,9 А соответственно), имеют меньшую интенсивность по сравнению с пиками кривых, отвечающих температурам кристаллической фазы и мезофазы. Кроме того, посредством визуального анализа (рис. 4), удалось установить, что при температуре мезофазы молекулы системы имеют колончатую гексагональную упаковку.
б
Рис. 4. Фрагмент системы молекул соединения I, при температуре мезофазы (а) и модельное представление надмолекулярной упаковки (б)
Так как исследуемая система является достаточно сложной, для того чтобы оценить ее пространственную структуру и фазовое состояние, функций радиального распределения недостаточно. Для более детального анализа полученных ФРР при помощи предложенной выше модели (формула (2)) была построена зависимость координационных чисел от расстояния (сплошная линия на рис. 5) для идеального кристалла с гексагональной упаковкой молекул со следующими параметрами: к = 3,5 А, В« = 30 А.
Рис. 5. ФРР, соответствующая температуре кристаллической фазы, полученная по результатам моделирования (%(г)_Сг), и зависимость, построенная при помощи формулы (2) (п(г)_Сг)
Анализируя зависимости, представленные на рис. 5, можно сделать вывод, что максимумы ФРР, полученной в результате МД-моделирования, соответствуют пикам кривой, построенной для идеального кристалла при помощи формулы (2), из чего можно сделать вывод о том, что при данной температуре система действительно находится в кристаллической фазе.
Теперь проанализируем ФРР в дальней области (от 20 до 30 А) (рис. 6), полученные в результате моделирования при различных температурах.
2,00 1,75 1,50 1,25 1,00 0,75 0,50
0,25 -I I
о.оо ..............1.1...1|1,.Н.....1Л1.....1|.........1,1....................
18 20 22 24 26 28 30 32
Г, А
Рис. 6. ФРР, полученные в результате МД-моделирования при различных температурах и зависимость, построенная по формуле (2), в области от 20 до 30 А
Можно заметить, что пики на кривой, соответствующей температуре изотропной жидкости, имеют более размытую форму и смещены вправо относительно пиков на кривых, соответствующих двум другим температурам и кривой, построенной для идеального кристалла. Данный факт свидетельствует об изменении подвижности системы, а, следовательно, - фазовом переходе (изотропная жидкость - мезофаза). Кроме того, в дальней области пики различимые на кривой, отвечающей температуре кристаллической фазы (отчетливо видны 6 пиков), становятся неразличимыми на кривой, отвечающей температуре мезофазы (на ней отчетливо видны только три пика с координационными числами 18, 28, и 34 соответственно), что обусловлено увеличением подвижности системы. Следовательно, здесь также имеет место фазовый переход (мезофаза - кристаллическая фаза).
Выводы
Посредством проведенных в настоящей работе вычислительных экспериментов удалось установить, что результаты молекулярно-динамического моделирования пространственной структуры и фазового состояния 2,3,6,7,10,11-гексагептилокси-трифенилена находятся в хорошем согласии с данными ранее выполненных физических исследований (поляризационной оптической микроскопии и
рентгеноструктурного анализа): отклонение пика, соответствующего первой координационной сфере, составляет 2,9 %, отклонение пика, соответствующего ширине колонки, - 1,4 %, кроме того, совпали температуры существования фаз. Следовательно, данный метод (модель) может быть применен для анализа соединений, представляющих схожие классы. Кроме того, предложена модель интерпретации результатов, полученных в ходе вычислительных экспериментов.
Работа поддержана грантом РФФИ№ 13-03-00481-а.
Список использованной литературы
1. Габдулсадыкова Г. Ф., Жарова М. А., СоцкийВ. В., Усолъцева Н. В. // Жидкие кристаллы и их практическое использование. 2011. Вып. 4. C. 32-40.
2. Акопова О. Б., Пестов С. М. // Жидкие кристаллы и их практическое использование. 2012. Вып. 4. С. 19-33.
3. Martinelli N., Olivier Y., Muccioli L., Minoia A., et al. // Functional Supramolecular Architectures for Organic Electronics and Nanotechnology / edited by F. Cacialli and P. Samori, Wiley-VCH, 2010. P. 3-37.
4. Москвин Д. О., СоцкийВ. В., Усолъцева H. В. // Изв. вузов. Сер. физ. 2012. Т. 55, № 10. С.38-40.
5. СоцкийВ. В., УсолъцеваН. В. Методы моделирования мезоморфных систем. Иваново, 2007. 58 с. Деп. в ВИНИТИ 03.12.2007, № 111б.
6. Alen M. P., Tildesley D. J. Computer simulations of liquids. Oxford : Clarendon press, 1991. 385 p.
7. Demus D., Gooldby J., Gray G.W., Spiess H.-W., Vill V. Hahdbook of liquid crystals. Willey-VCH, 1998. Vol. 1. 914 p.
8. Методы компьютерного моделирования для исследования полимеров и биополимеров / под ред. В. А. Иванова, А. Л. Рабиновича, А. Р. Хохлова. М. : КД Либрком, 2009. б9б с.
9. Аксёнова Е. В., Кшевецкий М. С. Вычислительные методы исследования молекулярной динамики. СПб. : СПбГУ, 2009. 50 с.
10. УсолъцеваН. В., Ясинский Ф. Н., СоцкийВ. В., Костин М. С. // Вести. ИГЭУ. 2009. Вып. 4. С. 45-48.
11. Ибрагимов И. М., Ковшов А. Н., Назаров Ю. Ф. Основы компьютерного моделирования наносистем : учеб. пособие. СПб. : Лань, 2010. 38б с.
12. Клочков С. А., Ясинский Ф. Н. Численные методы молекулярной динамики и нанотехнологии. Иваново, ИГЭУ, 2007. 8б с.
13. СоцкийВ. В. // Жидкие кристаллы и их практическое применение. Иваново : Иван. гос. ун-т, 2011. Вып. 3. С. 77-84.
14. Москвин Д. О., Соцкий В. В. // Тез. конф. фестиваля «Молодая наука в классическом университете». 2012. С. 19.
15. Москвин Д. О., Соцкий В. В. // Тез. конф. фестиваля «Молодая наука в классическом университете». 2013. С. 81-82.
16. Anderson J. A., Lorenz C. D., Travesset A. A. // Journal of Computational Physics. 2008. Vol. 227. P.5342-5359.
17. Боярченков А. С., Поташников С. И. // Вычислительные методы и программирование. 2009. Т. 10. С. 9-23.
18. Giupponi G., Harvey M. J., De Fabritis G. // Drug discovery today. 2008. Vol. 13. P. 10521058.
19. Москвин Д. О., Соцкий В. В., Усолъцева Н. В. // Сборник научных трудов S World по материалам конференции «Современные проблемы и пути их решения в науке, транспорте, производстве и образовании '2011». Одесса, 2011. Т. 38. С. 47-49.
20. Москвин Д. О., Соцкий В. В. // Иваново : Иван. гос. ун-т, 2012. 83 с. Деп. в ВИНИТИ 21.03.2012, № 101—В2012.
21. Москвин Д. О., Соцкий В. В., Усолъцева Н. В. // Вестн. ИГЭУ. Иваново, 2012. Вып. 2. С. 75-79.
22. Москвин Д. О., Соцкий В. В., Усолъцева Н. В. // Вестник ИГЭУ. Иваново, 2013. Вып. 2. С. 71-74.
23. Destrade C., Tinh N. H., Gasparoux H., Malthete J., Levelut A. M. // Mol. Cryst., Liq. Cryst. 1981. Vol. 71. P. 111-135.
24. Свидетельство № 2011613855 об официальной регистрации программы для ЭВМ «Система молекулярно-динамического моделирования MDsimGrid» / В. В. Соцкий. М., 2011.
25. Свидетельство № 2011613854 об официальной регистрации программы для ЭВМ «Программный модуль молекулярно-динамического моделирования MDsimGrid-GPU» / В. В. Соцкий. М., 2011.
Поступила в редакцию 2.12.2013 г.
