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

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

Уроки по работе с GROMACS находятся здесь.
Введение в скриптовании в Bash здесь.

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

Начнем с того, что подготовим файл координат и файл топологии.

a) Вам предоставлен gro файл(http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro) с координатами этана.

b) Построим файл топологии 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

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

   1 # загрузим модули
   2 import rdkit
   3 from rdkit import Chem
   4 from rdkit.Chem import AllChem, Draw
   5 from rdkit.Chem import Descriptors
   6 from rdkit import RDConfig
   7 from IPython.display import Image,display
   8 import numpy as np

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

   1 # и придумаем циклы 
   2 ## связи
   3 bonds = m3d.GetBonds()
   4 for i,b in enumerate(bonds):
   5     print b.GetBeginAtomIdx() , b.GetEndAtomIdx()
   6 
   7 ## углы    
   8 for 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            
  13 for b1 in m3d.GetBonds():
  14     for b2 in m3d.GetBonds():
  15         for b3 in ....
  16              .......

Вам даны 5 файлов с разными параметрами контроля температуры:
- be.mdp - метод Берендсена для контроля температуры.
- vr.mdp - метод "Velocity rescale" для контроля температуры.
- nh.mdp - метод Нуза-Хувера для контроля температуры.
- an.mdp - метод Андерсена для контроля температуры.
- sd.mdp - метод стохастической молекулярной динамики.

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

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

   1 grompp -f %s.mdp -c et.gro -p et.top -o et_%s.tpr
   2 # где i: be,vr,nh,an,sd  см. выше список mdp файлов
   3

Задать i вне скрипта можно командой export i="be".

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

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

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

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

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

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

   1 g_energy -f et_%s.edr -o et_%s_en.xvg -xvg none

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

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

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

   1 [ b ]
   2 1 2

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

   1 g_bond -f et_%s.trr -s et_%s.tpr -o bond_%s.xvg -n b.ndx -xvg none

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

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

Kodomo

Пользователь

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