In [1]:
from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys
import scipy as sp
from scipy import constants
from scipy.constants import codata
import numpy as np
In [2]:
import __main__
__main__.pymol_argv = [ 'pymol', '-cp' ]
import pymol
pymol.finish_launching()
from pymol import cmd
In [3]:
from IPython.display import Image
In [42]:
%%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
ЗАДАНИЕ 1¶
Оптимизируем структуру азулена с помощью MOPAC¶
In [5]:
%%bash
cd
cd Term6
cd Practice3
echo 'C1=CC=C2C=CC=C2C=C1 Azulene' > azulene.smi
obgen azulene.smi > azulene.mol
babel -imol azulene.mol -omop azulene_opt.mop -xk "PM6"
MOPAC2009.exe azulene_opt.mop
In [10]:
%%bash
cd
cd Term6
cd Practice3
babel -imopout azulene_opt.out -opdb azulene_opt.pdb
In [55]:
cmd.do('''
reinitialize
set ray_trace_mode, 0; red
set ray_opaque_background, off
set antialias, .5
set light_count, 8
set ambient, 0.5
set ray_trace_color, red
set cartoon_side_chain_helper, on
cd /home/students/y12/iltarn/Term6/Practice3
load azulene_opt.pdb
orient azulene_opt
as sticks, azulene_opt
zoom azulene_opt, 4
label all, " %s" % ID
set label_color, black
rotate x, -90
ray
png pic1.png
''')
In [56]:
Image(filename='pic1.png')
Out[56]:
Видим, что молекула не плоская. Займемся магией - попробуем разные силовые поля (MMFF94, MMFF94s, UFF) в obgen.
MMFF94¶
In [53]:
%%bash
obgen azulene.smi -ff MMFF94 > azulene_MMFF94.mol
babel -imol azulene_MMFF94.mol -omop azulene_MMFF94_opt.mop -xk "PM6"
MOPAC2009.exe azulene_MMFF94_opt.mop
In [54]:
%%bash
cd
cd Term6
cd Practice3
babel -imopout azulene_MMFF94_opt.out -opdb azulene_MMFF94_opt.pdb
In [64]:
cmd.do('''
delete all
load azulene_MMFF94_opt.pdb
orient azulene_MMFF94_opt
as sticks, azulene_MMFF94_opt
zoom azulene_MMFF94_opt, 4
label all, " %s" % ID
set label_color, black
rotate x, -90
ray
png pic3.png
''')
In [65]:
Image(filename='pic3.png')
Out[65]:
Ничего не изменилось
UFF¶
In []:
%%bash
obgen azulene.smi -ff UFF > azulene_UFF.mol
babel -imol azulene_UFF.mol -omop azulene_UFF_opt.mop -xk "PM6"
MOPAC2009.exe azulene_UFF_opt.mop
In [67]:
%%bash
cd
cd Term6
cd Practice3
babel -imopout azulene_UFF_opt.out -opdb azulene_UFF_opt.pdb
In [68]:
cmd.do('''
delete all
load azulene_UFF_opt.pdb
orient azulene_UFF_opt
as sticks, azulene_UFF_opt
zoom azulene_UFF_opt, 4
label all, " %s" % ID
set label_color, black
rotate x, -90
ray
png pic4.png
''')
In [69]:
Image(filename='pic4.png')
Out[69]:
Молекула стала плоской. Преобразуем ее в azu.mol
In [111]:
%%bash
babel -imopout azulene_UFF_opt.out -omol azu.mol
Оптимизируем структуру нафталена с помощью MOPAC.¶
In [43]:
%%bash
cd
cd Term6
cd Practice3
echo 'c1ccc2ccccc2c1 Naphthalene' > naphthalene.smi
obgen naphthalene.smi > naphthalene.mol
babel -imol naphthalene.mol -omop naphthalene_opt.mop -xk "PM6"
MOPAC2009.exe naphthalene_opt.mop
In [44]:
%%bash
cd
cd Term6
cd Practice3
babel -imopout naphthalene_opt.out -opdb naphthalene_opt.pdb
In [48]:
cmd.do('''
delete all
cd /home/students/y12/iltarn/Term6/Practice3
load naphthalene_opt.pdb
orient naphthalene_opt
as sticks, naphthalene_opt
zoom naphthalene_opt, 4
label all, " %s" % ID
set label_color, black
rotate x, -30
rotate y, 30
ray
png pic2.png
''')
In [49]:
Image(filename='pic2.png')
Out[49]:
In [50]:
cmd.do('''
rotate y, -30
rotate x, -60
ray
png pic2_90.png
''')
In [52]:
Image(filename='pic2_90.png')
Out[52]:
Молекула плоская. Делаем из нее nap.mol
In [112]:
%%bash
babel -imopout naphthalene_opt.out -omol nap.mol
ЗАДАНИЕ 2. Оптимизация геометрии средствами GAMESS¶
С помощью babel переформатируем координаты в gamin (или inp) формат:¶
In [115]:
%%bash
babel -imol azu.mol -ogamin azu_opt.inp
babel -imol nap.mol -ogamin nap_opt.inp
Добавим заголовки:¶
In [116]:
with open("azu_opt.inp") as f,open("azu_opt_h.inp",'w') as o:
data=f.read()
data=data.replace("$CONTRL COORD=CART UNITS=ANGS $END","$CONTRL COORD=CART UNITS=ANGS SCFTYP=RHF RUNTYP=OPTIMIZE $END\n $BASIS GBASIS=N31 NGAUSS=6 $end\n $system mwords=2 $end")
o.write(data)
In [117]:
with open("nap_opt.inp") as f,open("nap_opt_h.inp",'w') as o:
data=f.read()
data=data.replace("$CONTRL COORD=CART UNITS=ANGS $END","$CONTRL COORD=CART UNITS=ANGS SCFTYP=RHF RUNTYP=OPTIMIZE $END\n $BASIS GBASIS=N31 NGAUSS=6 $end\n $system mwords=2 $end")
o.write(data)
Запускаем GAMESS¶
In [119]:
%%bash
gms nap_opt_h.inp 1 >& nap_opt.log
gms azu_opt_h.inp 1 >& azu_opt.log
ЗАДАНИЕ 3. Рассчет энергий.¶
а) Расчёт по Хартри-Фоку¶
In [122]:
%%bash
# Переформатирование... Добавляем текст в файл .log
babel -igamout azu_opt.log -ogamin azu_hart.inp
babel -igamout nap_opt.log -ogamin nap_hart.inp
In [126]:
with open("azu_hart.inp") as f,open("azu_hart_n.inp",'w') as o:
data=f.read()
data=data.replace("$CONTRL COORD=CART UNITS=ANGS $END","$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")
o.write(data)
with open("nap_hart.inp") as f,open("nap_hart_n.inp",'w') as o:
data=f.read()
data=data.replace("$CONTRL COORD=CART UNITS=ANGS $END","$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")
o.write(data)
б) Расчет по теории функционала плотности¶
In [124]:
%%bash
# Переформатирование... Добавляем текст в файл .log
babel -igamout azu_opt.log -ogamin azu_th.inp
babel -igamout nap_opt.log -ogamin nap_th.inp
In [127]:
with open("nap_th.inp") as f,open("nap_th_n.inp",'w') as o:
data=f.read()
data=data.replace("$CONTRL COORD=CART UNITS=ANGS $END","$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")
o.write(data)
with open("azu_th.inp") as f,open("azu_th_n.inp",'w') as o:
data=f.read()
data=data.replace("$CONTRL COORD=CART UNITS=ANGS $END","$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")
o.write(data)
ЗАДАНИЕ 4. Энергии¶
In [128]:
%%bash
gms azu_hart_n.inp 1 >& azu_hart_n.log
gms azu_th_n.inp 1 >& azu_th_n.log
gms nap_hart_n.inp 1 >& nap_hart_n.log
gms nap_th_n.inp 1 >& nap_th_n.log
In [131]:
%%bash
awk '/TOTAL ENERGY =/{print $4}' *_n.log
# Вывод значений в алфавитном порядке по именам файлов.
Вещество |
E Naphthalene |
E Azulene |
ΔE , Hartree |
ΔE, kCal/mol |
Хартри-Фок |
-383.3549297881 |
-383.2824173049 |
0.072512483 |
45.50226983 |
DFT |
-385.6401629708 |
-385.5855568639 |
0.054606107 |
34.26584914 |
Из эксперимента известно, что энергия изомеризации нафталина в азулен составляет 35.3±2.2 kCal/mol. Значение, полученное по теории функционала плотности, ближе к экспериментальному значению. Следовательно, метод DFT лучше. Метод Хартри-Фока дает значительную ошибку.
Ссылки на файлы с результатами рассчетов:
azu_opt.log
azu_hart_n.log
azu_th_n.log
nap_opt.log
nap_hart_n.log
nap_th_n.log