Семестр VIII. Молекулярное моделирование.

Спецкурс по выбору

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

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

Рабочая директория - /home/students/y11/ janemoiseeva/Term8/Pr5.

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

Подготовительный этап

Подготовка файла координат и файла топологии

Имеется .gro файл с 38 молекулами этана.

first one
  304
    1ETH     C1    1   0.577   0.217   0.574
    1ETH     C2    2   0.680   0.252   0.467
    1ETH     H1    3   0.478   0.241   0.538
    1ETH     H2    4   0.597   0.274   0.664
    1ETH     H3    5   0.583   0.111   0.597
    1ETH     H4    6   0.676   0.358   0.445
    1ETH     H5    7   0.660   0.195   0.377
    1ETH     H6    8   0.780   0.228   0.504
    2ETH     C1    9   1.350   0.305   0.750
    2ETH     C2   10   1.410   0.333   0.888
    2ETH     H1   11   1.276   0.227   0.758
    2ETH     H2   12   1.428   0.275   0.681
	...

Создадим индекс файл котором будет группа из одной молекулы этана.

%%bash
## создание новой группы (из одной молекулы этана)
make_ndx -f box_38.gro -o 1.ndx
r 1
q
## запись новой группы в отдельный файл
editconf -f  box_38.gro -o et1.gro -n 1.ndx
#№ задание ячейки и расположение молекулы по центру ячейки
editconf -f et1.gro -o et.gro -d 2 -c

Теперь построим файл топологии et.top для этана, который выгляди примерно так et.top. Но этот вариант топологии не рабочий, надо изменить типы атомов, которые надо взять из файла на kodomo. Кроме того (как выяснилось уже при запуске mdrun) есть еще ошибка в разделе [ dihedrals ], а так же в описании связей в молекуле - не хватает строчки (исправим это сразу в файле et.top):

[ bonds ]
;  ai    aj funct  b0       kb
	...
	1   3   1
	...	 
[ angles ]

Для проверки правильности догадок используем команду:

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

Методом научного тыка, выяснилось, что подходит (grompp не ругается) несколько сочетаний разных типов С и Н атомов (см. файлы et_correct_*_*.top в рабочей директории). Методом здравого смысла из них было выбрано сочетание et_correct_7_11.top :

 opls_135   12.01100  ; alkane CH3 
 opls_140    1.00800  ; alkane H.
 

Итак, файл топологии: et_correct.top, файл координат молекулы этана et.gro

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

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

Скачаем файлы с параметрами:

for i in $(echo be vr nh an sd); do
    wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/$i.mdp -O $i.mdp
done

Входные файлы для mdrun (с помощью grompp)

%%bash
for i in $(echo be vr nh an sd); do
    grompp -f $i.mdp -c et.gro -p et_correct.top -o et_$i.tpr
    mdrun -deffnm et_$i -v -nt 1;
done

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

Визуальный анализ

Для визуального анализа необходима конвертация в .pdb

%%bash
for i in $(echo be vr nh an sd); do
trjconv -f et_${i}.trr -s et_${i}.tpr -o et_${i}.pdb
done

for i in $(echo be vr nh an sd); do
ffmpeg -i et_${i}_states/et_${i}_state0%03d.png -b:v 5000000 -c:v libx264 et_${i}_states.mp4
done

Все команды работа в PyMol тут.

Различия можно заметить между положениями молекулы в разное время (states в PyMol). Например, в методе Нуза-Хувера (nh)) в системе изменяется только длина одной C — H связи, а так же смена конформации (из заторможенной в заслоненную), в то время как другие методы изменяют как длины связей, так и углы молекулы, но не меняют конформацию (см. видео).

Видео 1. Движения молекулы этана, полученные методом Андерсена для контроля температуры.	Видео 2. Движения молекулы этана, полученные методом Берендсена для контроля температуры.
Видео 3. Движения молекулы этана, полученные методом Нуза-Хувера для контроля температуры.	Видео 4. Движения молекулы этана, полученные методом стохастической молекулярной динамики для контроля температуры.
Видео 5. Движения молекулы этана, полученные методом "Velocity rescale" для контроля температуры.

Сравнение потенциальной энергии связи и кинетической энергии всех систем

for i in $(echo be vr nh an sd); do
g_energy -f et_${i}.edr -o et_${i}_en.xvg
>10
>11
done

Графики изменения энергий посмтроим в Gnuplot в виде dot-plot.

%%gnuplot
set term png;
set xrange [-5:255];
set yrange [1:1000];
set logscale y 10;
set xlabel "time, ps";
set ylabel "log10 Energy";

set output "et_an_en_log.png";
plot "et_an_en.xvg" using 1:2 title "Potential energy, kJ/mol" lt rgb "#48D1CC" pointtype 7 pointsize 0.6, "et_an_en.xvg" using 1:3 title "Kinetic energy, kJ/mol" lt rgb "#BA55D3" pointtype 8 pointsize 0.6;

Как видно из графиков, изменение энергии в различных методах происходит очень по-разному.

Распределение длинны связи С-С за время моделирования

echo -e '[ b ]\n1 2' > b.ndx ##индекс файл с одной связью

for i in $(echo be vr nh an sd); do
g_bond -f et_${i}.trr -s et_${i}.tpr -o bond_${i}.xvg -n b.ndx ##утилита по анализу связей
done

Графики распределения длинн связей посмтроим в Gnuplot в виде histogram.

set term png;
set xlabel "bond length, nm";
set ylabel "quantity";
set xrange [0.14:0.17];
set yrange [0:400];

set output "et_an_bond.png";
plot "bond_an.xvg" using 1:2 title "An" with boxes lc rgb "#3CB371";

set output "et_be_bond.png";
plot "bond_be.xvg" using 1:2 title "Be" w boxes lc rgb "#3CB371";

set output "et_nh_bond.png";
plot "bond_nh.xvg" using 1:2 title "NH" w boxes lc rgb "#3CB371";

set output "et_sd_bond.png";
plot "bond_sd.xvg" using 1:2 title "SD" w boxes lc rgb "#3CB371";

set output "et_vr_bond.png";
plot "bond_vr.xvg" using 1:2 title "VR" w boxes lc rgb "#3CB371";
-->

Лучше всего распределением Больцмана описываюся методы "Velocity rescale" и стохастической молекулярной динамики, распределения для методов Нуза-Хувера и Берендсена имеют слишком острый пик, т.е. длина связи меняется слабо, метод Андерсена не описывается распределением Больцмана вообще.