In [3]:
# загрузим модули
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
In [2]:
# создадим этан
mol=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
Out[2]:
In [3]:
#скачиваем координаты этана
! wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
In [5]:
#скачиваем методы контроля температуры
! wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp #метод Берендсена для контроля температуры.
! wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp #метод "Velocity rescale" для контроля температуры.
! 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 #метод стохастической молекулярной динамики.
In [8]:
# и придумаем циклы
## связи
bonds = m3d.GetBonds()
for i,b in enumerate(bonds):
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
## 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
In [10]:
topology_table = """#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
"""
In [11]:
with open("et.top","w") as topology:
topology.write(topology_table)
Сравнение методов¶
In [1]:
methods=["be","vr","nh","an","sd"]
Строим файлы .tpr
In [13]:
for method in methods:
dofiles = 'grompp -f %s.mdp -c etane.gro -p et.top -o et_%s.tpr' %(method,method)
subprocess.Popen(dofiles,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
Запускаем mdrun
In [14]:
for method in methods:
start = time.time()
domdrun = 'mdrun -deffnm et_%s -v -nt 1' %method
subprocess.Popen(domdrun,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
Перевод в pdb
In [16]:
for method in methods:
dopdb = 'echo 0 | trjconv -f et_%s.trr -s et_%s.tpr -o et_%s.pdb' %(method,method,method)
subprocess.Popen(dopdb,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
Анализ просмотра pdb-шек в паймоле:¶
Метод стохастической молекулярной динамики показывает какую-то дискотеку этана.
Спокойнее всего молекула ведет себя при методе Андерсена: меняется длина связи С-С, без всяких вращений.
Velocity rescale и Нуза-Хувера похожи и создают вращение вокруг связи C-C.
В метода Берендсена также видно вращение вокруг связи C-C, но более активное.
Сравнение энергий¶
In [4]:
for method in methods:
dopotential = 'echo 10 | g_energy -f et_%s.edr -o et_%s_en_potential.xvg -xvg none' %(method,method)
subprocess.Popen(dopotential,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) #potential
dokinetic = 'echo 11 | g_energy -f et_%s.edr -o et_%s_en_kinetic.xvg -xvg none' %(method,method)
subprocess.Popen(dokinetic,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) #kinetic
In [14]:
for method in methods:
fig = plt.figure(figsize=(17, 5))
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='black')
plt.plot(t, e_kinetic, color='red')
plt.title('%s' %method)
pot_patch = mpatches.Patch(color='black', label='Potential energy')
kin_patch = mpatches.Patch(color='red', label='Kinetic energy')
plt.legend(handles=[kin_patch,pot_patch])
plt.show()
Видно, что значения кинетической и потенциальной энергий очень близки везде, кроме метода Нуза-Хувера - тут кинетическая энергия намного больше.
Распределение длины связи С-С за время моделирования¶
In [26]:
! echo -e "[ b ]\n1 2" > b.ndx
In [27]:
#И запустим утилиту по анализу связей g_bond:
for method in methods:
dosv = 'g_bond -f et_%s.trr -s et_%s.tpr -o bond_%s.xvg -n b.ndx -xvg none' %(method,method, method)
subprocess.Popen(dosv,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
In [36]:
for method in methods:
bonds = np.loadtxt('bond_%s.xvg' % method)
fig = plt.figure(figsize=(17, 5))
ax = fig.add_subplot(111)
plt.title('%s' % method)
ax.bar(np.arange(len(bonds[:,0])), bonds[:,1])
plt.show()
Берендсен, Нуза-Хувер и стохастик похожи на нормальные.
In [49]:
for method in methods:
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)))
Метод Берендсена кажется оптимальным, хоть и работает дольше всех