Полуэмпирические и Ab initio вычисления молекулярных систем
Цель: Сравнение структуры порфирина полученные pybel и двумя способами PM6 и AM1 в Mopac
In [1]:
#загружаем модуль
import pybel
In [2]:
#Создадем соединение и сразу сделаем его 3D:
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)
In [3]:
mol.write(format='pdb',filename='pro_1.pdb',overwrite=True)
PM6 и AM1
In [4]:
import subprocess
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)
#Смотрим, как оно может записывать формат программы сэмиэмпирических расчётов MOPAC и записываем в файл
mop_pm6=mol.write(format='mopin',opt={'k':'PM6 CHARGE=%d' % mol.charge})
mop_pm6=mol.write(format='mopin',filename='pm6.mop',opt={'k':'PM6 CHARGE=%d' % mol.charge},overwrite=True)
print mop_pm6
mol.write(format='pdb',filename='pm6_before.pdb',overwrite=True)
# запускаем MOPAC
cmb1 = "export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/'"
cmb2 = "/home/preps/golovin/progs/mopac/MOPAC2016.exe pm6.mop"
subprocess.call(cmb1, shell=True)
subprocess.call(cmb2, shell=True)
#Считаем геометрию
opt=pybel.readfile('mopout', 'pm6.out')
for i in opt:
print i
i.write(format='pdb',filename='pm6.pdb',overwrite=True)
In [5]:
import subprocess
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)
mop_am1=mol.write(format='mopin',opt={'k':'AM1 CHARGE=%d' % mol.charge})
mop_am1=mol.write(format='mopin',filename='am1.mop',opt={'k':'AM1 CHARGE=%d' % mol.charge},overwrite=True)
print mop_am1
mol.write(format='pdb',filename='am1_before.pdb',overwrite=True)
cmb1 = "export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/'"
cmb2 = "/home/preps/golovin/progs/mopac/MOPAC2016.exe am1.mop"
subprocess.call(cmb1, shell=True)
subprocess.call(cmb2, shell=True)
opt=pybel.readfile('mopout', 'am1.out')
for i in opt:
print i
i.write(format='pdb',filename='am1.pdb',overwrite=True)
In [6]:
from IPython.display import Image
Image(filename='all_done.png')
Out[6]:
На рисунке 1 изображена структура порфирина, полученная pybel, на рисунке 2 - AM1 и на рисунке 3 - PM6 Порфирин получается изогнутый, остальные методы подтверждают эту структуру.
Рассчитываем возбужденные состояния порфирина и на основе этих данных строим спектр поглощения молекулы.
In [7]:
%%bash
cp my.mop my_spectr.mop
echo "" >> my_spectr.mop
echo "cis c.i.=4 meci oldgeo" >> my_spectr.mop
echo "some description" >> my_spectr.mop
export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/'
/home/preps/golovin/progs/mopac/MOPAC2016.exe my_spectr.mop
In [8]:
opt=pybel.readfile('mopout','my_spectr.out')
for i in opt:
print i
i.write(format='pdb',filename='my_spectr.pdb',overwrite=True)
Получены значения энергий для электронных переходов
In [1]:
f = open('my_spectr.out', 'r')
lines = f.readlines()[-21:-13]
f.close()
energies = [float(line.split()[1]) for line in lines]
print energies
На основании этих значений и простой формулы были рассчитаны длины волн и их частоты, при которых происходят эти переходы.
In [12]:
wavelengths = [1239.84193/e for e in energies]
frequencies = [e/4.135667516 for e in energies]
print 'frequncy,THz\t wavelength,nm\tenergy,eV'
for i in range(len(energies)):
print '%.4f \t %.4f \t %.4f' % (frequencies[i]*1000, wavelengths[i], energies[i])