Учебный сайт

Бредихина Данилы

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

Запустим PyMol для работы в режиме сервера.

In [9]: 
from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys, time

# pymol launching
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]

# __main__.pymol_argv = [ 'pymol', '-cp' ]  # to disable GUI
import pymol
from pymol import cmd
pymol.finish_launching()
 
from IPython.display import Image

defaultImage = '/tmp/pymolimg.png'
def prepareImage(width=300, height=300, sleep=2, filename=defaultImage):
    ## To save the rendered image
    cmd.ray(width, height)
    cmd.png('/tmp/pymolimg.png')
    time.sleep(sleep)

# Define some shortcuts
def focus(x):
    cmd.center(x)
    cmd.zoom(x)

Построение пространственной структуры порфирина

На основе SMILES-аннотации порфирина создадим его пространственную структуру. Затем загрузим полученную структуру в PyMol.

In [21]: 
%%bash
echo "c1cc2cc3ccc(cc4ccc(cc5ccc(cc1n2)[nH]5)n4)[nH]3  porphyrin" > porphyrin.smi
obgen porphyrin.smi > porphyrin.mol
In [27]: 
cmd.do('''
load porphyrin.mol
show spheres, all
set sphere_scale, 0.2
set valence, on
''')

prepareImage(); Image(defaultImage)
Out[27]:

Удалим лишние водороды (в GUI, т. к. имена атомов одних элементов одинаковы) и сохраним файл в формате pdb.

In [13]: 
cmd.save('porphyrin.pdb', 'porphyrin')

Переформатируем координаты в porphyrin.pdb для использования как входного файла для Mopac. С помощью аргумента -xk "PM6" зададим соответствующий тип параметризации.

In [17]: 
%%bash
babel -ipdb porphyrin.pdb -omop porphyrin_1opt.mop -xk "PM6"

Запустим Mopac. Файл porphyrin_1opt.out переформатируем в pdb:

In [20]: 
%%bash
MOPAC2009.exe porphyrin_1opt.mop
babel -imopout porphyrin_1opt.out -opdb porphyrin_1opt.pdb

Загрузим структуру porphyrin_1opt.pdb в PyMol.

In [29]: 
cmd.do('''
reinit
load porphyrin_1opt.pdb
set valence, on
show spheres, all
set sphere_scale, 0.2
''')

prepareImage(); Image(defaultImage)
Out[29]:

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

In [30]: 
cmd.do('''
rotate x, 86
rotate z, -90
''')
prepareImage(); Image(defaultImage)
Out[30]:

Проведём оптимизацию с другой параметризацией – AM1.

In [24]: 
%%bash
babel -ipdb porphyrin.pdb -omop porphyrin_AM1.mop -xk "AM1"
MOPAC2009.exe porphyrin_AM1.mop
babel -imopout porphyrin_AM1.out -opdb porphyrin_AM1.pdb
In [32]: 
cmd.do('''
reinit
load porphyrin_AM1.pdb
set valence, on
show spheres, all
set sphere_scale, 0.2
''')

prepareImage(); Image(defaultImage)
Out[32]:

Полученная структура также является плоской и не имеет значительных отличий от предыдущей (porphyrin_1opt.pdb). Сравнение файлов porphyrin_AM1.mop и porphyrin_1opt.mop приводит к такому же выводу.

Структура, полученная с помощью obgen, визуально отличается нарушенной (неплоской) геометрией от этих структур.

In [33]: 
cmd.do('''
reinit
load porphyrin.mol
set valence, on
show spheres, all
set sphere_scale, 0.2
rotate y, 94
''')

prepareImage(); Image(defaultImage)
Out[33]:

Расчёт возбуждённых состояний порфирина

Рассчитаем возбуждённые состояния для порфирина (файл porphyrin_1opt.mop). Для указания программе Mopac необходимости расчёта возбуждённого состояния добавим в конец файла соответствуюшие строки.

In [45]: 
%%bash
cp porphyrin_1opt.mop 1_opt_spectr.mop
echo "
cis c.i.=4 meci oldgeo
For porphyrin spectrum" >> 1_opt_spectr.mop
MOPAC2009.exe 1_opt_spectr.mop

В конце файла можно найти значения энергий для электронных переходов.

In [1]: 
%%bash
tail -n19 1_opt_spectr.out | head -n12

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

    1+    0.000000    0.000000     1+  SINGLET     ????
    2     1.914650    1.914650     1   TRIPLET     ????
    3     2.267178    2.267178     2   SINGLET     ????
    4     2.464890    2.464890     2   TRIPLET     ????
    5     2.825739    2.825739     3   TRIPLET     ????
    6     3.364665    3.364665     4   TRIPLET     ????
    7     3.391936    3.391936     3   SINGLET     ????  0.2018  0.2373  0.0011
    8     3.669748    3.669748     4   SINGLET     ????  2.3998  2.0267  0.0129
    9     3.871564    3.871564     5   SINGLET     ????  1.5364  1.7970  0.0087

Рассчитаем длину волн, при которых происходят эти переходы. Сначала выведем формулу для расчёта, а затем применим её к полученным значениям энергии.

$$ E = h \cdot \nu = h \cdot \frac{c}{\lambda} \Rightarrow \lambda = \frac{h \cdot c}{E} $$

Если размерность энергии – эВ, то решение в численном виде можно записать так:

$$ \lambda = \frac{1.23984193}{E} $$

При этом длина волны λ будет иметь размерность микрометров (μm).

In [12]: 
%%bash
tail -n15 1_opt_spectr.out | head -n8 | awk '{print $2}' | \
    xargs -n1 echo "1.23984193 /" | bc -l > porphyrin_spectr.miM
cat porphyrin_spectr.miM

.64755539132478520878
.54686572029192238103
.50300091687661518363
.43876732069026898804
.36848896695510548598
.36552633363365346515
.33785478730419636443
.32024317046031009690

Для удобства восприятия выразим полученные значения в нанометрах, округлим до целых и отсортируем по возрастанию.

In [13]: 
cat porphyrin_spectr.miM | xargs -n1 echo "1000*" | bc | \
    cut -f1 -d'.' | sort -n | xargs -I 'q' echo 'q' nm

320 nm
337 nm
365 nm
368 nm
438 nm
503 nm
546 nm
647 nm

Пожалуй, можно сказать, что полученные значения согласуются с литературными данными. Так, в книге Macro to Nano Spectroscopy указывается, что в спектре поглощения порфирина отчётливо выделяются два участка с максимумами поглощения — в ближнем ультрафиолете и в видимой части спекта.

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

Для молекулы хинона определим геометрию с помощью obgen и Мopac.

In [121]: 
%%bash
echo "O=C1C=CC(=O)C=C1  quinone" > quinone.smi
obgen quinone.smi > quinone.mol
In [51]: 
cmd.do('''
reinit
bg white
set ray_trace_mode, 3
load quinone.mol
show spheres, all
set sphere_scale, 0.2
show sticks, all
set stick_radius, 0.1
zoom all
rotate y, -30
set valence, on
''')

prepareImage(); Image(defaultImage)

# Obtained with obgen
Out[51]:
In [123]: 
cmd.save('quinone.pdb', 'quinone')
In []: 
%%bash
babel -ipdb quinone.pdb -omop quinone_1opt.mop -xk "PM6"

MOPAC2009.exe quinone_1opt.mop
babel -imopout quinone_1opt.out -opdb quinone_1opt.pdb
In [50]: 
cmd.do('''
reinit
bg white
set ray_trace_mode, 3
load quinone_1opt.pdb
show spheres, all
set sphere_scale, 0.2
show sticks, all
set stick_radius, 0.1
rotate y, 30
zoom all
rotate y, -30
set valence, on
''')

prepareImage(400,350); Image(defaultImage)

# Obtained with Mopac
Out[50]:

Полученные структуры имеют плоскую геометрию и сходны между собой. Можно заметить изменение угла CC(=O)C: в полученной с помощью Mopac структуре он равен 116.2° (вместо 117.4° в структуре, полученной с помощью obgen).

Определим также геометрию дианиона этой молекулы.

Для начала в первую строчку файла .mop добавим CHARGE=-2. Затем явным способом укажем те атомы, на которых должен находиться отрицательный заряд, – атомы кислорода.

In [1]: 
%%bash
cp quinone_1opt.mop quinone_anion.mop
sed -i '' -e 's/PM6/PM6 CHARGE=-2/' -e 's/^O/O(-)/g' quinone_anion.mop
cat quinone_anion.mop

PM6 CHARGE=-2
quinone.pdb

C   2.54600 1  0.00200 1 -0.66900 1
C   3.08600 1  1.26700 1 -0.13200 1
C   4.03400 1  1.27200 1  0.81000 1
C   4.58200 1  0.01400 1  1.35500 1
C   4.04300 1 -1.25000 1  0.81700 1
C   3.09500 1 -1.25700 1 -0.12500 1
O(-)   1.67900 1 -0.00400 1 -1.53100 1
O(-)   5.44900 1  0.02000 1  2.21700 1
H   2.67400 1  2.18000 1 -0.54200 1
H   4.44100 1  2.19200 1  1.21500 1
H   4.45600 1 -2.16400 1  1.22700 1
H   2.68800 1 -2.17500 1 -0.52900 1
In [7]: 
%%bash
MOPAC2009.exe quinone_anion.mop
babel -imopout quinone_anion.out -opdb quinone_anion.pdb
In [40]: 
cmd.do('''
reinit
bg white
set ray_trace_mode, 3
load quinone_anion.pdb
show spheres, all
set sphere_scale, 0.2
show sticks, all
set stick_radius, 0.1
rotate y, 30
zoom all
rotate y, -30
set valence, on
''')

prepareImage(400,350); Image(defaultImage)

# For anion
Out[40]:

Геометрия структуры дианиона также плоская. Длина связей C=O увеличилась с 1.2 Å до 1.3 Å. Длины связей C-C в структуре дианона равны 1.4 Å, в то время как в структуре хинона они составляли 1.3 Å и 1.5 Å.

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

Известно, что ультрафиолет может превращать тимины в тиминовые димеры. Также известно, что ДНК-фотолиаза при облучении ультрафиолетом восстанавливает основания тиминов до нормального состояния.

Так как вычисления возбуждённых состояний в Mopac затруднены, сымитируем возбуждение, ионизируя оба кольца, т. е. указывая заряд системы +2. Полученное возбуждённое состояние снова оптимизируем при заряде 0.

In [6]: 
%%bash
# 1. Geometry optimization at 
babel -ipdb td.pdb -omop td_stage1.mop -xk "PM6"
MOPAC2009.exe td_stage1.mop
babel -imopout td_stage1.out -opdb td_stage1.pdb
In []: 
%%bash
# 2. Geometry optimization at +2
babel -ipdb td_stage1.pdb -omop td_stage2.mop -xk "PM6"

sed -i 's/PM6/PM6 CHARGE=+2/' td_stage2.mop

MOPAC2009.exe td_stage2.mop
babel -imopout td_stage2.out -opdb td_stage2.pdb
In []: 
%%bash
# 3. Geometry optimization at 0
babel -ipdb td_stage2.pdb -omop td_stage3.mop -xk "PM6"

MOPAC2009.exe td_stage3.mop
babel -imopout td_stage3.out -opdb td_stage3.pdb

Отобразим полученные структуры в PyMol.

In [12]: 
def showPDB(pdb):
    cmd.do('''
    reinit
    bg white
    set antialias, 2
    set ray_trace_mode, 3
    load ''' + pdb + '''
    show spheres, all
    set sphere_scale, 0.2
    show sticks, all
    set stick_radius, 0.1
    zoom all
    set valence, on
    ''')

showPDB('td_stage1.pdb')
cmd.rotate('y', '-30')
prepareImage(); Image(defaultImage)

# Initial thymine dimer (charge: 0)
Out[12]:
In [8]: 
showPDB('td_stage2.pdb')
cmd.rotate('y', '-30')
prepareImage(); Image(defaultImage)

# Excited state (charge: +2)
Out[8]:
In [13]: 
showPDB('td_stage3.pdb')
cmd.rotate('y', '-30')
prepareImage(); Image(defaultImage)

# Thymines (charge 0)
Out[13]:

Как видно, при возбуждении системы (заряд +2) переход димера в тимины, как это ожидалось, не наблюдается. Повторим этот этап анализа с использованием другой параметризации – AM1 (вместо PM6).

In []: 
%%bash
# 2. Geometry optimization at +2
babel -ipdb td_stage2.pdb -omop td_stage2x.mop -xk "AM1"

sed -i 's/PM6/PM6 CHARGE=+2/' td_stage2x.mop

MOPAC2009.exe td_stage2x.mop
babel -imopout td_stage2x.out -opdb td_stage2x.pdb
In [10]: 
showPDB('td_stage2x.pdb')
cmd.rotate('y', '-30')
prepareImage(); Image(defaultImage)

# Excited state (charge: +2); AM1 parameters
Out[10]:

Таким образом, удалось увидеть переход из димера в тимины при возбуждении системы.

Сравним энергии молекул в каждом из состояний.

In [6]: 
%%bash
for file in `echo td_stage*.out`; do
    echo -ne $file"\t"
    grep "TOTAL ENERGY" $file | awk '{print $4}';
done

td_stage1.out	-3273.58217
td_stage2.out	-3253.90834
td_stage2x.out	-3555.39717
td_stage3.out	-3273.69661

Различие величин энергий для систем с зарядом +2 (стадии 2 и 2x) свидетельствует о значительном вкладе выбранной параметризации в значение общей энергии.

Ссылки на файлы

△