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

#Подгружаю все утелиты:
import os
os.environ['PATH']=os.environ['PATH']+'/home/preps/golovin/progs/bin'
os.environ['MOPAC_LICENSE'] ='/home/preps/golovin/progs/bin'
os.environ['MOPAC_LICENSE']='/home/preps/golovin/progs/bin'

%%bash
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/home/preps/golovin/progs/lib

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
import __main__
__main__.pymol_argv = [ 'pymol', '-cp' ]
import pymol
pymol.finish_launching()
from pymol import cmd

%%bash
echo 'c1cc2cc3ccc(cc4ccc(cc5ccc(cc1n2)[nH]5)n4)[nH]3' >  porphyrin.smi
#Найденный для порфирина SMILES

Теперь из, грубо говоря, текстового файла хочется получить PDB

%%bash
obgen porphyrin.smi > porphyrin.mol

cmd.delete('all')
cmd.load('porphyrin.mol')

cmd.refresh()
cmd.do('''
bg_color black
as sticks, all
orient, all
set label_outline_color, white
set label_size, -0.2
label all, " %s" % ID
remove id 39+40
rotate x, 90
ray
png pic1.png
''')

Image(filename='pic1.png')

cmd.save('porphyrin.pdb')

Вот и получен PDB-файл со структурой порфирина

%%bash
babel -ipdb porphyrin.pdb -omop porphyrin_opt.mop -xk "PM6"

%%bash
#/home/preps/golovin/progs/bin/MOPAC2009.exe porphyrin_opt.mop
#На самом деле это так не работает, нужно было зайти в пути, вбить пути на вашу страницу и использовать MOPAC

%%bash 
babel -imopout porphyrin_opt.out -opdb porphyrin_opt.pdb

cmd.delete('all')
cmd.load('porphyrin_opt.pdb')

cmd.refresh()
cmd.do('''
bg_color black
as sticks, all
orient, all
set label_outline_color, white
set label_size, -0.2
label all, " %s" % ID
remove id 39+40
rotate x, 90
ray
png pic2.png
''')

Image(filename='pic2.png')

После удаления водородов, молекула оптимизировалась и стала плоской. УРА!

Рассчитайте возбужденные состояния порфирина и на основе этих данных прикиньте спектр поглощения молекулы. Для расчёта возбуждённых состояний сделайте копию mop файла из предыдущего задания

%%bash
cp porphyrin_opt.mop porphyrin_opt_spectr.mop
echo "
cis c.i.=4 meci oldgeo
For porphyrin spectrum" >> porphyrin_opt_spectr.mop

#Снова пришлось через Putty сделать расчет
#/home/preps/golovin/progs/bin/MOPAC2009.exe 1_opt_spectr.mop

Получили файл и нужно рассчитать энергию

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

Получается, порфирин поглащает в УФ и в видимом свете.

Для молекулы O=C1C=CC(=O)C=C1 определите геометрию как с помощью obgen так и Мopac (основные шаги см. выше). Определите геометрию дианиона этой молекулы. Для начала в первую строчку mop файла добавьте слово CHARGE=-2. Потом явным способом укажите на каких атомов по вашему мнению должен находиться отрицательный заряд

%%bash 
echo 'O=C1C=CC(=O)C=C1' >  mol.smi

%%bash 
obgen mol.smi > mol.mol

cmd.delete('all')
cmd.load('mol.mol')

cmd.refresh()
cmd.do('''
bg_color black
as sticks, all
orient, all
set valence, on 
set valence_mode, 3
set stick_ball, on
set stick_ball_ratio, 1
set label_outline_color, white
set label_size, -0.2
label all, " %s" % ID
remove id 39+40
rotate x, 90
ray
png pic3.png
''')

Image(filename='pic3.png')

Молекула плоская! Зделаес, как и в прошлый раз, PDB-файл

cmd.save('molecule.pdb')

%%bash
babel -ipdb molecule.pdb -omop molecule.mop -xk "PM6"

#/home/preps/golovin/progs/bin/MOPAC2009.exe molecule.mop

%%bash
babel -imopout molecule.out -opdb molecule_2.pdb

cmd.delete('all')
cmd.load('molecule_2.pdb')
cmd.refresh()
cmd.do('''
bg_color black
orient, all
set label_outline_color, white
set label_size, -0.5
label all, " %s" % ID
as sticks, all
set valence, on
set stick_ball, on
set stick_ball_ratio, 3
set stick_radius, 0.12
set valence_mode, 3
dist id 1, id 7
ray
png pic4.png

''')

Image(filename='pic4.png')

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

#Теперь поменяем заряд и глянем, что изменится

with open("molecule.mop") as f,open("molecule_newcharge.mop",'w') as o:
    data=f.read()
    data=data.replace("PM6\nmol.pdb","PM6 CHARGE=-2\nmol.pdb")
    o.write(data)
with open("molecule_newcharge.mop") as f,open("molecule_newcharge2.mop",'w') as o:
    data=f.read()
    data=data.replace("O","O(-)")
    o.write(data) 

#/home/preps/golovin/progs/bin/MOPAC2009.exe molecule_newcharge2.mop

%%bash
babel -imopout molecule_newcharge2.out -opdb molecule_newcharge2.pdb

cmd.delete('all')
cmd.load('molecule_newcharge2.pdb')
cmd.refresh()
cmd.do('''
bg_color black
orient, all
set label_outline_color, white
set label_size, -0.5
label all, " %s" % ID
as sticks, all
set valence, on
set stick_ball, on
set stick_ball_ratio, 3
set stick_radius, 0.12
set valence_mode, 3
dist id 1, id 7
ray
png pic5.png

''')

Image(filename='pic5.png')

Ничего не изменилось, хотя должно было. Видимо, где-то произошло ошибка при изменении заряда.

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

cmd.delete('all')
cmd.load('td.pdb')
cmd.refresh()
cmd.do('''
bg_color black
orient, all
set label_outline_color, white
set label_size, -0.5
label all, " %s" % ID
as sticks, all
set valence, on
set stick_ball, on
set stick_ball_ratio, 3
set stick_radius, 0.12
set valence_mode, 3
ray
png pic6.png

''')

Image(filename='pic6.png')

1) Проведите оптимизацию геометрии этого димера, при заряде системы 0. 2) Проведите оптимизацию результата из пункта 1, при заряде системы +2. 3) Проведите оптимизацию результата из пункта 2, при заряде системы 0.

%%bash 
babel -ipdb td.pdb -omop td.mop -xk "PM6"

#/home/preps/golovin/progs/bin/MOPAC2009.exe td.mop

%%bash 
babel -imopout td.out -opdb td_opt.pdb

cmd.delete('all')
cmd.load('td_opt.pdb')
cmd.refresh()
cmd.do('''
bg_color black
orient, all
set label_outline_color, white
set label_size, -0.5
label all, " %s" % ID
as sticks, all
set valence, on
set stick_ball, on
set stick_ball_ratio, 3
set stick_radius, 0.12
set valence_mode, 3
ray
png pic7.png

''')

Image(filename='pic7.png')

После выполнения этой оптимизации, геометрия немного поменялась

#Теперь сделаем заряд +2

%%bash
babel -ipdb td_opt.pdb -omop td_opt.mop -xk "PM6"

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)

#/home/preps/golovin/progs/bin/MOPAC2009.exe td2.mop

%%bash 
babel -imopout td2.out -opdb td2_opt.pdb

cmd.delete('all')
cmd.load('td2_opt.pdb')
cmd.refresh()
cmd.do('''
bg_color black
orient, all
set label_outline_color, white
set label_size, -0.5
label all, " %s" % ID
as sticks, all
set valence, on
set stick_ball, on
set stick_ball_ratio, 3
set stick_radius, 0.12
set valence_mode, 3
ray
png pic8.png

''')

Image(filename='pic8.png')

Внесения заряда +2 сильно изменило геометрию молекулы.

#Теперь к этой молекуле сделаем оптимизацию с зарядом +0

%%bash
babel -ipdb td2_opt.pdb -omop td2_opt.mop -xk "PM6"

with open("td2_opt.mop") as f,open("td3.mop",'w') as o:
    data=f.read()
    data=data.replace("PM6\ntd2_opt.pdb","PM6 CHARGE=+0\ntd_opt.pdb")
    o.write(data)

#/home/preps/golovin/progs/bin/MOPAC2009.exe td3.mop

%%bash 
babel -imopout td3.out -opdb td3_opt.pdb

cmd.delete('all')
cmd.load('td3_opt.pdb')
cmd.refresh()
cmd.do('''
bg_color black
orient, all
set label_outline_color, white
set label_size, -0.5
label all, " %s" % ID
as sticks, all
set valence, on
set stick_ball, on
set stick_ball_ratio, 3
set stick_radius, 0.12
set valence_mode, 3
ray
png pic9.png

''')

Image(filename='pic9.png')

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

Из файлов OUT сравнил энергии и понял, что энергия, как димера, так и отдельных нуклеотидов приблизительно одинаковая. Но энергию отдельных нуклетидов всё же совсем немного меньше, из-за чего при сдвиге заряда димер стал отдельными нуклеотидами