Домашняя работа 3. Abinitio вычисления для нафталина и азулена

Рабочая директория с файлами: H:83

Цель работы: научиться обращаться с пакетом квантово-механических рассчетов GAMESS.

Создаю .smi файлы и сразу пробую получить 3D-структуры с использованием разных иловых полей:

In [13]: 
%%bash
echo "C1=CC=C2C=CC=C2C=C1" > azulene.smi
echo "c1ccc2ccccc2c1" > naphthalene.smi

obgen azulene.smi > azulene.mol
obgen azulene.smi -ff MMFF94 > azulene_MMFF94.mol
obgen azulene.smi -ff MMFF94s > azulene_MMFF94s.mol
obgen azulene.smi -ff UFF > azulene_UFF.mol
babel --gen3D -ismi azulene.smi -omol azulene_gen3D.mol

obgen naphthalene.smi > naphthalene.mol
obgen naphthalene.smi -ff MMFF94 > naphthalene_MMFF94.mol
obgen naphthalene.smi -ff MMFF94s > naphthalene_MMFF94s.mol
obgen naphthalene.smi -ff UFF > naphthalene_UFF.mol
babel --gen3D -ismi naphthalene.smi -omol naphthalene_gen3D.mol

Смотрю, какой вариант мне больше всего нравится:

In [2]: 
from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys
import time
# pymol launching
import __main__

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

### Если вывод в графическое окно не нужен, то:
##__main__.pymol_argv = [ 'pymol', '-cp' ]


import pymol
 
pymol.finish_launching()
 
from pymol import cmd
In [58]: 
cmd.reinitialize()
cmd.load("/home/students/y11/agalicina/Term8/Practice3/azulene.mol")
cmd.do(
'''
ray 500,500
png images/a1.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/a1.png')
Out[58]:
In [59]: 
cmd.turn('x',60)
cmd.do(
'''
ray 500,500
png images/a11.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/a11.png')
Out[59]:
In [61]: 
cmd.reinitialize()
cmd.load("/home/students/y11/agalicina/Term8/Practice3/azulene_MMFF94.mol")
cmd.do(
'''
ray 500,500
png images/a2.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/a2.png')
Out[61]:
In [62]: 
cmd.turn('x',60)
cmd.do(
'''
ray 500,500
png images/a21.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/a21.png')
Out[62]:
In [63]: 
cmd.reinitialize()
cmd.load("/home/students/y11/agalicina/Term8/Practice3/azulene_MMFF94s.mol")
cmd.do(
'''
ray 500,500
png images/a3.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/a3.png')
Out[63]:
In [64]: 
cmd.turn('x',60)
cmd.do(
'''
ray 500,500
png images/a31.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/a31.png')
Out[64]:
In [67]: 
cmd.reinitialize()
cmd.load("/home/students/y11/agalicina/Term8/Practice3/azulene_UFF.mol")
cmd.do(
'''
ray 500,500
png images/a4.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/a4.png')
Out[67]:
In [68]: 
cmd.turn('x',65)
cmd.do(
'''
ray 500,500
png images/a41.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/a41.png')
Out[68]:
In [71]: 
cmd.reinitialize()
cmd.load("/home/students/y11/agalicina/Term8/Practice3/azulene_gen3D.mol")
cmd.do(
'''
ray 500,500
png images/a5.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/a5.png')
Out[71]:
In [72]: 
cmd.turn('x',40)
cmd.do(
'''
ray 500,500
png images/a51.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/a51.png')
Out[72]:

Вывод: для азулена и данного SMILES лучше применять модель UFF obgen. Беру файл azulene_UFF.mol для дальнейшей работы.

Повторяю процедуру для нафталина.

In [40]: 
cmd.reinitialize()
cmd.load("/home/students/y11/agalicina/Term8/Practice3/naphthalene.mol")
cmd.do(
'''
ray 500,500
png images/n1.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/n1.png')
Out[40]:
In [41]: 
cmd.turn('x',90)
cmd.do(
'''
ray 500,500
png images/n11.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/n11.png')
Out[41]:
In [42]: 
cmd.reinitialize()
cmd.load("/home/students/y11/agalicina/Term8/Practice3/naphthalene_MMFF94.mol")
cmd.do(
'''
ray 500,500
png images/n2.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/n2.png')
Out[42]:
In [43]: 
cmd.turn('x',90)
cmd.do(
'''
ray 500,500
png images/n21.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/n21.png')
Out[43]:
In [44]: 
cmd.reinitialize()
cmd.load("/home/students/y11/agalicina/Term8/Practice3/naphthalene_MMFF94s.mol")
cmd.do(
'''
ray 500,500
png images/n3.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/n3.png')
Out[44]:
In [45]: 
cmd.turn('x',90)
cmd.do(
'''
ray 500,500
png images/n31.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/n31.png')
Out[45]:
In [46]: 
cmd.reinitialize()
cmd.load("/home/students/y11/agalicina/Term8/Practice3/naphthalene_UFF.mol")
cmd.do(
'''
ray 500,500
png images/n4.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/n4.png')
Out[46]:
In [47]: 
cmd.turn('x',90)
cmd.do(
'''
ray 500,500
png images/n41.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/n41.png')
Out[47]:
In [48]: 
cmd.reinitialize()
cmd.load("/home/students/y11/agalicina/Term8/Practice3/naphthalene_gen3D.mol")
cmd.do(
'''
ray 500,500
png images/n5.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/n5.png')
Out[48]:
In [49]: 
cmd.turn('x',90)
cmd.do(
'''
ray 500,500
png images/n51.png
''')
time.sleep(2)
Image(filename='/home/students/y11/agalicina/Term8/Practice3/images/n51.png')
Out[49]:
In [50]: 
cmd.reinitialize()

Вывод. Нафталин "испортить" разными силовыми полями не получилось. Можно брать любую модель, но для удобства написания команд возьму ту же модель, что и для азулена: naphthalene_UFF.mol

Перезаписываю выбранные файлы во входные для GAMESS, добавляю шапку и меняю формат на .inp:

In [82]: 
%%bash
babel --imol azulene_UFF.mol -oinp azu_opt.inp
babel --imol naphthalene_UFF.mol -oinp nap_opt.inp
In [84]: 
%%bash
printf ' $CONTRL COORD=CART UNITS=ANGS   SCFTYP=RHF RUNTYP=OPTIMIZE $END\n $BASIS  GBASIS=N31 NGAUSS=6  $END\n $SYSTEM MWORDS=2 $END\n' > az_opt.inp
tail -n +3 azu_opt.inp >> az_opt.inp

printf ' $CONTRL COORD=CART UNITS=ANGS   SCFTYP=RHF RUNTYP=OPTIMIZE $END\n $BASIS  GBASIS=N31 NGAUSS=6  $END\n $SYSTEM MWORDS=2 $END\n' > na_opt.inp
tail -n +3 nap_opt.inp >> na_opt.inp
In [1]: 
%%bash
gms na_opt.inp  1 >& na_opt.log 
gms az_opt.inp  1 >& az_opt.log
In [1]: 
%%bash
babel --igamout na_opt.log --ogamin na_opt.gamin
babel --igamout az_opt.log --ogamin az_opt.gamin

1 molecule converted
17 audit log messages 
1 molecule converted
2 info messages 17 audit log messages

Делаю по два файла на каждую молекулу с замененными заголовками:

In [5]: 
%%bash
printf ' $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\n' > na_HF.inp
printf ' $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\n' > az_HF.inp
printf ' $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\n' > na_FT.inp
printf ' $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\n' > az_FT.inp

tail -n +4 na_opt.gamin >> na_HF.inp
tail -n +4 na_opt.gamin >> na_FT.inp
tail -n +4 az_opt.gamin >> az_HF.inp
tail -n +4 az_opt.gamin >> az_FT.inp
In []: 
gms na_HF.inp  1 >& na_HF.log 
gms az_HF.inp  1 >& az_HF.log
gms na_FT.inp  1 >& na_FT.log 
gms az_FT.inp  1 >& az_FT.log

Оформление результатов в виде таблицы:

In [6]: 
from IPython.display import HTML
In [19]: 
%%bash
printf "Na_HF\t" > resulting_file.txt
grep 'TOTAL ENERGY =' na_HF.log >> resulting_file.txt

printf "Na_FT\t" >> resulting_file.txt
grep 'TOTAL ENERGY =' na_FT.log >> resulting_file.txt

printf "Az_HF\t" >> resulting_file.txt
grep 'TOTAL ENERGY =' az_HF.log >> resulting_file.txt

printf "Az_HF\t" >> resulting_file.txt
grep 'TOTAL ENERGY =' az_FT.log >> resulting_file.txt
In [20]: 
with open('resulting_file.txt','r') as f:
    l = [ float( x.split()[4] ) for x in f.readlines() ]
In [36]: 
s = """<table>
<tr>
<th>Вещество</th>
<th>E Naptalene</th>
<th>E Azulene</th>
<th>delta E , Hartree</th>
<th>delta E, kCal/mol</th>
</tr>
<tr>
<td>Харти-Фок</td>
<td>%f</td>
<td>%f</td>
<td>%f</td>
<td>%f</td>
</tr>
<tr>
<td>DFT (теория функционала плотности) </td>
<td>%f</td>
<td>%f</td>
<td>%f</td>
<td>%f</td>
</tr>
</table>""" % (l[0], l[2], l[0]-l[2], 627.509*(l[0]-l[2]), l[1], l[3], l[1]-l[3], 627.509*(l[1]-l[3]) ) 
# 1 Hartee = 627.509 kcal mol-1
In [37]: 
h = HTML(s); h
Out[37]:
Вещество E Naptalene E Azulene delta E , Hartree delta E, kCal/mol
Харти-Фок -383.354921 -383.282266 -0.072655 -45.591907
DFT (теория функционала плотности) -385.640246 -385.585751 -0.054496 -34.196455

Сравним с экспериментальыми данными \(\Delta E = 35.3±2.2 \dfrac{kCal}{mol}\): лучшее согласование результатов применения теории функционала плотности.

In []: