Задание 1
Начнем с того, что подготовим файл координат и файл топологии. На прошлом занятии нам был предоставлен gro файл с 38 молекулами этана.Создадим индекс файл, в котором будет группа из одной молекулы этана.
make_ndx -f box_38.gro -o 1.ndx
После запуска команды у нас появилось приглашение к вводу. Сначала мы ознакомились с помощью, нажав "h" + enter.
Далее выбрали остаток номер 1. Нажали enter и увидели, что появилась новая группа. Создали gro файл с одной молекулой и задали ячейку . При запуске ediconf выбрали номер, соответствующий группе из одной молекулы.
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 из прошлого задания. В разделе [ molecules ] просто изменили количество молекул этана (исправили 38 на 1).
Задание 2
Нам были даны 5 файлов с разными параметрами контроля температуры:be.mdp - метод Берендсена для контроля температуры.
vr.mdp - метод "Velocity rescale" для контроля температуры.
nh.mdp - метод Нуза-Хувера для контроля температуры.
an.mdp - метод Андерсена для контроля температуры.
sd.mdp - метод стохастической молекулярной динамики.
Можно было написать скрипт по работе с 5ю системами, но мы предпочли делать все вручную.=)
Задание 3
Построим входные файлы для молекулярно-динамического движка mdrun с помощью grompp:
grompp -f ${i}.mdp -c et.gro -p et.top -o et_${i}.tpr
# где i: be,vr,nh,an,sd - см. выше список mdp файлов
Задание 4
У нас получилось 5 tpr файлов. Теперь для каждого из них запустим mdrun.mdrun -deffnm et_${i} -v -nt 1
Задание 5
Перешли к анализу результатов. Решили начать работу с визуального анализа. Для каждой из 5 систем провели конвертацию в pdb, посмотрев полученные файлы с помощью PyMol.trjconv -f et_${i}.trr -s et_${i}.tpr -o et_${i}.pdb
Полученные ролики для PyMOL:
Некоторые наблюдения:
В методе Андерсена видим небольшие колебания по длинам связей и валентным углам. Молекула при этом не вращается. Возможно, такое состояние наблюдается при достаточно низких температурах. Должно быть, такие молекулы образуют некую упорядоченную структуру, для которой не характерны существенные изменения положения составляющих ее частей.
В методе Берендсена молекула колеблется и вращается достаточно медленно, но постепенно скорость увеличивается, молекула начинает все больше вращаться, амплитуда колебаний уменьшается. Это, конечно, интересный подход, но вот существование систем из подобных единиц порождает вопросы.
В методе Нуза-Хувера длина связи С-С не изменяется. Наблюдаем лишь активное вращение. В основном, молекула находится в более выгодной заторможенной конформации, но можно увидеть и заслоненную.
В методе стохастической молекулярной динамики молекула как только не крутится и не перемещается. Проследить какой-то алгоритм ее движения невозможно, по крайней мере, пока она одна. Может быть, это неплохо будет моделировать модель случайных столкновений.
Метод "Velocity rescale" похож на метод Нуза-Хувера. Отличие состоит в том, что здесь амплитуды колебаний больше, а вращения по связи меньше. Невыгодные конформации, близкие к заслоненной, видим достаточно редко.
Задание 6
Сравнили потенциальную энергию связи и кинетическую энергию для каждой из 5 систем.g_energy -f et_${i}.edr -o et_${i}_en.xvg
Построили графики изменения энергий. Рекомендуемый вид - dot-plot.
При использовании Gnuplot для построения графиков перешли в рабочую директорию и запустили Gnuplot:
set datafile commentschars "#@&"
plot "./et_be_en.xvg" using 1:2, "./et_be_en.xvg" using 1:3
....
plot "./et_sd_en.xvg" using 1:2, "./et_sd_en.xvg" using 1:3
Получили следущие графики (кинетическая энергия - зеленая; потенциальная энергия связей - красная):
в методе Андерсена значения потенциальной и кинетической энергий примерно одинаковые и очень маленькие. Как-то такой термостат не очень подходит для того, чтобы сообщить молекуле необходимую энергию для того, чтобы она попала в глобальный минимум;
в методе Берендсена потенциальная энергия достаточно быстро падает, а кинетическая выходит на плато. Значения кинетической энергии выше, чем в методе Андерсена, но все равно еще достаточно маленькие;
в методе Нуза-Хувера можно увидеть единичные достаточно высокие значения кинетической энергии. Это нам уже больше подходит;
в методе стохастической молекулярной динамики видим средние, примерно постоянные, значения кинетической энергии;
в методе "Velocity rescale" ситуация примерно такая же, как и в методе стохастической молекулярной динамики.
Задание 7
Рассмотрели распределение длины связи С-С за время моделирования. Сначала создали для этого индекс-файл с одной связью. В текстовом редакторе создали файл 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.
Получили следущие графики:
В методе Андерсена
В методе Берендсена
В методе Нуза-Хувера
В методе стохастической молекулярной динамики
В методе "Velocity rescale"
Задание 8
Форма распределения Больцмана:
Этому распределению совершенно не соответствуют графики методов Андерсена и Берендсена, неплохо соответствуют графики методов стохастической молекулярной динамики и "Velocity rescale".
Исходя из всех наблюдений, можно сказать, что наиболее вероятны 2 метода: стохастической молекулярной динамики и "Velocity rescale". На мой взгляд, метод "Velocity rescale" позволяет наиболее реалистично поддерживать температуру в системе, сложно что-то сказать о методе стохастической молекулярной динамики, поскольку в Pymol мы видели что-то не вполне понятное.