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
from collections import defaultdict
In [4]:
class M3dTop(object):
def __init__(self, smiles, atom_types = {'C' : "opls_076", 'H' : 'opls_140'}, name = "et",
nrexcl=3):
mol=Chem.MolFromSmiles(smiles)
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200)
AllChem.ComputeGasteigerCharges(m3d)
self.m3d = m3d
self.atom_types = atom_types
self.name = name
self.nrexcl = nrexcl
def write(self, out_file_name, ft=1, torsft=3):
self.out_file = open(out_file_name, "w")
self.out_file.write("#include \"/usr/share/gromacs/top/oplsaa.ff/forcefield.itp\"\n")
self.out_file.write("[ moleculetype ]\n{}\t{}\n".format(self.name, self.nrexcl))
self.get_atoms(ft)
self.get_bonds(ft)
self.get_angles(ft)
self.get_torsions(torsft)
self.out_file.write("[ System ]\n; any text here\nfirst one\n[ molecules ]\n;Name count\n{} 1".
format(self.name))
self.out_file.close()
def get_bonds(self, ft=1):
self.out_file.write("[ bonds ]\n")
self.out_file.write("; ai aj funct b0 kb\n")
for b in self.m3d.GetBonds():
self.out_file.write("{} {} {}\n".format(b.GetBeginAtomIdx() + 1,
b.GetEndAtomIdx() + 1,
ft))
def get_angles(self, ft=1):
self.out_file.write("[ angles ]\n")
self.out_file.write("; ai aj ak funct phi0 kphi\n")
for b1 in self.m3d.GetBonds():
for b2 in self.m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetIdx() != b2.GetIdx():
self.out_file.write("{}\t{}\t{}\t{}\n".format(b1.GetEndAtomIdx() + 1,
b1.GetBeginAtomIdx() + 1,
b2.GetEndAtomIdx() + 1,
ft))
def get_torsions(self, ft=3):
self.out_file.write("[ dihedrals ]\n")
self.out_file.write("; ai aj ak al funct\n")
atom14 = list()
for b1 in self.m3d.GetBonds():# a - b
for b2 in self.m3d.GetBonds():
for b3 in self.m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetIdx() != b2.GetIdx() and\
b2.GetEndAtomIdx() == b3.GetBeginAtomIdx() and b2.GetIdx() != b3.GetIdx():
ind1, ind2, ind3, ind4 = b1.GetEndAtomIdx() + 1, b1.GetBeginAtomIdx() + 1,\
b2.GetEndAtomIdx() + 1, b3.GetEndAtomIdx() + 1
self.out_file.write("{}\t{}\t{}\t{}\t{}\n".format(
ind1, ind2, ind3, ind4, ft))
atom14.append((ind1, ind4) )
self.out_file.write("[ pairs ]\n")
self.out_file.write("; ai aj\n")
for i, j in atom14:
self.out_file.write("{}\t{}\n".format(i, j))
def get_atoms(self, atom_types, ft=1, resnr=1, name="etan"):
self.out_file.write("[ atoms ]\n")
counter = defaultdict(lambda : 0)
for a in self.m3d.GetAtoms():
sym = a.GetSymbol()
counter[sym] += 1
ind = a.GetIdx() + 1
self.out_file.write("{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\n".format(
ind, self.atom_types[sym], resnr, name, sym + str(counter[sym]), ind,\
a.GetFormalCharge(), a.GetMass()))
In [13]:
top_file_name = "et.top"
a = M3dTop("CC")
a.write(top_file_name)
In [8]:
from subprocess import Popen, PIPE
import re
import os
In [11]:
def run_field(gro_file_name, top_file_name, mdp):
print "grompp -f {0}.mdp -c et.gro -p {1} -o et_{0}.tpr".format(mdp, top_file_name)
proc = Popen(
"grompp -f {0}.mdp -c et.gro -p {1} -o et_{0}.tpr".format(mdp, top_file_name),
shell=True,
stdout=PIPE, stderr=PIPE
)
_ = proc.communicate()
print "done"
print "mdrun -deffnm et_{0}.tpr -v -nt 1".format(mdp)
proc = Popen(
"mdrun -deffnm et_{0}.tpr -v -nt 1".format(mdp),
shell=True,
stdout=PIPE, stderr=PIPE
)
_ = proc.communicate()
print "done"
print "echo '2\\n' | trjconv -f et_{0}.tpr.trr -s et_{0}.tpr -o et_{0}.pdb".format(mdp)
#visual
proc = Popen(
"echo '2\\n' | trjconv -f et_{0}.tpr.trr -s et_{0}.tpr -o et_{0}.pdb".format(mdp),
shell=True,
stdout=PIPE, stderr=PIPE
)
_ = proc.communicate()
print "done"
print 'echo -e "10\\n11\\n\\n" | g_energy -f et_{0}.tpr.edr -o et_{0}_en.xvg -xvg none'.format(mdp)
#potential and kinetic
proc = Popen(
'echo -e "10\\n11\\n\\n" | g_energy -f et_{0}.tpr.edr -o et_{0}_en.xvg -xvg none'.format(mdp),
shell=True,
stdout=PIPE, stderr=PIPE
)
info = proc.communicate()[0]
result = []
for line in info.splitlines():
if re.match("Potential", line):
result.append(line)
elif re.match("Kinetic En.", line):
result.append(line)
return result
Запустим молекулярную динамику с 5 параметрами контроля температуры
be.mdp - метод Берендсена для контроля температуры.
vr.mdp - метод "Velocity rescale" для контроля температуры.
nh.mdp - метод Нуза-Хувера для контроля температуры.
an.mdp - метод Андерсена для контроля температуры.
sd.mdp - метод стохастической молекулярной динамики
In [14]:
run_field("etane.gro", "et.top", "be")
Out[14]:
In [15]:
run_field("etane.gro", "et.top", "vr")
Out[15]:
In [16]:
run_field("etane.gro", "et.top", "nh")
Out[16]:
In [17]:
run_field("etane.gro", "et.top", "an")
Out[17]:
In [19]:
run_field("etane.gro", "et.top", "sd")
Out[19]:
In [22]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
for filename in ['be','vr','nh','an','sd']:
energy_data = np.loadtxt('et_{}_en.xvg'.format(filename))
x = energy_data[:, 0]
pot = energy_data[:, 1]
kin = energy_data[:, 2]
plt.plot(x, pot, "b-", x, kin, "r-" )
red_patch = mpatches.Patch(color='red', label='Kinetic')
blue_patch = mpatches.Patch(color='blue', label='Potential')
plt.legend(handles=[red_patch,blue_patch])
plt.title('Energies of {} system'.format(filename))
plt.show()
Рассмотрим распределение длинны связи С-С за время моделирования
Сначала создадим индекс файл с одной связью. В текстовом редакторе создадим файл b.ndx со следующим содержимым:
[ b ]
1 2
In [24]:
for meth in ['be','vr','nh','an','sd']:
proc = Popen(
'g_bond -f et_{0}.tpr.trr -s et_{0}.tpr -o bond_{0}.xvg -n b.ndx -xvg none'.format(meth),
shell=True,
stdout=PIPE, stderr=PIPE
)
proc.communicate()[0]
In [29]:
for meth in ['be','vr','nh','an','sd']:
bond_data = np.loadtxt('bond_{}.xvg'.format(meth))
x =bond_data[:,0]
y =bond_data[:,1]
width = 0.0001
plt.bar(x,y,width, color="blue")
plt.title('Length of C-C bond in {} system'.format(meth))
plt.show()
Сравним наши распределения с распределением Максвелла
In [32]:
from IPython.display import Image
Image(filename='mack.jpg')
Out[32]:
Видим, что более-менее похоже на распределение Максвелла распределение для метода Берендсена, и распределение для метода Нуза-Хувера
У распределения Максвелла правый хвост тяжелее, но не одно из распределений этого не отражает
Сравним время работы методов
In [35]:
from timeit import default_timer as timer
In [42]:
def test_time(gro_file_name, top_file_name, mdp):
start = timer()
#print "grompp -f {0}.mdp -c et.gro -p {1} -o et_{0}.tpr".format(mdp, top_file_name)
proc = Popen(
"grompp -f {0}.mdp -c et.gro -p {1} -o et_{0}.tpr".format(mdp, top_file_name),
shell=True,
stdout=PIPE, stderr=PIPE
)
_ = proc.communicate()
#print "done"
#print "mdrun -deffnm et_{0}.tpr -v -nt 1".format(mdp)
proc = Popen(
"mdrun -deffnm et_{0}.tpr -v -nt 1".format(mdp),
shell=True,
stdout=PIPE, stderr=PIPE
)
print mdp, timer() - start
In [43]:
for meth in ['be','vr','nh','an','sd']:
test_time("etane.gro", "et.top", meth)
Видим, что быстрее всего отработал sd, медленнее всего - an