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

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

• Уроки по работе с GROMACS находятся https://tutorials.gromacs.org/ http://www.mdtutorials.com/gmx/.
• Введение в скриптование в Bash http://www.opennet.ru/docs/RUS/bash_scripting_guide .
• Colab заготовка находится тут https://colab.research.google.com/drive/1ztk3rEzEJWLd3mgHS45ElsM-Vj5zLROf?usp=sharing .
• Важно установить Gromacs если вы работаете локально.

Сегодня мы будем изучать как реализован контроль температуры в молекулярной динамике на примере 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

Этот вариант топологии работать не будет, надо изменить типы атомов.

%%bash
cat /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp

или намек как это сделать программно:

# загрузим модули
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

# создадим этан
mol=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )

 # и придумаем циклы 
 ## связи
 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 ....
             .......

• 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/sd.mdp - метод стохастической молекулярной динамики.

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

• 3 Очень краткое описание программ и типов файлов вы можете найти здесь Сначала надо построить входные файлы для молекулярно-динамического движка mdrun с помощью grompp:

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". *

• 4 У Вас должно получиться 5 tpr файлов. Теперь для каждого из них запустим mdrun.
  

  
  gmx mdrun -deffnm et_%s -v -nt 1
  

• 5 Теперь переходим к анализу результатов. Начнем с визуального анализа. Для каждой из 5 систем проведите конвертацию в pdb и просмотрите в PyMol.
  

  gmx trjconv -f et_%s.trr -s et_%s.tpr -o et_%s.pdb
  

В отчёт занесите ваши наблюдения и предварительные выводы.

• 6 Сравним потенциальную энергию и кинетическую энергию для каждой из 5 систем.

echo  -e "12\n13" | gmx energy -f et_%s.edr -o et_%s_en.xvg -xvg none

Постройте графики изменения энергий. Рекомендуемый вид это lines. Графики добавьте в отчёт.
Удобно загружать данные с numpy.loadtxt

import numpy as np
a= np.loadtxt('%s.xvg' % name)
t=a[:,0]
y=... догадайтесь

• 7 Рассмотрим распределение длинны связи С-С за время моделирования. Сначала создадим индекс файл с одной связью. В текстовом редакторе создайте файл b.ndx со следующим содержимым:
  

   
  [ b ]
  1 2 
  

И запустим утилиту по анализу связей g_bond:

gmx distance -f et_%s.trr -s et_%s.tpr -o bond_%s.xvg -n b.ndx -xvg none

Постройте графики распределения длинн связей. Рекомендуемый вид это гистограмма. Графики добавьте в отчёт.

• 8 Учитывая форму распределения Больцмана и все Ваши наблюдения сделайте вывод о том какой из методов позволяет наиболее реалистично поддерживать температуру в системе. Опишите найденные Вами недостатки предложенных алгоритмов. Постройте таблицу зависимости быстродействия от алгоритма.