Семи-эмпирические вычисления: MOPAC
Все расчеты и конвертирование файлов проводились на сервере kodomo через putty. Для начала установили переменные:
export PATH=${PATH}:/home/preps/golovin/progs/bin
export MOPAC_LICENSE=/home/preps/golovin/progs/bin
Цель работы - поэтапное освоение возможностей MOPAC как пакета для оптимизации структуры молекул и расчета некоторых свойств.
1. Оптимизация структуры порфирина
Был создан файл 1.smi , содержащий SMILE-аннотацию молекулы порфирина и название молекулы. Затем на основе этого описания с помощью программы из openbabel была построена трехмерная структура порфирина:
obgen 1.smi > 1.mol
Из молекулы были удалены лишние атомы водорода (с помощью PyMol, рисунок 1).
Структура была сохранена в файле 1.pdb .
Рисунок 1. Структура порфирина, полученная с помощью программы из openbabel (слева), та же молекула после удаления ненужных атомов водорода (справа).
Затем координаты в .mol были переформатированы во входной файл для MOPAC.
babel -ipdb myfile.pdb -omop 1_opt.mop -xk "PM6"
Запустили MOPAC, результат переформатировали в .pdb-файл .
MOPAC2009.exe 1_opt.mop babel -imopout 1_opt.out -opdb 1_opt.pdb
Сравнение полученной структуры с исходной показано на рисунке 2. Отметим, что оптимизированная структура - плоская (как ей и положено быть), в отличие от исходной.
Рисунок 2. Исходная структура порфирина (оранжевым цветом) и оптимизированная плоская структура (голубым цветом).
Далее те же самые операции провели для оптимизации с параметризацией am1 (до этого была параметризация pm6).
babel -ipdb myfile.pdb -omop 1_opt.mop -xk "AM1".pdb-файл .
На рисунке 3 производится сравнение новой оптимизированной структуры и исходной структуры порфирина.
Рисунок 3. Исходная структура порфириа (оранжевым цветом) и оптимизировавнная структура с параметризацией am6 (зеленый цвет).
Сравним оптимизированные структуры между собой (рисунок 4).
Рисунок 4. Сравнение оптимизированных структур: с параметризацией pm6 (голубым цветом) и am1 (зеленым цветом).
Видим, что структуры немного отличаются (голубые области сдвинуты ближе к периферии молекулы), но данные различия в оптимизации несущественны.
2. Расчет возбужденный состояний порфирина
Чтобы оценить спектр поглощения порфирина, надо рассчитать возбужденные состояния молекулы. Для этого сделаем копию .mop файла из предыдущего задания, и добавим в него в конце следующие строчки (файл - 1_opt_spectr.mop
cis c.i.=4 meci oldgeo For porphyrin spectrum
Для того, чтобы оценить спектр поглощения порфирина, рассчитаем возбужденные состояния молекулы. Для расчёта возбуждённых состояний в копию .mop файла из предыдущего задания 1_spec.mop добавили строку в конец файла:
Затем запуск MOPAC как в предыдущем задании. В выходном файле 1_opt_spectr.out содержатся значения энергий для электронных переходов:
STATE ENERGY (EV) Q.N. SPIN SYMMETRY POLARIZATION
ABSOLUTE RELATIVE X Y Z
1+ 0.000000 0.000000 1+ SINGLET ????
2 1.914025 1.914025 1 TRIPLET ????
3 2.266622 2.266622 2 SINGLET ????
4 2.463768 2.463768 2 TRIPLET ????
5 2.824944 2.824944 3 TRIPLET ????
6 3.363247 3.363247 4 TRIPLET ????
7 3.390559 3.390559 3 SINGLET ???? 0.2381 0.2007 0.0013
8 3.669361 3.669361 4 SINGLET ???? 2.0374 2.3896 0.0128
9 3.871992 3.871992 5 SINGLET ???? 1.8075 1.5323 0.0101
Из полученных данных об энергии мы можем посчитать длины волн переходов. Используем следующую формулу:
E = h*ν => λ = c/v = c*h/E = (2.99792458*1014мкм/с * 4.135667517*10-15 эВ*с) / Е = 1,23984193 / Е мкм
Список волн получился следующим:
338 нм; 338 нм; 366 нм; 366 нм; 369 нм; 439 нм; 503 нм; 547 нм; 648 нм
3. Определение геометрии молекулы
Для молекулы O=C1C=CC(=O)C=C1 (хинон) были повторены операции задания 1 для получения структур с использованием obgen и MOPAC. Изображения структур на рисунке 4.
Рисунок 4. Структуры хинона - полученная с помощью obgen (оранжевым цветом, слева) и оптимизированная с помощью MOPAC с параметризацией PM6 (голубая, справа)
Обе структуры плоские, по геометрии сходны. MOPAC изображает ароматику как систему сопряженных связей.
Далее нужно определить геометрию дианиона молекулы хинона. Для этого сначала в файл .mop была добавлена первая строчка CHARGE=-2
и помечены отрицательные заряды у атомов кислорода:
PM6 CHARGE=-2 quinone.pdb ... C 2.98000 1 -1.31200 1 0.02700 1 O(-) 1.05100 1 -0.05900 1 -0.48300 1 O(-) 6.19000 1 -0.03500 1 0.87800 1 H 2.40700 1 2.12500 1 -0.12400 1 ...Сравним полученную структуру с оптимизированной MOPAC (рисунок 5).
Рисунок 5. Структура хинона, оптимизированного MOPAC (слева, голубой цвет) и дианиона хинона, также оптимизированного MOPAC (справа, розовый цвет). Подписаны длины связей.
Изминились длины связей в ароматическом кольце между атомами углерода (в незаряженной молекуле были 1.5 и 1.3, в дианионе стали все 1.4), а также длины связей C=O (в незаряженной молекуле были 1.2, в дианионе стали 1.3).
4. Тиминовые димеры
Известно, что ультрафиолет может превращать тимины в тиминовые димеры, также известно, что ДНК фотолиаза при облучении ультрафиолетом
востановливает основания тиминов до нормального. Нам дана структура тимидинового димера (рисунок 6).
Рисунок 6. Структура тимидинового димера
Цель работы - увидеть переход из димера в тимины при возбуждении системы. Чтобы облегчить вычисления MOPAC, ионизируем оба кольца, т.е. укажем заряд системы +2.
Сначала проведем оптимизацию геометрии этой структуры при заряде системы 0:
babel -ipdb td.pdb -omop td1.mop -xk "PM6" MOPAC2009.exe td1.mop babel -imopout td1.out -opdb td1.pdb
Рисунок 7. Структура оптимизированного димера с зарядом системы 0 (слева, голубой цвет) и сравнение оптимизированной (голубой цвет) и неоптимизированной
структуры (оранжевый цвет) с зарядом системы 0 (справа).
Затем присвоим системе из предыдущего пункта заряд +2 и проведем оптимизацию:
babel -ipdb td1.pdb -omop td2.mop -xk "PM6" sed -i 's/PM6/PM6 CHARGE=+2/' td2.mop MOPAC2009.exe td2.mop babel -imopout td2.out -opdb td2.pdb
Рисунок 8. Оптимизированная структура димера при заряде +2.
Следом оптимизируем результат из предыдущего пункта при заряде системы 0:
babel -ipdb td2.pdb -omop td3.mop -xk "PM6" MOPAC2009.exe td3.mop babel -imopout td3.out -opdb td3.pdb
Рисунок 9. Оптимизированная структура при заряде системы 0. Произошел переход димера в тимины.
Сравним энергии всех состояний:
td1.out -3273.58217 EV
td2.out -3253.90834 EV
td3.out -3273.69661 EV
Энергии конечного и исходного состояния практически равны. При добавлении электронов к возбужденному состоянию
происходит переход димера в тимины, при этом происходит образование двух ароматических структур.
Состояние двух димеров энергетически должно быть более выгодным состояния димеров (еще и с напряженной пространственной структурой).
Однако разница энергий этих состояний невелика. Этим можно объяснить возможность образования тиминовых димеров при облучении УФ.