Часть 1. Cеми-эмпирические вычисления: Mopac

In [33]:
import re
import pybel
import subprocess, os
from IPython.display import Image
import pandas as pd
In [16]:
mol=pybel.readstring('smi','C1=CC2=CC3=CC=C(N3)C=C4C=CC(=N4)C=C5C=CC(=N5)C=C1N2')
mol.addh()
mol.make3D(steps=1000)
mol.write(format='pdb',filename='mol.pdb',overwrite=True)
mol
Out[16]:
H N N H H H H H H H H H H H H H N N
In [4]:
# for mopac
def creat_file (mol, file_name, param, charge):
    mop = mol.write(format='mopin', filename='%s.mop' % file_name, opt={'k':'%s CHARGE=%s XYZ' % (param, charge)}, overwrite=True)
    os.system("export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/'")
    os.system("echo | /home/preps/golovin/progs/mopac/MOPAC2016.exe %s.mop" % file_name)
    # считаем геометрию
    opt=pybel.readfile('mopout','%s.out' % file_name)
    for i in opt:
        i.write(format='pdb',filename='%s.pdb' % file_name, overwrite=True)
In [16]:
creat_file (mol, "mol_pm6", "PM6", mol.charge)
creat_file (mol, "mol_am1", "AM1", mol.charge)

c1c/c/2=C/c3cc/c(=C/c4ccc([n]4)/C=C\4/C=CC(=N4)/C=c/1\[nH]2)/[nH]3	mol_pm6.out

c1c/c/2=C/c3cc/c(=C/c4ccc([n]4)/C=C\4/C=CC(=N4)/C=c/1\[nH]2)/[nH]3	mol_am1.out
In [2]:
Image(filename='mol.png')
Out[2]:
In [4]:
Image(filename='mol_am1.png')
Out[4]:
In [6]:
Image(filename='mol_pm6.png')
Out[6]:

Выше представлены структуры порфирина, полученные pybel и двумя способами PM6 и AM1 в Mopac. Структура, полученная при помощи PM6 является наиболее корректной, потому что кольцо порфирина в ней плоское.

Рассчитаем возбужденные состояния порфирина
In [27]:
os.system("cp mol_pm6.mop mol_spectr.mop")
os.system('echo "" >> mol_spectr.mop')
os.system('echo "cis c.i.=4 meci oldgeo" >> mol_spectr.mop')
os.system('echo "some description" >> mol_spectr.mop')
os.system("export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/'")
os.system("echo | /home/preps/golovin/progs/mopac/MOPAC2016.exe mol_spectr.mop")
Out[27]:
0
In [11]:
f = open('mol_spectr.out', 'r')
lines = f.readlines()[1515:1523]
f.close()
energies = []
for line in lines:
    energies.append(float(line.split()[1]))
print(energies)

[1.765672, 2.170938, 2.327822, 2.693383, 3.091472, 3.143967, 3.756628, 3.7652]

На основании этих значений и простой формулы были расчитаны длины волн при которых происходят эти переходы.

In [13]:
wavelength = list(map(lambda energy: 1239.84193 / energy, energies))
print(wavelength)

[702.1926665881319, 571.1088616994128, 532.6188729206959, 460.32886150985587, 401.05229159442496, 394.3558981376077, 330.0411778861255, 329.2897933708701]

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

In [34]:
creat_file (pybel.readfile("pdb", "td.pdb").next(), "td_charge_0", "PM6", 0)
creat_file (pybel.readfile("pdb", "td_charge_0.pdb").next(), "td_charge_2", "PM6", 2)
creat_file (pybel.readfile("pdb", "td_charge_2.pdb").next(), "td_charge_again_0", "PM6", 0)
In [35]:
%%bash
echo dimer
grep 'TOTAL ENERGY' td_charge_0.out
echo semi-dimer
grep 'TOTAL ENERGY' td_charge_2.out
echo two thymins
grep 'TOTAL ENERGY' td_charge_again_0.out

dimer
          TOTAL ENERGY            =      -3273.57701 EV
semi-dimer
          TOTAL ENERGY            =      -3253.90535 EV
two thymins
          TOTAL ENERGY            =      -3273.68570 EV

Тимидиновый димер

In [37]:
Image(filename='td.png')
Out[37]:

Тимидиновый димер после оптимизации при заряде 0

In [39]:
Image(filename='td_charge_0.png')
Out[39]:

Тимидиновый димер после оптимизации при заряде 2

In [40]:
Image(filename='td_charge_2.png')
Out[40]:

Тимины после повторной оптимизации при заряде 0

In [42]:
Image(filename='td_charge_again_0.png')
Out[42]:

При возбуждении системы тиминовый димер перешел в более энергетически невыгодную конформацию, которая стала переходным состоянием между тимининовым димером и двумя тиминами. Два тимина обладают немного меньшей энергией, по сравнению с тимидиновым димером, поэтому возбуждённое состояние при добавлении электронов не перешло обратно в исходно состояние.

Часть 2. Ab inito вычисления: ORCA

Найдём оптимальную геометрию для нафталена и азулена и рассчитаем теплоты образования этих молекул разными подходами квантовой механики.

In [10]:
nap=pybel.readstring('smi','C1=CC=C2C=CC=CC2=C1')
nap.addh()
nap.make3D(steps=1000)
nap.write(format='pdb',filename='nap.pdb', overwrite=True)
nap
Out[10]:
H H H H H H H H
In [11]:
azu=pybel.readstring('smi','C1=CC=C2C=CC=C2C=C1')
azu.addh()
azu.make3D(steps=1000)
azu.write(format='pdb',filename='azu.pdb', overwrite=True)
azu
Out[11]:
H H H H H H H H
In [12]:
creat_file(nap, "nap_opt", "PM6", nap.charge)
creat_file(azu, "azu_opt", "PM6", azu.charge)

Оптимизируем с помощью MOPAC структуры нафталена и азулена

Нафталин до оптимизации

In [8]:
Image(filename='nap.png')
Out[8]:

Нафталин после оптимизации

In [9]:
Image(filename='nap_opt.png')
Out[9]:

Видно, что pybel с задачей построения нафталина не справился. А азулен после оптимизации практически не изменился.

Азулен до оптимизации

In [3]:
Image(filename='azu.png')
Out[3]:

Азулен после оптимизации

In [4]:
Image(filename='azu_opt.png')
Out[4]:
In [35]:
os.system('echo "C1=CC=C2C=CC=CC2=C1" > nap.smi')
os.system('obgen nap.smi > nap.mol')
Out[35]:
0
In [8]:
nap = pybel.readfile('sdf', 'nap.mol').next()
nap.write(format='pdb',filename='nap_.pdb', overwrite=True)
nap
Out[8]:
H H H H H H H H
In [9]:
creat_file(nap, "nap_opt", "PM6", nap.charge)

Исправленный нафталин

In [10]:
Image(filename='nap_opt_new.png')
Out[10]:
In [25]:
def orca(infile, file_name, param):
    mol = pybel.readfile('pdb', '%s.pdb' % infile).next()
    mol.write(format='orcainp',filename='nap.inp', overwrite=True)
    orca = mol.write(format='orcainp',opt={'k': param}, overwrite=True)
    orca_out = orca.replace('! insert inline commands here ', param)
    outfile = open('%s.inp' % file_name, 'w')
    outfile.write(orca_out)
    outfile.close()
    os.system('/srv/databases/orca/orca %s.inp | tee %s.log' % (file_name, file_name))
In [26]:
orca('nap_opt', 'nap_hf', '!HF RHF 6-31G')
In [27]:
orca('nap_opt', 'nap_dft', '!DFT RHF 6-31G')
In [28]:
orca('azu_opt', 'azu_hf', '!HF RHF 6-31G')
In [29]:
orca('azu_opt', 'azu_dft', '!DFT RHF 6-31G')
In [31]:
%%bash
grep "Total Energy" *.log

azu_dft.log:Total Energy       :         -382.08969469 Eh          -10397.18918 eV
azu_hf.log:Total Energy       :         -383.03184550 Eh          -10422.82640 eV
nap_dft.log:Total Energy       :         -382.25559313 Eh          -10401.70350 eV
nap_hf.log:Total Energy       :         -383.22018089 Eh          -10427.95127 eV
In [36]:
Image(filename='table.png')
Out[36]:

Если из эксперимента известно, что энергия изомеризации нафталина в азулен составляет 35.3±2.2 kCal/mol, то метод DFT справился немного лучше, потому что значение, полученное с помощью этого метода, немного ближе к 35.3±2.2 kCal/mol. Однако, полученные значения довольно сильно отличются от экспериментальных, что ставит под сомнение корректность вывода.