inp_init = '''!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 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
inp=inp_init
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]
out=out.split('Total Energy')[1]
out=out.splitlines()
out=out[0].split(" ")
matches = [x for x in out if (x!="" and x!=":")]
return(float(matches[0]))
def exchange(inp,n,pos,step):
a=inp.splitlines()[pos]
b=a.split(" ")
b[n]=str(float(b[n])+step)
a=" ".join(b)
c=inp.splitlines()
c[pos]=a
inp="\n".join(c) + '\n'
return (inp)
def Energy(inp,n,pos,step,init,N):
E=[]
inp=exchange(inp,n,pos,-init)
E.append(run_orca(inp))
for i in range(1,N):
inp=exchange(inp,n,pos,step)
E.append(run_orca(inp))
return(E)
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
def x_generation(init,step,r,N):
x=[init-r]
for i in range (N):
x.append(x[-1]+step)
return(x)
Посчитаем энергию для различных длин связи и построим график
x_o=x_generation(1.52986,0.02,0.2,20)
y_o=Energy(inp,4,3,0.02,0.2,21)
#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
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()
Энергия у нас в Eh - энергия Хартри.
Константа ковалентного взаимодействия:
Kr=0.61699 Eh/Aˆ2
Переведем в ккал/моль/Аˆ2.
Kr=387.161225 ккал/моль/Аˆ2
Исходя из статьи по разработке поля gaff, можно посчитать эту константу через молекулярную механику. Формула, полученная МНК на экспериментальных данных:
lnKr=-4.9575*ln(r)+7.8187, для атомов С-С, где r - равновесная длина связи
Kr=311 ккал/моль/Аˆ2
Т.о. мы видим, что константы, посчитанные разными методами, в общем хоть как-то сопоставимы. Возможно, различие связано с тем, что в основе Orca квантовая химия, а в gaff - молекулярная динамика, которая лучше описывает большие молекулы (белки, днк). Возможно, в данном случае не стоит пренебрегать квантовым эффектом, приближения МД не очень хорошо подходят.
Посчитаем энергию для разных валентных углов.
E_val_angle=Energy(inp_init,5,4,0.2,2.,21)
x_o=x_generation(111.2,0.2,2.0,20)
y_o=E_val_angle
#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(109,114)
plt.show()
Здесь аппроксимация лучше. Проведем те же манипуляции, что и для константы взаимодействия Kr. Получили вот что:
$\begin{equation} \\K^\theta = 83.72 \; ккал/моль/рад^2 \end{equation}$
Уравнение для константы из молекулярной механики:
$\begin{equation}
\\K^\theta_{ijk} = 143.9 Z_iC_jZ_k\left(r_{ij}^{eq}+r_{jk}^{eq}\right)^{-1}\cdot\left(\theta_{ijk}^{eq}\right)^{-2}exp\left({-2}\cdot\left(\frac{r_{ij}^{eq}-r_{jk}^{eq}}{r_{ij}^{eq}+r_{jk}^{eq}}\right)^2\right)
\end{equation}$
ijk - это С-С-Н
Параметры были взяты из таблиц в статье: $\begin{equation} \\K^\theta_{ijk} = 16.98\;ккал/моль/рад^2 \end{equation}$
Полученные константы схожи по порядку величины, но все же сильно различаются. В статье сказано, что, обычно, если один из крайних атомов угла является C, то константа приблизительно равна 50 ккал/моль/радˆ2. При этом при рассчете по формуле из этой статьи почему-то получилось около 17. Ошибку найти не удалось. Причины расхождения значений могут быть те же, что и в предыдущем пункте.
Расчет энергии для торсионного угла
E_torsion=Energy(inp_init,6,7,12,360,31)
x_o=x_generation(180,12,360,30)
y_o=E_torsion
fitfunc = lambda p, x: p[0]/2*(1+numpy.cos(p[1]*x))+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(-200,200)
plt.show()
На заданном диапазоне торсионных углов функция имеет 3 минимума (-180 и 180 - это один минимум)
Энергия для связи С-С с шагом 0.1
E_bind2=Energy(inp_init,4,3,0.1,1.,21)
x_o=x_generation(1.52986,0.1,1.,20)
y_o=E_bind2
fitfunc = lambda p, x: numpy.power((1-numpy.exp(-p[0]*(x-p[1]))),2)+ p[2]#p[0]/(x-p[1])+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(0,2)
plt.show()
На заданном интервале расстояний между С-С эту зависимость можно было бы аппроксимировать следующей функцией:
$\begin{equation}
\\f\left(x\right) = \left(1-e^{-a(x-b)}\right)^2+c
\end{equation}$
параметры: [ 1.34834595, 1.43021128, -79.38876095]