Task8a
Version 5 (Andrey Golovin, 14.10.2021 18:21) → Version 6/8 (Andrey Golovin, 20.10.2021 18:57)
h1. Изучение работы методов контроля температуры в молекулярной динамике
Традиционные ссылки на полезные ресурсы:
* Уроки по работе с GROMACS находятся http://www.gromacs.org/Documentation/Tutorials .
* Введение в скриптовании в 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/an.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 g_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 g_bond -f et_%s.trr -s et_%s.tpr -o bond_%s.xvg -n b.ndx -xvg none
</code></pre>
Постройте графики распределения длинн связей. Рекомендуемый вид это гистограмма. Графики добавьте в отчёт.
* 8 Учитывая форму распределения Больцмана и все Ваши наблюдения сделайте вывод о том какой из методов позволяет наиболее реалистично поддерживать температуру в системе. Опишите найденные Вами недостатки предложенных алгоритмов. Постройте таблицу зависимости быстродействия от алгоритма.
Традиционные ссылки на полезные ресурсы:
* Уроки по работе с GROMACS находятся http://www.gromacs.org/Documentation/Tutorials .
* Введение в скриптовании в 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/an.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 g_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 g_bond -f et_%s.trr -s et_%s.tpr -o bond_%s.xvg -n b.ndx -xvg none
</code></pre>
Постройте графики распределения длинн связей. Рекомендуемый вид это гистограмма. Графики добавьте в отчёт.
* 8 Учитывая форму распределения Больцмана и все Ваши наблюдения сделайте вывод о том какой из методов позволяет наиболее реалистично поддерживать температуру в системе. Опишите найденные Вами недостатки предложенных алгоритмов. Постройте таблицу зависимости быстродействия от алгоритма.