ShadTask6
Version 1 (Andrey Golovin, 16.10.2017 21:34) → Version 2/3 (Andrey Golovin, 16.10.2017 21:34)
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.
*
1. Вам предоставлена оптимизированная структура этана в виде 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 ангстрема при расчёте связи. Постройте зависимость. Укажите какой функцией можно было бы аппроксимировать наблюдаемую зависимость.
Необходимые сведения о работе с 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.
*
1. Вам предоставлена оптимизированная структура этана в виде 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 ангстрема при расчёте связи. Постройте зависимость. Укажите какой функцией можно было бы аппроксимировать наблюдаемую зависимость.