Structural Bioinformatics. Homework 8¶
Elen Tevanyan
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
import matplotlib
from matplotlib import pyplot as plt
%matplotlib inline
Загрузка этана и подготовка топологии этана¶
#etane.pro - имя файла
%%bash
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
%%bash
cat /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp
Этан - алкан, H3CCH3. Из файла выше для определения типа атомов внимание привлекают следующие строчки:
opls_135 12.01100 ; alkane CH3
opls_136 12.01100 ; alkane CH2
opls_137 12.01100 ; alkane CH
opls_138 12.01100 ; alkane CH4
opls_139 12.01100 ; alkane C
opls_140 1.00800 ; alkane H.
В этане два углерода, связанных с тремя водородами, поэтому беру 135 для описания типа молекулы углерода
ettop ='''
#include "/usr/share/gromacs/top/oplsaa.ff/forcefield.itp"
[ moleculetype ]
; Name nrexcl
et 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 opls_135 1 ETH C1 1 -0.189 12.01
2 opls_135 1 ETH C2 2 -0.155 12.01
3 opls_140 1 ETH H1 3 0.0059 1.008
4 opls_140 1 ETH H2 4 0.0059 1.008
5 opls_140 1 ETH H3 5 0.0059 1.008
6 opls_140 1 ETH H4 6 0.0056 1.008
7 opls_140 1 ETH H5 7 0.0056 1.008
8 opls_140 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
1 2 6 1
1 2 7 1
1 2 8 1
6 2 7 1
6 2 8 1
7 2 8 1
[ 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
4 1 2 6 1
4 1 2 7 1
4 1 2 8 1
5 1 2 6 1
5 1 2 7 1
5 1 2 8 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
'''
#запись файла
with open ('et.top','w') as e:
e.write(ettop)
Применение методов контроля температуры¶
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/an.mdp - метод Андерсена для контроля температуры. http://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp - метод стохастической молекулярной динамики.
%%bash
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
%%bash
for s in be vr nh an sd
do
grompp -f $s.mdp -c etane.gro -p et.top -o et_$s.tpr
done
%%bash
for s in be vr nh an sd
do
mdrun -deffnm et_$s -v -nt 1
done
Конвертацию проводила через терминал:
for s in be vr nh an sd
do
trjconv -f et$s.trr -s et$s.tpr -o et_$s.pdb
done
Визуализация¶
Видео созданы на десктопной версии PyMol на родной машине, поэтому сюда только грузим:
import io
from base64 import b64encode
from IPython.core.display import HTML, display
for i in ["be","vr","nh","an","sd"]:
video = io.open('%s.mp4'%i, 'r+b').read()
video_encoded = b64encode(video)
html_code = '<h1>{0}</h1><video controls alt="PyMol Movie" src="data:video/mp4;base64,{1}" type="video/mp4">'.format(i, video_encoded)
display(HTML(html_code))
Если располагать методы по принципу от самого "статичного" до самого "неустойчивого", то рейтинг такой:
1) метод Андерсена: ощущение, что колеблется только длина связи С-С
2) метод "Velocity rescale": во все стороны, аккуратно
3) метод Нуза-Хувера: водородные связи
4) метод Берендсена: похож по характеру на VS, но движение связи С-С более аккуратное
5) метод стохастической молекулярной динамики: молекула крутится во все стороны и ведет себя безумно.
Потенциальная и кинетические энергии¶
Энергию считала через терминал вручную типом команд, выбирая опции 10 (потенциальная) и 11 (кинетическая):
g_energy -f et_be.edr -o et_be_en.xvg -xvg none¶
for i in ['be','vr', 'nh', 'an', 'sd']:
a = np.loadtxt('et_%s_en.xvg'%i)
t = a[:,0]
ep = a[:,1]
ek = a[:,2]
plt.plot(t,ep, color='blue',label='Potential Energy')
plt.plot(t,ek, color='Red', label='Kinetic Energy')
plt.legend()
plt.title('Method %s'%i)
plt.show()
Вижу два вывода из графиков:
1) Характер изменения потенциальной и кинетической энергии совпадает
2) Сильно разнится от метода к методу амплитуда изменений обеих энергий
Распределение длин связи С-С¶
%%writefile b.ndx
[ b ]
1 2
%%bash
for s in be vr nh an sd
do
g_bond -f et_$s.trr -s et_$s.tpr -o bond_$s.xvg -n b.ndx -xvg none
done
for i in ['be','vr', 'nh', 'an', 'sd']:
a = np.loadtxt('bond_%s.xvg'%i)
t = a[:,0]
l = a[:,1]
plt.bar(t,l, 0.0001,color='green', edgecolor = 'black')
plt.title('C-C length distribution: %s'%i)
plt.show()
Время работы¶
Замечу, что результаты искажаются, если повторно замерять время, т.к. перезаписываются имена файлов на прошлом шаге
import subprocess
import time
import pandas as pd
times = []
for i in ["be","vr","nh","an","sd"]:
start = time.time()
do = 'mdrun -deffnm et_%s -v -nt 1' %i
subprocess.call(do,shell=True)
times.append(time.time() - start)
methods = ['Berendsen Method','Velocity Rescale','Nose-Hoover MEthod','Andersen Method','Stochastic MD']
table = pd.DataFrame({"Method:":methods, "Time":times})
table
Выводы¶
По времени, самый быстрый метод - метод Андерсена, самый длительный - метод стохастической молекулярной динамики.
- Метод Берендсена: кинетическая и потенциальные энергии сокращают амплитуду колебаний и стабилизируются на схожих уровнях; длина связи С-С колеблется, на распределение Больцмана похоже
- Метод Velocity Rescale: сильные колебания энергий, небольшой разброс в колебании длины С-С, на распределение Больцмана похоже
- Метод Нуза-Хувера: самые высокие уровни колебаний энергии, средние колебания длины С-С и стабильные, на распределение Больцмана немного похоже
- Метод Андерсена: самый маленький по абсолютной величине диапазон колебании энергий, связь С-С растянута (максимальные значения), на распределение Больцмана не сильно похоже
- Метод стохастической молекулярной динамики: энергии колеблятся на среднем уровне, длина связи С-С изменяется скромно, работает дольше всех и самое безумное визуальное поведение молекулы
Мне кажется удачным метод Берендсена: стабилизируется энергия, длина С-С изменяется средне, работает быстрее почти всех методов
import subprocess
converted = subprocess.call(["jupyter-nbconvert", '--to', 'html',"SB. HW8. Tevanyan.ipynb"], shell=False)
if converted==0:
print 'Your ipynb-file was sucessfully converted!'
else:
print 'Smth went wrong. for instance, check the filename...'
show_link = '<a href="%s" target="_blank">You can download it here!</a>'%('SB. HW8. Tevanyan.html')
display(HTML(show_link))