Virtual Structural Biology Group » Учебные курсы » Структурная биоинформатика в ВШЭ

h1.  Вычисление параметров для молекулярной механики 

Необходимые сведения о работе с ORCA http://kodomo.fbb.msu.ru/~golovin/pdf/OrcaManual_2_9.pdf  или https://sites.google.com/site/orcainputlibrary/home 

Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов. Пример подобной работы представлен в  http://ambermd.org/antechamber/gaff.pdf это статье о разработке поля gaff.

* Вам предоставлена оптимизированная структура этана в виде z-matrix :

<pre>

inp = '''!HF RHF 6-31G

* int 0 1

C 0 0 0 0 0 0 

C 1 0 0 1.52986 0 0 

H 1 2 0 1.08439 111.200 0

H 1 2 3 1.08439 111.200 120

H 1 2 3 1.08439 111.200 -120

H 2 1 3 1.08439 111.200 180

H 2 1 5 1.08439 111.200 120

H 2 1 5 1.08439 111.200 -120

*

'''

</pre>

Как вы видите, вместо значений длин и углов связей можно вставить значение. Наша цель состоит в том, что бы рассчитать 20 разных длинн связи с шагом 0.02 ангстрема для расчёта энергии в ORCA с разными значениями по длине '''одной''' из связей (C-C).

* Конечно можно сделать эти файлы вручную, но ожидается, что вы сделаете это средствами  '''python''' . 

* Предлагаю следующую функцию по запуску ORCA 

<pre>

def run_orca(inp):

    with open('orca.inp', 'w') as outfile:

        outfile.write(inp)

    p = subprocess.Popen("/home/shad/progs/bin/orca orca.inp", 

                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)

    out=p.communicate()[0]

    out.splitlines()

    # extract energy: FINAL SINGLE POINT ENERGY'

    # and return it as float

    return ....

</pre>

* У вас есть зависимость энергии молекулы от длины одной связи. Эту зависимость можно построить в Excel (...). Вам предлагается сделать это в  matplotlib в Jupyter Notebook. 

Для matplotlib:

<pre>

%matplotlib inline

import numpy as np

import matplotlib.pyplot as plt

from scipy import optimize 

</pre>

<pre>

# fake x array , replace with real one

x_o=np.arange(1,2,0.1)

# fake y array, replace with energies

y_o=-70+2*(x_o -1.58)**2 

#function is  f(x)=k(b-x)^2 + a

fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] # Target function

errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [1,1, -79] # Initial guess for the parameters

p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))

print "Optimized params:", p1

#Plot it

plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)

plt.xlim(1,2)

plt.show()

</pre>

* Проделайте аналогичные операции для валентного угла HCС, его значения должны изменяться от 109.2 до 113.2. Сохраните полученные коэффициенты в отчёт. Примечание: не перезаписывайте файл со значениями энергий они Вам нужны для отчёта. Сравните полученные константы с данными из статьи о GAFF. Укажите возможные причины расхождений ваших результатов и публикацией.

* Проделайте аналогичные операции для торсионного угла CC, его значения должны изменяться от -180 до 180 c шагом 12. Сохраните изображение и укажите в отчёте количество минимумов функции.

* Увеличьте шаг до 0.1 ангстрема при расчёте связи. Постройте зависимость. Укажите какой функцией можно было бы аппроксимировать наблюдаемую зависимость.