« Previous -
Version 3/8
(diff) -
Next » -
Current version
Andrey Golovin, 14.12.2017 01:35
Изучение работы методов контроля температуры в молекулярной динамике¶
Традиционные ссылки на полезные ресурсы:
- Уроки по работе с GROMACS находятся http://www.gromacs.org/Documentation/Tutorials .
- Введение в скриптовании в Bash http://www.opennet.ru/docs/RUS/bash_scripting_guide .
Сегодня мы будем изучать как реализован контроль температуры в молекулярной динамике на примере GROMACS. Объект исследования это одна молекула этана.
- 1. Начнем с того, что подготовим файл координат и файл топологии.
- Вам предоставлен gro файл(http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro) с координатами этана.
- Построим файл топологии 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
Этот вариант топологии работать не будет, надо изменить типы атомов.
Вам придется угадать типы атомов из файла на shadbox:
%%bash cat /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp
или намек как это сделать программно:
1# загрузим модули
2import rdkit
3from rdkit import Chem
4from rdkit.Chem import AllChem, Draw
5from rdkit.Chem import Descriptors
6from rdkit import RDConfig
7from IPython.display import Image,display
8import numpy as np
1# создадим этан
2mol=Chem.MolFromSmiles('CC')
3AllChem.Compute2DCoords(mol)
4m3d=Chem.AddHs(mol)
5Chem.AllChem.EmbedMolecule(m3d)
6AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
1# и придумаем циклы
2## связи
3bonds = m3d.GetBonds()
4for i,b in enumerate(bonds):
5 print b.GetBeginAtomIdx() , b.GetEndAtomIdx()
6
7## углы
8for b1 in m3d.GetBonds():
9 for b2 in m3d.GetBonds():
10 if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetIdx() != b2.GetIdx():
11 print b1.GetEndAtomIdx() , b1.GetBeginAtomIdx(), b2.GetEndAtomIdx()
12## dihedrals
13for b1 in m3d.GetBonds():
14 for b2 in m3d.GetBonds():
15 for b3 in ....
16 .......
- 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:
1grompp -f %s.mdp -c et.gro -p et.top -o et_%s.tpr
2# где i: be,vr,nh,an,sd см. выше список mdp файлов
*Задать i вне скрипта можно командой export i="be". *
- 4 У Вас должно получиться 5 tpr файлов. Теперь для каждого из них запустим mdrun.
1 2mdrun -deffnm et_%s -v -nt 1
- 5 Теперь переходим к анализу результатов. Начнем с визуального анализа. Для каждой из 5 систем проведите конвертацию в pdb и просмотрите в PyMol.
1trjconv -f et_%s.trr -s et_%s.tpr -o et_%s.pdb
- 6 Сравним потенциальную энергию и кинетическую энергию для каждой из 5 систем.
1g_energy -f et_%s.edr -o et_%s_en.xvg -xvg none
Постройте графики изменения энергий. Рекомендуемый вид это lines. Графики добавьте в отчёт.
Удобно загружать данные с numpy.loadtxt
1import numpy as np
2a= np.loadtxt('%s.xvg' % name)
3t=a[:,0]
4y=... догадайтесь
- 7 Рассмотрим распределение длинны связи С-С за время моделирования. Сначала создадим индекс файл с одной связью. В текстовом редакторе создайте файл b.ndx со следующим содержимым:
[ b ] 1 2
И запустим утилиту по анализу связей g_bond:
1g_bond -f et_%s.trr -s et_%s.tpr -o bond_%s.xvg -n b.ndx -xvg none
Постройте графики распределения длинн связей. Рекомендуемый вид это гистограмма. Графики добавьте в отчёт.
- 8 Учитывая форму распределения Больцмана и все Ваши наблюдения сделайте вывод о том какой из методов позволяет наиболее реалистично поддерживать температуру в системе. Опишите найденные Вами недостатки предложенных алгоритмов. Постройте таблицу зависимости быстродействия от алгоритма.