Полуэмпирические и Ab initio вычисления молекулярных систем
Часть 1. Полуэмпирические вычисления в Mopac
Задача - сравнить структуры порфирина, полученные с помощью pybel и двумя способами (PM6 и AM1) в Mopac.
import pybel
from os import system
from IPython.display import display, Image
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=500)
mol.write(format='pdb', filename='porph_pybel.pdb', overwrite=True)
mol
Pybel не смог понять, как на самом деле должна выглядеть структура порфирина (до плоскости и ароматичности всех 4 циклов далековато), поэтому ее необходимо оптимизировать с помощью Mopac.
def do_mopac(mol, name, ham="PM6", charge=False, option=''):
command = "export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/'\n"
if type(charge) == int or type(charge) == str:
ch = charge
else:
ch = mol.charge
mop = mol.write(format='mopin', filename='%s.mop' % name,
opt={'k':'%s PRECISE EF CHARGE=%s %s' % (ham, ch, option)},
overwrite=True)
system("echo | /home/preps/golovin/progs/mopac/MOPAC2016.exe %s.mop" % name)
opt = pybel.readfile('mopout','%s.out' % name)
for c,i in enumerate(opt):
i.write(format='pdb', filename='%s_%d.pdb' % (name,c), overwrite=True)
Сначала оптимизируем с использованием гамильтонианов PM6.
do_mopac(mol, "mopac_PM6_2")
 |  |
Рис. 2. Сравнение структур порфирина, полученных pybel (цветом salmon) и Mopac с использованием метода PM6 (оранжевым).
Почему-то Морас не справился с уплощением порфирина, однако все 4 кольца стали ароматичными - в отличие от структуры, сгенерированной pybel.
do_mopac(mol, "mopac_AM1", ham="AM1")
 |  |
Рис. 3. Слева: сравнение структур порфирина, полученных pybel (цветом salmon) и Mopac с использованием метода AM1 (голубым). Справа: сравнение структур порфирина, полученных Mopac методами PM6 (оранжевым цветом) и АМ1 (голубым).
С помощью метода PM6 удалось добиться идеальной структуры порфирина - плоской, с 4 ароматическими циклами. Морас с гамильтонианами АМ1 тоже справился с задачей неплохо, хотя структура получилась не совсем плоской (Рис. 3).
Далее нужно было рассчитать возбужденные состояния порфирина и оценить его спектр поглощения.
Для этого была сделана копия mop файла, и в его конец для расчета возбужденных состояний были добавлены 2 строчки:
%%bash
cp mopac_PM6_2.mop mopac_PM6_exited.mop
echo '' >> mopac_PM6_exited.mop
echo 'cis c.i.=4 meci oldgeo' >> mopac_PM6_exited.mop
echo 'some description' >> mopac_PM6_exited.mop
%%bash
export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/'
echo | /home/preps/golovin/progs/mopac/MOPAC2016.exe mopac_PM6_exited.mop
por_energ = [1.763459, 2.168394, 2.323063, 2.690530, 3.087475, 3.139846, 3.755972, 3.765160]
def wavelen(x):
wl = 1239.842 / x #1239.842 = постоянная Планка, эВ*с, * скорость света, м/с - для перевода длины волны в нм
return (wl)
por_wvlen = list(map(wavelen, por_energ))
print 'State\tWavelength, nm'
for c in range(len(por_wvlen)):
print '%s\t%.2f' % (c+2,por_wvlen[c])
В полученном списке для 8 возбужденных состояний с известными энергиями приведены длины волн, при которых происходят соответствующие переходы.
В следующем задании нужно было пронаблюдать переход тиминов из димерного состояния в обычное.
! wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/td.pdb
Исходный тиминовый димер выглядит как на Рис.4.
Далее проведем последовательную оптимизацию системы с разными зарядами: 0, +2 и снова 0.
td_0 = pybel.readfile("pdb", "td.pdb").next()
do_mopac(td_0, "thymine_1", ham='PM6', charge='0', option='XYZ')
td_1 = pybel.readfile("pdb", "thymine_1_0.pdb").next()
do_mopac(td_1, "thymine_2", ham='PM6', charge='+2', option='XYZ')
td_2 = pybel.readfile("pdb", "thymine_2_0.pdb").next()
do_mopac(td_2, "thymine_3", ham='PM6', charge='0', option='XYZ')
В возбужденном состоянии при заряде системы +2 рвется одна С-С связь, и ставшие неспаренными атомы С становятся карбокатионами -[CH+]-. При оптимизации такой системы при заряде 0 димер распадается на два тимина вместо того, чтобы вернуться в димерную форму с двумя связями.
%%bash
grep "TOTAL ENERGY" thymine_*.out >> td.log
with open("td.log", "r") as f:
for num, val in enumerate(f):
print "State %s " % str(num+1), val.split()[-2], "eV"
Невозбужденные состояния (1 и 3) близки по энергиям, но в состоянии 3 - отдельные мономеры тиминов - энергия немного более выгодна, и из-за этого в это состояние сваливается система из возбужденного состояния при сообщении достаточного количества энергии.
Часть 2. Ab initio вычисления: Orca
В этой части практикума надо найти оптимальную геометрию нафталина и азулена и теплоты их образования.