In [ ]:
Занятие 6. Вычисление параметров для молекулярной механики
In [1]:
#Definig orca function
def run_orca(inp):
with open('orca.inp', 'w') as outfile:
outfile.write(inp)
p = subprocess.Popen("/srv/databases/orca/orca orca.inp",
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out=p.communicate()[0]
outlines=out.splitlines()
# extract energy: FINAL SINGLE POINT ENERGY'
# and return it as float
for i in outlines:
if "FINAL SINGLE POINT ENERGY" in i:
return float(i.split()[4])
In [ ]:
Теперь рассчитаем энергии для 20 разных длинн связи с шагом 0.02 ангстрема для расчёта энергии в ORCA с разными значениями по длине одной из связей (C-C).
In [2]:
import subprocess, os
bond=1.52986
x_axis=[]
y_axis=[]
print "bond\tenergy"
for i in range(-10,10):
bond = 1.52986 + i*0.02
x_axis.append(bond)
inp = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 %f 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
*
''' % bond
energy = run_orca(inp)
print str(bond)+"\t"+str(energy)
y_axis.append(energy)
In [ ]:
Построим график и линию тренда
In [3]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
x_o=np.array(x_axis)
y_o=np.array(y_axis)
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, "bo", x_o,fitfunc(p1,x_o),"g-",alpha=0.5)
plt.xlim(1.3,1.73)
plt.ylabel('Energy, eV')
plt.xlabel('C-C bond length, angstrem')
plt.show()
In [ ]:
Сделаем то же самое, но для валентного угла HCС. Его значения должны изменяться от 109.2 до 113.2
In [8]:
x_axis=[]
y_axis=[]
print "angle\tenergy"
for i in range(-10,10):
angle = 111.2 + i*0.2
x_axis.append(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 %f 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
*
''' % angle
energy = run_orca(inp)
print str(angle)+"\t"+str(energy)
y_axis.append(energy)
In [9]:
x_o=np.array(x_axis)
y_o=np.array(y_axis)
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, "bo", x_o,fitfunc(p1,x_o),"g-",alpha=0.5)
plt.xlim(108,115)
plt.ylabel('Energy, eV')
plt.xlabel('Angle')
plt.show()
In [ ]:
И тепеь то же самое, но для торсионного угла CC, его значения должны изменяться от -180 до 180 c шагом 12.
In [15]:
x_axis=[]
y_axis=[]
print "dihedral\tenergy"
for i in range(-15,16):
dihedral = 0 + i*12
x_axis.append(dihedral)
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 %f
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
''' % dihedral
energy = run_orca(inp)
print str(dihedral)+"\t"+str(energy)
y_axis.append(energy)
In [21]:
x_o=np.array(x_axis)
y_o=np.array(y_axis)
#I was forced to think about a new function, the old one looked like ... a power, but here is definetly a periodic one
#It was strange that sin didn't work, but cos did.
fitfunc = lambda p, x: p[0]*np.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, "bo", x_o,fitfunc(p1,x_o),"g-",alpha=0.5)
plt.xlim(-200,200)
plt.ylabel('Energy, eV')
plt.xlabel('C-C bond length, angstrem')
plt.show()
In [ ]:
Теперь увеличим шаг до 0.1 при расчёте энергии для длин связи.
In [22]:
print "bond\tenergy"
for i in range(-20,20):
bond = 1.52986 + i*0.01
x_axis.append(bond)
inp = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 %f 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
*
''' % bond
energy = run_orca(inp)
print str(bond)+"\t"+str(energy)
y_axis.append(energy)
In [27]:
x_o=np.array(x_axis)
y_o=np.array(y_axis)
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, "bo", x_o,fitfunc(p1,x_o),"g-",alpha=0.5)
plt.xlim(1.3,1.73)
plt.ylabel('Energy, eV')
plt.xlabel('C-C bond length, angstrem')
plt.show()
x_o_1=np.array(x_axis)
y_o_1=np.array(y_axis)
In [ ]:
При построении линии тренда выходит что-то странное. Я построила в Excel то же самое и там точки прекрасно аппроксимируется линией тренда y = 6E-05x2 - 0,0029x - 79,165. Но вдруг получится аппроксимировать потенциалом Морзе? Потенциал Морзе — функция потенциальной энергии элетростатического поля, предложенная американским физиком Филипом Морзе, как приближение для энергии двухатомной молекулы.
In [26]:
x_o=np.array(x_axis)
y_o=np.array(y_axis)
#function is f(x)=(1-e^(-k(x-a)))^2
fitfunc = lambda p, x: p[0]*pow(1-np.e**(-p[1]*(x-p[2])),2) # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [-79,45,0.5] # 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.title('Energy(bond)')
plt.show()
In [29]:
print "bond\tenergy"
for i in range(-10,10):
bond = 1.52986 + i*0.01
x_axis.append(bond)
inp = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 %f 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
*
''' % bond
energy = run_orca(inp)
print str(bond)+"\t"+str(energy)
y_axis.append(energy)
In [30]:
x_o=np.array(x_axis)
y_o=np.array(y_axis)
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,0.029, -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, "bo", x_o,fitfunc(p1,x_o),"g-",alpha=0.5)
plt.xlim(1.3,1.73)
plt.ylabel('Energy, eV')
plt.xlabel('C-C bond length, angstrem')
plt.show()
In [ ]: