Домашнее задание 8¶
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
from tabulate import tabulate
import io
import base64
from IPython.display import HTML
# создадим этан
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):
#нумерование вроде с 1 у атомов
print b.GetBeginAtomIdx()+1 , b.GetEndAtomIdx()+1,1
## углы
for b1 in m3d.GetBonds():
for b2 in m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetIdx() != b2.GetIdx():
print b1.GetEndAtomIdx()+1 , b1.GetBeginAtomIdx()+1, b2.GetEndAtomIdx()+1,1
Удалим дупликаты углов: например, для 2 1 3 1 (с1) можем удалить 3 1 2 1 и тому подобное.
## dihedrals
for b1 in m3d.GetBonds():
for b2 in m3d.GetBonds():
for b3 in m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b2.GetEndAtomIdx() == b3.GetBeginAtomIdx() and b1.GetIdx() != b2.GetIdx():
print b1.GetEndAtomIdx()+1 , b1.GetBeginAtomIdx()+1, b2.GetEndAtomIdx()+1, b3.GetEndAtomIdx()+1
Тогда топологический файл будет выглядеть следующим образом:
ethan="""
#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
2 6 1
2 7 1
2 8 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
7 2 8 1
8 2 6 1
7 2 6 1
1 2 6 1
1 2 7 1
1 2 8 1
[ dihedrals ]
; ai aj ak al funct
3 1 2 6 3
3 1 2 7 3
3 1 2 8 3
4 1 2 6 3
4 1 2 7 3
4 1 2 8 3
5 1 2 6 3
5 1 2 7 3
5 1 2 8 3
[ 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
"""
Скачиваем файлы:
!wget "http://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp"
!wget "http://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp"
!wget "http://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp"
!wget "http://kodomo.fbb.msu.ru/FBB/year_08/term6/an.mdp"
!wget "http://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp"
!wget "http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro"
Дальше я делал все вручную(всего 5 файлов - на все это ушло не так много времени):
1) С помощью команды grompp создал 5 файлов
2) Затем применил к каждому файлу mdrun
3) Затем произвел конвертацию в pdb.
На выходе мы получили pdb файл со множеством моделирований. Из него удобно сделать видео в PYMOL.
i=['be','vr','nh','an','sd']
method=['Berendsen','Velocity rescale','Nose-Hoover','Andersen','Stochastic Molecular Dynamic']
for j,name in enumerate(i):
print (method[j])
video = io.open('{0}.mp4'.format(name), 'r+b').read()
encoded = base64.b64encode(video)
display(HTML(data='<video alt="test" controls><source src="data:video/mp4;base64,{0}" type="video/mp4" /></video>'.format(encoded)))
Вывод по видео¶
Ну что ж. Метод Андерсона в целом демонстрирует молекулу почти в статическом положении (разве что, видимо, изменения происходят только со связью С-С). Стохастический метод "вселяет" беса в этан и тот бьется в диких конвульсиях. Относительно остальных трех методов имеет место быть вращение водородных связей, которое наиболее ярко выражено в методе Берендсена. Там скорость этого вращения ,судя по всему, наибольшая.
Время работы алгоритмов¶
for j,name in enumerate(i):
print (method[j])
display(Image(filename='time/et_{0}.png'.format(name)))
Потенциальная и кинетическая энергии¶
C помощью genergy извлечем данные о потенциальной и кинетической энергий из et%s.edr файлов. В итоге: второй столбец в файле - потенциальная, третий - кинетическая энергии.
Оранжевый - потенциальная энергия, синий - кинетическая энергия.
import matplotlib.pyplot as plt
import matplotlib.patches as ptch
for j,name in enumerate(i):
print (method[j])
a= np.loadtxt('energy/et_{0}_en.xvg'.format(name))
#moment of time
t=a[:,0]
#potential
plt.plot(t,a[:,1],'orange',label='Potential energy')
#kinetic
plt.plot(t,a[:,2],'b',label='Kinetic energy')
pot=ptch.Patch(color='orange',label='Potential energy')
kin=ptch.Patch(color='blue',label='Kinetic energy')
plt.legend(handles=[pot,kin])
plt.show()
#y=... догадайтесь
Распределение длин связей С-С за время моделирования.¶
for j,name in enumerate(i):
print (method[j])
a= np.loadtxt('C-C_distr/bond_{0}.xvg'.format(name))
t=a[:,0]
#бины выбрал по алгоритму Скотта
plt.hist(a[:,1],bins='scott')
plt.show()
#y=... догадайтес
Табличка со временем алгоритмов:
time=[(3.332,3.674,3.193,3.776,3.974)]
print tabulate(tabular_data=time,headers=method,showindex=[('Real time')])
Вывод¶
Судя по построенным распределениям можно сделать вывод, что больше всего для наших целей могут подходить метод Берендсена и метод Нуса-Хувера. Их распределния больше всего походят на Больцмана. Такой же вывод напрашивается и из видео. Кроме того, по времени эти два алгоритма показали наилучший результат, что однозначно можно также занести им в плюс)