Запустим PyMol для работы в режиме сервера.
from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys, time
# pymol launching
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
# __main__.pymol_argv = [ 'pymol', '-cp' ] # to disable GUI
import pymol
from pymol import cmd
pymol.finish_launching()
from IPython.display import Image
defaultImage = '/tmp/pymolimg.png'
def prepareImage(width=300, height=300, sleep=2, filename=defaultImage):
## To save the rendered image
cmd.ray(width, height)
cmd.png('/tmp/pymolimg.png')
time.sleep(sleep)
# Define some shortcuts
def focus(x):
cmd.center(x)
cmd.zoom(x)
def showPDB(pdb):
cmd.do('''
reinit
bg white
set antialias, 2
set ray_trace_mode, 3
load ''' + pdb + '''
show spheres, all
set sphere_scale, 0.2
show sticks, all
set stick_radius, 0.1
set valence, on
zoom all
''')
GAMESS — стандартный квантово-химический пакет. С его помощью найдём оптимальную геометрию для нафталена и азулена и рассчитаем теплоты образования этих молекул разными подходами квантовой механики.
Построение и оптимизация структур нафталена и азулена
На основе SMILES-аннотаций нафталена и азулена создадим их пространственные структуры. Затем оптимизируем с помощью MOPAC и загрузим полученные структуры в PyMol.
%%bash
echo "C1=CC=C2C=CC=C2C=C1 azulene" > azulene.smi
obgen azulene.smi > azulene.mol
echo "c1ccc2ccccc2c1 naphthalene" > naphthalene.smi
obgen naphthalene.smi > naphthalene.mol
%%bash
babel -imol azulene.mol -omop azu.mop -xk "PM6"
babel -imol naphthalene.mol -omop nap.mop -xk "PM6"
MOPAC2009.exe azu.mop
MOPAC2009.exe nap.mop
babel -imopout azu.out -opdb azu.pdb
babel -imopout nap.out -opdb nap.pdb
showPDB('azu.pdb')
cmd.rotate('y', '-30')
prepareImage(); Image(defaultImage)
showPDB('nap.pdb')
cmd.rotate('y', '-30')
prepareImage(); Image(defaultImage)
Как это и ожидалось, ароматическая система нафталена имеет плоскую геометрию.
Молекула азулена, однако, не плоская. Попробуем использовать, например, силовое поле UFF в obgen.
%%bash
obgen azulene.smi -ff UFF > azulene_UFF.mol
babel -imol azulene_UFF.mol -omop azu_UFF.mop -xk "PM6"
MOPAC2009.exe azu_UFF.mop
babel -imopout azu_UFF.out -opdb azu_UFF.pdb
showPDB('azu_UFF.pdb')
cmd.rotate('y', '-60')
prepareImage(); Image(defaultImage)
Использование силового поля UFF в obgen позволило получить плоскую ароматическую систему структуры нафталена.
Оптимизация геометрии средствами GAMESS
Чтобы получить входные файлы для GAMESS, переформатируем координаты, записанные в файлах nap.out и azu.out в gamin формат.
%%bash
babel -imopout azu_UFF.out -ogamin azu.gamin
babel -imopout nap.out -ogamin nap.gamin
%%bash
HEADER=" \$CONTRL COORD=CART UNITS=ANGS SCFTYP=RHF RUNTYP=OPTIMIZE \$END\n\
\$BASIS GBASIS=N31 NGAUSS=6 \$end\n\
\$system mwords=2 \$end\n\
\$DATA"
echo -en " " > azu_opt.inp
echo -e $HEADER >> azu_opt.inp
tail -n+4 azu.gamin >> azu_opt.inp
echo -en " " >> nap_opt.inp
echo -e $HEADER > nap_opt.inp
tail -n+4 nap.gamin >> nap_opt.inp
Проведём оптимизацию геометрии для обоих молекул средствами GAMESS.
%%bash
gms nap_opt.inp 1 >& nap_opt.log
gms azu_opt.inp 1 >& azu_opt.log
# NOTE: `1` is a number of cores to use
Отобразим полученные структуры с помощью PyMol.
%%bash
babel -igamout azu_opt.log -opdb azu_gamess.pdb
babel -igamout nap_opt.log -opdb nap_gamess.pdb
Структуры азулена и нафталена имеют плоскую геометрию.
# Azulene, en face
showPDB("azu_gamess.pdb")
prepareImage(); Image(defaultImage)
# Azulene, profile
cmd.rotate('y', '-90')
prepareImage(); Image(defaultImage)
# Naphthalene, en face
showPDB("nap_gamess.pdb")
prepareImage(); Image(defaultImage)
# Naphthalene, profile
cmd.rotate('y', '-90')
prepareImage(); Image(defaultImage)
В результате оптимизации в структуре азулена незначительно изменились длины двух связей между углеродами пятичленного кольца (с 1.5 Å до 1.4 Å). В структуре нафталена длины связей после выполнения оптимизации не изменились.
Как следует из приведённого выше кода, оптимизация геометрии проводилась с базисом 6-31G. В документации к GAMESS он описан как Pople's N-31G split valence basis set, а N — число гауссианов — было выбрано равным 6. Использование такого базиса означает, что орбитали остова (core) представлены шестью гауссианами, а валентные орбитали представлены двумя составляющими — из трёх гауссианов и из одного гауссиана.
Расчёт по Хартри-Фоку и DFT
На основе полученных координат составим новые входные файлы для расчёта энергии: для расчёта методом Хартри-Фока и для использования теории функционала плотности.
Расчёт по Хартри-Фоку
Переформатируем файлы для расчёта по Хартри-Фоку (Hartree-Fock).
%%bash
babel -igamout azu_opt.log -ogamin azu_gamout.inp
babel -igamout nap_opt.log -ogamin nap_gamout.inp
HF_HEADER="\$CONTRL COORD=CART UNITS=ANGS SCFTYP=RHF RUNTYP=ENERGY \$END\n\
\$BASIS GBASIS=N31 NGAUSS=6\n\
POLAR=POPN31 NDFUNC=1 \$END\n\
\$GUESS GUESS=HUCKEL \$END\n\
\$system mwords=2 \$end\n\
\$DATA"
echo -en " " > azu_HF.inp
echo -e $HF_HEADER >> azu_HF.inp
tail -n+4 azu_gamout.inp >> azu_HF.inp
echo -en " " > nap_HF.inp
echo -e $HF_HEADER >> nap_HF.inp
tail -n+4 nap_gamout.inp >> nap_HF.inp
Запустим GAMESS для расчёта энергии.
%%bash
gms nap_HF.inp 1 >& nap_HF.log
gms azu_HF.inp 1 >& azu_HF.log
Выпишем значения энергии.
%%bash
grep "TOTAL ENERGY =" azu_HF.log nap_HF.log
Расчёт с использованием теории функционала плотности (DFT)
Аналогичным образом проведём расчёт энергии с использованием теории функционала плотности (Density functional theory).
%%bash
DFT_HEADER="\$CONTRL COORD=CART UNITS=ANGS dfttyp=b3lyp RUNTYP=ENERGY \$END\n\
\$BASIS GBASIS=N31 NGAUSS=6\n\
POLAR=POPN31 NDFUNC=1 \$END\n\
\$GUESS GUESS=HUCKEL \$END\n\
\$system mwords=2 \$end\n\
\$DATA"
echo -en " " > azu_DFT.inp
echo -e $DFT_HEADER >> azu_DFT.inp
tail -n+4 azu_gamout.inp >> azu_DFT.inp
echo -en " " > nap_DFT.inp
echo -e $DFT_HEADER >> nap_DFT.inp
tail -n+4 nap_gamout.inp >> nap_DFT.inp
%%bash
gms nap_DFT.inp 1 >& nap_DFT.log
gms azu_DFT.inp 1 >& azu_DFT.log
grep "TOTAL ENERGY =" azu_DFT.log nap_DFT.log
Сравнение методов
На основе результатов расчётов составим таблицу.
| Вещество | ENaphthalene | EAzulene | ΔE, Hartree | ΔE, kCal/mol |
|---|---|---|---|---|
| Хартри-Фок | -383.35 | -383.28 | 0.07 | 45.62 |
| DFT | -385.64 | -385.59 | 0.05 | 34.32 |
Сравнение с экспериментальным значением энергии изомеризации нафталина в азулен, которое равно 35.3±2.2 kCal/mol, показывает, что более точного результата позволила добиться теория функционала плотности.
Визуализация орбиталей в Gabedit
Для азулена и нафталена с помощью программы Gabedit выполним визуализацию орбиталей HOMO (highest occupied molecular orbital) и LUMO (lowest unoccupied molecular orbital).
Image('azu_homo.png')
Image('azu_lumo.png')
Image('nap_homo.png')
Image('nap_lumo.png')
Ссылки на файлы
Для азулена:
- SMILES,
- результат
obgen(.mol), - результат
babel(.mop), - результат
MOPAC(.out), результат
MOPAC(.pdb);//
результат
obgen(.mol, силовое поле UFF),- результат
babel(.mop, силовое поле UFF), - результат
MOPAC(.out, силовое поле UFF), результат
MOPAC(.pdb, силовое поле UFF);//
результат
GAMESS(.log),результат
GAMESS(.pdb);расчёт энергии по HF (.log),
- расчёт энергии по DFT (.log).
Для нафталина:
- SMILES,
- результат
obgen(.mol), - результат
babel(.mop), - результат
MOPAC(.out), результат
MOPAC(.pdb);//
результат
GAMESS(.log),результат
GAMESS(.pdb);расчёт энергии по HF (.log),
- расчёт энергии по DFT (.log).