Task8a

Version 7 (Andrey Golovin, 07.11.2022 18:39)

1 1 Andrey Golovin
h1. Изучение работы методов контроля температуры в молекулярной динамике
2 1 Andrey Golovin
3 1 Andrey Golovin
4 1 Andrey Golovin
Традиционные ссылки на полезные ресурсы: 
5 1 Andrey Golovin
6 3 Andrey Golovin
* Уроки по работе с GROMACS находятся http://www.gromacs.org/Documentation/Tutorials .
7 3 Andrey Golovin
* Введение в скриптовании в Bash http://www.opennet.ru/docs/RUS/bash_scripting_guide .
8 5 Andrey Golovin
* Colab заготовка находится тут https://colab.research.google.com/drive/1ztk3rEzEJWLd3mgHS45ElsM-Vj5zLROf?usp=sharing .
9 5 Andrey Golovin
* Важно установить Gromacs если вы работаете локально. 
10 1 Andrey Golovin
11 1 Andrey Golovin
12 1 Andrey Golovin
*Сегодня мы будем изучать как реализован контроль температуры в молекулярной динамике на примере GROMACS. Объект исследования это одна молекула этана.*
13 1 Andrey Golovin
14 1 Andrey Golovin
15 1 Andrey Golovin
16 2 Andrey Golovin
* 1. Начнем с того, что подготовим файл координат и файл топологии. 
17 1 Andrey Golovin
18 1 Andrey Golovin
> #  Вам предоставлен gro файл(http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro) с координатами этана.
19 1 Andrey Golovin
> #  Построим  файл топологии et.top для этана, который выгляди примерно так: 
20 1 Andrey Golovin
 
21 1 Andrey Golovin
<pre>
22 1 Andrey Golovin
23 1 Andrey Golovin
#include "/usr/share/gromacs/top/oplsaa.ff/forcefield.itp"
24 1 Andrey Golovin
25 1 Andrey Golovin
[ moleculetype ]
26 1 Andrey Golovin
; Name            nrexcl
27 1 Andrey Golovin
et            3
28 1 Andrey Golovin
29 1 Andrey Golovin
[ atoms ]
30 1 Andrey Golovin
;   nr  type  resnr  residue  atom   cgnr     charge       mass
31 1 Andrey Golovin
    1   CX      1    ETH      C1      1    -0.189      12.01
32 1 Andrey Golovin
    2   CX      1    ETH      C2      2    -0.155      12.01
33 1 Andrey Golovin
    3   HX      1    ETH      H1      3     0.0059       1.008
34 1 Andrey Golovin
    4   HX      1    ETH      H2      4     0.0059       1.008
35 1 Andrey Golovin
    5   HX      1    ETH      H3      5     0.0059       1.008
36 1 Andrey Golovin
    6   HX      1    ETH      H4      6     0.0056       1.008
37 1 Andrey Golovin
    7   HX      1    ETH      H5      7     0.0056       1.008
38 1 Andrey Golovin
    8   HX      1    ETH      H6      8     0.0056       1.008
39 1 Andrey Golovin
    
40 1 Andrey Golovin
[ bonds ]
41 1 Andrey Golovin
;  ai    aj funct  b0       kb
42 1 Andrey Golovin
     1   2   1  
43 1 Andrey Golovin
     1   3   1
44 1 Andrey Golovin
     1   4   1  
45 1 Andrey Golovin
     1   5   1  
46 1 Andrey Golovin
........... 
47 1 Andrey Golovin
надо дописать
48 1 Andrey Golovin
49 1 Andrey Golovin
[ angles ]
50 1 Andrey Golovin
;  ai    aj    ak funct  phi0   kphi
51 1 Andrey Golovin
;around c1
52 1 Andrey Golovin
    3     1     4     1  
53 1 Andrey Golovin
    4     1     5     1  
54 1 Andrey Golovin
    3     1     5     1  
55 1 Andrey Golovin
    2     1     3     1  
56 1 Andrey Golovin
    2     1     4     1  
57 1 Andrey Golovin
    2     1     5     1  
58 1 Andrey Golovin
;around c2
59 1 Andrey Golovin
........... 
60 1 Andrey Golovin
надо дописать 
61 1 Andrey Golovin
62 1 Andrey Golovin
[ dihedrals ]
63 1 Andrey Golovin
;  ai    aj    ak    al funct  
64 1 Andrey Golovin
    3    1     2     6      1  
65 1 Andrey Golovin
    3    1     2     7      1 
66 1 Andrey Golovin
    3    1     2     8      1  
67 1 Andrey Golovin
    4    1     2     6      1  
68 1 Andrey Golovin
 ........... 
69 1 Andrey Golovin
надо дописать 
70 1 Andrey Golovin
 
71 1 Andrey Golovin
[ pairs ]
72 1 Andrey Golovin
; список атомов 1-4
73 1 Andrey Golovin
;  ai    aj funct
74 1 Andrey Golovin
   3  6
75 1 Andrey Golovin
   3  7
76 1 Andrey Golovin
   3  8
77 1 Andrey Golovin
   4  6
78 1 Andrey Golovin
   4  7
79 1 Andrey Golovin
   4  8
80 1 Andrey Golovin
   5  6
81 1 Andrey Golovin
   5  7
82 1 Andrey Golovin
   5  8
83 1 Andrey Golovin
84 1 Andrey Golovin
[ System ]
85 1 Andrey Golovin
; any text here
86 1 Andrey Golovin
first one
87 1 Andrey Golovin
[ molecules ]
88 1 Andrey Golovin
;Name count
89 1 Andrey Golovin
 et    1
90 1 Andrey Golovin
91 1 Andrey Golovin
</pre>
92 1 Andrey Golovin
93 1 Andrey Golovin
Этот вариант топологии работать не будет, надо изменить типы атомов.
94 1 Andrey Golovin
95 1 Andrey Golovin
<pre>
96 1 Andrey Golovin
%%bash
97 1 Andrey Golovin
cat /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp
98 1 Andrey Golovin
</pre>
99 1 Andrey Golovin
100 1 Andrey Golovin
или намек как это сделать программно:
101 1 Andrey Golovin
<pre><code class="python">
102 1 Andrey Golovin
# загрузим модули
103 1 Andrey Golovin
import rdkit
104 1 Andrey Golovin
from rdkit import Chem
105 1 Andrey Golovin
from rdkit.Chem import AllChem, Draw
106 1 Andrey Golovin
from rdkit.Chem import Descriptors
107 1 Andrey Golovin
from rdkit import RDConfig
108 1 Andrey Golovin
from IPython.display import Image,display
109 1 Andrey Golovin
import numpy as np
110 1 Andrey Golovin
</code></pre>
111 1 Andrey Golovin
112 1 Andrey Golovin
<pre><code class="python">
113 1 Andrey Golovin
# создадим этан
114 1 Andrey Golovin
mol=Chem.MolFromSmiles('CC')
115 1 Andrey Golovin
AllChem.Compute2DCoords(mol)
116 1 Andrey Golovin
m3d=Chem.AddHs(mol)
117 1 Andrey Golovin
Chem.AllChem.EmbedMolecule(m3d)
118 1 Andrey Golovin
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
119 1 Andrey Golovin
</code></pre>
120 1 Andrey Golovin
121 1 Andrey Golovin
<pre><code class="python">
122 1 Andrey Golovin
# и придумаем циклы 
123 1 Andrey Golovin
## связи
124 1 Andrey Golovin
bonds = m3d.GetBonds()
125 1 Andrey Golovin
for i,b in enumerate(bonds):
126 1 Andrey Golovin
    print b.GetBeginAtomIdx() , b.GetEndAtomIdx()
127 1 Andrey Golovin
128 1 Andrey Golovin
## углы    
129 1 Andrey Golovin
for b1 in m3d.GetBonds():
130 1 Andrey Golovin
    for b2 in m3d.GetBonds():
131 1 Andrey Golovin
        if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetIdx() != b2.GetIdx():
132 1 Andrey Golovin
            print b1.GetEndAtomIdx() , b1.GetBeginAtomIdx(), b2.GetEndAtomIdx()
133 1 Andrey Golovin
## dihedrals            
134 1 Andrey Golovin
for b1 in m3d.GetBonds():
135 1 Andrey Golovin
    for b2 in m3d.GetBonds():
136 1 Andrey Golovin
        for b3 in ....
137 1 Andrey Golovin
             .......
138 1 Andrey Golovin
</code></pre>            
139 1 Andrey Golovin
140 1 Andrey Golovin
141 2 Andrey Golovin
* 2 Вам даны 5 файлов с разными параметрами контроля температуры:
142 2 Andrey Golovin
> * http://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp - метод Берендсена для контроля температуры.
143 2 Andrey Golovin
> * http://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp - метод "Velocity rescale" для контроля температуры.
144 2 Andrey Golovin
> * http://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp - метод Нуза-Хувера для контроля температуры.
145 2 Andrey Golovin
> * http://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp - метод стохастической молекулярной динамики.
146 1 Andrey Golovin
147 1 Andrey Golovin
Начиная с этого момента вы можете написать функции, а можете делать всё вручную. 
148 2 Andrey Golovin
  
149 4 Student HSE
* 3 Очень краткое описание программ и типов файлов вы можете найти "здесь":https://kodomo.fbb.msu.ru/wiki/Main/Modelling/BioNanoPrac Сначала надо построить входные файлы для молекулярно-динамического движка mdrun с помощью grompp: 
150 1 Andrey Golovin
151 1 Andrey Golovin
<pre><code class="python">
152 6 Andrey Golovin
gmx grompp -f %s.mdp -c et.gro -p et.top -o et_%s.tpr
153 1 Andrey Golovin
# где i: be,vr,nh,an,sd  см. выше список mdp файлов
154 1 Andrey Golovin
</code></pre>
155 1 Andrey Golovin
156 2 Andrey Golovin
*Задать i вне скрипта можно командой export i="be". * 
157 1 Andrey Golovin
158 2 Andrey Golovin
* 4 У Вас должно получиться 5 tpr файлов. Теперь для каждого из них запустим mdrun.
159 1 Andrey Golovin
<pre><code class="python"> 
160 6 Andrey Golovin
gmx mdrun -deffnm et_%s -v -nt 1
161 1 Andrey Golovin
</code></pre>
162 1 Andrey Golovin
163 2 Andrey Golovin
* 5 Теперь переходим к анализу результатов. Начнем с визуального анализа. Для каждой из 5 систем проведите конвертацию в pdb и просмотрите в PyMol.
164 1 Andrey Golovin
<pre><code class="python">
165 6 Andrey Golovin
gmx trjconv -f et_%s.trr -s et_%s.tpr -o et_%s.pdb
166 1 Andrey Golovin
</code></pre>
167 1 Andrey Golovin
168 1 Andrey Golovin
В отчёт занесите ваши наблюдения и предварительные выводы. 
169 2 Andrey Golovin
* 6 Сравним потенциальную энергию и кинетическую энергию для каждой из 5 систем.
170 1 Andrey Golovin
171 1 Andrey Golovin
<pre><code class="python">
172 6 Andrey Golovin
echo  -e "12\n13" | gmx energy -f et_%s.edr -o et_%s_en.xvg -xvg none
173 1 Andrey Golovin
</code></pre>
174 1 Andrey Golovin
175 1 Andrey Golovin
Постройте графики изменения энергий. Рекомендуемый вид это lines. Графики добавьте в отчёт.
176 1 Andrey Golovin
Удобно загружать данные с numpy.loadtxt
177 1 Andrey Golovin
178 1 Andrey Golovin
<pre><code class="python">
179 1 Andrey Golovin
import numpy as np
180 1 Andrey Golovin
a= np.loadtxt('%s.xvg' % name)
181 1 Andrey Golovin
t=a[:,0]
182 1 Andrey Golovin
y=... догадайтесь
183 1 Andrey Golovin
</code></pre>
184 1 Andrey Golovin
185 1 Andrey Golovin
186 2 Andrey Golovin
* 7 Рассмотрим распределение длинны связи С-С за время моделирования. Сначала создадим индекс файл с одной связью. В текстовом редакторе создайте файл b.ndx со следующим содержимым:
187 1 Andrey Golovin
<pre> 
188 1 Andrey Golovin
[ b ]
189 1 Andrey Golovin
1 2 
190 1 Andrey Golovin
</code></pre>
191 1 Andrey Golovin
192 1 Andrey Golovin
И запустим утилиту по анализу связей g_bond:
193 1 Andrey Golovin
194 1 Andrey Golovin
<pre><code class="python">
195 6 Andrey Golovin
gmx distance -f et_%s.trr -s et_%s.tpr -o bond_%s.xvg -n b.ndx -xvg none
196 1 Andrey Golovin
</code></pre>
197 1 Andrey Golovin
Постройте графики распределения длинн связей. Рекомендуемый вид это гистограмма. Графики добавьте в отчёт.
198 1 Andrey Golovin
199 2 Andrey Golovin
* 8 Учитывая форму распределения Больцмана и все Ваши наблюдения сделайте вывод о том какой из методов позволяет наиболее реалистично поддерживать температуру в системе. Опишите найденные Вами недостатки предложенных алгоритмов. Постройте таблицу зависимости быстродействия от алгоритма.