Virtual Structural Biology Group » Учебные курсы » Структурная биоинформатика в ВШЭ

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

Традиционные ссылки на полезные ресурсы: 

* Уроки по работе с GROMACS находятся https://tutorials.gromacs.org/ http://www.mdtutorials.com/gmx/.

* Введение в скриптование в Bash http://www.opennet.ru/docs/RUS/bash_scripting_guide .

* Colab заготовка находится тут https://colab.research.google.com/drive/1ztk3rEzEJWLd3mgHS45ElsM-Vj5zLROf?usp=sharing .

* Важно установить Gromacs если вы работаете локально. 

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

* 1. Начнем с того, что подготовим файл координат и файл топологии. 

> #  Вам предоставлен gro файл(http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro) с координатами этана.

> #  Построим  файл топологии et.top для этана, который выгляди примерно так: 

<pre>

#include "/usr/share/gromacs/top/oplsaa.ff/forcefield.itp"

[ moleculetype ]

; Name            nrexcl

et            3

[ atoms ]

;   nr  type  resnr  residue  atom   cgnr     charge       mass

    1   CX      1    ETH      C1      1    -0.189      12.01

    2   CX      1    ETH      C2      2    -0.155      12.01

    3   HX      1    ETH      H1      3     0.0059       1.008

    4   HX      1    ETH      H2      4     0.0059       1.008

    5   HX      1    ETH      H3      5     0.0059       1.008

    6   HX      1    ETH      H4      6     0.0056       1.008

    7   HX      1    ETH      H5      7     0.0056       1.008

    8   HX      1    ETH      H6      8     0.0056       1.008

[ bonds ]

;  ai    aj funct  b0       kb

     1   2   1  

     1   3   1

     1   4   1  

     1   5   1  

........... 

надо дописать

[ angles ]

;  ai    aj    ak funct  phi0   kphi

;around c1

    3     1     4     1  

    4     1     5     1  

    3     1     5     1  

    2     1     3     1  

    2     1     4     1  

    2     1     5     1  

;around c2

........... 

надо дописать 

[ dihedrals ]

;  ai    aj    ak    al funct  

    3    1     2     6      1  

    3    1     2     7      1 

    3    1     2     8      1  

    4    1     2     6      1  

 ........... 

надо дописать 

[ pairs ]

; список атомов 1-4

;  ai    aj funct

   3  6

   3  7

   3  8

   4  6

   4  7

   4  8

   5  6

   5  7

   5  8

[ System ]

; any text here

first one

[ molecules ]

;Name count

 et    1

</pre>

Этот вариант топологии работать не будет, надо изменить типы атомов.

<pre>

%%bash

cat /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp

</pre>

или намек как это сделать программно:

<pre><code class="python">

# загрузим модули

import rdkit

from rdkit import Chem

from rdkit.Chem import AllChem, Draw

from rdkit.Chem import Descriptors

from rdkit import RDConfig

from IPython.display import Image,display

import numpy as np

</code></pre>

<pre><code class="python">

# создадим этан

mol=Chem.MolFromSmiles('CC')

AllChem.Compute2DCoords(mol)

m3d=Chem.AddHs(mol)

Chem.AllChem.EmbedMolecule(m3d)

AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )

</code></pre>

<pre><code class="python">

# и придумаем циклы 

## связи

bonds = m3d.GetBonds()

for i,b in enumerate(bonds):

    print b.GetBeginAtomIdx() , b.GetEndAtomIdx()

## углы    

for b1 in m3d.GetBonds():

    for b2 in m3d.GetBonds():

        if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetIdx() != b2.GetIdx():

            print b1.GetEndAtomIdx() , b1.GetBeginAtomIdx(), b2.GetEndAtomIdx()

## dihedrals            

for b1 in m3d.GetBonds():

    for b2 in m3d.GetBonds():

        for b3 in ....

             .......

</code></pre>            

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

> * http://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp - метод Берендсена для контроля температуры.

> * http://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp - метод "Velocity rescale" для контроля температуры.

> * http://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp - метод Нуза-Хувера для контроля температуры.

> * http://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp - метод стохастической молекулярной динамики.

Начиная с этого момента вы можете написать функции, а можете делать всё вручную. 

* 3 Очень краткое описание программ и типов файлов вы можете найти "здесь":https://kodomo.fbb.msu.ru/wiki/Main/Modelling/BioNanoPrac Сначала надо построить входные файлы для молекулярно-динамического движка mdrun с помощью grompp: 

<pre><code class="python">

gmx grompp -f %s.mdp -c et.gro -p et.top -o et_%s.tpr

# где i: be,vr,nh,an,sd  см. выше список mdp файлов

</code></pre>

*Задать i вне скрипта можно командой export i="be". * 

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

<pre><code class="python"> 

gmx mdrun -deffnm et_%s -v -nt 1

</code></pre>

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

<pre><code class="python">

gmx trjconv -f et_%s.trr -s et_%s.tpr -o et_%s.pdb

</code></pre>

В отчёт занесите ваши наблюдения и предварительные выводы. 

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

<pre><code class="python">

echo  -e "12\n13" | gmx energy -f et_%s.edr -o et_%s_en.xvg -xvg none

</code></pre>

Постройте графики изменения энергий. Рекомендуемый вид это lines. Графики добавьте в отчёт.

Удобно загружать данные с numpy.loadtxt

<pre><code class="python">

import numpy as np

a= np.loadtxt('%s.xvg' % name)

t=a[:,0]

y=... догадайтесь

</code></pre>

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

<pre> 

[ b ]

1 2 

</code></pre>

И запустим утилиту по анализу связей g_bond:

<pre><code class="python">

gmx distance -f et_%s.trr -s et_%s.tpr -o bond_%s.xvg -n b.ndx -xvg none

</code></pre>

Постройте графики распределения длинн связей. Рекомендуемый вид это гистограмма. Графики добавьте в отчёт.

* 8 Учитывая форму распределения Больцмана и все Ваши наблюдения сделайте вывод о том какой из методов позволяет наиболее реалистично поддерживать температуру в системе. Опишите найденные Вами недостатки предложенных алгоритмов. Постройте таблицу зависимости быстродействия от алгоритма.