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

Суть задания состоит в поэтапном освоении возможностей Mopac, как пакета для оптимизации структуры молекул и расчёта некоторых свойств.

Вся работа по расчётам и конвертированию файлов будет проходить на сервере kodomo через putty. Начнём работу с установки переменных:

Оптимизация структуры порфирина

In [1]: 
%%bash
export PATH=${PATH}:/home/preps/golovin/progs/bin
export MOPAC_LICENSE=/home/preps/golovin/progs/bin
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/home/preps/golovin/progs/lib
In [2]: 
from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys
import numpy as np

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


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

import pymol
from pymol import cmd
pymol.finish_launching()

#
from IPython.display import Image

from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys, time

# images
defaultImage = 'pymolimg.png'
def prepareImage(width=300, height=300, sleep=5, filename=defaultImage):
    ## To save the rendered image
    cmd.ray(width, height)
    cmd.png('pymolimg.png')
    time.sleep(sleep)
    
# Define some shortcuts
def focus(x):
    cmd.center(x)
    cmd.zoom(x)
In [3]: 
os.environ['PATH']=os.environ['PATH']+'/home/preps/golovin/progs/bin'
os.environ['MOPAC_LICENSE'] ='/home/preps/golovin/progs/bin'

Создадим SMILES порфирина и построим его структуру с помощью babel.

In [5]: 
%%bash 
echo 'c1cc2cc3ccc(cc4ccc(cc5ccc(cc1n2)[nH]5)n4)[nH]3' >  porphyrin.smi
In [6]: 
%%bash 
obgen porphyrin.smi > porphyrin.mol
In [8]: 
cmd.delete('all')
cmd.load('porphyrin.mol')

prepareImage()
Image(defaultImage)
Out[8]:
In [9]: 
cmd.save('porphyrin.pdb')

Удалили возможные ненужные водороды из pdb.

In [12]: 
cmd.load('porphyrin-h.pdb')

prepareImage()
Image(defaultImage)
Out[12]:

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

In [29]: 
%%bash 
babel -ipdb porphyrin-h.pdb -omop porphyrin_opt.mop -xk "PM6"

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

In [39]: 
%%bash 
/home/preps/golovin/progs/bin/MOPAC2009.exe porphyrin_opt.mop
#выдачу porphyrin_opt.out переконверитуем в pdb далее
In [41]: 
%%bash 
babel -imopout porphyrin_opt.out -opdb porphyrin_opt.pdb

Было

In [38]: 
cmd.delete('all')
cmd.load('porphyrin-h.pdb')
cmd.refresh()

prepareImage()
Image(defaultImage)
Out[38]:

Стало

In [42]: 
cmd.delete('all')
cmd.load('porphyrin_opt.pdb')
cmd.refresh()

prepareImage()
Image(defaultImage)
Out[42]:

Полученная структура имеет плоскую геометрию.

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

Рассчитаем возбуждённые состояния для порфирина (файл porphyrin_1opt.mop).

Чтобы указать программе Mopac необходимость расчёта возбуждённого состояния, добавим в конец файла соответствуюшие строки:

In [3]: 
%%bash
cp porphyrin_opt.mop 1_opt_spectr.mop
echo "
cis c.i.=4 meci oldgeo
For porphyrin spectrum" >> 1_opt_spectr.mop
In [16]: 
%%bash
/home/preps/golovin/progs/bin/MOPAC2009.exe 1_opt_spectr.mop
In [15]: 
%%bash
MOPAC2009.exe 1_opt_spectr.mop #почему-то будто работает и это, и выше, но файл out не появляется
In [4]: 
fl=open('1_opt_spectr.out','r')
lines=fl.readlines()
for l in lines[-30:-1]:
    print l.split()

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

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

Запишем эту часть в файл:

In [6]: 
%%bash
for i in 'SINGLET' 'DOUBLET' 'TRIPLET'
    do
    grep ${i} '1_opt_spectr.out' >> states.txt
done

Файл states.txt (по каким-то причинам ipython notebook перестал показывать output).

Длины волн, при которых происходят данные электронные переходы, рассчитывается по формуле: длина волны (нм) = h*c/E, где h - постоянная Планка, c - скорость света, E - энергия перехода. Они подсчитаны далее:

In [8]: 
%%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

Файл porphyrin_spectr.miM (по каким-то причинам ipython notebook перестал показывать output). В нем находятся значения в мкМ:

.64743737725189412631

.54687199182057488867

.50310194997149401780

.43869280123019482920

.36848644807482892916

.36548635780900726260

.33791942919222053907

.32024631372745527215

Таким образом, порфирин может поглощать в УФ и видимом спектре.

Определение геометрии молекулы

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

In [9]: 
%%bash 
echo 'O=C1C=CC(=O)C=C1' >  mol.smi
In [10]: 
%%bash 
obgen mol.smi > mol.mol
In [12]: 
cmd.delete('all')
cmd.load('mol.mol')
cmd.refresh()
cmd.do('''
set valence, on
rotate y, -90
''')

prepareImage()
Image(defaultImage)
Out[12]:

Молекула (это хинон) плоская, судя по выдаче obgen. Посмотрим на оптимизацию Mopac.

In [18]: 
cmd.save('quinone_1opt.pdb')
In [21]: 
%%bash
babel -ipdb quinone_1opt.pdb -omop quinone_1opt.mop -xk "PM6"

START HERE

In [22]: 
%%bash
/home/preps/golovin/progs/bin/MOPAC2009.exe quinone_1opt.mop
In [24]: 
%%bash
babel -imopout quinone_1opt.out -opdb quinone_mop.pdb
In [30]: 
cmd.delete('all')
cmd.load('quinone_mop.pdb')
cmd.refresh()
cmd.do('''
set valence, on
rotate y, -30
''')

prepareImage()
Image(defaultImage)
Out[30]:

По виду полученные геометрии молекулы сходны: они плоские.

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

In [31]: 
%%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
In []: 
%%bash
/home/preps/golovin/progs/bin/MOPAC2009.exe quinone_anion.mop
In [32]: 
%%bash
babel -imopout quinone_anion.out -opdb quinone_anion.pdb
In [33]: 
cmd.delete('all')
cmd.load('quinone_anion.pdb')
cmd.refresh()
cmd.do('''
set valence, on
rotate y, -30
''')

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

Геометрия структуры дианиона также плоская.

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

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

In [34]: 
cmd.delete('all')
cmd.load('td.pdb')
cmd.refresh()

prepareImage()
Image(defaultImage)
Out[34]:

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

In [35]: 
%%bash 
babel -ipdb td.pdb -omop td.mop -xk "PM6"
In []: 
%%bash 
/home/preps/golovin/progs/bin/MOPAC2009.exe td.mop
In [36]: 
%%bash 
babel -imopout td.out -opdb td_opt.pdb
In [37]: 
cmd.delete('all')
cmd.load('td_opt.pdb')
cmd.refresh()

prepareImage()
Image(defaultImage)
Out[37]:

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

In [38]: 
%%bash 
babel -ipdb td_opt.pdb -omop td_opt.mop -xk "PM6"
In [39]: 
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 []: 
%%bash 
/home/preps/golovin/progs/bin/MOPAC2009.exe td2.mop
In [40]: 
%%bash 
babel -imopout td2.out -opdb td2_opt.pdb
In [42]: 
cmd.delete('all')
cmd.load('td2_opt.pdb')
cmd.refresh()

prepareImage()
Image(defaultImage)
Out[42]:

Проведем оптимизацию полученного результата при заряде системы 0.

In [43]: 
%%bash 
babel -ipdb td2_opt.pdb -omop td2_opt.mop -xk "PM6"
In [45]: 
with open("td2_opt.mop") as f,open("td3.mop",'w') as o:
    data=f.read()
    data=data.replace("PM6\ntd_opt2.pdb","PM6 CHARGE=+2\ntd_opt2.pdb")
    o.write(data)
In []: 
%%bash 
/home/preps/golovin/progs/bin/MOPAC2009.exe td3.mop
In [46]: 
%%bash 
babel -imopout td3.out -opdb td3_opt.pdb
In [48]: 
cmd.delete('all')
cmd.load('td3_opt.pdb')
cmd.refresh()

prepareImage()
Image(defaultImage)
Out[48]:

Сначала при возбуждении системы (заряд +2) произошел переход димера в тимины. Однако после оптимизации полученной системы при заряде 0 переход не произошел. Почему возбуждённое состояние при добавлении электронов не перешло обратно в исходно состояние? Попроубем сравнить энерегии всех состояний.

Значения энергий:

td.out:

TOTAL ENERGY = -3273.58217 EV

ELECTRONIC ENERGY = -22320.64805 EV

td2.out:

TOTAL ENERGY = -3253.90834 EV

ELECTRONIC ENERGY = -21567.92464 EV

td3.out:

TOTAL ENERGY = -3273.69661 EV

ELECTRONIC ENERGY = -20274.13483 EV

Полная энергия вернулась к исходному состянию, а энергия электронов повысилась в третьем состоянии. Не очень ясно, как связаны энергии и невозвращение тиминов в состояние димеров.