Вычисление параметров молекулы этана для молекулярной механики¶
Для начала импортируем нобходимые библиотеки и инициализируем функцию, которая будет обращаться к процессу ORCA
import subprocess
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from numpy import cos, exp
def killer_whale(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]
return out
Создадим входной файл для программы ORCA и посчитаем энергию С-С связи для разных длин связи (шаг 0.02)
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 180
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
''' # оптимизированная молекула этана
x=[]
y=[]
for i in range(-10,10):
ini=inp.splitlines()
ind=ini[3].split()
ind[4]='{:.5f}'.format(float(ind[4])+0.02*i) # цикл в котором мы подсчитываем энергию связей разной длины
ini[3]=' '.join(ind)
inp1='\n'.join(ini)
s=killer_whale(inp1)
k=s.decode('utf-8').index('Total Energy')
x.append(float(ind[4])) # длина связи
y.append(float(s[k:k+100].split()[5])) # энергия в миллиэлектрон-вольтах
График зависимости энергии молекулы от длины С-С связи с аппроксимацией параболой
x_0=x
y_0=y
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] #параметры аппроксимации f(x)=k(b-x)^2 + a
p0 = [1,1, -79]
errfunc = lambda p, x, y: fitfunc(p, x) - y
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_0, y_0))
print ("Optimized params:", p1)
plt.plot(x_0, y_0, "ro", x_0,fitfunc(p1,x_0),"r-",c='blue',alpha=0.5)
plt.show()
Полученный в ходе аппроксимации оптимальное значение длины связи близко к литературному, но сама линия на точки ложится плохо. Возможно, что на данном масштабе аппроксимация параболой еще адекватна, а на большем - нет.
Проведем такую же итерацию по разным значениям угла связи Н-С-С
x=[]
y=[]
for i in range(-10,10):
ini=inp.splitlines()
ind=ini[4].split()
ind[5]='{:.3f}'.format(float(ind[5])+0.2*i)
ini[4]=' '.join(ind)
inp1='\n'.join(ini)
s=killer_whale(inp1)
k=s.decode('utf-8').index('Total Energy')
x.append(float(ind[5]))
y.append(float(s[k:k+100].split()[5]))
x_0=x
y_0=y
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2]
p0 = [1,1, -79]
errfunc = lambda p, x, y: fitfunc(p, x) - y
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_0, y_0))
print ("Optimized params:", p1)
plt.plot(x_0, y_0, "ro", x_0,fitfunc(p1,x_0),"r-",c='green',alpha=0.5)
plt.show()
Видно, что аппроксимация вписывается лучше, оптимальное значение угла хорошо соответствует научному консенсусу.
Производим аналогичные расчеты для торсионного угла. Аппроксимация будет тригонометрической
x=[]
y=[]
for i in range(0,361,12):
ini=inp.splitlines()
ind=ini[7].split()
ind[6]=str(int(ind[6])-i)
ini[7]=' '.join(ind)
inp1='\n'.join(ini)
s=killer_whale(inp1)
k=s.decode('utf-8').index('Total Energy')
x.append(float(ind[6]))
y.append(float(s[k:k+100].split()[5]))
x_0=x
y_0=y
fitfunc = lambda p, x: p[0]*cos((x-p[1])*p[2]) + p[3]
errfunc = lambda p, x, y: fitfunc(p, x) - y
p0 = [0.1,0,1,-2155]
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_0, y_0))
print ("Optimized params:", p1)
plt.plot(x_0, y_0, "ro", x_0,fitfunc(p1,x_0),"r-",c='blue',alpha=0.5)
plt.show()
Молекула этана обладает трехсторонней симметрией, что видно из графика энергии. Минимумы энергии у заторможенных конформаций, а максимумы - у заслоненных
Дополним мой практикум расчетом энергии связи в более широких пределах
x=[]
y=[]
for i in range(-10,20):
ini=inp.splitlines()
ind=ini[3].split()
ind[4]='{:.5f}'.format(float(ind[4])+0.1*i)
ini[3]=' '.join(ind)
inp1='\n'.join(ini)
s=killer_whale(inp1)
k=s.decode('utf-8').index('Total Energy')
x.append(float(ind[4]))
y.append(float(s[k:k+100].split()[5]))
x_0=x
y_0=y
fitfunc = lambda p, x: np.multiply(p[0],pow(x,p[1])) + np.multiply(p[2], pow(x,p[3])) - p[4]
p0 = [1, -12, -10, -2, 79]
errfunc = lambda p, x, y: fitfunc(p, x) - y
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_0, y_0), maxfev=15000)
print("Optimized params:", p1)
plt.plot(x_0, y_0, "ro", x_0,fitfunc(p1,x_0),"r-",c='orange',alpha=0.5)
plt.show()
На более широком масштабе аппрокисмировать надо полиномом