Учебная страница курса биоинформатики,
год поступления 2013

Семестры Студенты Преподаватели

• 2024
• 2023
• 2022
• 2021
• 2020
• 2019
• 2018
• 2017

• ...

Занятие 4

Часть 1. Cеми-эмпирические вычисления: Mopac

Сегодня ключевым модулем на сегодня будет openbabel https://openbabel.org/docs/dev/index.html, это небольшой фрейворк для конвертации молекул между разными форматами.

Итак, как обычно создаем рабочую директорию и в ней запускаем jypiter notebook как мы это делали раньше .

• Первым шагом проверим в какой мы директории.
• Загрузка модуля

   1 import pybel

посмотрим какие форматы он может читать и писать:

   1 print pybel.informats
   2 print pybel.outformats

• Создадим соединение и сразу сделаем его 3D:

   1 mol=pybel.readstring('smi','CC(=O)C')
   2 mol.addh()
   3 mol.make3D()

• Посмотрим как оно может записывать формат программы сэмиэмпирических расчётов MOPAC, и вывод покажем на экран:

   1 mop=mol.write(format='mopin',opt={'k':'PM6 CHARGE=%d' % mol.charge})

Заряд правильный? можно и сохранить как файл:

   1 mop=mol.write(format='mopin',filename='my.mop',opt={'k':'PM6 CHARGE=%d' % mol.charge})

• Осталось запустить MOPAC, как это сделать можно узнать http://stackoverflow.com/questions/2231227/python-subprocess-popen-with-a-modified-environment тут или внутри notebook, обратите внимание на magic %%bash:

   1 %%bash
   2 export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/'
   3 /home/preps/golovin/progs/mopac/MOPAC2016.exe my.mop

Подсказака:

   1 import subprocess, os
   2 my_env = os.environ.copy()
   3 my_env["MOPAC_LICENSE"] = '/home/preps/golovin/progs/mopac/'
   4 p = subprocess.Popen("/home/preps/golovin/progs/mopac/MOPAC2016.exe %s" % 'my.mop', 
   5                      shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=my_env)
   6 p.communicate()

• Обратите внимание, что просмотреть файл вы можете в FAR или в терминале запущенном в jypiter notebook (new->terminal)

• Считать геометрию можно так:

   1 opt=pybel.readfile('mopout','my.out')
   2 for i in opt:
   3     print i
   4     i.write(format='pdb',filename='my.pdb')

Обратите внимание, что объект чтения это итератор, так как геометрий там может быть много.

• Сравните структуры порфирина полученные pybel и двумя способами PM6 и AM1 в Mopac, выводы, картинки и наблюдения занесите в отчёт.
• Рассчитайте возбужденные состояния порфирина и на основе этих данных прикиньте спектр поглощения молекулы. Для расчёта возбуждённых состояний сделайте копию mop файла из предыдущего задания. Например opt_spectr.mop. Для указания Mopac о необходимости расчёта возбуждённого состояния добавьте в конец файла:

cis c.i.=4 meci oldgeo
some description

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

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

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

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

Сравните энерегии всех состояний. Опишите причину по которой по Вашему мнению возбуждённое состояние при добавлении электронов не перешло обратно в исходно состояние.

Часть 2. Ab inito вычисления: ORCA

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

• Постройте и оптимизируйте с помощью MOPAC структуры нафталена и азулена
• Создадим входные файлы для Gamess.

format='orcainp'
header='''!HF RHF OPT 6-31G ''' 

• Вероятно pybel не изменит заголовок тогда надо:

mop=mol.write(format='orcainp',opt={'k':''})
out = mop.replace('! insert inline commands here ', 'My commands')
save out to file orca.inp

• Проведите оптимизацию геометрии для обоих молекул. Для этого запустите ORCA следующим образом:

   1 /srv/databases/orca/orca/orca orca.inp | tee orca-opt.log 

• На основе полученных координат составьте новые входные файлы для расчёта энергии. Теперь Вам надо будет построить по два файла на каждую молекулу, в первом случае Вы будете вести расчёт методом Хартри-Фока, а во втором используя теорию функционала плотности. Для расчёта по Хартри-Фоку надо составить файл с таким заголовком:

   1  !HF RHF 6-31G

• В случае теории функционала плотности заголовок будет таким:

   1  !DFT RHF 6-31G

• Рассчитайте четыре системы: два способа на каждую молекулу.
  • Найдите в log файлах расчёта энергии строчку с "TOTAL ENERGY = " и выпишите значения этой энергии в таблицу:

• Определите какой метод лучше если из эксперимента известно, что энергия изомеризации нафталина в азулен составляет 35.3±2.2 kCal/mol

* Допонительно: ORCA может сохранять cube файлы с плотностью примерный синтаксис в input:

   1 %plots Format Cube
   2  MO("xxx.cube",n,0);
   3 end

где n орбитали. Узнать номер можно если найти в output:

----------------
ORBITAL ENERGIES
----------------

Тогда HOMO будет орбиталь где макс энергия при OCC =2 , а LUMO следующая