Вычисление параметров для молекулярной механики¶
Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов. Считать будем для молекулы этана.
Это оптимизированная структура этана в виде z-matrix. Нужно сделать 20 значений связи С-С с шагом +-0.02 (по 10 в разные стороны). Запишем в 20 файлов ethaneXX.txt, где ХХ - номер от 0 до 20.
def gen_inp():
for i in range(-10,10):
bond = 1.52986 + i*0.02
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
output_name = 'ethane_%i.txt' %(i+10)
output_file = open(output_name,"w")
output_file.writelines(inp)
output_file.close()
Будем запускать ORCA с этими параметрами. На выходе - значение энергии структуры (FINAL SINGLE POINT ENERGY).
import subprocess
def run_orca(num):
cmd1 = '/srv/databases/orca/orca ethane_%d.txt > orca_temp.txt' %num
cmd2 = 'grep "FINAL SINGLE POINT ENERGY" orca_temp.txt'
subprocess.call(cmd1, shell=True)
p = subprocess.Popen(cmd2, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
p_status = p.wait()
out=iter(p.stdout.readline, b'')
for line in out:
energy = line.strip(None).strip('FINAL SINGLE POINT ENERGY').strip()
return float(energy)
Получим значения энергии для каждой длины С-С связи:
for i in range(0,21):
print run_orca(i)
Нарисуем полученную зависимость энергии от длины С-С связи.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from scipy.fftpack import fft, ifft
import math
# x array , lengths of C-C bonds
x_o=[]
# y array, energies
y_o=[]
for i in range(-10,10):
bond = 1.52986 + i*0.02
x_o.append(bond)
energy = run_orca(i+10)
y_o.append(energy)
x_o=np.array(x_o)
y_o=np.array(y_o)
print x_o
print 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()
Смотрим на строку Optimized params: [ 0.66824637 1.55389273 -79.08233322]
Первое число - это коэффициент жесткости связи как пружинки: 0.66824637 эВ = 15.41 ккал/моль = 64.48 кДж/моль. Для пересчета я использовала конвертер: http://www.colby.edu/chemistry/PChem/Hartree.html
Парабола не точно совпадает с точками, т.к. межатомный потенциал на самом деле не является параболическим, а описывается более сложными функциями.
Теперь проведем аналогичные операции для валентного угла HCС. Его значения должны изменяться от 109.2 до 113.2. Аналогично, запишем в 20 файлов ethaneXX.txt, где ХХ - номер от 20 до 40.
def gen_inp2():
for i in range(0,21):
angle = 109.2 + i*0.2
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
output_name = 'ethane_%i.txt' %(i+20)
output_file = open(output_name,"w")
output_file.writelines(inp)
output_file.close()
gen_inp2()
ls
for i in range(20,41):
print run_orca(i)
# x array , HCC angles
x_o=[]
# y array, energies
y_o=[]
for i in range(0,21):
angle = 109.2 + i*0.2
x_o.append(angle)
energy = run_orca(i+20)
y_o.append(energy)
x_o=np.array(x_o)
y_o=np.array(y_o)
print x_o
print 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.title('Energy(angle)')
plt.xlim(108,115)
plt.show()
Здесь коэффициент жесткости равен 0.00397 кДж/моль.
Проделаем аналогичные операции для торсионного угла HCCH, его значения должны изменяться от -180 до 180 c шагом 12. Аналогично, запишем в 30 файлов ethaneXX.txt, где ХХ - номер от 41 до 72.
def gen_inp3():
x_o=[]
for i in range(0,31):
dihedral = -180 + i*12
x_o.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
output_name = 'ethane_%i.txt' %(i+41)
output_file = open(output_name,"w")
output_file.writelines(inp)
output_file.close()
return x_o
print gen_inp3()
x_o = gen_inp3()
y_o=[]
for i in range(41,72):
energy = run_orca(i)
print energy
y_o.append(energy)
x_o=np.array(x_o)
y_o=np.array(y_o)
print x_o
print y_o
#former 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
#and what if function is k*cos(nx) + a?
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))
yfft = fft(y_o)
print "Optimized params:", p1
#Plot it
#plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",x_o, yfft,"r-",c='blue',alpha=0.5)
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.title('Energy(dihedral)')
plt.xlim(-200,200)
plt.show()
Таким образом, полученную зависимость энергии от торсионного угла можно хорошо аппроксимировать функцией 0,002*cos(0,99x)-79.
Наблюдаем три минимума энергии при изменении торсионного угла HCCH от -180 до 180. Это соответствует трем заторможенным конформациям молекулы этана.
Увеличим шаг до 0,1 при расчете связи. Полученная зависимость плохо аппроксимируется параболой и хорошо - потенциалом Морзе.
def gen_inp4():
x_o=[]
for i in range(-10,10):
bond = 1.52986 + i*0.1
x_o.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
output_name = 'ethane_%i.txt' %(i+100)
output_file = open(output_name,"w")
output_file.writelines(inp)
output_file.close()
return x_o
print gen_inp4()
x_o = gen_inp4()
y_o=[]
for i in range(90,110):
energy = run_orca(i)
print energy
y_o.append(energy)
x_o=np.array(x_o)
y_o=np.array(y_o)
print x_o
print 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.title('Energy(bond)')
#plt.xlim(108,115)
plt.show()
x_o=np.array(x_o)
y_o=np.array(y_o)
print x_o
print y_o
#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.xlim(108,115)
plt.show()
