HW 8. Didkovskaya¶
Загрузка модулей
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 subprocess
import time
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
Часть 1. Подготовка файлов¶
Скачаем все стартовые файлы
%%bash
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
for method in ["be","vr","nh","an","sd"]:
get='wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/%s.mdp' %method
subprocess.Popen(get,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
Создадим файл топологии
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 opls_139 1 ETH C1 1 -0.189 12.01
2 opls_139 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
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
"""
with open("et.top","w") as topology:
topology.write(top)
Cоздадим этан
mol=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
Циклы¶
1) Связи:
bonds = m3d.GetBonds()
for i,b in enumerate(bonds):
print b.GetBeginAtomIdx()+1 , b.GetEndAtomIdx()+1, 1
2) Углы:
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
3) 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, 1
Часть 2. Сравнение методов контроля¶
Оценим работу 5 методов контроля температуры
%%bash
ls
Строим tpr файлы
for method in ["be","vr","nh","an","sd"]:
#строим файлы .tpr для запуска молдинамики
do1 = 'grompp -f %s.mdp -c etane.gro -p et.top -o et_%s.tpr' %(method,method)
subprocess.Popen(do1,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
%%bash
ls
Запускаем для каждого tpr mdrun
for method in ["be","vr","nh","an","sd"]:
start = time.time()
do = 'mdrun -deffnm et_%s -v -nt 1' %method
subprocess.Popen(do,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
print("Runtime for %s-method: %s seconds" % (method,round(time.time() - start, 3)))
Переводим в pdb формат
for method in ["be","vr","nh","an","sd"]:
do = 'echo 0| trjconv -f et_%s.trr -s et_%s.tpr -o et_%s.pdb' %(method,method,method)
subprocess.Popen(do,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
Анализ результатов:¶
Pdb файлы были скопированы на локальный компьютер и визуализированы в Pymol.
Визуально оценка:
1) наиболее хаотично движется молекула этана из метода стохастической молекулярной динамики по неопознанным осям :)
2) самая статичная молекула при методе Андерсена(только колебание длин связей СС, без вращения самой молекулы)
3) вращение молекулы вокруг связи СС происходит для метода Нуза-Хувера
4) для метода Берендсена и метода "Velocity rescale" – вращение по другим осям в пространстве
Часть 3. Сравнение энергий¶
for method in ["be","vr","nh","an","sd"]:
#создаем файл с потенциальной (do1) энергией системы и кинетической (do2)
do1 = 'echo 10 | g_energy -f et_%s.edr -o et_%s_en_potential.xvg -xvg none' %(method,method)
subprocess.Popen(do1,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
do2 = 'echo 11 | g_energy -f et_%s.edr -o et_%s_en_kinetic.xvg -xvg none' %(method,method)
subprocess.Popen(do2,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
Получили для каждого метода два файла .xvg для потенциальной и кинетической энергий
Визуализируем энергии для каждого метода:
for method in ["be","vr","nh","an","sd"]:
a_potential= np.loadtxt('et_%s_en_potential.xvg' %method)
a_kinetic= np.loadtxt('et_%s_en_kinetic.xvg' %method)
t=a_kinetic[:,0]
e_potential=a_potential[:,1]
e_kinetic=a_kinetic[:,1]
plt.plot(t, e_potential, color='blue')
plt.plot(t, e_kinetic, color='yellow')
plt.title('Method %s' %method)
kin_patch = mpatches.Patch(color='yellow', label='Kinetic energy')
pot_patch = mpatches.Patch(color='blue', label='Potential energy')
plt.legend(handles=[kin_patch,pot_patch])
plt.show()
Анализ графиков:¶
1) значения кинеической и потенциальной энергий достаточно близки для большинства случаев
2) в методе Нуза-Хувера кинетическая энергия заметно больше потенциальной
3) на всех графиках видны колебания
4) наименьшая амплитуда у метода Андерсена и затухающие колебания в методе Берендсена
Часть 4. Распределение длинны связи С-С¶
Сначала создадим индекс файл с одной связью. В текстовом редакторе создайте файл b.ndx со следующим содержимым:
%%bash
echo -e "[ b ]\n1 2" > b.ndx
И запустим утилиту по анализу связей g_bond:
for method 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' %(method,method, method)
subprocess.Popen(do,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
Получили для каждого метода файл с информацией о связи
Визуализация:
for method in ["be","vr","nh","an","sd"]:
bonds = np.loadtxt('bond_%s.xvg' % method)
x = bonds[:,0]
y = bonds[:,1]
plt.bar(x, y, 0.001, color="blue")
plt.title('Length of C-C bond (%s method)' % method)
plt.xlim(0.14, 0.17)
plt.show()
Анализ графиков распределений:¶
1) распределение из метода Андерсена демонстрирует маленькую дисперсию, что не соответсвует наблюдениям из Pymol
2) на нормальное распределение похожи: распределения Берендсена, Нуза-Хувера и стохастическое
Выводы:¶
Показатели времени работы (из части 2):
Runtime for be-method: 0.012 seconds
Runtime for vr-method: 0.01 seconds
Runtime for nh-method: 0.008 seconds
Runtime for an-method: 0.007 seconds
Runtime for sd-method: 0.007 seconds
1) Метод Андерсена:
имеет очень маленький разброс длины связи C-C
молекула не движется
маленькая энергия молекулы
работает быстро
2) Метод Нуза-Хувера:
сильные выбросы значений энергии системы
работает быстро
молекула вращается вокруг СС
амплитуда длины связи СС не очень маленькая
3) Метода Берендсена:
амплитуда C-C связи изменяется не сильно
работает не сильно дольше остальных быстро
распределение наиболее схоже с распределением Максвела-Больцмана
4) Методы стохастической молекулярной динамики и метод Velocity rescale:
схожие распределения с распределением Максвелла-Больцмана
стохастический метод придаёт молекуле хаотичные и быстрые движения
скорости работы примерно одинаковые
На мой взгляд оптимальными являются методы Velocity rescale и Берендсена