Structural bionformatics. Homework №5¶
Elen Tevanyan
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
import math
import subprocess
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]
# extract energy: FINAL SINGLE POINT ENERGY'
for i in out.splitlines():
if i.startswith("FINAL SINGLE POINT ENERGY"):
return(float(i[len("FINAL SINGLE POINT ENERGY"):]))
# and return it as float
return out.splitlines()
def dependency_plot (x,y,plt_title,lim1,lim2):
x_o = x
y_o = y
#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(lim1,lim2)
plt.title(plt_title)
plt.show()
1. Длины связей
Текущая длина связи "С-С" 1.52986. Генерируем 20 возможных длин связи с шагом в 0.02 ангстрема.
length = np.arange(-10,11)*0.02 + 1.52986
print (length)
Записываем функцию для подготовки файлов с разными длинами связей
def inp_file(length):
inp = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 {} 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
*
'''.format(length)
return inp
Считаем энергию
energy = []
for i in length:
energy.append(run_orca(inp_file(i)))
print('Energy levels')
print(energy)
name = 'Bond Length and Energy Dependency'
dependency_plot(length,energy, name,1,2)
Аппроксимация не идеальная, но вполне аккуратная и улавливающая суть колебаний
2. Валентные углы
Значения валентного угла должны изменяться от 109.2 до 113.2
Cчитаем, что все еще 20 различных значений ожидаются, поэтому шаг $$|\frac{113.2-109.2}{20}|=0.2$$
angle = np.arange(-10,11)*0.2+111.2
print(angle)
Записываем функцию для генерации файлов с разными валентными углами.
def inp_file_angle(angle):
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 {0} 0
H 1 2 3 1.08439 {0} 120
H 1 2 3 1.08439 {0} -120
H 2 1 3 1.08439 {0} 180
H 2 1 5 1.08439 {0} 120
H 2 1 5 1.08439 {0} -120
*
'''.format(angle)
return inp
energy_angle = []
for i in angle:
energy_angle.append(run_orca(inp_file_angle(i)))
print('Energy levels')
print(energy_angle)
name = 'Bond Angle and Energy Dependency'
dependency_plot(angle,energy_angle, name, 108,114)
Подозрительно удачная аппроксимация
3. Торсионные углы
Торсионный угол CС изменяется от -180 до 180 c шагом 12
torsion = np.arange(-15,16)*12
print(torsion)
def inp_file_torsion(torsion):
torsion_1 = torsion - 60
torsion_2 = torsion +60
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 {0}
H 2 1 5 1.08439 111.200 {1}
H 2 1 5 1.08439 111.200 {2}
*
'''.format(torsion, torsion_1, torsion_2)
return inp
energy_torsion = []
for i in torsion:
energy_torsion.append(run_orca(inp_file_torsion(i)))
print('Energy levels')
print(energy_torsion)
name = 'Torsion angle and Energy Dependency'
dependency_plot(torsion,energy_torsion, name, -190,190)
У функции три минимума.
Аппроксимируется квадратичной функцией из рук вон плохо; визуально похожа на циклическую функцию (синус, косинус).
Попробуем такую функцию
$$ f(x) = a\cdot sin(b\cdot x+c)+d$$
x_o = torsion
y_o = energy_torsion
fitfunc = lambda p, x: p[0]*np.sin(p[1]*x+p[2]) + p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1,1, -79,1] # 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(-190,190)
plt.title(name)
plt.show()
А вот теперь аппроксимация выглядит значительно лучше
4. Длины связи с увеличинным шагом
Текущая длина связи "С-С" 1.52986. Генерируем 20 возможных длин связи с шагом в 0.1 ангстрема.
length_2 = np.arange(-10,11)*0.1 + 1.52986
print (length_2)
energy_2 = []
for i in length_2:
energy_2.append(run_orca(inp_file(i)))
print('Energy levels')
print(energy_2)
name = 'Bond Length and Energy Dependency witn Increased Step'
dependency_plot(length_2,energy_2, name,0,3)
Квадратичная аппроксимация неудачна. Из стандартных функций можно попробовать что-то из семейства $e^{-x}$: $$f(x)=a\cdot e^{-b\cdot x+c}+d$$
lim1 = 0
lim2 = 3
x_o = length_2
y_o = energy_2
fitfunc = lambda p, x: p[0]*np.exp(p[1]*x+p[2]) + p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1,1, 1,1] # 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(lim1,lim2)
plt.title(name)
plt.show()
Экспонентой аппроксимируется куда лучше.
Сравнение констант со статьей
1) Длина связи
Эмпирическое уравнение для связей С-С (Figure_2) из статьи:
$$ ln(K_r) = -4.9575\cdot r +7.8187 $$
Оптимальная длиная связи (по статье):
$$ r = 1.526 $$
Тогда
$$K_r^1 = exp(-4.9526\cdot 1.526 +7.8187) = 306.582\,Ккал/моль$$
У меня:
$$K_r^2 = 0.63991948 E_h$$
Переведем в Ккал/моль:
$$1 \;Хартри = 627,51\,ккал/моль$$
$$K_r^2=0.63991948\cdot 627,51=401.56\,Ккал/моль$$
2) Валентные углы
В статье пишут, что когда один из крайних атомов - углерод, константа получается около 50 Ккал/моль/рад^-2 (стр.1165)