Для запуска ipython notebook надо запустить специальную версию python в терминале, лучше это делать из той директории где будут лежать файлы практикума.
Форматирование задаётся Markdown: http://daringfireball.net/projects/markdown/basics

сведения о Морас можно получить здесь: http://openmopac.net/manual/index.html

Работа проводится с помощью putty и обработки pymol.

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


Задание 1. Оптимизация 3D структур.


установим переменные: 

export PATH=${PATH}:/home/preps/golovin/progs/bin
export MOPAC_LICENSE=/home/preps/golovin/progs/bin
In [1]: 
#Для работы с Pymol: прокси и загрузка картинок:
from xmlrpclib import ServerProxy
from IPython.display import Image
In [4]: 
cmd = ServerProxy(uri="http://localhost:9123/RPC2")
import time

нашел smiles порфирина и записал в файл 1.smi строку "c1cc2cc3ccc(cc4ccc(cc5ccc(cc1n2)[nH]5)n4)[nH]3 porphyrin"

далее пользуясь пакетом openbabel можно сгенерировать файл со структурой молекулы

obgen 1.smi > 1.mol

Получившаяся молекула ппредставляет собой просто изображение в пространстве того, что закодировано в smiles, в частности она несет лишние два водорода, а ее форма более-менее произвольна (не плоская).

удалив водороды и сохранив файл как pdb-файл, который далее надо переформатировать, используя модель PM6.

babel -ipdb 1.pdb -omop 1_opt.mop -xk "PM6"

где 1.pdb - входной файл. запустим MOPAC для расчета структуры в пространстве:

MOPAC2009.exe 1_opt.mop

переформатируем файл обратно для визуализации

babel -imopout 1_opt.out -opdb 1_opt.pdb

и повторим действия для модели AM1.

babel -ipdb myfile.pdb -omop  1_opt.mop -xk "AM1"

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

In [5]: 
cmd.ray(600,400)
cmd.png('tmp\mypng1.png')
time.sleep(2)

Image(filename='tmp\mypng1.png')
Out[5]:

Рис. 1. Изображения порфирина, полученного с помощью MOPAC с использованием моделей PM6 (зеленым) и AM1 (голубым).
Можно заметить, что хотя структуры получены из одного и того же файла, и они очень похожи, но все же, они не идеально плоские, и видна небольшая разница в положении атомов азотсодержащих колец. Тем не менее, визуально отличие получилось весьма незначительным.

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


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

cis c.i.=4 meci oldgeo
some description


получаем выходной файл, где в конце приведено следующее:

...
STATE       ENERGY (EV)        Q.N.  SPIN   SYMMETRY     POLARIZATION        
     ABSOLUTE     RELATIVE                            X       Y       Z    
                                                                           
1+    0.000000    0.000000     1+  SINGLET     ????                        
2     1.913361    1.913361     1   TRIPLET     ????                        
3     2.266459    2.266459     2   SINGLET     ????                        
4     2.463841    2.463841     2   TRIPLET     ????                        
5     2.824794    2.824794     3   TRIPLET     ????                        
6     3.362879    3.362879     4   TRIPLET     ????                        
7     3.390595    3.390595     3   SINGLET     ????  0.2337  0.2042  0.0005
8     3.669119    3.669119     4   SINGLET     ????  2.0838  2.3506  0.0070
9     3.871059    3.871059     5   SINGLET     ????  1.7797  1.5674  0.0025
The "+" symbol indicates the root used.

что по формулам \[E = h*\nu = h * \frac{c}{\lambda} {\Longrightarrow} \lambda = \frac{h*c}{E}\]

Следовательно, при энергии измеряемой в эВ, длину волны в микрометрах считаем так:

\[\lambda = \frac{1.23984193}{E}\]

в результате энергии, соответственно, составляют 
    0.64755539132478520878      (в микрометрах)        648  (nm)
    0.54686572029192238103                             547   
    0.50300091687661518363                             503
    0.43876732069026898804                             439
    0.36848896695510548598            или              368
    0.36552633363365346515                             366
    0.33785478730419636443                             338
    0.32024317046031009690                             320

то есть должен поглощаться и УФ, и видимый свет.

Задание 3. геометрия хинона и его дианиона.


повторим наши действия для хинона (из "O=C1C=CC(=O)C=C1" получим .mop , а затем .pdb файлы),
в результате obgen получаем

In [7]: 
cmd.ray(600,400)
cmd.png('tmp\mypng2.png')
time.sleep(2)

Image(filename='tmp\mypng2.png')
Out[7]:
In [12]: 
#снова
#babel -ipdb quinone.pdb -omop quinone_1opt.mop -xk "PM6"
#MOPAC2009.exe quinone_1opt.mop
#babel -imopout quinone_1opt.out -opdb quinone_1opt.pdb

cmd.do('''
reinit
bg white
load Practice2\quinone_1opt.pdb
set sphere_scale, 0.2
set stick_radius, 0.1
show spheres, all
show sticks, all
set valence, on
''')

cmd.ray(600,400)
cmd.png('tmp\mypng3.png')
time.sleep(2)

Image(filename='tmp\mypng3.png')
Out[12]:

Различие между результатами внешне незначительное, хотя мы видим, что результат MOPAC дает представление в виде пи-системы,
то есть, возможно, слегка различаются углы при атомах углерода в кольце.



Для рассмотрения дианиона этой молекулы впишем в первую строку файла 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

далее снова MOPAC2009.exe q_anion.mop babel -imopout q_anion.out -opdb q_anion.pdb

In [14]: 
cmd.do('''
reinit
bg white
load Practice2\q_anion.pdb
set sphere_scale, 0.2
set stick_radius, 0.1
show spheres, all
show sticks, all
set valence, on
''')

cmd.ray(600,400)
cmd.png('tmp\mypng4.png')
time.sleep(2)

Image(filename='tmp\mypng4.png')
Out[14]:

Получившийся анион тоже плоский, связи при кислородах двойные, так как заряд записать не получится.
Зато можно записать длину связей, которая действительно немного изменилась - с 1.2 Ангстрем до 1.3 Ангстрем.

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

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

In [3]: 
# Cперва оптимизируем имеющийся димер
babel -ipdb td.pdb -omop td_1.mop -xk "PM6"
MOPAC2009.exe td_1.mop
babel -imopout td_1.out -opdb td_1.pdb
In [3]: 
# затем ионизируем димер, записав в файл заряд (куда конкретно ставить заряды мы не знаем, потому только общий заряд)
babel -ipdb td_1.pdb -omop td_2.mop -xk "PM6"

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

MOPAC2009.exe td_2.mop
babel -imopout td_2.out -opdb td_2.pdb
In [3]: 
# наконец обратная конвертация (сразу снимает заряд, так как в pdb он не пишется) и оптимизируем при заряде 0.
babel -ipdb td_2.pdb -omop td_3.mop -xk "PM6"

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

Что получается:

In [18]: 
cmd.do('''
reinit
bg white
load Practice2/td_1.pdb
set sphere_scale, 0.2
set stick_radius, 0.1
show spheres, all
show sticks, all
set valence, on
''')

cmd.ray(600,400)
cmd.png('tmp\mypng5.png')
time.sleep(2)

Image(filename='tmp\mypng5.png')
Out[18]:
In [19]: 
cmd.do('''
reinit
bg white
load Practice2/td_2.pdb
set sphere_scale, 0.2
set stick_radius, 0.1
show spheres, all
show sticks, all
set valence, on
''')

cmd.ray(600,400)
cmd.png('tmp\mypng6.png')
time.sleep(2)

Image(filename='tmp\mypng6.png')
Out[19]:
In [22]: 
cmd.do('''
reinit
bg white
load Practice2/td_3.pdb
set sphere_scale, 0.2
set stick_radius, 0.1
show spheres, all
show sticks, all
set valence, on
rotate y, -30
''')

cmd.ray(600,400)
cmd.png('tmp\mypng7.png')
time.sleep(2)

Image(filename='tmp\mypng7.png')
Out[22]:

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

Энергии состояний можно посмотреть в файлах .out MOPAC, где есть запись "TOTAL ENERGY" 
Для состояний 1-3
энергия составляет соответственно

    -3273.58217
    -3253.90834
    -3273.69661 (Эв)
    
То есть, ионизация значительно повысила энергию системы, которая затем скатилась к почти исходному состоянию. То, что состояния
обладают все же разной энергией и в реальности очевидно разные, можно понимать, как попадание в локальный минимум.
В норме образованию тиминовых димеров способствует их близкая укладка, то есть фактически стэкинг.
В нашей же системе этот фактор при возврате из состояния 2 в состояние 3 отсутствовал - тимины не были прижаты друг к другу, а их
отталкивание в пространстве, как противоположно заряженных частей одной молекулы, вызвало увеличение расстояния между кольцами.
В этих условиях возврат в состояние 1 уже не происходит, оставшаяся связь рвется, так как можно образовать две pi-системы, 
что оказывается более выгодно.