History

Изучение работы методов контроля температуры в молекулярной динамике

Традиционные ссылки на полезные ресурсы:

Сегодня мы будем изучать как реализован контроль температуры в молекулярной динамике на примере GROMACS. Объект исследования это одна молекула этана.

  • 1. Начнем с того, что подготовим файл координат и файл топологии.
  1. Вам предоставлен gro файл(http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro) с координатами этана.
  2. Построим файл топологии 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 файлов с разными параметрами контроля температуры:

Начиная с этого момента вы можете написать функции, а можете делать всё вручную.

  • 3 Очень краткое описание программ и типов файлов вы можете найти здесь Сначала надо построить входные файлы для молекулярно-динамического движка 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 Учитывая форму распределения Больцмана и все Ваши наблюдения сделайте вывод о том какой из методов позволяет наиболее реалистично поддерживать температуру в системе. Опишите найденные Вами недостатки предложенных алгоритмов. Постройте таблицу зависимости быстродействия от алгоритма.