Cеми-эмпирические вычисления: MOPAC

1. Построение структуры порфирина

In [1]: 
from xmlrpclib import ServerProxy
from IPython.display import Image
import time
In [6]: 
cmd = ServerProxy(uri="http://localhost:9123/RPC2")
In [7]: 
cmd.bg_color('white')
Dir = '/Users/aleksandrabezmenova/Documents/fbb/term8/Golovin/Pr3/'

ImPath = '/tmp/pymolpng.png'
def MakeImage():
    cmd.ray(500,400)
    cmd.png(ImPath)
    time.sleep(2)
    
def focus(sele):
    cmd.center(sele)
    cmd.zoom(sele)

Загрузим SMILES аннотацию порфирина и построим по ней 3d структуру с помощью obgen:

$ export PATH=${PATH}:/home/preps/golovin/progs/bin
$ export MOPAC_LICENSE=/home/preps/golovin/progs/bin
$ echo 'c1cc2cc3ccc(cc4ccc(cc5ccc(cc1n2)[nH]5)n4)[nH]3   porphyrin' > 1.smi
$ obgen 1.smi > 1.mol
$ babel -imol 1.mol -opdb 1.pdb

Загрузим получившуюся структуру в Pymol:

In [16]: 
cmd.load(Dir+'1.pdb','porph')
cmd.set('stick_radius','0.1')
cmd.set('sphere_scale','0.1')
cmd.show('sticks','porph')
cmd.show('sphere','porph')
MakeImage()
Image(ImPath)
Out[16]:

и удалим лишние водороды:

In [18]: 
cmd.extract('hydr','id 39+40')
cmd.delete('hydr')
focus('porph')
MakeImage()
Image(ImPath)
Out[18]:

Экспортируем структуру, преобразуем формат и оптимизируем с помощью программы MOPAC. Используем два разных типа параметризации - PM6 и AM1.

$ babel -ipdb 1_1.pdb -omop 1_opt_pm6.mop -xk "PM6"
$ MOPAC2009.exe 1_opt_pm6.mop
$ babel -imopout 1_opt_pm6.out -opdb 1_opt_pm6.pdb

$ babel -ipdb 1_1.pdb -omop 1_opt_am1.mop -xk "AM1"
$ MOPAC2009.exe 1_opt_am1.mop
$ babel -imopout 1_opt_am1.out -opdb 1_opt_am1.pdb

Сравним получившиеся структуры в PyMol:

In [23]: 
cmd.load(Dir+'1_opt_pm6.pdb','porph_pm6')
cmd.load(Dir+'1_opt_am1.pdb','porph_am1')
cmd.color('grey60','porph')
cmd.color('forest','porph_pm6')
cmd.color('tv_red','porph_am1')
cmd.show('sticks','all')
cmd.show('sphere','all')
MakeImage()
Image(ImPath)
Out[23]:
In [24]: 
MakeImage()
Image(ImPath)
Out[24]:
In [26]: 
cmd.hide('everything','porph')
MakeImage()
Image(ImPath)
Out[26]:

Структура, полученная с помощью Obgen, помимо того, что содержит лишние водороды, явно недостаточно плоская для ароматической структуры. Структуры, полученные с помощью MOPAC, более плоские. При разных типах параметризации структуры получаюся немного разными: если смотреть на структуру, приведенную на рисунках, то в структуре Obgen кольца с водородами "опущены", а кольца без водородов - "подняты". В структурах MOPAC на месте колец с водородами мы видим красную структуру (AM1), а на месте колец без водородав - зеленую (PM6). Таким образом, структура AM1 уплостилась чуть-чуть сильнее.

2. Возбужденное состояние порфирина

Создадим для MOPAC файл .mop с указанием того, что надо рассчитать возбужденное состояние:

$ cp 1_opt_pm6.mop 1_opt_pm6_spectr.mop
$ echo >> 1_opt_pm6_spectr.mop   # empty line is necessary 
$ echo "cis c.i.=4 meci oldgeo" >> 1_opt_pm6_spectr.mop
$ echo "Porfirine spectrum" >> 1_opt_pm6_spectr.mop

Запустим MOPAC:

$ MOPAC2009.exe 1_opt_pm6_spectr.mop

В конце файла 1_opt_pm6_spectr.out есть информация об энергиях электронных переходов:

$ tail -19 1_opt_pm6_spectr.out 

  STATE       ENERGY (EV)        Q.N.  SPIN   SYMMETRY     POLARIZATION
         ABSOLUTE     RELATIVE                            X       Y       Z

    1+    0.000000    0.000000     1+  SINGLET     ????
    2     1.914982    1.914982     1   TRIPLET     ????
    3     2.267356    2.267356     2   SINGLET     ????
    4     2.464879    2.464879     2   TRIPLET     ????
    5     2.826510    2.826510     3   TRIPLET     ????
    6     3.365178    3.365178     4   TRIPLET     ????
    7     3.392490    3.392490     3   SINGLET     ????  0.2399  0.1995  0.0013
    8     3.669238    3.669238     4   SINGLET     ????  2.0363  2.3850  0.0129
    9     3.872141    3.872141     5   SINGLET     ????  1.8024  1.5312  0.0100
 The "+" symbol indicates the root used.



 TOTAL CPU TIME:            40.29 SECONDS

 == MOPAC DONE ==

Рассчитаем длины волн, при которых происходят эти переходы.

\[E = h \nu = \frac{h c}{\lambda};\quad \lambda [nm] = \frac{h c}{E} = \frac{1239.84}{E [eV]}\]

Длины волн в нм, отсортированные по возрастанию:

320.19
337.9
365.47
368.43
438.65
503
546.82
647.44

Math(r'F(k) = _{-}^{} f(x) e^{2i k} dx')

3. Определение геометрии хинона и его дианиона

Определим геометрию хинона с помощью obgen:

$ echo 'O=C1C=CC(=O)C=C1   quinone' > qui.smi
$ obgen qui.smi > qui.mol
$ babel -imol qui.mol -opdb qui.pdb

и MOPAC:

$ babel -imol qui.mol -omop qui_opt_pm6.mop -xk "PM6"
$ MOPAC2009.exe qui_opt_pm6.mop
$ babel -imopout qui_opt_pm6.out -opdb qui_opt_pm6.pdb
In [30]: 
cmd.reinitialize()
cmd.bg_color('white')
cmd.load(Dir+'qui.pdb','qui')              # green
cmd.load(Dir+'qui_opt_pm6.pdb','qui_opt')  # cyan
cmd.set('stick_radius','0.1')
cmd.set('sphere_scale','0.2')
cmd.show('sticks','qui*')
cmd.show('sphere','qui*')
cmd.select('o','qui_opt and n. O')
cmd.color('orange','o')
In [31]: 
MakeImage()
Image(ImPath)
Out[31]:
In [32]: 
MakeImage()
Image(ImPath)
Out[32]:

Обе программы выдали очень похожие вполне плоские структуры.

Для определения геометрии дианиона сначала создадим .mop файл

$ cp qui_opt_pm6.mop qui-2_opt_pm6.mop
$ vim qui-2_opt_pm6.mop
$ cat qui-2_opt_pm6.mop

PM6 CHARGE=-2
quinone

O(-)   0.92680 1  0.07740 1 -0.40960 1
C   2.10940 1  0.08290 1 -0.09680 1
C   2.84530 1  1.34830 1  0.09810 1
C   4.13720 1  1.35350 1  0.44030 1
C   4.88440 1  0.09590 1  0.63800 1
O(-)   6.06670 1  0.10130 1  0.95150 1
C   4.14840 1 -1.16880 1  0.44240 1
C   2.85650 1 -1.17550 1  0.10030 1
H   2.28300 1  2.26140 1 -0.05060 1
H   4.69190 1  2.27310 1  0.58730 1
H   4.71120 1 -2.08300 1  0.59100 1
H   2.30230 1 -2.09390 1 -0.04690 1

$ MOPAC2009.exe qui-2_opt_pm6.mop
$ babel -imopout qui-2_opt_pm6.out -opdb qui-2_opt_pm6.pdb
In [37]: 
cmd.reinitialize()
cmd.bg_color('white')
cmd.load(Dir+'qui_opt_pm6.pdb','qui_opt')      # green
cmd.load(Dir+'qui-2_opt_pm6.pdb','qui-2_opt')  # cyan
cmd.set('stick_radius','0.1')
cmd.set('sphere_scale','0.1')
cmd.show('sticks','qui*')
cmd.show('sphere','qui*')
cmd.select('o','qui-2_opt and n. O')
cmd.color('orange','o')
In [38]: 
MakeImage()
Image(ImPath)
Out[38]:
In [40]: 
MakeImage()
Image(ImPath)
Out[40]:

Структура по-прежнему плоская, но немного увеличились длины C-O связей.

4. Тиминовые димеры

  1. Оптимизируем геометрию тиминового димера с зарядом 0:
In [22]: 
%%bash
export PATH=${PATH}:/home/preps/golovin/progs/bin
export MOPAC_LICENSE=/home/preps/golovin/progs/bin
babel -ipdb td.pdb -omop td.mop -xk "PM6"  # pdb файл пришлось пересохранить через pymol
MOPAC2009.exe td.mop
babel -imopout td.out -opdb td0.pdb

1 molecule converted
14 audit log messages 
1 molecule converted
23 audit log messages
  1. Оптимизируем геометрию из предыдущего пункта с зарядом +2:
In [27]: 
%%bash
export PATH=${PATH}:/home/preps/golovin/progs/bin
export MOPAC_LICENSE=/home/preps/golovin/progs/bin
babel -ipdb td0.pdb -omop td1.mop -xk "PM6"
sed -i 's/PM6/PM6 CHARGE=+2/' td1.mop
MOPAC2009.exe td1.mop
babel -imopout td1.out -opdb td1.pdb -xk "PM6"

1 molecule converted
14 audit log messages 
1 molecule converted
23 audit log messages
  1. Оптимизируем геометрию из предыдущего пункта с зарядом 0:
In [28]: 
%%bash
export PATH=${PATH}:/home/preps/golovin/progs/bin
export MOPAC_LICENSE=/home/preps/golovin/progs/bin
babel -ipdb td1.pdb -omop td2.mop -xk "PM6"
sed -i 's/PM6 CHARGE=+2/PM6/' td2.mop
MOPAC2009.exe td2.mop
babel -imopout td2.out -opdb td2.pdb -xk "PM6"

1 molecule converted
14 audit log messages 
1 molecule converted
23 audit log messages

Структура i:

In [8]: 
MakeImage()
Image(ImPath)
Out[8]:

Структура ii:

In [10]: 
MakeImage()
Image(ImPath)
Out[10]:

Структура iii:

In [11]: 
MakeImage()
Image(ImPath)
Out[11]:

Энергии состояний:

In [30]: 
%%bash
grep "TOTAL ENERGY" td.out
grep "TOTAL ENERGY" td1.out
grep "TOTAL ENERGY" td2.out

          TOTAL ENERGY            =      -3273.58208 EV
          TOTAL ENERGY            =      -3253.90802 EV
          TOTAL ENERGY            =      -3273.68999 EV

При добавлении обратно электронов не образуются опять тиминовый димер, так как энергия двух тиминув чуть ниже, чем энергия димера.