Task8a
Version 3 (Andrey Golovin, 14.12.2017 01:35)
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 | 1 | Andrey Golovin | |
9 | 1 | Andrey Golovin | |
10 | 1 | Andrey Golovin | |
11 | 1 | Andrey Golovin | *Сегодня мы будем изучать как реализован контроль температуры в молекулярной динамике на примере GROMACS. Объект исследования это одна молекула этана.* |
12 | 1 | Andrey Golovin | |
13 | 1 | Andrey Golovin | |
14 | 1 | Andrey Golovin | |
15 | 2 | Andrey Golovin | * 1. Начнем с того, что подготовим файл координат и файл топологии. |
16 | 1 | Andrey Golovin | |
17 | 1 | Andrey Golovin | > # Вам предоставлен gro файл(http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro) с координатами этана. |
18 | 1 | Andrey Golovin | > # Построим файл топологии et.top для этана, который выгляди примерно так: |
19 | 1 | Andrey Golovin | |
20 | 1 | Andrey Golovin | <pre> |
21 | 1 | Andrey Golovin | |
22 | 1 | Andrey Golovin | #include "/usr/share/gromacs/top/oplsaa.ff/forcefield.itp" |
23 | 1 | Andrey Golovin | |
24 | 1 | Andrey Golovin | [ moleculetype ] |
25 | 1 | Andrey Golovin | ; Name nrexcl |
26 | 1 | Andrey Golovin | et 3 |
27 | 1 | Andrey Golovin | |
28 | 1 | Andrey Golovin | [ atoms ] |
29 | 1 | Andrey Golovin | ; nr type resnr residue atom cgnr charge mass |
30 | 1 | Andrey Golovin | 1 CX 1 ETH C1 1 -0.189 12.01 |
31 | 1 | Andrey Golovin | 2 CX 1 ETH C2 2 -0.155 12.01 |
32 | 1 | Andrey Golovin | 3 HX 1 ETH H1 3 0.0059 1.008 |
33 | 1 | Andrey Golovin | 4 HX 1 ETH H2 4 0.0059 1.008 |
34 | 1 | Andrey Golovin | 5 HX 1 ETH H3 5 0.0059 1.008 |
35 | 1 | Andrey Golovin | 6 HX 1 ETH H4 6 0.0056 1.008 |
36 | 1 | Andrey Golovin | 7 HX 1 ETH H5 7 0.0056 1.008 |
37 | 1 | Andrey Golovin | 8 HX 1 ETH H6 8 0.0056 1.008 |
38 | 1 | Andrey Golovin | |
39 | 1 | Andrey Golovin | [ bonds ] |
40 | 1 | Andrey Golovin | ; ai aj funct b0 kb |
41 | 1 | Andrey Golovin | 1 2 1 |
42 | 1 | Andrey Golovin | 1 3 1 |
43 | 1 | Andrey Golovin | 1 4 1 |
44 | 1 | Andrey Golovin | 1 5 1 |
45 | 1 | Andrey Golovin | ........... |
46 | 1 | Andrey Golovin | надо дописать |
47 | 1 | Andrey Golovin | |
48 | 1 | Andrey Golovin | [ angles ] |
49 | 1 | Andrey Golovin | ; ai aj ak funct phi0 kphi |
50 | 1 | Andrey Golovin | ;around c1 |
51 | 1 | Andrey Golovin | 3 1 4 1 |
52 | 1 | Andrey Golovin | 4 1 5 1 |
53 | 1 | Andrey Golovin | 3 1 5 1 |
54 | 1 | Andrey Golovin | 2 1 3 1 |
55 | 1 | Andrey Golovin | 2 1 4 1 |
56 | 1 | Andrey Golovin | 2 1 5 1 |
57 | 1 | Andrey Golovin | ;around c2 |
58 | 1 | Andrey Golovin | ........... |
59 | 1 | Andrey Golovin | надо дописать |
60 | 1 | Andrey Golovin | |
61 | 1 | Andrey Golovin | [ dihedrals ] |
62 | 1 | Andrey Golovin | ; ai aj ak al funct |
63 | 1 | Andrey Golovin | 3 1 2 6 1 |
64 | 1 | Andrey Golovin | 3 1 2 7 1 |
65 | 1 | Andrey Golovin | 3 1 2 8 1 |
66 | 1 | Andrey Golovin | 4 1 2 6 1 |
67 | 1 | Andrey Golovin | ........... |
68 | 1 | Andrey Golovin | надо дописать |
69 | 1 | Andrey Golovin | |
70 | 1 | Andrey Golovin | [ pairs ] |
71 | 1 | Andrey Golovin | ; список атомов 1-4 |
72 | 1 | Andrey Golovin | ; ai aj funct |
73 | 1 | Andrey Golovin | 3 6 |
74 | 1 | Andrey Golovin | 3 7 |
75 | 1 | Andrey Golovin | 3 8 |
76 | 1 | Andrey Golovin | 4 6 |
77 | 1 | Andrey Golovin | 4 7 |
78 | 1 | Andrey Golovin | 4 8 |
79 | 1 | Andrey Golovin | 5 6 |
80 | 1 | Andrey Golovin | 5 7 |
81 | 1 | Andrey Golovin | 5 8 |
82 | 1 | Andrey Golovin | |
83 | 1 | Andrey Golovin | [ System ] |
84 | 1 | Andrey Golovin | ; any text here |
85 | 1 | Andrey Golovin | first one |
86 | 1 | Andrey Golovin | [ molecules ] |
87 | 1 | Andrey Golovin | ;Name count |
88 | 1 | Andrey Golovin | et 1 |
89 | 1 | Andrey Golovin | |
90 | 1 | Andrey Golovin | </pre> |
91 | 1 | Andrey Golovin | |
92 | 1 | Andrey Golovin | Этот вариант топологии работать не будет, надо изменить типы атомов. |
93 | 1 | Andrey Golovin | |
94 | 1 | Andrey Golovin | Вам придется угадать типы атомов из файла на shadbox: |
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/an.mdp - метод Андерсена для контроля температуры. |
146 | 2 | Andrey Golovin | > * http://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp - метод стохастической молекулярной динамики. |
147 | 1 | Andrey Golovin | |
148 | 1 | Andrey Golovin | Начиная с этого момента вы можете написать функции, а можете делать всё вручную. |
149 | 2 | Andrey Golovin | |
150 | 2 | Andrey Golovin | * 3 Очень краткое описание программ и типов файлов вы можете найти [[https://kodomo.fbb.msu.ru/wiki/Main/Modelling/BioNanoPrac|здесь]] Сначала надо построить входные файлы для молекулярно-динамического движка mdrun с помощью grompp: |
151 | 1 | Andrey Golovin | |
152 | 1 | Andrey Golovin | <pre><code class="python"> |
153 | 1 | Andrey Golovin | grompp -f %s.mdp -c et.gro -p et.top -o et_%s.tpr |
154 | 1 | Andrey Golovin | # где i: be,vr,nh,an,sd см. выше список mdp файлов |
155 | 1 | Andrey Golovin | </code></pre> |
156 | 1 | Andrey Golovin | |
157 | 2 | Andrey Golovin | *Задать i вне скрипта можно командой export i="be". * |
158 | 1 | Andrey Golovin | |
159 | 2 | Andrey Golovin | * 4 У Вас должно получиться 5 tpr файлов. Теперь для каждого из них запустим mdrun. |
160 | 1 | Andrey Golovin | <pre><code class="python"> |
161 | 1 | Andrey Golovin | mdrun -deffnm et_%s -v -nt 1 |
162 | 1 | Andrey Golovin | </code></pre> |
163 | 1 | Andrey Golovin | |
164 | 2 | Andrey Golovin | * 5 Теперь переходим к анализу результатов. Начнем с визуального анализа. Для каждой из 5 систем проведите конвертацию в pdb и просмотрите в PyMol. |
165 | 1 | Andrey Golovin | <pre><code class="python"> |
166 | 1 | Andrey Golovin | trjconv -f et_%s.trr -s et_%s.tpr -o et_%s.pdb |
167 | 1 | Andrey Golovin | </code></pre> |
168 | 1 | Andrey Golovin | |
169 | 1 | Andrey Golovin | В отчёт занесите ваши наблюдения и предварительные выводы. |
170 | 2 | Andrey Golovin | * 6 Сравним потенциальную энергию и кинетическую энергию для каждой из 5 систем. |
171 | 1 | Andrey Golovin | |
172 | 1 | Andrey Golovin | <pre><code class="python"> |
173 | 1 | Andrey Golovin | g_energy -f et_%s.edr -o et_%s_en.xvg -xvg none |
174 | 1 | Andrey Golovin | </code></pre> |
175 | 1 | Andrey Golovin | |
176 | 1 | Andrey Golovin | Постройте графики изменения энергий. Рекомендуемый вид это lines. Графики добавьте в отчёт. |
177 | 1 | Andrey Golovin | Удобно загружать данные с numpy.loadtxt |
178 | 1 | Andrey Golovin | |
179 | 1 | Andrey Golovin | <pre><code class="python"> |
180 | 1 | Andrey Golovin | import numpy as np |
181 | 1 | Andrey Golovin | a= np.loadtxt('%s.xvg' % name) |
182 | 1 | Andrey Golovin | t=a[:,0] |
183 | 1 | Andrey Golovin | y=... догадайтесь |
184 | 1 | Andrey Golovin | </code></pre> |
185 | 1 | Andrey Golovin | |
186 | 1 | Andrey Golovin | |
187 | 2 | Andrey Golovin | * 7 Рассмотрим распределение длинны связи С-С за время моделирования. Сначала создадим индекс файл с одной связью. В текстовом редакторе создайте файл b.ndx со следующим содержимым: |
188 | 1 | Andrey Golovin | <pre> |
189 | 1 | Andrey Golovin | [ b ] |
190 | 1 | Andrey Golovin | 1 2 |
191 | 1 | Andrey Golovin | </code></pre> |
192 | 1 | Andrey Golovin | |
193 | 1 | Andrey Golovin | И запустим утилиту по анализу связей g_bond: |
194 | 1 | Andrey Golovin | |
195 | 1 | Andrey Golovin | <pre><code class="python"> |
196 | 1 | Andrey Golovin | g_bond -f et_%s.trr -s et_%s.tpr -o bond_%s.xvg -n b.ndx -xvg none |
197 | 1 | Andrey Golovin | </code></pre> |
198 | 1 | Andrey Golovin | Постройте графики распределения длинн связей. Рекомендуемый вид это гистограмма. Графики добавьте в отчёт. |
199 | 1 | Andrey Golovin | |
200 | 2 | Andrey Golovin | * 8 Учитывая форму распределения Больцмана и все Ваши наблюдения сделайте вывод о том какой из методов позволяет наиболее реалистично поддерживать температуру в системе. Опишите найденные Вами недостатки предложенных алгоритмов. Постройте таблицу зависимости быстродействия от алгоритма. |