Ab initio вычисления для нафталина и азулена¶
Запустим PyMol:
from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys, time
# pymol launching
import __main__
__main__.pymol_argv = [ 'pymol', '-cp' ]
# __main__.pymol_argv = [ 'pymol', '-x' ] # for GUI
import pymol
from pymol import cmd
pymol.finish_launching()
defaultImage = 'pymolimg.png'
def prepareImage(width=300, height=300, sleep=5, filename=defaultImage):
## To save the rendered image
cmd.ray(width, height)
cmd.png('pymolimg.png')
time.sleep(sleep)
# Define some shortcuts
def focus(x):
cmd.center(x)
cmd.zoom(x)
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/bin'
%%bash
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/home/preps/golovin/progs/lib
%%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
%%bash
echo C1=CC=C2C=CC=C2C=C1 Azulene >> azulene.smi
echo "c1ccc2ccccc2c1 naphthalene" > naphthalene.smi
%%bash
obgen azulene.smi -ff UFF > azulene.mol
obgen naphthalene.smi > naphthalene.mol
Азулен¶
cmd.load('azulene.mol')
cmd.save('azulene.pdb')
prepareImage()
Image(defaultImage)
В поле UFF азулен плоский. Далее проведем оптимизацию с помощью MOPAC.
%%bash
babel -ipdb azulene.pdb -omop azulene.mop -xk "PM6"
/home/preps/golovin/progs/bin/MOPAC2009.exe azulene.mop
%%bash
babel -imopout azulene.out -omol azulene_new.mol
Нафтален¶
cmd.delete('all')
cmd.load('naphthalene.mol')
cmd.save('naphthalene.pdb')
prepareImage()
Image(defaultImage)
В поле UFF нафтален также плоский. Далее проведем оптимизацию с помощью MOPAC.
%%bash
babel -ipdb naphthalene.pdb -omop naphthalene.mop -xk "PM6"
/home/preps/golovin/progs/bin/MOPAC2009.exe naphthalene.mop
%%bash
babel -imopout naphthalene.out -omol naphthalene_new.mol
Оптимизация геометрии средствами GAMESS¶
Переформатируем полученные mol-файлы для обеих молекул для подачи на оптимизацию средствами GAMESS.
%%bash
babel -imol azulene_new.mol -ogamin azulene_opt.inp
babel -imol naphthalene_new.mol -ogamin naphthalene_opt.inp
with open("azulene_opt.inp") as f,open("azu_opt.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)
with open("naphthalene_opt.inp") as f,open("nap_opt.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 следующим образом (1 - это количество ядер для расчёта; на текущий момент на kodomo параллельное использование GAMESS не реализовано):
%%bash
gms nap_opt.inp 1 >& nap_opt.log
gms azu_opt.inp 1 >& azu_opt.log
На основе полученных координат составим новые входные файлы для расчёта энергии.
Цель: построить по два файла на каждую молекулу, в первом случае проводя расчёт методом Хартри-Фока, а во втором используя теорию функционала плотности.
%%bash
babel -igamout azu_opt.log -ogamin azu_hart.inp
babel -igamout nap_opt.log -ogamin nap_hart.inp
Для расчёта по Хартри-Фоку изменим заголовок следующим образом (для обеих молекул):
with open("azu_hart.inp") as f,open("azu_hart_hf.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_hf.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("azu_hart.inp") as f,open("azu_hart_tfp.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("nap_hart.inp") as f,open("nap_hart_tfp.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)
Рассчет по Хартри-Фоку:
%%bash
gms azu_hart_hf.inp 1 >& azu_hart_hf_opt.log
gms nap_hart_hf.inp 1 >& nap_hart_hf_opt.log
Рассчет с применением теории функционала плотности:
%%bash
gms azu_hart_tfp.inp 1 >& azu_hart_tfp_opt.log
gms nap_hart_tfp.inp 1 >& nap_hart_tfp_opt.log
Найдем в log файлах расчёта энергии строчку с "TOTAL ENERGY = ".
Цель: понять, какой метод лучше если из эксперимента известно, что энергия изомеризации нафталина в азулен составляет 35.3±2.2 kCal/mol
Для рассчета по Хартри-Фоку:
cat azu_hart_hf_opt.log | grep "TOTAL ENERGY =" > 1.log
%%bash
grep "TOTAL ENERGY =" azu_hart_hf_opt.log
grep "TOTAL ENERGY =" nap_hart_hf_opt.log
Для рассчета с применением т. функционала плотности:
%%bash
grep "TOTAL ENERGY =" azu_hart_tfp_opt.log
grep "TOTAL ENERGY =" nap_hart_tfp_opt.log
Для рассчета по Хартри-Фоку ΔE, kCal/mol (=ΔE , Hartree (= total1 - total2) -> kCal/mol):
0.0725179177 ->
45.50569011112
Для рассчета с применением т. функционала плотности:
0.0546314356 ->
34.2817507395
Таким образом, более точного результата позволила добиться теория функционала плотности.