In [2]: 
import os
os.environ['PATH']=os.environ['PATH']+'/home/preps/golovin/progs/bin'
os.environ['MOPAC_LICENSE']='/home/preps/golovin/progs/bin'
os.environ['LD_LIBRARY_PATH']='/home/preps/golovin/progs/lib'
In [4]: 
%%bash
export PATH=${PATH}:/home/preps/golovin/progs/bin
export MOPAC_LICENSE=/home/preps/golovin/progs/bin
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/home/preps/golovin/progs/lib
In [23]: 
%%bash
echo C1=CC=C2C=CC=C2C=C1 Azulene >> azulene.smi
In [24]: 
%%bash
echo c1ccc2ccccc2c1 Naphthalene >> naphtalene.smi
In []: 
##Создали файлы smiles для нафталина и азулена
In [25]: 
from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys
In [26]: 
import __main__

__main__.pymol_argv = [ 'pymol', '-cp' ]



import pymol
 
pymol.finish_launching()
 
from pymol import cmd
In [27]: 
from IPython.display import Image
In [28]: 
%%bash
obgen azulene.smi -ff UFF > azu.mol
In [29]: 
%%bash
obgen naphtalene.smi -ff UFF > naph.mol
In [30]: 
cmd.delete('all')
cmd.load('azu.mol')

cmd.do('''
ray
png azu1.png
''')
In [31]: 
Image(filename='azu1.png')
Out[31]:
In [32]: 
cmd.save("azul_start.pdb")
In [33]: 
cmd.delete('all')
cmd.load('naph.mol')
cmd.do('''
rotate x, 90
ray
png nap1.png
''')
In [34]: 
Image(filename='nap1.png')
Out[34]:
In [35]: 
cmd.save("naph_start.pdb")
In [19]: 
%%bash
babel -ipdb azul_start.pdb -omop azul_opt.mop -xk "PM6"
/home/preps/golovin/progs/bin/MOPAC2016.exe  azul_opt.mop
In [1]: 
%%bash
babel -imopout azul_opt.out -omol azul_opt.mol

1 molecule converted
2 info messages 22 audit log messages
In [82]: 
cmd.delete('all')
cmd.load('azul_opt.mol')
cmd.do('''
ray
png azul2.png
''')
In [84]: 
Image(filename='azul2.png')
Out[84]:
In [7]: 
%%bash
export PATH=${PATH}:/home/preps/golovin/progs/bin
export MOPAC_LICENSE=/home/preps/golovin/progs/bin
In [4]: 
%%bash
babel -ipdb naph_start.pdb -omop naph_opt.mop -xk "PM6"
/home/preps/golovin/progs/bin/MOPAC2016.exe  naph_opt.mop
##В итоге использовался MOPAC2009, запущенный через консоль

1 molecule converted
42 info messages 14 audit log messages 


          MOPAC Job: "naph_opt.mop" ended normally on Mar 18, 2016, at 14:16.
In [2]: 
%%bash
babel -imopout naph_opt.out -omol naph_opt.mol

1 molecule converted
22 audit log messages
In []: 
##Таким образом, провели преоптимизацию для обеих молекул с помощью MOPAC.
In [3]: 
%%bash
babel -imol azul_opt.mol -ogamin azul_opt_new.inp
babel -imol naph_opt.mol -ogamin naph_opt_new.inp

1 molecule converted
18 audit log messages 
1 molecule converted
18 audit log messages
In [6]: 
%%bash
gms naph_opt_new.inp  1 >& nap_opt.log 
gms azul_opt_new.inp  1 >& azu_opt.log
In [8]: 
%%bash
babel -igamout azu_opt.log -ogamin azu_hart.inp
babel -igamout nap_opt.log -ogamin nap_hart.inp

1 molecule converted
2 info messages 17 audit log messages 
1 molecule converted
17 audit log messages
In [9]: 
%%bash
babel -igamout azu_opt.log -ogamin azu_dft.inp
babel -igamout nap_opt.log -ogamin nap_dft.inp

1 molecule converted
2 info messages 17 audit log messages 
1 molecule converted
17 audit log messages
In [11]: 
%%bash
gms nap_hart.inp  1 >& nap_hart_opt.log 
gms azu_hart.inp  1 >& azu_hart_opt.log
gms nap_dft.inp  1 >& nap_dft_opt.log 
gms azu_dft.inp  1 >& azu_dft_opt.log
In []: 
##В полученных файлах почему-то не было строки TOTAL ENERGY, поэтому пришлось складывать кинетическую энергию и потенциальную.
In [15]: 
%%bash 
grep 'TOTAL POTENTIAL ENERGY' nap_hart_opt.log >> energy_nap_hart.txt
grep 'TOTAL KINETIC ENERGY' nap_hart_opt.log >> energy_nap_hart.txt
In [49]: 
%%bash 
grep 'TOTAL POTENTIAL ENERGY' nap_dft_opt.log >> energy_nap_dft.txt
grep 'TOTAL KINETIC ENERGY' nap_dft_opt.log >> energy_nap_dft.txt
grep 'TOTAL POTENTIAL ENERGY' azu_hart_opt.log >> energy_azu_hart.txt
grep 'TOTAL KINETIC ENERGY' azu_hart_opt.log >> energy_azu_hart.txt
grep 'TOTAL POTENTIAL ENERGY' azu_dft_opt.log >> energy_azu_dft.txt
grep 'TOTAL KINETIC ENERGY' azu_dft_opt.log >> energy_azu_dft.txt
In [4]: 
file = open('energy_azu_dft.txt', 'r')
a = []
for line in file: 
    a.append(float(line.split()[-1]))
print(a)
total_azu_dft=sum(a)
print(total_azu_dft)

file = open('energy_nap_dft.txt', 'r')
b = []
for line in file: 
    b.append(float(line.split()[-1]))
print(b)
total_nap_dft=sum(b)
print(total_nap_dft)

delta_dft_h = total_azu_dft - total_nap_dft 
print(delta_dft_h)
##энергия изомеризации по dft в единицах Хартри
delta_kcal=delta_dft_h*627.509469
print(delta_kcal)
##энергия изомеризации по dft в единицах kcal/mol

[-768.0135367173, 382.428225882]
-385.585310835
[-768.112155235, 382.4724979851]
-385.63965725
0.0543464146
34.1028897677
In [3]: 
file = open('energy_nap_hart.txt', 'r')
a = []
for line in file: 
    a.append(float(line.split()[-1]))
print(a)
sum(a)
total_nap_hart=sum(a)
print(total_nap_hart)

file = open('energy_azu_hart.txt', 'r')
b = []
for line in file: 
    b.append(float(line.split()[-1]))
print(b)
sum(b)
total_azu_hart=sum(b)
print(total_azu_hart)

delta_hart_h = total_azu_hart - total_nap_hart
print(delta_hart_h)
##энергия изомеризации по Хартри-Фоку в единицах Харти
delta_kcal=delta_hart_h*627.509469
print(delta_kcal)
##энергия изомеризации по Хартри-Фоку в единицах kcal/mol

[-766.3700375516, 383.0155220424]
-383.354515509
[-766.2296430979, 382.9473501296]
-383.282292968
0.0722225409001
45.3203282901
In []: 
##Рассчитали суммарную энергию молекул по Хартри-Фоку и DFT, посчитали энергию изомеризации в каждом из случаев
##Больше похожа на экспериментальную (35.3±2.2 kCal/mol) получилась при рассчете по DFT

E Naphthalene

E Azulene

?E , Hartree

?E, kCal/mol
-

-383.282292968

-383.354515509

0.0722225409001

45.3203282901
DFT

-385.585310835

-385.63965725

0.0543464146

34.1028897677


DFT DFT