Домашняя работа 2

Основная ссылка: H:\Term8\Practice2

После некоторых интеллектуальных усилий справляюсь с запуском IPython на kodomo и "перетуннелированием" его порта на порт локали. При ближайшем рассмотрени процесс очень логичный и простой.

1. Захожу на сервер в режиме X11 trusty forwarding (кстати, стоит проверить, будет ли работать опция -X - без trusty). Этот режим предоставляет возможность перенаправления графического окна на порт.

ssh -Y agalicina@kodomo.fbb.msu.ru

2. Запускаю IPython с опцией "без браузера". Выход графического окна уже перенаправлен на какой-то порт сервера. Запоминаю номер этого порта

ipython notebook --no-browser

Внимание! К порту сервера имеет доступ любой, кто знает номер этого порта! Таким образом, IPython notebook - очень ненадежный метод, если не запаролить его. Для этого получаем хэш-код желаемого пароля через IPython:

In []: 
from IPython.lib import passwd
passwd()

Получаю хэш, который копирую и вставляю в файл ~/.ipython/profile_USERNAME/ipython_notebook_config.py

# Password to use for web authentication
c = get_config()
c.NotebookApp.password = PASSWORD_HASH_HERE

Для того, чтобы найти файл ipython_notebook_config.py нужно сначала создать пользователя на сервере. Для этого достаточно запустить команду в командной строке:

ipython profile create USERNAME

Теперь запускать ноубук нужно с указанием пользователя:

ipython notebook --no-browser --profile=USERNAME

И да, теперь вход в ноутбук запрашивает тот пароль, для которого изначально был получен hash.

3. Перенаправляю выход порта с сервера, номер которого я запомнила (XXXX), на порт локального компьютера (8888):

ssh -L 8888:localhost:XXXX agalicina@kodomo.fbb.msu.ru

4. Теперь я могу смотреть выход своего порта через браузер. И более того, запускать X11 и перенаправлять их к себе.

http://localhost:8888/

5. Запускаю pymol прямо из ноутбука в режиме RPC-сервера как модуль.

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

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

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


import pymol
 
pymol.finish_launching()
 
from pymol import cmd

Если все прошло успешно, то должно появиться окно с pymol, запущенным на сервере.

Ну, и последнее, что нужно сделать - установить переменные для дальнейшей работы с программами из практикума (по каким-то причинам это не работает, ipython notebook не видит новые переменные и не запускает mopac с ошибкой "Command not found"):

In [11]: 
%%bash
export PATH=${PATH}:/home/preps/golovin/progs/bin
export MOPAC_LICENSE=/home/preps/golovin/progs/bin

Задание 1

Нахожу SMILES-нотацию порфирина, записываю в файл на kodomo, перевожу в 3D-структуру и отображаю в PyMol

In [4]: 
%%bash
echo "C1=CC2=NC1=CC3=NC(=CC4=NC(=CC5=NC(=C2)C=C5)C=C4)C=C3" > porphyrin.smi
obgen porphyrin.smi > porphyrin.mol
In [5]: 
cmd.load("/home/students/y11/agalicina/Term8/Practice2/porphyrin.mol")
cmd.do(
'''
ray 500,500
png images/p1.png
''')
Image(filename='/home/students/y11/agalicina/Term8/Practice2/images/p1.png')
Out[5]:

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

In [6]: 
cmd.delete('all')
In [7]: 
%%bash
echo "c1c/2[nH]c(c1)/cc/3\nc(/cc/4\[nH]/c(c\c5n/c(c2)/C=C5)/cc4)C=C3" > porphyrin.smi
obgen porphyrin.smi > porphyrin.mol
In [8]: 
cmd.load("/home/students/y11/agalicina/Term8/Practice2/porphyrin.mol")
cmd.do(
'''
set valence,1
set valence_mode, 0
ray 500,500
png images/p2.png
''')

Image(filename='/home/students/y11/agalicina/Term8/Practice2/images/p2.png')
Out[8]:
In [9]: 
# Вид сбоку
cmd.set_view("(-0.7655993103981018,\
 0.43343475461006165,\
 -0.47538262605667114,\
 -0.4332069754600525,\
 0.1989416778087616,\
 0.8790627717971802,\
 0.475589781999588,\
 0.8789519667625427,\
 0.03545708209276199,\
 0.0,\
 0.0,\
 -27.551538467407227,\
 1.7244553565979004,\
 5.390379428863525,\
 0.08007892966270447,\
 21.721843719482422,\
 33.38123321533203,\
 -20.0)")
cmd.do(
'''
ray 500,500
png images/p3.png
''')


Image(filename='/home/students/y11/agalicina/Term8/Practice2/images/p3.png')
Out[9]:
In [11]: 
# Ненужных водородов нет, просто сохраняю структуру в .PDB:
cmd.save("porphyrin.pdb")

Создаю входной файл для Mopac из координат pdb:

In [12]: 
%%bash
babel -ipdb porphyrin.pdb -omop 1_opt.mop -xk "PM6"

Запускаю MOPAC:

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

Использование другого параметра:

In [25]: 
%%bash
babel -ipdb porphyrin.pdb -omop 2_opt.mop -xk "AM1"
MOPAC2009.exe 2_opt.mop
babel -imopout 2_opt.out -opdb 2_opt.pdb

Смотрю, что получилось с этими двумя опциями и пытаюсь сравнить:

In [19]: 
# Вид на кольцо
cmd.set_view("(0.03581596538424492,\
 0.9812870621681213,\
 -0.18917885422706604,\
 -0.999190628528595,\
 0.03171350061893463,\
 -0.024669542908668518,\
 -0.018209027126431465,\
 0.18991151452064514,\
 0.9816315770149231,\
 0.0,\
 0.0,\
 -27.551538467407227,\
 1.7244553565979004,\
 5.390379428863525,\
 0.08007892966270447,\
 21.721843719482422,\
 33.38123321533203,\
 -20.0)")
cmd.delete('all')
cmd.load("porphyrin.mol")
cmd.load("1_opt.pdb")
cmd.load("2_opt.pdb")
cmd.do(
'''
ray 500,500
png images/p_mopac.png
''')

Image(filename='/home/students/y11/agalicina/Term8/Practice2/images/p_mopac.png')
Out[19]:
In [24]: 
cmd.do(
'''
hide everything, porphyrin
ray 500,500
png images/p_mopac1.png
''')

Image(filename='/home/students/y11/agalicina/Term8/Practice2/images/p_mopac1.png')
Out[24]:
In [25]: 
cmd.set_view("(-0.7655993103981018,\
 0.43343475461006165,\
 -0.47538262605667114,\
 -0.4332069754600525,\
 0.1989416778087616,\
 0.8790627717971802,\
 0.475589781999588,\
 0.8789519667625427,\
 0.03545708209276199,\
 0.0,\
 0.0,\
 -27.551538467407227,\
 1.7244553565979004,\
 5.390379428863525,\
 0.08007892966270447,\
 21.721843719482422,\
 33.38123321533203,\
 -20.0)")
cmd.do(
'''
show lines, porphyrin
ray 500,500
png images/p_mopac2.png
''')


Image(filename='/home/students/y11/agalicina/Term8/Practice2/images/p_mopac2.png')
Out[25]:

Как можно увидеть, MOPAC проводит оптимизацию молекулы, причем зависимость от метода есть, но малозаметная. Правда, то, что молекула порфирина не совсем плоская, программа не исправила.

Задание 2

In [33]: 
%%bash
cp 1_opt.mop 1_opt_spectr.mop
printf "\ncis c.i.=4 meci oldgeo\nsome description" >> 1_opt_spectr.mop
MOPAC2009.exe 1_opt_spectr.mop
tail -n 20 1_opt_spectr.out | head -n 13

Выдача:

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

    1+    0.000000    0.000000     1+  SINGLET     ????
    2     1.914250    1.914250     1   TRIPLET     ????
    3     2.267103    2.267103     2   SINGLET     ????
    4     2.465049    2.465049     2   TRIPLET     ????
    5     2.825994    2.825994     3   TRIPLET     ????
    6     3.364647    3.364647     4   TRIPLET     ????
    7     3.391904    3.391904     3   SINGLET     ????  0.1657  0.2746  0.0003
    8     3.669233    3.669233     4   SINGLET     ????  2.7405  1.6954  0.0004
    9     3.871653    3.871653     5   SINGLET     ????  1.2768  2.0618  0.0035

Считаем по формулам: \(E=h\cdot\nu \rightarrow \lambda=\dfrac{c}{\nu}=\dfrac{c\cdot h}{E}\)

Переведем в атомные единицы измерения:

\(a_{0}=\dfrac{c\cdot \hbar}{E_{h}}\)

\([bohr]=\dfrac{[ат.ед.скорости]}{[hartree]}\)

\(a_{0}=\dfrac{137\cdot1}{E_{h}}\), где:

\(E_{h}[hartree]=27,211 3845(23)\cdot E[eV]\),

\(a_{0}[bohr]=0.052917721092(17)\cdot \lambda[nm]\)

\([E]=[eV]\)

In [20]: 
def compute_lambda(E):
    if not E: 
        print "Zero enegry"
        raise Exception
    Eh = 27.211*E
    a0 = 137/Eh
    l = a0/0.0529
    return l

energy_list = [1.914250,2.267103,2.465049,2.825994,3.364647,3.391904,3.669233,3.871653]
state_transitions = [energy_list[i+1]-energy_list[i] for i in range(0,7)]

lambda_full_transitions = [compute_lambda(x) for x in energy_list]
lambda_state_transitions = [compute_lambda(x) for x in state_transitions]

print "\nWavelengths in nm for transitions from zero energy state:\n"+str(lambda_full_transitions)
print "\nWavelengths in nm for transitions between states:\n"+str(lambda_state_transitions)


Wavelengths in nm for transitions from zero energy state:
[49.71892580991721, 41.980648312685396, 38.6095585652188, 33.67822215179297, 28.286608886945352, 28.05930053787903, 25.938514597365174, 24.58238218446592]

Wavelengths in nm for transitions between states:
[269.7283393697488, 480.8101893023049, 263.68131912516867, 176.68973111007276, 3491.743542269335, 343.18247904703037, 470.1830537083006]

Задание 3

Провести всю процедуру, что и для порфирина, но для заряженной молекулы O=C1C=CC(=O)C=C1 - дианиона хинона.

In [34]: 
%%bash
echo "O=C1C=CC(=O)C=C1" > compound.smi
obgen compound.smi > compound.mol
In [137]: 
cmd.delete("all")
cmd.load("compound.mol")

cmd.set_view("(-0.042833708226680756,\
 0.9922202229499817,\
 -0.11688100546598434,\
 -0.9975786209106445,\
 -0.036062151193618774,\
 0.05944821238517761,\
 0.05477065593004227,\
 0.11914627254009247,\
 0.9913651347160339,\
 0.0,\
 0.0,\
 -14.575573921203613,\
 3.622666835784912,\
 -0.03835010528564453,\
 0.1896916776895523,\
 11.491493225097656,\
 17.65965461730957,\
 -20.0)")

cmd.do(
'''
set valence,1
set valence_mode, 0
ray 500,500
png images/c1.png
''')
cmd.save("compound.pdb")
Image(filename='images/c1.png')
Out[137]:

Получили обычный незаряженный хинон. Попробую получить заряженную форму с помощью MOPAC.

In [20]: 
%%bash
babel -ipdb compound.pdb -omop 3_opt1.mop -xk "PM6"
printf "PM6 CHARGE=-2\ngg\n\n" > 3_opt.mop
grep O 3_opt1.mop >> 3_opt.mop
printf "\n" >> 3_opt.mop
less 3_opt1.mop >> 3_opt.mop
In [24]: 
%%bash
MOPAC2009.exe 3_opt.mop
babel -imopout 3_opt.out -opdb 3_opt.pdb
In [26]: 
cmd.delete("all")
cmd.load("3_opt.pdb")
cmd.do(
'''
ray 500,500
png images/c2.png
''')
Image(filename='images/c2.png')
Out[26]:

Двойная связь у кислорода сохраняется, хотя это неверно. Вообще, двойные связи строит Pymol, и решает, где отображать делокализацию заряда. Попробую при составлении SMILES надо прописать анион кислорода явно:

In [33]: 
%%bash
echo "[O-]C1=CC=C([O-])C=C1" > compound1.smi
obgen compound1.smi > compound1.mol
In [38]: 
cmd.delete("all")
cmd.load("compound1.mol")
cmd.do(
'''
ray 500,500
png images/c3.png
''')
cmd.save("compound1.pdb")
Image(filename='images/c3.png')
Out[38]:
In [40]: 
%%bash
babel -ipdb compound1.pdb -omop 4_opt1.mop -xk "PM6"
printf "PM6 CHARGE=-2\ngg\n\n" > 4_opt.mop
grep O 4_opt1.mop >> 4_opt.mop
printf "\n" >> 4_opt.mop
less 4_opt1.mop >> 4_opt.mop
MOPAC2009.exe 4_opt.mop
babel -imopout 4_opt.out -opdb 4_opt.pdb
In [73]: 
cmd.delete("all")
cmd.load("4_opt.pdb")
cmd.do(
'''
ray 500,500
png images/c4.png
''')
Image(filename='images/c4.png')
Out[73]:

Вижу, что такая особенность положения двойных связей - происходит от самого python. Буду считать, что MOPAC умнее меня и нарисую ball & stick модель.

In [27]: 
cmd.do('bg_color white')
cmd.set('ambient','0')
cmd.set('direct','0.7')
cmd.set('specular','off')
cmd.set("ray_shadow","off")
In [81]: 
cmd.do(
'''
set valence,0
show spheres
set sphere_scale, 0.1
ray 500,500
png images/c5.png
''')
Image(filename='images/c5.png')
Out[81]:

Наложу молекулы с зарядом и без заряда:

In [230]: 
cmd.delete("all")
cmd.set("valence",0)
cmd.load("compound.pdb")
cmd.do("util.cbag compound") # Зеленый цвет для pdb выдачи obgen для SMILES с двойными связями около O
cmd.load("compound1.pdb")
cmd.do("util.cbay compound1") # Желтый цвет для pdb выдачи obgen для SMILES с одинарными связями около O
cmd.load("3_opt.pdb")
cmd.do("util.cbam 3_opt") # Розовый цвет для pdb выдачи MOPAC для SMILES с двойными связями около O
cmd.load("4_opt.pdb")
cmd.do("util.cbab 4_opt") # Синий цвет для pdb выдачи MOPAC для SMILES с одинарными связями около O
cmd.do(
'''
center
zoom
move z,-3
hide lines
show sticks
set stick_radius=0.02
show spheres
set sphere_scale, 0.04
ray 500,500
png images/c6.png
''')
In [231]: 
Image(filename='images/c6.png')
Out[231]:

Как видно, исходная формула SMILES играет большую роль в геометрии построенных молекул.

In [232]: 
cmd.pair_fit("3_opt","compound1")
cmd.pair_fit("4_opt","compound1")
cmd.pair_fit("compound","compound1")
cmd.do("color green, compound") # Зеленый цвет для pdb выдачи obgen для SMILES с двойными связями около O
cmd.do("color yellow, compound1") # Желтый цвет для pdb выдачи obgen для SMILES с одинарными связями около O
cmd.do("color magenta, 3_opt") # Розовый цвет для pdb выдачи MOPAC для SMILES с двойными связями около O
cmd.do("color blue,4_opt") # Синий цвет для pdb выдачи MOPAC для SMILES с одинарными связями около O
cmd.do(
'''
center
zoom
move z,-3

ray 500,500
png images/c7.png
''')
In [234]: 
Image(filename='images/c7.png')
Out[234]:

В этой структуре атомы O окрашены по тому, в какой модели они находятся, нужно помнить, что они по-прежнему сверху и снизу.

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

Следующая отличающаяся модель - полученная из SMILES с двойными связями выдачи MOPAC. Можно заметить, что CO-связи в ней немного короче, чем в оставшихся двух моделях. При этом безнзольное кольцо слегка вытянуто в направлении к O-атомам.

Как ни странно, модель, полученная из SMILES с двойной связью obgen и из SMILES с одинарной MOPAC - практически идентичны.

Задание 4

Дан файл td.pdb тимидинового димера, нужно с помощью MOPAC оценить его энергию в заряженном и нейтральном состоянии.

In [5]: 
cmd.reinitialize()
cmd.load("td.pdb")
In [28]: 
cmd.set_view("(-0.701326310634613,\
 -0.7124418020248413,\
 -0.023816648870706558,\
 -0.3272649645805359,\
 0.3514811098575592,\
 -0.8771309852600098,\
 0.6332755088806152,\
 -0.607361376285553,\
 -0.47966060042381287,\
 0.0,\
 0.0,\
 -20.163301467895508,\
 29.799999237060547,\
 34.61433410644531,\
 60.675331115722656,\
 15.896900177001953,\
 24.429702758789062,\
 -20.0)"
)
cmd.do(
'''
show spheres
set sphere_scale, 0.05
ray 500,500
png images/td1.png
''')
In [9]: 
Image(filename='images/td1.png')
Out[9]:

Видно, что молекула "неровная", хотелось бы её оптимизировать.

I модель: оптимизация исходной структуры

In [21]: 
%%bash
babel -ipdb td.pdb -omop td.mop -xk "PM6"
In [25]: 
%%bash
MOPAC2009.exe td.mop
babel -imopout td.out -opdb td1.pdb # td1 - первая оптимизация, модель без заряда
In [32]: 
cmd.delete("all")
cmd.load("td1.pdb")
In [33]: 
cmd.do(
'''
show spheres
set sphere_scale, 0.05
ray 500,500
png images/td2.png
''')
In [34]: 
Image(filename='images/td2.png')
Out[34]:
In [48]: 
cmd.delete("all")
cmd.load("td1.pdb")
cmd.load("td.pdb")
cmd.do("util.cbag td") # Зеленый для неоптимизированной молекулы
cmd.do("util.cbay td1") # Желтый для оптимизированной молекулы
time.sleep(2)
cmd.do(
'''
show spheres
set sphere_scale, 0.05
ray 500,500
png images/td3.png
''')
In [50]: 
Image(filename='images/td3.png')
Out[50]:

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

Подсчитаю энергию:

In [79]: 
%%bash
babel -imopout td.out -omop td1.mop -xk "PM6"
cp td1.mop td1_spectr.mop
printf "\ncis c.i.=4 meci oldgeo\nsome description">> td1_spectr.mop
MOPAC2009.exe td1_spectr.mop
In [83]: 
%%bash
tail -n 20 td1_spectr.out | head -n 13

Выдача:


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

    1+    0.000000    0.000000     1+  SINGLET     A
    2     4.955172    4.955172     1   TRIPLET     A
    3     5.241099    5.241099     2   TRIPLET     A
    4     5.273208    5.273208     2   SINGLET     A     0.0071  0.0544  0.0020
    5     5.550135    5.550135     3   SINGLET     A     0.0257  0.0960  0.0505
    6     6.162894    6.162894     3   TRIPLET     A
    7     6.305476    6.305476     4   SINGLET     A     0.0137  0.0304  0.0068
    8     6.395284    6.395284     4   TRIPLET     A
    9     6.476016    6.476016     5   SINGLET     A     0.0009  0.0152  0.0064

II модель: оптимизация заряженной структуры I

Теперь можно добавить по положительному заряду в тимидиновые кольца и произвести оптимизацию такой модели возбужденного состояния димера.

In [53]: 
%%bash
babel -ipdb td1.pdb -omop td1.mop -xk "PM6"
printf "PM6 CHARGE=+2\ngg\n\n" > td1_charged.mop # td1_charged - результат первой оптимизации с зарядом
grep N td1.mop | head -n 2 | tail -n 1 >> td1_charged.mop  # заряды располагаем на N3 обоих колец 
grep N td1.mop | tail -n 1 >> td1_charged.mop
printf "\n" >> td1_charged.mop
less td1.mop >> td1_charged.mop
MOPAC2009.exe td1_charged.mop
babel -imopout td1_charged.out -opdb td2_charged.pdb # td2_charged - оптимизация результата предыдущей оптимизации, но с зарядом
In [66]: 
cmd.delete("all")
cmd.load("td2_charged.pdb")
cmd.load("td1.pdb")
cmd.do("util.cbay td1") # Желтый для незаряженной молекулы
cmd.do("util.cbac td2_charged") # Синий для заряженной молекулы
time.sleep(2)
cmd.do(
'''
show spheres
set sphere_scale, 0.05
ray 500,500
png images/td4.png
''')
cmd.rms("td1","td2_charged") # Видим разницу, но небольшую
Out[66]:
0.0022289834450930357
In [63]: 
Image(filename='images/td4.png')
Out[63]:

Теперь подсчитаю энергию:

In [75]: 
%%bash
babel -imopout td1_charged.out -omop td2_charged.mop -xk "PM6"
cp td2_charged.mop td2_charged_spectr.mop
printf "\ncis c.i.=4 meci oldgeo\nsome description">> td2_charged_spectr.mop
MOPAC2009.exe td2_charged_spectr.mop
In [77]: 
%%bash
tail -n 20 td2_charged_spectr.out | head -n 13

Выдача:

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

    1+    0.000000    0.000000     1+  SINGLET     A
    2     4.957464    4.957464     1   TRIPLET     A
    3     5.239384    5.239384     2   TRIPLET     A
    4     5.274602    5.274602     2   SINGLET     A     0.0069  0.0540  0.0021
    5     5.549203    5.549203     3   SINGLET     A     0.0258  0.0970  0.0502
    6     6.164134    6.164134     3   TRIPLET     A
    7     6.306077    6.306077     4   SINGLET     A     0.0138  0.0299  0.0067
    8     6.393083    6.393083     4   TRIPLET     A
    9     6.474037    6.474037     5   SINGLET     A     0.0009  0.0150  0.0066

III модель: оптимизация модели II, нет заряда

Теперь еще раз оптимизируем систему без заряда (по итогам предыдущей итерации):

In [84]: 
%%bash
MOPAC2009.exe td2_charged.mop
babel -imopout td2_charged.out -opdb td3.pdb # td3 - оптимизация результата предыдущей оптимизации, нет заряда
In [99]: 
cmd.delete("all")
cmd.load("td3.pdb")
cmd.load("td2_charged.pdb")
cmd.do("util.cbac td2_charged") # Синий для заряженной молекулы
cmd.do("util.cbag td2_charged") # Зеленый для незаряженной молекулы, после второй оптимизации
time.sleep(2)
cmd.do(
'''
show spheres
set sphere_scale, 0.05
ray 500,500
png images/td5.png
''')
cmd.rms("td3","td2_charged") # Видим разницу, но небольшую
Out[99]:
0.002505919197574258
In [100]: 
Image(filename='images/td5.png')
Out[100]:
In [101]: 
cmd.delete("all")
cmd.load("td3.pdb")
cmd.load("td1.pdb")
cmd.do("util.cbay td1") # Желтый для незаряженной молекулы
cmd.do("util.cbag td2_charged") # Зеленый для незаряженной молекулы после второй оптимизации
time.sleep(2)
cmd.do(
'''
show spheres
set sphere_scale, 0.05
ray 500,500
png images/td6.png
''')
cmd.rms("td1","td3") # Видим разницу, но небольшую
Out[101]:
0.004314545076340437
In [102]: 
Image(filename='images/td6.png')
Out[102]:

Теперь подсчитаю энергию:

In [103]: 
%%bash
babel -imopout td2_charged.out -omop td3.mop -xk "PM6"
cp td3.mop td3_spectr.mop
printf "\ncis c.i.=4 meci oldgeo\nsome description">> td3_spectr.mop
MOPAC2009.exe td3_spectr.mop
In [104]: 
%%bash
tail -n 20 td3_spectr.out | head -n 13

Выдача:

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

    1+    0.000000    0.000000     1+  SINGLET     A
    2     4.957136    4.957136     1   TRIPLET     A
    3     5.240092    5.240092     2   TRIPLET     A
    4     5.274514    5.274514     2   SINGLET     A     0.0070  0.0541  0.0020
    5     5.549612    5.549612     3   SINGLET     A     0.0258  0.0967  0.0502
    6     6.163186    6.163186     3   TRIPLET     A
    7     6.305374    6.305374     4   SINGLET     A     0.0138  0.0301  0.0067
    8     6.393410    6.393410     4   TRIPLET     A
    9     6.474350    6.474350     5   SINGLET     A     0.0009  0.0150  0.0066

Сравнение энергий трех состояний

`STATE ENERGY I II III 1+ 0.000000 0.000000 0.000000 2 4.955172 4.957464 4.957136 3 5.241099 5.239384 5.240092 4 5.273208 5.274602 5.274514 5 5.550135 5.549203 5.549612 6 6.162894 6.164134 6.163186 7 6.305476 6.306077 6.305374 8 6.395284 6.393083 6.393410 9 6.476016 6.474037 6.474350

Третье состояние больше похоже на первое по большинству энергий. Видимо, это свидетельствует о том, что молекула не возвращается после возбуждения с исходное состояние.