Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов. Будем использовать молекулу этана.
Это оптимизированная структура этана в виде z-matrix. Сделаем 40 значений связи С-С с шагом +-0.02 (20 в сторону уменьшения и 20 в сторону увеличения). Запишем в 40 файлов ethaneXX.txt, где ХХ - номер от 0 до 40.
def gen_inp():
for i in range(-20,20):
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+20)
output_file = open(output_name,"w")
output_file.writelines(inp)
output_file.close()
gen_inp()
Запускаем 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,40):
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(-20,20):
bond = 1.52986 + i*0.02
x_o.append(bond)
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(x-b)^2 + a
fitfunc = lambda p, x: p[0]*pow(x-p[1],2) + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
#function is f(x)=k*(ln(сx-b))^2 + a
my_fitfunc = lambda p, x: p[0]*pow(np.log(p[3]*x-p[1]),2) + p[2]
my_errfunc = lambda p, x, y: my_fitfunc(p, x) - y
#function is f(x)=с*(1-e^(-k(x-a)))^2
mors_fitfunc = lambda p, x: p[2]*pow(1-np.e**(-p[1]*(x-p[0])),2)
mors_errfunc = lambda p, x, y: mors_fitfunc(p, x) - y
p0 = [1, 1, -80, 1] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:-1], args=(x_o, y_o))
p2, success = optimize.leastsq(my_errfunc, p0[:], args=(x_o, y_o))
p3, success = optimize.leastsq(mors_errfunc, p0[:-1], args=(x_o, y_o))
p4 = [1.30722816, 0.30222351, -79.26679609, 0.80124823]
print "Optimized params for x^2:", p1
print "Optimized params for ln^2:", p2
print "Optimized params for Morse:", p3
#Plot it
plt.plot(x_o, y_o, "bo", x_o,fitfunc(p1,x_o),"r-", x_o, my_fitfunc(p2,x_o), "b-", x_o, mors_fitfunc(p3,x_o), "g-", x_o, my_fitfunc(p4,x_o), 'k-', alpha=0.5)
plt.xlim(1,2)
plt.show()
Разберемся с подбранными параметрами оптимизации (строки Optimized params for ...): [0.77924885 1.59808014 -79.20515722] для x^2 [0.27906518 1.14337602 -79.19734782 1.39222098] для (ln x)^2
Для параболы первое число - это коэффициент жесткости связи как пружинки: 0.77924881 эВ = 17.97 ккал/моль = 75.19 кДж/моль для x^2 Коэффициенты логарифма интерпретировать трудно, зато если их подобрать, то можно учесть, что связь - не пружинка. Для пересчета использовался конвертер: http://www.colby.edu/chemistry/PChem/Hartree.html
Парабола (показана красной линией) не точно совпадает с точками, т.к. межатомный потенциал на самом деле не является параболическим, а описывается более сложными функциями. Логарифм в квадрате (показан синей линией) выглядит гораздо убедительнее, то есть математически логарифм подходит лучше. Потенциал Морзе на таких маленьких расстояниях не работает (показан зелёным). Черным показан график логарифма в квадрате, оптимизированный под большой шаг (см. дальше)
Увеличим шаг до 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+110)
output_file = open(output_name,"w")
output_file.writelines(inp)
output_file.close()
return x_o
x_o = gen_inp4()
y_o=[]
for i in range(100,120):
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(x-b)^2 + a
fitfunc = lambda p, x: p[0]*pow(x-p[1],2) + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
#function is f(x)=k*(ln(сx-b))^2 + a
my_fitfunc = lambda p, x: p[0]*pow(np.log(p[3]*x-p[1]),2) + p[2]
my_errfunc = lambda p, x, y: my_fitfunc(p, x) - y
#function is f(x)=(1-e^(-k(x-a)))^2
mors_fitfunc = lambda p, x: p[2]*pow(1-np.e**(-p[1]*(x-p[0])),2)
mors_errfunc = lambda p, x, y: mors_fitfunc(p, x) - y
p0 = [0.5, 0.5, -80, 1] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:-1], args=(x_o, y_o))
p2, success = optimize.leastsq(my_errfunc, p0[:], args=(x_o, y_o))
p3, success = optimize.leastsq(mors_errfunc, p0[:-1], args=(x_o, y_o))
print "Optimized params for x^2:", p1
print "Optimized params for ln^2:", p2
print "Optimized params for Morse:", p3
p4 = [0.27906518, 1.14337602, -79.19734782, 1.39222098]
#Plot it
plt.plot(x_o, y_o, "bo", x_o,fitfunc(p1,x_o),"r-", x_o, my_fitfunc(p2,x_o), "b-", x_o, mors_fitfunc(p3,x_o), "g-", x_o, my_fitfunc(p4,x_o), 'k-', alpha=0.5)
plt.xlim(0.5,2.5)
plt.show()
Ошибка возникает из-за ограниченности области определения логарифма, оптимизтрованного под маленький шаг. Интересно, что логарифм в квадрате хорошо апроксимирует зависимости энергии связи от расстояния между атомами, но, к сожалению, для каждого шага разные параметры.
Теперь проведем аналогичные операции для валентного угла HCС. Его значения должны изменяться от 109.2 до 113.2 с шагом 0.2. Аналогично, запишем в 20 файлов ethaneXX.txt, где ХХ - номер от 40 до 60.
def gen_inp2():
for i in range(0,20):
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+40)
output_file = open(output_name,"w")
output_file.writelines(inp)
output_file.close()
gen_inp2()
for i in range(40,60):
print run_orca(i)
# x array , HCC angles
x_o=[]
# y array, energies
y_o=[]
for i in range(0,20):
angle = 109.2 + i*0.2
x_o.append(angle)
energy = run_orca(i+40)
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, -80] # 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.00392 кДж/моль.
Проделаем аналогичные операции для торсионного угла HCCH, его значения должны изменяться от -180 до 180 c шагом 12. Аналогично, запишем в 31 файл ethaneXX.txt, где ХХ - номер от 60 по 90.
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+60)
output_file = open(output_name,"w")
output_file.writelines(inp)
output_file.close()
return x_o
x_o = gen_inp3()
print x_o
y_o=[]
for i in range(60,91):
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 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, -80] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
yfft = fft(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.title('Energy(dihedral)')
plt.xlim(-200,200)
plt.show()
Таким образом, полученную зависимость энергии от торсионного угла можно хорошо аппроксимировать функцией 0,002*cos(0,99x)-79.
Наблюдаем три минимума энергии при изменении торсионного угла HCCH от -180 до 180. Это соответствует трем заторможенным конформациям молекулы этана.