Изучение работы методов контроля температуры в GROMACS

В процессе выполнения этого задания мы будем изучать, как реализован контроль температуры в молекулярной динамике на примере GROMACS. Объект исследования - одна молекула этана.

1. Подготовка файлов.

Сначала подготовим файл координат и файл топологии. Дан .gro-файл с 38 молекулами этана. Создадим индекс файл, в котором будет группа из одной молекулы этана: 

make_ndx -f box_38.gro -o 1.ndx
r 1
q

В результате получился файл 1.ndx . Теперь создадим .gro-файл с одной молекулой, зададим ячейку и расположим молекулу в центре этой ячейки: 

editconf -f box_38.gro -o et1.gro -n 1.ndx
3

editconf -f et1.gro -o et.gro -d 2 -c

Нам дал примерный файл топологии этана et_example.top . Однако этот вариант работать не будет, нужно изменить типы атомов. В этом нам поможет файл atomtypes.atp. Выбираем такие типы атомов: 

opls_135   12.01100  ; alkane CH3 
opls_140    1.00800  ; alkane H.

Также в файле в качестве примера не хватает одной связи. Добавим ее. Получается файл et.top. Теперь нужно проверить, что это будет работать (для этого нам понадобится дополнительный файл be.mdp). 

grompp -f be.mdp -c et.gro -p et.top -o et_test.tpr

Все работает, ура!)

Параметры контроля температуры

Даны 5 файлов с разными параметрами контроля температуры:

  • be.mdp - метод Берендсена для контроля температуры;
  • vr.mdp - метод "Velocity rescale" для контроля температуры;
  • nh.mdp - метод Нуза-Хувера для контроля температуры;
  • an.mdp - метод Андерсена для контроля температуры;
  • sd.mdp - метод стохастической молекулярной динамики.

    Сначала строим входные файлы для молекулярно-динамического движка mdrun с помощью grompp: 

     grompp -f ${i}.mdp -c et.gro -p et.top -o et_${i}.tpr
# где i: be,vr,nh,an,sd

    В результате получилось 5 .tpr-файлов. Теперь для каждого из них запустим mdrun: 

    mdrun -deffnm et_${i} -v -nt 1

    2. Анализ результатов

    Сначала проводим визуальный анализ. Для этого для каждой из 5 систем проводим конвертацию в pdb: 

    trjconv -f et_${i}.trr -s et_${i}.tpr -o et_${i}.pdb
2

    В первом состоянии все молекулы полностью накладываются друг на друга. Молекулы этана находятся в заторможенной конформации.
    Однако на видео уже заметны различия в различных параметрах контроля температуры.

    Видео 1. Движения молекулы этана, полученные методом Берендсена для контроля температуры

    Видео 2. Движения молекулы этана, полученные методом "Velocity rescale" для контроля температуры

    Видео 3. Движения молекулы этана, полученные методом Нуза-Хувера для контроля температуры

    Видео 4. Движения молекулы этана, полученные методом Андерсена для контроля температуры

    Видео 5. Движения молекулы этана, полученные методом стохастической молекулярной динамики для контроля температуры.

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

    Рисунок 1. Молекула этана в заслоненной конформации.

    Теперь сравним потенциальную энергию связи и кинетическую энергию для каждой из 5 систем: 

    g_energy -f et_${i}.edr -o et_${i}_en.xvg
10
11

    С помощью Gnuplot построим графики изменения энергий в виде dot-plot (рисунки 2-6):

    Рисунок 2. График изменения энергий в системе с методом Берендсена для контроля температуры.

    Рисунок 3. График изменения энергий в системе с методом "Velocity rescale" для контроля температуры.

    Рисунок 4. График изменения энергий в системе с методом Нуза-Хувера для контроля температуры.

    Рисунок 5. График изменения энергий в системе с методом Андерсена для контроля температуры.

    Рисунок 6. График изменения энергий в системе с методом стохастической молекулярной динамики.

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

    Далее рассмотрим распределение длинны связи С-С за время моделирования. Сначала создадим индекс файл с одной связью b.ndx. Затем запустим утилиту по анализу связей g_bond: 

    g_bond -f et_${i}.trr -s et_${i}.tpr -o bond_${i}.xvg -n b.ndx

    В Gnuplot построим гистограммы распределения длинн связей (рисунки 7-11):

    Рисунок 7. Распределение длинн связей в системе с методом Берендсена для контроля температуры.

    Рисунок 8. Распределение длинн связей в системе с методом "Velocity rescale" для контроля температуры.

    Рисунок 9. Распределение длинн связей в системе с методом Нуза-Хувера для контроля температуры.

    Рисунок 10. Распределение длинн связей в системе с методом Андерсена для контроля температуры.

    Рисунок 11. Распределение длинн связей в системе с методом стохастической молекулярной динамики.

    На распределение Больцмана больше всего похожи графики для методов стохастической молекулярной динамики и "Velocity rescale". Графики для методов Берендсена и Нуза-Хувера больше похожи на нормальное распределение. Для графика метода Андерсена сложно сходу подобрать распределение.

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