Работа 7

Материал из md
Перейти к: навигация, поиск
Моделирование и анализ жидких С6-алканов

Скачать файлы для выполнения работы: 07_hexane.tar.gz

Скачать описание работы в pdf-формате: 07_hexane.pdf

Цель

Научиться работать с файлом топологии для органической молекулы. Получить МД-модель жидкого гексана. Получить МД-модель жидкого циклогексана и/или 2,3-диметилбутана, предварительно создав файл топологии.

Задания

  1. Получить МД-модель жидкого гексана (создать исходную конфигурацию, провести релаксацию, провести основное моделирование).
  2. На основе известной химической структуры молекулы создать файл топологии для циклогексана и/или 2,3-диметилбутана. Используя полученную топологию, получить МД-модели указанных жидкостей.
  3. Рассчитать полную и парциальные ФРР для полученных жидкостей (с учетом и без учета внутримолекулярных расстояний). Рассчитать ФРР для центров масс молекул.
  4. Рассчитать плотности и коэффициенты самодиффузии полученных жидкостей.

Вопросы

  1. Объяснить пики на функциях радиального распределения для гексана и циклогексана (2,3-диметилбутана) для:
  • полноатомной g(r), учитывая все атомы системы, включая внутримолекулярные,
  • парциальных gС-С(r) (для атомов углерода), с учетом и без учета внутримолекулярных расстояний,
  • gЦТ(r) для центров тяжести молекул.
  1. Объяснить причину различия (сходства) в плотностях и коэффициентах самодиффузии данных жидкостей.

Общие замечания

Молекулы алканов представляют собой насыщенные цепочки метильных групп (линейные, циклические или ветвистые), которые могут служить фрагментами макромолекул. Умение моделировать жидкости, состоящие из сравнительно простых молекул, является шагом к работе с большими молекулами, полимерами, белками.

Несмотря на химическое сходство молекул алканов, их физические свойства могут существенно различаться. Считается, что это связано со структурой молекул и их конформациями, из-за которых структура твердой и жидкой фазы, а также температуры фазовых переходов, получаются различными.

Однако вопросы вызывает не только различие, но и сходство некоторых свойств. Например, на Рис. 7.1 показана плотность С6-алканов в зависимости от температуры кипения. Видно, что все изомеры гексана имеют практически одинаковую плотность (различие в пределах 3%), хотя форма молекул и температуры кипения различаются существенно. При этом, плотность циклогексана отличается от них почти на 20%. Все это кажется странным. Форма этих молекул изменяется от линейной для нормального гексана (nh) до практически сферической для 2,3-дибутилметана (2,3-dmb) и циклогексана (ch), см. Рис. 7.2. Кроме формы у молекул нет других существенных различий, межмолекулярное взаимодействие у них практически одинаковое.

Fig7.1.png
Рис. 7.1. Плотность С6-алканов в зависимости от температуры кипения (экспериментальные данные). Справа показаны структурные формулы молекул (нормальный гексан, 3-метилпентан, 2-метилпентан, 2,2-диметилбутан, 2,3-диметилбутан, циклогексан).
Fig7.2.png
Рис. 7.2. Ван-дер-Ваальсовские (верхний ряд) и молекулярные (нижний ряд, Rprobe =1.9A) поверхности молекул гексана, 2,3-диметилбутана и циклогексана (слева на право).

Инертные молекулы с формой, близкой к сферической образуют, так называемую, пластическую кристаллическую фазу. Известно, что циклогексан кристаллизуется при 6,5°С в плотнейшую упаковку (ГЦК), характерную для сферически-симметричных атомов и твёрдых шаров. Почему это происходит, если молекула циклогексана все же не шар? (Заметим, что устойчивая кристаллическая фаза циклогексана возникает при – 93°С и имеет моноклинную структуру). Пластические кристаллы характерны тем, что в них молекулы продолжают вращаться. Можно думать, что благодаря этому они располагаются в пространстве подобно шарам. Проявляется ли это свойство компактных молекул в жидкости? Если это так, то должно быть нечто общее в структуре жидкого циклогексана и простых жидкостей. Но то же самое можно предположить и про 2,3-диметлбутан, несмотря на то, что его плотность существенно отличается от плотности циклогексана.

Указания к выполнению работы

Гексан

Моделирование гексана выполняется в директории 07_hexane/nhex/. Подготовка исходных данных и релаксация проводится в поддиректориях preparing, equi1, equi2.

Подготовка начальной конфигурации модели

Для подготовки начальной конфигурации перейдите в поддиректорию preparing/. В этой директории находится файл координат атомов одной молекулы гексана nhex.pdb. Тут же находится файл топологии молекулы гексана nhex.top.

В нём используется поле сил OPLS. В разделе [atoms] указано соответствие между названиями атомов (CA, HA1 и т.д.), используемые в файле молекулярной структуры, с типами атомов поля сил OPLS. Обратим внимание на разделы ковалентно-связанных атомов:

  • [bonds] – пары атомов, связанные ковалентной связью,
  • [angels] – тройки связанных атомов, участвующих в деформационных движениях,
  • [dihedrals] – четвёрки связанных атомов, участвующих в торсионных движениях,
  • [pairs] – пары атомов, разделённых тремя связями, для которых частично включается межатомное Леннард-Джонсовское взаимодействие.

Убедимся, что топология молекулы, представленная в файле, соответствует нормальному гексану, атомы которого пронумерованы как показано на Рис. 7.38.

Fig7.3.jpg
Рис. 7.3. Нумерация атомов молекулы гексана, используемая в файле топологии

Если топология молекулы записана правильно, приступим к дальнейшей работе. Создадим вокруг молекулы кубический модельный бокс такого размера, чтобы его объём примерно соответствовал объёму, приходящемуся на одну молекулу гексана в жидкости при нормальных условиях. Таким боксом является куб со стороной примерно 0.6 нм. Создадим его программой editconf:

editconf -f nhex.pdb -box 0.6 -o nhex_box.gro

После этого с помощью команды genconf создадим модель, содержащую 216 молекул гексана. Для этого размножим наш бокс с одной молекулой, транслируя его по шесть раз вдоль каждой оси (6х6х6):

genconf -f nhex_box.gro -nbox 6 6 6 -o nhex_216.gro

Получившуюся конфигурацию модели можно посмотреть с помощью VMD:

vmd nhex_216.gro

Теперь следует исправить файл топологии. Для этого достаточно установить (в разделе [molecules]) правильное число молекул: 216. После этого скопируем файл топологии в основную директорию для гексана 07_hexane/nhex/ под именем topol.top . Как в предыдущих наших работах, так и в дальнейшем, будем держать файл топологии в основной директории.

Релаксация

Первичную релаксацию проведём в папке equi1/ с мелким шагом, 0,001 ps, в течение 20 пс. Мелкий шаг необходим, поскольку молекула достаточно сложная, и возможны некоторые перекрывания молекул. Обратим внимание, что здесь не нужно делать минимизацию энергии, как было в работе №6 при создании модели раствора воды и метанола. Там бокс с молекулами налагается на бокс с другими молекулами и исключаются только явные перекрывания.

Создадим ссылки на нужные файлы, чтобы использовать их по умолчанию:

ln -s ../preparing/nhex_216.gro conf.gro
ln -s ../topol.top .

Теперь запускаем препроцессинг и моделирование:

grompp
mdrun -nt 4 -nice 11 -v

Обратите внимание, что в данной работе параметр constraints в файле параметров моделирования grompp.mdp определён как all-bonds, что означает, что длины всех связей не меняются, т.е. мы не будем моделировать валентные колебания. Как правило это и не нужно, так как их частоты очень высоки по сравнению с деформационными и торсионными колебаниями.

После завершения расчёта посмотрим изменения потенциальной энергии и плотности системы (для этого используется программа g_energy). Должно получиться нечто похожее на то, что на Рис. 7.4. Обратим внимание, что энергия выходит на стационар намного быстрее, чем плотность.

Fig7.4.png
Рис. 7.4. Изменение потенциальной энергии (слева) и плотности (справа) модели жидкого гексана на первом этапе релаксации.

Посмотрим траекторию с помощью VMD:

vmd conf.gro traj.xtc

На первой конфигурации видно исходное регулярное расположение молекул внутри бокса. Затем расположение молекул становится неупорядоченным. Однако при этом появляются линии, пересекающие весь бокс, Рис. 7.5 (слева). Это вызвано периодическими граничными условиями. Атомы молекул гексана, выходящие за границу бокса, «возникают» с противоположной стороны, и связи к этими атомами тянутся через весь бокс. Такая ситуация типична для молекул, содержащих большое число атомов.

Fig7.5.jpg
Рис. 7.5. Изображение модели гексана с помощью VMD. Слева: одна из конфигураций, где видно, что некоторые молекулы «разорваны» периодическими граничными условиями. Справа: молекулы восстановлены при помощи команды trjconv.

Поэтому прежде, чем визуализировать траекторию молекулярной системы с помощью VMD, рекомендуется восстановить молекулы, выходящие за границу бокса так, чтобы молекула рисовалась целиком со всеми своими атомами с той стороны бокса, к которому ближе её центр масс. Это делается командой trjconv с ключом -pbc mol (см. подсказку по команде trjconv -h):

trjconv -pbc mol

По умолчанию берётся файл траектории traj.xtc, а результат записывается в файл с именем trajout.xtc. Таким образом, теперь мы используем для VMD этот файл:

vmd conf.gro trajout.xtc

Можно убедиться, что теперь молекулы остаются целыми, см. Рис. 7.5. (справа).

Для дальнейшей релаксации модели перейдем в папку equi2/. Там находится файл grompp.mdp с параметрами основного моделирования. В качестве начальной конфигурации используется последняя конфигурация предыдущей релаксации – файл confout.gro из директории equi1. Создадим новые ссылки:

ln -s ../equi1/confout.gro conf.gro
ln -s ../topol.top .

и запустим моделирование:

grompp
mdrun -nt 4 -nice 11 -v

Полученный файл конечной конфигурации confout.gro будет начальной конфигурацией для равновесной МД-траектории.

Основное моделирование

Возвращаемся в основную директорию nhex/, делаем ссылку на только что полученный файл confout.gro:

ln -s equi2/confout.gro conf.gro

Проверяем параметры моделировании (давление Р=1бар, температура Т=300К, число шагов моделирования должно соответствовать 100 пс). Запускаем препроцессинг и моделирование:

grompp
mdrun -nt 8 -nice 11 -v

После выполнения программы в директории nhex/ появится файл с равновесной МД-траекторией traj.xtc и сопутствующие файлы. Чтобы молекулы не выглядели разорванными, сразу применим к траектории команду trjconv:

trjconv -pbc mol

2,3-диметилбутан

Основное моделирование жидкого 2,3-диметилбутана проводится в директории 07_hexane/dmb/. Для подготовки исходных данных и релаксации используются её поддиректории preparing, equi1, equi2.

В директории dmb/preparing находится файл координат атомов этой молекулы в pdb-формате. Файл топологии молекулы 2,3-диметилбутана предлагается создать самим на основе файла для гексана nhex.top, который нужно скопировать из директории nhex/ в директорию dmb/preparing/ под именем dmb.top на Рис. 7.41 показана молекула 2,3-диметилбутана, где каждому атому присвоен номер. Используя эту нумерацию, выписать группы атомов для разделов [bonds], [angels], [dihedrals] и [pairs] и заменить на них соответствующие разделы в файле dmb.top.

Fig7.6.jpg
Рис. 7.6. Нумерация атомов молекулы 2,3-диметилбутана для создания файла топологии.

Следуя процедуре, изложенной для моделирования гексана, получить равновесную молекулярно-динамическую модель 2,3-диметилбутана при Р=1бар, Т=300К на интервале 100 пс. В случае, если первый этап релаксации не пойдёт, попробуйте уменьшить для него шаг интегрирования. Для результирующей траектории не забудьте "склеить" разорванные молекулы на границах периодического бокса.

Циклогексан

Основное моделирование жидкого циклогексана проводится в директории 07_hexane/chex/. Аналогично описанным выше вариантам с н-гексаном и 2,3-диметилбутаном, здесь для подготовки исходных данных и релаксации используются поддиректории preparing, equi1, equi2.

В директории chex/preparing находится файл координат атомов этой молекулы в pdb–формате (координаты подобных молекул обычно можно скачать из какой-нибудь химической базы данных). Файл топологии молекулы циклогексана предлагается создать самим на основе файла для гексана nhex.top, который нужно скопировать из директории nhex/ в директорию chex/preparing под именем chex.top.

На Рис. 7.42 показана молекула циклогексана, где каждому атому присвоен номер. Используя эту нумерацию, выписать группы атомов для разделов [bonds], [angels], [dihedrals] и [pairs] и заменить на них соответствующие разделы в файле chex.top. Обратим внимание, что циклогексан имеет на два атома водорода меньше, чем гексан и 2,3-диметилбутан. Поэтому здесь ещё нужно убрать из файла топологии два атома H.

Fig7.7.jpg
Рис. 7.7. Нумерация атомов молекулы циклогексана для создания файла топологии

Следуя процедуре, изложенной для моделирования гексана, получите в директории chex/ равновесную молекулярно-динамическую модель циклогексана при Р=1бар, Т=300К на интервале 100 пс.

Анализ моделей

Анализ каждой жидкости проводится в своей директории, т.е. в nhex/, dmb/ или chex/.

Расчёт функций радиального распределения:

По умолчанию команда g_rdf с файлом conf.gro после ключа -s приводит к расчёту всех возможных расстояний между всеми атомами системы, включая внутримолекулярные расстояния, и строит по ним полноатомную парную корреляционную функцию. Плотность системы, используемая для нормировки функции g(r), рассчитывается как отношение полного числа атомов в модели к объёму модельного бокса.

Напомним, что производится усреднение по всем конфигурациям модели. Если же файл траектории длинный, то нет смысла использовать все конфигурации, достаточно выбрать «независимые», например, через каждую пикосекунду с помощью ключа -dt 1 (см. подсказку g_rdf -h). Итак, команда

g_rdf -s conf.gro -dt 1 -o rdf-all-intramol.xvg -w

выдаёт парную корреляционную функцию для всех атомов системы, состоящей из молекул, определённых в файле (conf.gro), используя конфигурации молекулярно-динамической траектории traj.xtc (по умолчанию) через одну пикосекунду (-dt 1), записывает полученную функцию под именем rdf-all-intramol.xvg и, благодаря ключу -w сразу изображает кривую на экране монитора, используя стандартную графическую программу xmgrace.

На Рис. 7.8 для иллюстрации показаны такие парные корреляционные функции для гексана, 2,3-демитилбутана и циклогексана. Видно много острых пиков, которые соответствуют внутримолекулярным расстояниям. Некоторые из них очень узкие, они соответствуют фиксированным длинам связей между атомами молекулы. Некоторые пики имеют заметную ширину. Почему?

Обычно атомы водорода исключают при расчёте функций радиального распределения. В этом случае говорят, что учитываются только «тяжелые» атомы. У нас такими атомами являются углероды. Для расчёта такой радиальной функции необходимо сначала создать группу углеродов в индекс-файле при помощи команды make_ndx:

make_ndx -f conf.gro

В данном случае надо выбрать (в интерактивном режиме) атомы с именем, начинающимся на С*. Также ввод информации в интерактивный режим можно осуществить через echo:

echo -e "a C*\n q" | make_ndx -f
Fig7.8.png
Рис. 7.8. Полноатомные парные корреляционные функции g(r) для жидкого гексана, 2,3-демитилбутана и циклогексана (с учетом всех, как межмолекулярных, так и внутримолекулярных расстояний).

Итоговый индекс-файл можно подать на вход g_rdf посредством ключа -n. Тогда функция радиального распределения для углеродов gС-С(r) может быть рассчитана так:

echo "C* C*" | g_rdf -n -dt 1 -o rdf-cc-intramol.xvg -w

На Рис. 7.9 собраны функции радиального распределения углеродов для всех рассматриваемых жидкостей.

Fig7.9.png
Рис. 7.9. Парциальные парные корреляционные функции углеродов gС-С(r) для жидкого гексана, 2,3-демитилбутана и циклогексана (с учетом внутримолекулярных расстояний).

Для исключения внутримолекулярных расстояний необходимо использовать информацию, содержащуюся в файле топологии. Напомним, что параметр nrexcl в разделе [moleculetype] файла topol.top означает число связей между атомами, взаимодействие между которыми считается особым образом. Обычно оно равно трем, что соответствует атомам, удалённым друг от друга не более чем на три связи, для которых рассчитываются валентное, деформационное и торсионное взаимодействия. Говорят, что это число задаёт «исключения» для леннард-джонсовского и кулоновского взаимодействий. Взаимодействие между всеми остальными атомами одной и той же молекулы рассчитывается так же, как между атомами разных молекул. При расчёте функции радиального распределения с файлом топологии параметр nrexcl используется для выявления межмолекулярных расстояний. Поэтому по умолчанию к ним будут отнесены только те атомы которые разделенным не более чем тремя связями. Для малых молекул (типа метанола) этого достаточно, а для наших молекул – мало, так как часть атомов одной молекулы может быть расположена на расстоянии большем, чем три связи. Чтобы исключить все возможные внутримолекулярные расстояния, можно задать параметр nrexcl равным максимальному числу связей, разделяющих атомы внутри данной молекулы (для гексана это 7). Поэтому, специально для расчёта функций радиального распределения, сделаем копию файл топологии под именем topol_rdf.top:

cp topol.top topol_rdf.top

и пропишем в нем в разделе [moleculetype] необходимое количество связей:

[ moleculetype ]
 ; Name            nrexcl
 nHexane              7

Программа g_rdf использует файл topol.tpr, т.е. нужно заново провести препроцессинг с новым файлом топологии:

grompp -p topol_rdf.top -o topol_rdf.tpr

Здесь мы называем новый tpr-файл как topol_rdf.tpr. Теперь можно запустить расчёт радиальной функции распределения без внутримолекулярных расстояний, указав новый tpr-файл после ключа -s.

echo 0 0 | g_rdf -s topol_rdf.tpr -dt 1 -o rdf-all.xvg -w

Соответствующие функции показаны на Рис. 7.10. Получившиеся функции кажутся весьма невыразительными, что типично для молекулярных жидкостей. Даже первый максимум оказывается слабо выраженным. Это понято, так как ближайшие молекулы могут быть ориентированы по-разному, и не дают четкого межмолекулярного расстояния между атомами углерода

Fig7.10.png
Рис. 7.10. Парциальные межмолекулярные парные корреляционные функции gС-С(r) для жидкого гексана, 2,3-демитилбутана и циклогексана.

Представляет интерес расчёт радиального распределения для центров масс молекул. Это делается командой:

g_rdf -f -s -rdf mol_com -o rdf-com.xvg -w

Ключ -rdf mol_com означает, что будут использованы центры масс молекул.

Функции gЦМ(r) для центров тяжестей наших жидкостей показаны на Рис. 7.11. Видно, что кривые для 2,3-демитилбутана и циклогексана имеют вид, характерный для простых жидкостей: высокий первый максимум и четкие, плавно затухающие последующие осцилляции.

Fig7.11.png
Рис. 7.11. Функции радиального распределения для центров масс молекул жидкого н-гексана, циклогексана и 2,3-диметилбутана.

Это показывает, что компактные молекулы (близкие по форме к сферическим), такие как циклогексан и 2,3-диметилбутан, располагаются в жидкости подобно сферическим атомам простой жидкости. Для н-гексана это не так. Вопрос: почему первый максимум у гексана такой широкий?

Расчёт плотности

Плотность можно извлечь из файла ener.edr c помощью команды:

echo "dens" | g_enegry -o density.xvg -w > dens.output

Для разных жидкостей программа выдаст следующую информацию:

nhex:

Statistics over 50001 steps [ 0.0000 through 100.0000 ps ], 1 da-ta sets
All statistics are over 10001 points

Energy                      Average   Err.Est.       RMSD  Tot-Drift
----------------------------------------------------------------
Density                     656.069          2    8.88654    1.17983  (kg/m^3)

dmb:

Energy                      Average   Err.Est.       RMSD  Tot-Drift
----------------------------------------------------------------
Density                     672.488        3.1    10.4125    17.6811  (kg/m^3)

chex:

Energy                      Average   Err.Est.       RMSD  Tot-Drift
----------------------------------------------------------------
Density                     772.072        3.2    9.66881    12.6439  (kg/m^3)

Заметим, что плотности гексана и 2,3-демитилбутана близкие, а у циклогексана заметно выше. Почему?

Расчёт среднеквадратичного смещения центров масс молекул

Для расчёта среднеквадратичного смещения центров масс больших молекул используется команда g_msd с ключом -mol:

echo "0 0" | g_msd -mol -w > self_diffusion.output

nhex:

Selected 0: 'System'
<D> = 5.8292 Std. Dev. = 5.7195 Error = 0.3892
Fitting from 10 to 90 ps
D[    System] 4.0108 (+/– 0.7875) 1e-5 cm^2/s

dmb:

<D> = 5.7621 Std. Dev. = 5.5546 Error = 0.3779
D[    System] 4.3902 (+/– 2.4034) 1e-5 cm^2/s

chex:

<D> = 3.3371 Std. Dev. = 4.5299 Error = 0.3082
D[    System] 2.6148 (+/– 0.1289) 1e-5 cm^2/s
Fig7.12.png
Рис. 7.12. Среднеквадратичное смещение молекул в жидких гексане, циклогексане и 2,3-диметилбутане.

Заметим, что самодиффузия молекул циклогексана заметно ниже, чем у н-гексана и 2,3-демитилбутана. Почему?