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

Семестры Студенты Преподаватели

• 2024
• 2023
• 2022
• 2021
• 2020
• 2019
• 2018
• 2016

• ...

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

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

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

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

2. Вам даны 5 конфигурационных файлов с разными алгоритмами контроля температуры:

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

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

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

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

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

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

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

   1 gmx 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".

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

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

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

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

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

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

   1 gmx 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=... догадайтесь

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

   1 [ b ]
   2 1 2 

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

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

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

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