Kodomo

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

Учебная страница курса биоинформатики,
год поступления 2015

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

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

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

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

  3. Построим файл топологии 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 после номеров атомов это тип потенциала. Неожиданно в поле OPLS тип функционала для торсионных углов не "1". Узнайте какой это тип. Посмотрите файл /usr/share/gromacs/top/oplsaa.ff/forcefield.itp там будет ссылка на файл с парметрами для ковалентных взаимодействий, в том файле и надо найти нужное значение для типа dihedrals.

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

   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              .......
  1. Вам даны 5 конфигурационных файлов с разными алгоритмами контроля температуры:
    • be.mdp - метод Берендсена для контроля температуры.

    • vr.mdp - метод "Velocity rescale" для контроля температуры.

    • nh.mdp - метод Нуза-Хувера для контроля температуры.

    • an.mdp - метод Андерсена для контроля температуры.

    • sd.mdp - метод стохастической молекулярной динамики.

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

  1. Очень краткое описание программ и типов файлов вы можете найти здесь Сначала надо построить входные файлы для молекулярно-динамического движка 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".

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

   1 mdrun -deffnm et_%s -v -nt 1
  1. Теперь переходим к анализу результатов. Начнем с визуального анализа. Для каждой из 5 систем проведите конвертацию в pdb и просмотрите в PyMol.

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

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

  1. Сравним потенциальную энергию и кинетическую энергию для каждой из 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=... догадайтесь
  1. Рассмотрим распределение длинны связи С-С за время моделирования. Сначала создадим индекс файл с одной связью. В текстовом редакторе создайте файл 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

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

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