Вычисление параметров для молекулярной механики¶
Импортировать необходимые модули
import subprocess
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from numpy import cos, exp
Задать функцию, запускающую Orca для заданной молекулы и выдающую результат
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]
return out
Задать цикл, в котором для разных длин С-С связи c шагом 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
*
''' # optimized ethane molecule
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)# change the CC bond length
ini[3]=' '.join(ind)
inp1='\n'.join(ini)
s=run_orca(inp1)
k=s.decode('utf-8').index('Total Energy')
x.append(float(ind[4])) # bond length
y.append(float(s[k:k+100].split()[5])) # energy in meV
Построить график зависимости энергии молекулы от длины связи С-С, аппроксимировать квадратичной функцией
x_o=x
y_o=y
#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 = [16,1.52986, -2151] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))# optimization of function
print ("Optimized params:", p1)
#Plot
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()
Как можно заметить, аппроксимирующая линия плохо ложится на точки, но формой более менее им соответствует. Минимум энергии находится на расстоянии 1,55 ангстрем, тогда как по данным GAFF оптимальная длина 1.526, причина в плохой аппроксимации.
Аналогично рассчитать энергии и пoстроить график для угла HCC.
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=run_orca(inp1)
k=s.decode('utf-8').index('Total Energy')
x.append(float(ind[5]))
y.append(float(s[k:k+100].split()[5]))
Построить график зависимости энергии молекулы от угла, аппроксимировать квадратичной функцией
x_o=x
y_o=y
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2]
errfunc = lambda p, x, y: fitfunc(p, x) - y
p0 = [0.001,111.2, -2151]
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, 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.xlim(105,115)
plt.show()
Здесь аппроксимирующая кривая вписалась лучше, чем в предыдущем случае. Оптимальное значение угла составило 111,2, что соответствует оптимальному значению.
Аналогично рассчитать энергии и построить график для поворота по торсионному углу CC
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=run_orca(inp1)
k=s.decode('utf-8').index('Total Energy')
x.append(float(ind[6]))
y.append(float(s[k:k+100].split()[5]))
В отличие от предыдущих случаев аппроксимирующая функция синусоидальная
x_o=x
y_o=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_o, 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.xlim(-180,180)
plt.show()
Аппроксимирующая функция имеет три минимума энергии, соответствующих трем заторможенным конформациям.
Рассчитать зависимость энергии от длины связи с шагом 0,1 ангстрема (таким образом интервал длин связи увеличивается в 5 раз). Построить график
x=[]
y=[]
for i in range(-10,10):
ini=inp.splitlines()
ind=ini[3].split()
ind[4]='{:.5f}'.format(float(ind[4])+0.1*i)# change the CC bond length
ini[3]=' '.join(ind)
inp1='\n'.join(ini)
s=run_orca(inp1)
k=s.decode('utf-8').index('Total Energy')
x.append(float(ind[4]))
y.append(float(s[k:k+100].split()[5]))
x_o=x
y_o=y
plt.plot(x_o,y_o,"ro",c='blue',alpha=0.5)
plt.xlim(0,3)
plt.show()
Как можно заметить, зависимость энергии от длины связи далека от квадратичной. При длине связи менее оптимальной для аппроксимации можно использовать экспоненциальную функцию, а при длине связи больше оптимальной - линейную.
Рассчитать зависимость энергии молекулы от величины угла НСС, также в увеличенном в 5 раз интервале
x=[]
y=[]
for i in range(-10,10):
ini=inp.splitlines()
ind=ini[4].split()
ind[5]='{:.3f}'.format(float(ind[5])+i)
ini[4]=' '.join(ind)
inp1='\n'.join(ini)
s=run_orca(inp1)
k=s.decode('utf-8').index('Total Energy')
x.append(float(ind[5]))
y.append(float(s[k:k+100].split()[5]))
x_o=x
y_o=y
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2]
errfunc = lambda p, x, y: fitfunc(p, x) - y
p0 = [0.001,111.2, -2151]
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, 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.xlim(100,125)
plt.show()
Что касается зависимости энергии от угла, то его можно аппроксимировать квадратичной функцией в широких пределах