Будем изучать как реализован контроль температуры в молекулярной динамике на примере GROMACS. Объект исследования это одна молекула этана. А начнем с того, что подготовим файл координат и файл топологии.
make_ndx -f box_38.gro -o 1.ndx editconf -f box_38.gro -o et1.gro -n 1.ndx editconf -f et1.gro -o et.gro -d 2 -cВ итоге получился следующий файл: et.gro. А файл топологии был исправлен: et.top. После этого были скачены для работы следующие файлы:
be.mdp - метод Берендсена для контроля температуры. vr.mdp - метод "Velocity rescale" для контроля температуры. nh.mdp - метод Нуза-Хувера для контроля температуры. an.mdp - метод Андерсена для контроля температуры. sd.mdp - метод стохастической молекулярной динамики.Для получения файлов на вход молекулярно-динамическому движку mdrun с помощью grompp был использован следующий скрипт:
#!/bin/bash
for i in {be,vr,nh,an,sd}; do
grompp -f ${i}.mdp -c et.gro -p et.top -o et_${i}.tpr
done
Получилось, как и ожидалось, 5 trp файлов. Подаем их на вход mdrun, добавив в скрипт следующую строчку:
mdrun -deffnm et_${i} -v -nt 1
Теперь переходим к анализу результатов. Начнем с визуального анализа. Для каждой из 5 систем проведем конвертацию в pdb и просмотрим в PyMol.
Добавим в скрипт такую строчку:
trjconv -f et_${i}.trr -s et_${i}.tpr -o et_${i}.pdb
А теперь опишем полученные pdb файлы!
et_an.pdb - молекула колеблется возле одного положения, нет вращений, есть слабые колебания длин связей и углов. et_be.pdb - в самом начале колеблятся только длины связей и углов. et_nh.pdb - тут наблюдается вращение метильных групп относительно С-С связи, а также небольшие колебания валентных углов. et_sd.pdb - молекула меняет конформацию и положение в пространстве (довольно хаотично). et_vr.pdb - наблюдаются смены конформаций и колебания длин связей и валентных углов; молекула большую часть времени находится в заторможенной конформации.Далее, надо сравнить потенциальную энергию связи и кинетическую энергию для каждой из 5 систем. Добавим в скрипт следующую строчку:
g_energy -f et_${i}.edr -o et_${i}_en.xvg
А теперь построим графики изменения энергий с помощью gnuplot. Для этого был создан файл (для каждой из систем) такого содержания:
set terminal png set output "et_vr_en.png" plot "./et_vr_en.xvg" using 1:2, "./et_vr_en.xvg" using 1:3Рассмотрим распределение длинны связи С-С за время моделирования. Сначала создадим индекс файл с одной связью. Для этого в текстовом редакторе создадим файл b.ndx со следующим содержимым:
[ b ] 1 2Запустим утилиту по анализу связей g_bond (добавив соответствующую строку в скрипт):
g_bond -f et_${i}.trr -s et_${i}.tpr -o bond_${i}.xvg -n b.ndx
Ну и построим графики распределения длинн связей (boxes в gnuplot).
Привожу все картинки (что за файлы и откуда, написано на графиках):
