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

In [1]: 
from xmlrpclib import ServerProxy
from IPython.display import Image
import time
In [2]: 
cmd = ServerProxy(uri="http://localhost:9123/RPC2")
In [3]: 
cmd.cd('/Users/aleksandrabezmenova/Documents/fbb/term8/Golovin/Pr4')
cmd.bg_color('white')
In [6]: 
ImPath = '/tmp/pymolpng.png'
def MakeImage():
    cmd.ray(600,400)
    cmd.png(ImPath)
    time.sleep(2)
    
def focus(sele):
    cmd.center(sele)
    cmd.zoom(sele,'1')

1. Построение и оптимизация структуры нафталена и азулена с помощью MOPAC

$ echo 'C1=CC=C2C=CC=C2C=C1' > Azulene.smi
$ echo 'c1ccc2ccccc2c1' > Napthalene.smi

Подберем такие силовые поля в obgen, чтобы структуры получились плоские.

$ obgen Napthalene.smi -ff MMFF94 > Napthalene.mol
$ obgen Azulene.smi -ff UFF > Azulene.mol
In [7]: 
MakeImage()
Image(ImPath)
Out[7]:
$ babel -imol Azulene.mol -omop Azulene.mop -xk "PM6"
$ babel -imol Napthalene.mol -omop Napthalene.mop -xk "PM6"

$ MOPAC2009.exe Azulene.mop
$ MOPAC2009.exe Napthalene.mop

2. Оптимизация геометрии обеих молекул с помощью GAMESS

In []: 
%%bash
# Первформатируем вывод MOPAC во вход для GAMESS
babel -imopout Azulene.out -ogamin azu_opt.inp
babel -imopout Napthalene.out -ogamin nap_opt.inp

# Добавим заголовки
sed -i 's/$CONTRL COORD=CART UNITS=ANGS $END/\
\$CONTRL COORD=CART UNITS=ANGS   SCFTYP=RHF RUNTYP=OPTIMIZE $END\
\$BASIS  GBASIS=N31 NGAUSS=6 $end\
\$system mwords=2 $end/' nap_opt.inp

sed -i 's/$CONTRL COORD=CART UNITS=ANGS $END/\
\$CONTRL COORD=CART UNITS=ANGS   SCFTYP=RHF RUNTYP=OPTIMIZE $END\
\$BASIS  GBASIS=N31 NGAUSS=6 $end\
\$system mwords=2 $end/' azu_opt.inp

# Запустим оптимизацию
gms nap_opt.inp  1 >& nap_opt.log 
gms azu_opt.inp  1 >& azu_opt.log

3. Подготовка входных файлов для рассчета энергии

In []: 
%%bash
babel -igamout nap_opt.log -ogamin nap_opt.inp
babel -igamout azu_opt.log -ogamin azu_opt.inp

# Input for Hartree-Fock
sed 's/$CONTRL COORD=CART UNITS=ANGS $END/\
\$CONTRL COORD=CART UNITS=ANGS   SCFTYP=RHF RUNTYP=ENERGY $END\
\$BASIS  GBASIS=N31 NGAUSS=6\
\ POLAR=POPN31 NDFUNC=1 $END\
\$GUESS  GUESS=HUCKEL $END\
\$system mwords=2 $end/' nap_opt.inp > nap_HF.inp

sed 's/$CONTRL COORD=CART UNITS=ANGS $END/\
\$CONTRL COORD=CART UNITS=ANGS   SCFTYP=RHF RUNTYP=ENERGY $END\
\$BASIS  GBASIS=N31 NGAUSS=6\
\ POLAR=POPN31 NDFUNC=1 $END\
\$GUESS  GUESS=HUCKEL $END\
\$system mwords=2 $end/' azu_opt.inp > azu_HF.inp

# Input for functional density
sed 's/$CONTRL COORD=CART UNITS=ANGS $END/\
\$CONTRL COORD=CART UNITS=ANGS   dfttyp=b3lyp RUNTYP=ENERGY $END\
\$BASIS  GBASIS=N31 NGAUSS=6\
\ POLAR=POPN31 NDFUNC=1 $END\
\$GUESS  GUESS=HUCKEL $END\
\$system mwords=2 $end/' nap_opt.inp > nap_DF.inp

sed 's/$CONTRL COORD=CART UNITS=ANGS $END/\
\$CONTRL COORD=CART UNITS=ANGS   dfttyp=b3lyp RUNTYP=ENERGY $END\
\$BASIS  GBASIS=N31 NGAUSS=6\
\ POLAR=POPN31 NDFUNC=1 $END\
\$GUESS  GUESS=HUCKEL $END\
\$system mwords=2 $end/' azu_opt.inp > azu_DF.inp

4. Рассчет энергии структур нафталина и азулена двумя способами

Рассчет методом Хартри-Фока:

In []: 
%%bash
gms nap_HF.inp  1 >& nap_HF.log 
gms azu_HF.inp  1 >& azu_HF.log
In [3]: 
%%bash
grep "TOTAL ENERGY =" nap_HF.log azu_HF.log

nap_HF.log:                       TOTAL ENERGY =    -383.3549297881
azu_HF.log:                       TOTAL ENERGY =    -383.2823803802

Рассчет c использованием теории функционала плотности:

In []: 
%%bash
gms nap_DF.inp  1 >& nap_DF.log 
gms azu_DG.inp  1 >& azu_DF.log
In [4]: 
%%bash
grep "TOTAL ENERGY =" nap_DF.log azu_DF.log

nap_DF.log:                       TOTAL ENERGY =    -385.6401629708
azu_DF.log:                       TOTAL ENERGY =    -385.5855055121

Сравнение энергий, полученных разными методами:

Вещество ENaphthalene EAzulene ΔE, Hartre ΔE, kCal/mol
Хартри-Фок -383.35 -383.28 0.07 45.62
DFT -385.64 -385.59 0.05 34.32

К экспериментальным результатам (ΔE изомеризации нафталина в азулен равна 35.3±2.2 kCal/mol) ближе результат, полученный с использованием теории функционала плотности.