Домашнее задание №5¶
#imports
import subprocess
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
Молекула этана:
\begin{equation*} C_{2}H_{6} \end{equation*}Сначала зададим 20 значений длины связи C-C с шагом 0.02. Для каждого из случаев рассчитаем величину значения энергии с помощью ORCA. Затем на этих значениях восстановим функцию зависимости энергии от длины связи.
На вход ORCA будет подаваться шаблон, где будет изменено значение длины связи C-C.
inp_length = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 {length} 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
*
'''
Функция по расчету энергии для входных файлов,где будет меняться значение длины связи|величины торсионного угла|величины валентного угла. Грубо говоря, функция возвращает иксы и игрики(величины энергии):
def energy(template,variable,begin,step,n):
value=begin
e=[]
x=[]
if (variable=="length"):
for i in range(0,n):
value=begin+i*step
dupl=template.format(length=str(value))
e.append(run_orca(dupl))
x.append(value)
elif (variable=="val_angle"):
for i in range(0,n):
value=begin+i*step
dupl=template.format(val_angle=str(value))
e.append(run_orca(dupl))
x.append(value)
elif (variable=="tors_angle"):
for i in range(0,n):
value=begin+i*step
dupl=template.format(tors_angle_1=str(value),tors_angle_2=str(value-60),tors_angle_3=str(value+60))
e.append(run_orca(dupl))
x.append(value)
return e,x
Далее идет функция по заупуску ORCA, где для каждого входного файла в формате строки будет найдено значение энергии:
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'
# and return it as float
for s in out.splitlines():
if s.startswith("FINAL SINGLE POINT ENERGY"):
return(float(s[len("FINAL SINGLE POINT ENERGY"):]))
Рассчитаем величины энергий для входных файлов по длине связи С-С
#path="/home/shad/hse/mrPuk/hw5/length/"
e_length,x_length=energy(inp_length,variable="length",begin=1.52986-10*0.02,step=0.02,n=20)
Функция построения зависимостей (целевая функция - квадратичная):
def plot(x_o,y_o):
#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()
plot(x_length,e_length)
Как видно из вышеприведенного графика, аппроксимация оказалась вполне удачной. Теперь сделаем аналогичные действия для валентного угла:
inp_val_angle = '''!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 {val_angle} 0
H 1 2 3 1.08439 {val_angle} 120
H 1 2 3 1.08439 {val_angle} -120
H 2 1 3 1.08439 {val_angle} 180
H 2 1 5 1.08439 {val_angle} 120
H 2 1 5 1.08439 {val_angle} -120
*
'''
e_val_angle,x_val_angle=energy(inp_val_angle,variable="val_angle",begin=111.2-10*0.2,step=0.2,n=21)
plot(x_val_angle,e_val_angle)
Как видим, аппроксимиция удалась. Теперь про торсионный угол:
inp_tors_angle = '''!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.2 0
H 1 2 3 1.08439 111.2 120
H 1 2 3 1.08439 111.2 -120
H 2 1 3 1.08439 111.2 {tors_angle_1}
H 2 1 5 1.08439 111.2 {tors_angle_2}
H 2 1 5 1.08439 111.2 {tors_angle_3}
*
'''
e_tors_angle,x_tors_angle=energy(inp_tors_angle,variable="tors_angle",begin=-180,step=12,n=31)
plot(x_tors_angle,e_tors_angle)
Как видим, на графике сверху не получилось аппроксимировать функцию. Скорее всего эта функция имеет зависимость от синуса и косинуса. Относительно минимумов функции - их три.
def plot_s_c(x_o,y_o):
#function is f(x)=a*sin(b*x)+c*cos(d*x)
fitfunc = lambda p, x: p[0]*np.sin(p[1]*x) + p[2]*np.cos(p[3]*x)+p[4] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1,1,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()
plot_s_c(np.array(x_tors_angle),np.array(e_tors_angle))
Ну вот, все очень даже неплохо заапроксимировалось.
Теперь увеличим шаг до 0.1 ангстрема при расчете связи и построим зависимость:
e_length,x_length=energy(inp_length,variable="length",begin=1.52986-12*0.1,step=0.1,n=41)
plot(x_length,e_length)
Не совсем хорошо. Попробуем степенную функцию...
def plot_2(x_o,y_o):
#function is f(x)=a*(1/x)^b+c
fitfunc = lambda p, x: p[0]*pow(1/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(1,2)
plt.show()
plot_2(np.array(x_length),e_length)
Как видно, данную зависимость можно аппроксимировать степенной функцией. Кстати, кол-во точек в промежутке (0;1), которые берутся для построения зависимости играет важную роль: чем их больше, тем точнее кривая "ложится" на заданные значения. Но с другой стороны тенденция возрастания значений функции при большем "х" говорит о том, что тут может пригодится и другая функция. Попробуем Morse potential. Его формула следующая:
Используем ее при построении зависимости:
def plot_morse(x_o,y_o):
#morse function
fitfunc = lambda p, x: p[0]*pow(1-np.exp(-p[1]*(x-p[2])),2) + p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1,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()
plot_morse(np.array(x_length),e_length)
Как видно, аппроксимация стала еще лучше.