In [1]: 
from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys
import scipy as sp
from scipy import constants
from scipy.constants import codata
import numpy as np
In [2]: 
import __main__

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


import pymol
 
pymol.finish_launching()
 
from pymol import cmd
In [3]: 
from IPython.display import Image
In [4]: 
os.environ['PATH']=os.environ['PATH']+'/home/preps/golovin/progs/bin'
os.environ['MOPAC_LICENSE'] ='/home/preps/golovin/progs/bin'

Исправление порфирина

In [6]: 
%%bash 
echo 'c1cc2cc3ccc(cc4ccc(cc5ccc(cc1n2)[nH]5)n4)[nH]3' >  porphyrin.smi
In [7]: 
%%bash 
obgen porphyrin.smi > porphyrin.mol
In [8]: 
cmd.delete('all')
cmd.load('porphyrin.mol')
In [9]: 
cmd.refresh()
cmd.do('''
bg_color white
as sticks, all
orient, all
set label_outline_color, black
set label_size, -0.5
label all, " %s" % ID
remove id 39+40
rotate x, 90
ray
png pic1.png
''')
In [10]: 
Image(filename='pic1.png')
Out[10]:
In [11]: 
cmd.save('porphyrin.pdb')
In [12]: 
%%bash 
babel -ipdb porphyrin.pdb -omop porphyrin_opt.mop -xk "PM6"
In [13]: 
%%bash 
/home/preps/golovin/progs/bin/MOPAC2009.exe porphyrin_opt.mop
In [14]: 
%%bash 
babel -imopout porphyrin_opt.out -opdb porphyrin_opt.pdb
In [15]: 
cmd.delete('all')
cmd.load('porphyrin_opt.pdb')
In [16]: 
cmd.refresh()
cmd.do('''
bg_color white
as sticks, all
orient, all
set label_outline_color, black
set label_size, -0.5
label all, " %s" % ID
ray
png pic2.png
''')
In [17]: 
Image(filename='pic2.png')
Out[17]:

После удаления лишних водородов молекула порфирина стала плоской.

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

In [18]: 
%%bash
cp porphyrin_opt.mop porphyrin_opt_spectr.mop
echo "
cis c.i.=4 meci oldgeo
For porphyrin spectrum" >> porphyrin_opt_spectr.mop
In [19]: 
%%bash
/home/preps/golovin/progs/bin/MOPAC2009.exe porphyrin_opt_spectr.mop
In [26]: 
fl=open('porphyrin_opt_spectr.out','r')
lines=fl.readlines()
for l in lines[-30:-1]:
    print l.split()

         8         0   1   1   0
      2.4974       1   1   0   0

         9         0   1   0   1
      3.4825       1   1   0   0

DIAGONAL OF C.I. MATRIX 0.000000 3.241984 2.658889 2.497420 3.482493 3.241984 2.658889 2.497420 3.482493

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

1+    0.000000    0.000000     1+  SINGLET     ????
2     1.913804    1.913804     1   TRIPLET     ????
3     2.266288    2.266288     2   SINGLET     ????
4     2.463311    2.463311     2   TRIPLET     ????
5     2.824593    2.824593     3   TRIPLET     ????
6     3.362738    3.362738     4   TRIPLET     ????
7     3.390622    3.390622     3   SINGLET     ????  0.2634  0.1676  0.0065
8     3.669215    3.669215     4   SINGLET     ????  0.0155  0.0639  4.3621
9     3.871001    3.871001     5   SINGLET     ????  2.0305  1.2716  0.0493

The "+" symbol indicates the root used.

TOTAL CPU TIME: 14.20 SECONDS

Мы получили энергии электронных переходов. Посчитаем длины волн, при которых они происходят. \[ \lambda=(h*c)/E \]

In [27]: 
for i in lines[-15:-7]:
     l=i.split()
     k=float(l[2])
     print l[0], l[4], (constants.c*constants.h*(10**9)/codata.value('electron volt-joule relationship'))/k, 'nm'
fl.close()

2 TRIPLET 647.841643763 nm
3 SINGLET 547.080481033 nm
4 TRIPLET 503.323343744 nm
5 TRIPLET 438.945338036 nm
6 TRIPLET 368.700127456 nm
7 SINGLET 365.667989295 nm
8 SINGLET 337.903864778 nm
9 SINGLET 320.289746554 nm

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

Определить геометрию молекулы O=C1C=CC(=O)C=C1 и её дианиона с помощью obgen и Мopac

In [28]: 
%%bash 
echo 'O=C1C=CC(=O)C=C1' >  mol.smi
In [29]: 
%%bash 
obgen mol.smi > mol.mol
In [29]: 
cmd.delete('all')
cmd.load('mol.mol')
In [30]: 
cmd.refresh()
cmd.do('''
bg_color white
orient, all
set label_outline_color, black
set label_size, -0.5
label all, " %s" % ID
show sticks
set valence, on
set stick_ball, on
set stick_ball_ratio, 3
set stick_radius, 0.12
set valence_mode, 2
ray
png pic3.png
rotate x, 90
ray
png pic4.png
''')
In [31]: 
Image(filename='pic3.png')
Out[31]:
In [32]: 
Image(filename='pic4.png')
Out[32]:
In [33]: 
cmd.save('mol.pdb')

С помощью obgen определили, что молекула плоская. Похоже, что это хинон. Теперь попробуем с помощью MOPAC.

In [49]: 
%%bash 
babel -ipdb mol.pdb -omop mol.mop -xk "PM6"
In [50]: 
%%bash 
/home/preps/golovin/progs/bin/MOPAC2009.exe mol.mop
babel -imopout mol.out -opdb mol_new.pdb
In [108]: 
cmd.delete('all')
cmd.load('mol_new.pdb')
cmd.refresh()
cmd.do('''
bg_color white
orient, all
set label_outline_color, black
set label_size, -0.4
label all, " %s" % ID
show sticks
set valence, on
set stick_ball, on
set stick_ball_ratio, 3
set stick_radius, 0.12
set valence_mode, 2
dist id 1, id 7
dist id 4, id 8
ray
png pic5.png
rotate x, 90
ray
png pic6.png
''')
In [109]: 
Image(filename='pic5.png')
Out[109]:
In [110]: 
Image(filename='pic6.png')
Out[110]:

MOPAC тоже сделал плоскую молекулу, только ещё показывает делокализованный заряд, в отличие от первого случая, хотя в насройках pymol стоит set valence_mode, 2, и он вроде бы тоже должен его показывать.

In [56]: 
with open("mol.mop") as f,open("mol_charge.mop",'w') as o:
    data=f.read()
    data=data.replace("PM6\nmol.pdb","PM6 CHARGE=-2\nmol.pdb")
    o.write(data)
with open("mol_charge.mop") as f,open("mol_charge2.mop",'w') as o:
    data=f.read()
    data=data.replace("O","O(-)")
    o.write(data)
In [65]: 
%%bash 
/home/preps/golovin/progs/bin/MOPAC2009.exe mol_charge2.mop
babel -imopout mol_charge2.out -opdb mol_charge2.pdb
In [105]: 
cmd.delete('all')
cmd.load('mol_charge2.pdb')
cmd.refresh()
cmd.do('''
bg_color white
orient, all
set label_outline_color, black
set label_size, -0.4
label all, " %s" % ID
show sticks
set valence, on
set stick_ball, on
set stick_ball_ratio, 3
set stick_radius, 0.12
set valence_mode, 2
dist id 1, id 7
dist id 4, id 8
ray
png pic7.png
rotate x, 90
ray
png pic8.png
''')
In [106]: 
Image(filename='pic7.png')
Out[106]:
In [107]: 
Image(filename='pic8.png')
Out[107]:

Длина связи увеличилась на 0.1 Å в молекуле с указанным зарядом.

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

In [44]: 
cmd.delete('all')
cmd.load('td.pdb')
cmd.refresh()
cmd.do('''
bg_color white
orient, all
show sticks
set valence, on
set stick_ball, on
set stick_ball_ratio, 3
set stick_radius, 0.12
set valence_mode, 2
set ray_trace_mode, 3
ray
png pic9.png
''')
In [45]: 
Image(filename='pic9.png')
Out[45]:

Проведите оптимизацию геометрии этого димера, при заряде системы 0.

In [5]: 
%%bash 
babel -ipdb td.pdb -omop td.mop -xk "PM6"
/home/preps/golovin/progs/bin/MOPAC2009.exe td.mop
In [6]: 
%%bash 
babel -imopout td.out -opdb td_opt.pdb
In [46]: 
cmd.delete('all')
cmd.load('td_opt.pdb')
cmd.refresh()
cmd.do('''
bg_color white
orient, all
show sticks
set valence, on
set stick_ball, on
set stick_ball_ratio, 3
set stick_radius, 0.12
set valence_mode, 2
set ray_trace_mode, 3
ray
png pic10.png
''')
In [47]: 
Image(filename='pic10.png')
Out[47]:

При сравнении с прошлой картинкой, можно заметить, что геометрия димера немного изменилась

Проведите оптимизацию результата при заряде системы +2.

In [10]: 
%%bash 
babel -ipdb td_opt.pdb -omop td_opt.mop -xk "PM6"
In [11]: 
with open("td_opt.mop") as f,open("td2.mop",'w') as o:
    data=f.read()
    data=data.replace("PM6\ntd_opt.pdb","PM6 CHARGE=+2\ntd_opt.pdb")
    o.write(data)
In [6]: 
%%bash 
/home/preps/golovin/progs/bin/MOPAC2009.exe td2.mop

Не работает MOPAC =( Хотя сначала его удалось запустить
Пишет: "Your MOPAC executable has expired. Please go to web-site: http://openmopac.net/Download_MOPAC_Executable_Step2.html to get a new version of MOPAC."

In []: 
%%bash 
babel -imopout td2.out -opdb td2_opt.pdb