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

  1. Подготовим файл координат и файл топологии. В прошлом занятии был предоставлен gro файл с 38 молекулами этана. Создадим индекс файл в котором будет группа из одной молекулы этана. 

    make_ndx -f box_38.gro -o 1.ndx
    Выбераем остаток номер 1. Появиляется новая группа.
    Теперь создадим gro файл с одной молекулой и зададим ячейку . При запуске editconf выбераем номер, соответствующей группе из одной молекулы. 
    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. Создаём отдельную папку для каждого метода. 

    mkdir ${i}
    # где i: be,vr,nh,an,sd
    Сначала надо построить входные файлы для молекулярно-динамического движка mdrun с помощью grompp: 
    grompp -f mdp/${i}.mdp -c et.gro -p et.top -o ${i}/et_${i}.tpr

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

    mdrun_d -deffnm ${i}/et_${i} -v -nt 1

  5. Теперь переходим к анализу результатов. Начнем с визуального анализа. Для каждой из 5 систем проведите конвертацию в pdb и просмотрите в PyMol. 

    trjconv -f ${i}/et_${i}.trr -s ${i}/et_${i}.tpr -o ${i}/et_${i}.pdb
    Полученные ролики для PyMOL:
    et_an.pdb et_be.pdb et_nh.pdb et_sd.pdb et_vr.pdb

      Некоторые наблюдения:
    • В методе Андерсена видны небольшие небольшие колебания по длинам связей и валентным углам и совершенно нет вращения. Такое состояние, наверное, характерно для низких температур, то есть для кристалла.
    • В методе Берендсена вначале видны всевозможные амплитудные колебания и вращения, но постепенно молекула начинает всё больше вращаться, а амплитуда колебаний уменьшается. Это выглядит совсем нереалистично.
    • В методе Нуза-Хувера движение кажется реалистичным. Небольшие колебания и вращение. Молекула в основном находится в заторможенной конформации, довольно часто видна и заслонённая.
    • В методе стохастической молекулярной динамике молекула очень быстро меняет своё положение. Похоже на вращение, но трудно проследить какое-то направление, будто она случайным образом меняет свою ориентацию или сталкивается с другими молекулами.
    • Метод "Velocity rescale" похож на метод Нуза-Хувера. Отличие в том, что здесь более амплитудные колебания и меньше вращения по связи. Конформации, близкие к заслонённой, прослеживаются редко.

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

    echo -e "Bond\nKinetic-En.\n0" | g_energy -f ${i}/et_${i}.edr -o ${i}/et_${i}_en.xvg > ${i}/et_${i}.txt
    Построим графики изменения энергий.
    Создадим файлы скриптов для Gnuplot и выполним их (сохраним графики в формате PNG): 
    echo "set datafile commentschars '#@&' 
set term 'png'
set output '${i}/en_${i}.png'
set pointsize 0.5
plot '${i}/et_${i}_en.xvg' using 1:2,  '${i}/et_${i}_en.xvg' using 1:3
" > ${i}/gnu_${i}.gnu

#запуск gnuplot
cat ${i}/gnu_${i}.gnu | gnuplot

  7. Рассмотрим распределение длинны связи С-С за время моделирования. Сначала создадим индекс файл с одной связью. В текстовом редакторе создаём файл b.ndx со следующим содержимым: 

    [ b ]
1 2
    И запустим утилиту по анализу связей g_bond: 
    g_bond -f ${i}/et_${i}.trr -s ${i}/et_${i}.tpr -o ${i}/bond_${i}.xvg -n b.ndx
    Построим графики распределения длинн связей. Сделаем отображение в Gnuplot в виде boxes.
    echo "set datafile commentschars '#@&' 
set term 'png'
set output '${i}/bond_${i}.png'
plot '${i}/bond_${i}.xvg'with boxes
" > ${i}/gnu_bond_${i}.gnu

cat ${i}/gnu_bond_${i}.gnu | gnuplot\
    Создадим список наших изображений и после работы цикла сольём их в одно. 
    ls -r ${i}/*.png >> image.list
...
montage +frame +shadow +label -tile 2x5 -geometry 640x480 `cat image.list` images.png

    Полученные изображения - слева энергия для каждого состояния (красным цветом потенциальная энергия связей, зелёным - кинетическая энергия), справа распределения длин связей.

    AN
    BE
    NH
    SD
    VR

  8. Распределение Максвелла-Больцмана имеет следующую форму:

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

    Что ж, чёткий пик плотности в небольшом отдалении от нуля не наблюдается, распределению Максвелла-Больцмана не соответствует.

    Исходя из всех наблюдений, можно сказать, что реалистичны 2 метода: стохастической молекулярной динамики и "Velocity rescale". Но на роликах в PyMOL первый кажется немного странным, поэтому более реалистичен, наверное, метод "Velocity rescale".

Семестр 6
На главную

© Айдарханов Руслан 2008