#pragma css /css/2017.css
<<BI>>
= Изучение работы методов контроля температуры в молекулярной динамике =
Традиционные ссылки на полезные ресурсы:
* Уроки по работе с GROMACS находятся [[http://www.gromacs.org/Documentation/Tutorials|здесь]].
'''Сегодня мы будем изучать как реализован контроль температуры в молекулярной динамике на примере GROMACS. Объект исследования это одна молекула этана.'''
<<BR>>
1. Начнем с того, что подготовим файл координат и файл топологии.
a. Вам предоставлен gro файл(http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro) с координатами этана.
a. Построим файл топологии et.top для этана, который выгляди примерно так:
{{{
#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
}}}
Этот вариант тпологии работать не будет, надо изменить типы атомов.
Вам придется угадать типы атомов из файла на kodomo:
{{{
%%bash
cat /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp
}}}
{{{#!wiki green/solid
Есть ещё один квест, 1 после номеров атомов это тип потенциала. Неожиданно в поле OPLS тип функционала для торсионных углов не "1". Узнайте какой это тип. Посмотрите файл /usr/share/gromacs/top/oplsaa.ff/forcefield.itp там будет ссылка на файл с парметрами для ковалентных взаимодействий, в том файле и надо найти нужное значение для типа dihedrals.
}}}
или намек как это сделать программно:
{{{#!highlight 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
}}}
{{{#!highlight python
# создадим этан
mol=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
}}}
{{{#!highlight 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 ....
.......
}}}
1.#2 Вам даны 5 конфигурационных файлов с разными алгоритмами контроля температуры:
* [[http://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp|be.mdp]] - метод Берендсена для контроля температуры.
* [[http://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp|vr.mdp]] - метод "Velocity rescale" для контроля температуры.
* [[http://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp|nh.mdp]] - метод Нуза-Хувера для контроля температуры.
* [[http://kodomo.fbb.msu.ru/FBB/year_08/term6/an.mdp|an.mdp]] - метод Андерсена для контроля температуры.
* [[http://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp|sd.mdp]] - метод стохастической молекулярной динамики.
Начиная с этого момента вы можете написать функции, а можете делать всё вручную.
1.#3 Очень краткое описание программ и типов файлов вы можете найти [[https://kodomo.fbb.msu.ru/wiki/Main/Modelling/BioNanoPrac|здесь]] Сначала надо построить входные файлы для молекулярно-динамического движка mdrun с помощью grompp:
{{{#!highlight bash
gmx grompp -f %s.mdp -c et.gro -p et.top -o et_%s.tpr
# где i: be,vr,nh,an,sd см. выше список mdp файлов
}}}
'''Задать i вне скрипта можно командой export i="be". '''
1.#4 У Вас должно получиться 5 tpr файлов. Теперь для каждого из них запустим mdrun.
{{{#!highlight bash
gmx mdrun -deffnm et_%s -v -nt 1
}}}
1.#5 Теперь переходим к анализу результатов. Начнем с визуального анализа. Для каждой из 5 систем проведите конвертацию в pdb и просмотрите в PyMol.
{{{#!highlight bash
gmx trjconv -f et_%s.trr -s et_%s.tpr -o et_%s.pdb
}}}
В отчёт занесите ваши наблюдения и предварительные выводы.
1.#6 Сравним потенциальную энергию и кинетическую энергию для каждой из 5 систем.
{{{#!highlight bash
gmx energy -f et_%s.edr -o et_%s_en.xvg -xvg none
}}}
Постройте графики изменения энергий. Рекомендуемый вид это lines. Графики добавьте в отчёт.
Удобно загружать данные с numpy.loadtxt
{{{#!highlight python
import numpy as np
a= np.loadtxt('%s.xvg' % name)
t=a[:,0]
y=... догадайтесь
}}}
1.#7 Рассмотрим распределение длинны связи С-С за время моделирования. Сначала создадим индекс файл с одной связью. В текстовом редакторе создайте файл b.ndx со следующим содержимым:
{{{#!highlight bash
[ b ]
1 2
}}}
И запустим утилиту по анализу связей distance :
{{{#!highlight bash
gmx distance -f et_%s.trr -s et_%s.tpr -oh bond_%s.xvg -n b.ndx -xvg none
}}}
Постройте графики распределения длинн связей. Рекомендуемый вид это гистограмма. Графики добавьте в отчёт.
1.#8 Учитывая форму распределения Больцмана и все Ваши наблюдения сделайте вывод о том какой из методов позволяет наиболее реалистично поддерживать температуру в системе. Опишите найденные Вами недостатки предложенных алгоритмов. Постройте таблицу зависимости быстродействия от алгоритма.