In [1]:
import pybel
from subprocess import Popen, PIPE
In [2]:
import threading
import nglview
In [3]:
def create_3d(form, meth = "PM6", mop_name = "mop.mop", pdf_file_name = "temp.pdb"):
    proc = Popen(
        "export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/'",
        shell=True,
        stdout=PIPE, stderr=PIPE
    )
    proc.wait()
    
    mol=pybel.readstring('smi', form)
    mol.addh()
    mol.make3D(steps=10000)
    mop=mol.write(format='mopin',filename=mop_name, opt={'k':'{} CHARGE={}'.format(meth, mol.charge)}, overwrite=True)
    proc = Popen(
        "/home/preps/golovin/progs/mopac/MOPAC2016.exe {}".format(mop_name),
        shell=True,
        stdout=PIPE, stderr=PIPE
    )
    _ = proc.communicate()
    
    opt=pybel.readfile('mopout', mop_name[:-3] + "out")
    mol = opt.next()
    print pdf_file_name
    mol.write(format='pdb',filename= pdf_file_name, overwrite=True)
    
    
    return mol
In [4]:
class MolThread(threading.Thread):   
    def __init__(self, **kwargs):
        self.kwargs = kwargs
        threading.Thread.__init__(self)

    def run(self):  
        self.mol = create_3d(**self.kwargs)
In [5]:
from IPython import display
In [53]:
processes = list()
for meth in ("PM6", "AM1"):
    e = MolThread(**{"form" : "C1=CC2=CC3=CC=C(N3)C=C4C=CC(=N4)C=C5C=CC(=N5)C=C1N2",
                                   "meth" : meth, 
                                   "mop_name" : "hem_{}_mop.mop".format(meth),
                                   "pdf_file_name" : "hem_{}.pdb".format(meth)
                                }
                    )
    e.start()
    processes.append(e)

    
    
for proc in processes:
    proc.join()

hem_PM6.pdb
hem_AM1.pdb
In [54]:
for proc in processes:
    display.display(proc.mol)
H N N H H H H H H H H H H H H H N . N
H N N H H H H H H H H H H H H H N . N
In [1]:
from IPython.display import Image

Сравним изображение, полученное при использовании AM1:

In [3]:
Image(filename="hem_AM1.png")
Out[3]:

И использовании PM6

In [4]:
Image(filename="hem_PM6.png")
Out[4]:

Видим, что в случае PM6 структура получилось плоской, что корректно.

Построим приблизительный спектр поглощения молекулы, используя формулу $E=hv$

In [6]:
%%bash
cp hem_PM6_mop.mop hem_PM6_mop_mod.mop
echo "" >> hem_PM6_mop_mod.mop
echo "cis c.i.=4 meci oldgeo" >> hem_PM6_mop_mod.mop
echo "some description" >> hem_PM6_mop_mod.mop
In [9]:
%%bash
export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/'
/home/preps/golovin/progs/mopac/MOPAC2016.exe hem_PM6_mop_mod.mop


          MOPAC Job: "hem_PM6_mop_mod.mop" ended normally on May 24, 2017, at 14:47.
In [19]:
with open('hem_PM6_mop_mod.out', 'r') as in_file:
    start = False
    end = False
    energies = []
    for line in in_file:
        if "ABSOLUTE     RELATIVE                            X       Y       Z" in line:
            start = True
        elif "The \"+\" symbol indicates the root used" in line:
            end = True
        elif start and not end and line.strip():
            energies.append(line)
    energies = [float(energy_info.split()[1]) for energy_info in energies]
In [26]:
energies = energies[1:]
In [33]:
lengths = [1239.84193/e for e in energies]
import matplotlib.pyplot as plt
plt.figure(figsize=(10,10))
plt.bar(lengths, energies, align='center', alpha=0.5)
plt.xticks(lengths, ["{:.00f}".format(x) for x in lengths])
plt.ylabel('Energy')
plt.title('Spectrum')
plt.xticks(rotation=90)
plt.show()

Попробуем теперь промоделировать переход тиминового димера в отдельные тимины при возбуждении

In [38]:
td_file_name= "td.pdb"
init_file_name = "td_init.mop"
middle_file_name = "td_middle.mop"
end_file_name = "td_end.mop"
init_out_file_name = "td_init.out"
middle_out_file_name = "td_middle.out"
end_out_file_name = "td_end.out"
init_pdb_file = "td_init.pdb"
middle_pdb_file = "td_middle.pdb"
end_pdb_file = "td_end.pdb"

Изначальная структура

In [41]:
Image(filename="td.png")
Out[41]:

Оптимизируем изначальную структуру

In [39]:
td = pybel.readfile("pdb", td_file_name).next()
td.write(format='mop',filename=init_file_name,opt={"k":"PM6 CHARGE={}".format(0)},overwrite=True)
In [27]:
proc = Popen(
    "/home/preps/golovin/progs/mopac/MOPAC2016.exe {}".format(init_file_name),
    shell=True,
    stdout=PIPE, stderr=PIPE
)

_ = proc.communicate()
td = pybel.readfile("mopout", init_out_file_name).next()
td.write(format="pdb", filename=init_pdb_file, overwrite=True)

Ничего особо не изменилось

In [43]:
Image("td_init.png")
Out[43]:

Имитируем возбуждение системы, ставя ее заряд равным $+2$

In [29]:
td.write(format="mop",filename=middle_file_name,opt={"k":"PM6 CHARGE={}".format(2)},overwrite=True)
proc = Popen(
    "/home/preps/golovin/progs/mopac/MOPAC2016.exe {}".format(middle_file_name),
    shell=True,
    stdout=PIPE, stderr=PIPE
)

_ = proc.communicate()

td = pybel.readfile("mopout", middle_out_file_name).next()
td.write(format="pdb", filename=middle_pdb_file, overwrite=True)

После оптимизации получаем тимины, соединенные уже одной связью

In [45]:
Image("td_middle.png")
Out[45]:

Ставим заряд системы опять 0 и оптимизируем

In [31]:
td.write(format='mop',filename=end_file_name,opt={"k":"PM6 CHARGE={}".format(0)},overwrite=True)
proc = Popen(
    "/home/preps/golovin/progs/mopac/MOPAC2016.exe {}".format(end_file_name),
    shell=True,
    stdout=PIPE, stderr=PIPE
)

td = pybel.readfile("mopout", end_out_file_name).next()
td.write(format="pdb", filename=end_pdb_file, overwrite=True)

Получаем отдельные тимины

In [46]:
Image("td_end.png")
Out[46]:

Смотрим на энергию каждого состояния, указанную в файлах

In [51]:
%%bash
echo "Init state"
grep "TOTAL ENERGY" "td_init.out"
echo "Middle state"
grep "TOTAL ENERGY" "td_middle.out"
echo "End state"
grep "TOTAL ENERGY" "td_end.out"

Init state
          TOTAL ENERGY            =      -3273.57488 EV
Middle state
          TOTAL ENERGY            =      -3253.90535 EV
End state
          TOTAL ENERGY            =      -3273.69848 EV

Видим, что состояния почти эквивалентны и после возбуждения наша молекула может перейти как в состояние димера, так и в состояниие двух отдельных тиминов