запускаем 'ipython notebook' и 'pymol -R' из рабочей директории.
Форматирование задаётся Markdown: http://daringfireball.net/projects/markdown/basics

сведения о Морас можно получить здесь: http://openmopac.net/manual/index.html

Работа проводится с помощью putty и обработки pymol. Команды для bash выделены табуляцией.

Abinitio вычисления для нафталина и азулена

Задание 1. Оптимизация 3D структур нафталина и азулена.

сперва установим переменные:

export PATH=${PATH}:/home/preps/golovin/progs/bin
export MOPAC_LICENSE=/home/preps/golovin/progs/bin

#Для работы с Pymol: прокси и загрузка картинок:
from xmlrpclib import ServerProxy
from IPython.display import Image
import time

cmd = ServerProxy(uri="http://localhost:9123/RPC2")

smiles азулена и нафталина запишем каждый в свой файл:

echo "C1=CC=C2C=CC=C2C=C1   Azulene" > Azulene.smi

echo "c1ccc2ccccc2c1     Napthalene" > Napthalene.smi

далее пользуясь пакетом openbabel можно сгенерировать файл со структурой молекулы

obgen Azulene.smi > azu.mol

Откроем файл. Получилась более-менее плоская структура, но все же, не идеал, потому сохраним файл как pdb-файл, который далее надо переформатировать, используя модель например PM6.

babel -ipdb azu.pdb -omop azu.mop -xk "PM6"

где az.pdb - входной файл. запустим MOPAC для расчета структуры в пространстве:

MOPAC2009.exe azu.mop

переформатируем файл обратно для визуализации, чтобы убедиться в плоскости геометрии

babel -imopout azu.out -opdb azu.pdb

и повторим действия для нафталина, получим файл nap.out после работы MOPAC. Интересно, что нафталин уже после первой команды получился плоским, но я все же проделал аналогичные операции и для него.

cmd.do('''
load azu.mol
show spheres, all
set sphere_scale, 0.2
set valence, on
''')

cmd.ray(600,400)
cmd.png('tmp\mypng1.png')
time.sleep(2)

Image(filename='tmp\mypng1.png')

Рис. 1. Изображения азулена, полученное с помощью MOPAC с использованием модели PM6.

Задание 2. Оптимизация геометрии с помощью GAMESS

GAMESS — квантово-химический пакет, которым мы найдём оптимальную геометрию для нафталена и азулена и рассчитаем теплоты образования этих молекул используя квантово-механические подходы.
Чтобы собственно использовать GAMESS, надо переформатировать данные в gamin формат и добавить заголовок попробуем все в putty:

babel -imopout azu.out -ogamin azu.gamin
babel -imopout nap.out -ogamin nap.gamin

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

далее проведем непосредственно расчет, что занимает в случае азулена немало времени, и сразу переконвертируем в pdb:

gms nap_opt.inp  1 >& nap_opt.log 
gms azu_opt.inp  1 >& azu_opt.log

babel -igamout azu_opt.log -opdb azu_gamess.pdb
babel -igamout nap_opt.log -opdb nap_gamess.pdb

cmd.do('''
load azu_gamess.pdb
show spheres, all
set sphere_scale, 0.2
set valence, on
''')

cmd.ray(600,400)
cmd.png('tmp\mypng2.png')
time.sleep(2)

Image(filename='tmp\mypng2.png')

Рис. 2. Изображения азулена, полученное с помощью GAMESS. Вид с плоскости колец.

cmd.do('rotate y, -90')

cmd.ray(600,400)
cmd.png('tmp\mypng3.png')
time.sleep(2)

Image(filename='tmp\mypng3.png')

Рис. 3. Изображения азулена, полученное с помощью GAMESS. Вид сбоку.
Можно видеть, что структура совершенно плоская, впрочем, как и полученная MOPAC.

Ниже сравним эти структуры: видно, что углеродный скелет почти не изменился, а длина связей С-Н чуть-чуть варьирует и меньше в результате GAMESS.

cmd.do('''
load azu_gamess.pdb
load azu.pdb
fit azu, azu_gamess
show spheres, all
set sphere_scale, 0.2
set valence, on
zoom all
''')

cmd.ray(600,400)
cmd.png('tmp\mypng5.png')
time.sleep(2)

Image(filename='tmp\mypng5.png')

Рис. 4. Изображения азулена, полученные с помощью GAMESS (С показан голубым) и MOPAC (С показан зеленым).
Аналогично полученная GAMESS структура для нафталина не отличается от результата MOPAC.

cmd.do('''
load nap_gamess.pdb
load nap.pdb
fit nap, nap_gamess
show spheres, all
set sphere_scale, 0.2
set valence, on
zoom all
''')

cmd.ray(600,400)
cmd.png('tmp\mypng4.png')
time.sleep(2)

Image(filename='tmp\mypng4.png')

Рис. 5. Изображения нафталина, полученные с помощью GAMESS (С показан голубым) и MOPAC (С показан зеленым).
Базис для оптимизации геометрии мы задали в тот момент, когда прописывали заголовок. Посмотрев на него еще раз и уточнив мысли по документации, убедимся, что использовался базис 6-31G, где N — число гауссианов — было выбрано равным 6. Такой базис означает, что орбитали остова представлены шестью гауссианами, а валентные орбитали представлены двумя составляющими: из трёх гауссианов и из одного гауссиана.

Подготовка к расчёту по Чартри-Фоку и DFT (теории функционала плотности)

На основе полученных GAMESS координат составим новые входные файлы для расчёта энергии: 1) для расчёта методом Хартри-Фока и 2) для использования теории функционала плотности.

1) Расчёт по Хартри-Фоку

Переформатируем файлы для расчёта по Hartree-Fock: сперва из gamout сделаем gamin, затем добавим заголовок.

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

2) То же, но с другим заголовком надо прописать для расчетов с использованием терии функционала плотности:

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

Непосредственно расчет и сравнение.

Расчет при этом аналогичен предыдущим вычислениям с GAMESS - все параметры уже заложены в заголовке.

1) Расчёт по Хартри-Фоку

gms nap_HF.inp  1 >& nap_HF.log 
gms azu_HF.inp  1 >& azu_HF.log

2) Расчёт по DFT

gms nap_DFT.inp  1 >& nap_DFT.log 
gms azu_DFT.inp  1 >& azu_DFT.log

Сравним результаты

Вытащим из log файлов расчёта энергии строчку с "TOTAL ENERGY = " и выпишите значения этой энергии в таблицу ниже.


 grep "TOTAL ENERGY =" azu_HF.log nap_HF.log azu_DFT.log nap_DFT.log
 
 
azu_HF.log:                       TOTAL ENERGY =    -383.2818684841
nap_HF.log:                       TOTAL ENERGY =    -383.3523474531
azu_DFT.log:                       TOTAL ENERGY =    -385.5852140593
nap_DFT.log:                       TOTAL ENERGY =    -385.6417295785

Известно, что энергия изомеризации нафталина в азулен составляет 35.3±2.2 kCal/mol.

1 атомная единица энергии (хартри)= 4.3597482*10-18 Дж = 627.5095 ккал/моль = 27.2116 эВ

Вещество	E _Naptalene	E _Azulene	ΔE , Hartree	ΔE, kCal/mol
Хартри-Фок	-383.35	-383.28	0.0705	44.23
DFT	-385.64	-385.59	0.0565	35.46

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

---------------------------------------------------------------------------------------

Ссылки с файлами:

SMILES --> Azulene.smi

GAMESS --> azu_opt.log

GAMESS --> azu_gamess.pdb

К сожалению, .log файлы тяжелые, поэтому добавить их все на свою страницу у меня не получилось.

SMILES --> Napthalene.smi

GAMESS --> nap_gamess.pdb