Вычисление параметров для молекулярной механики
In [1]:
import subprocess
In [2]:
%matplotlib inline
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [43]:
template = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 %s 0 0
H 1 2 0 1.08439 %s 0
H 1 2 3 1.08439 111.200 %s
H 1 2 3 1.08439 111.200 -120
H 2 1 3 1.08439 111.200 180
H 2 1 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
In [50]:
refvalues = ['1.52986', '111.200', '180']
In [54]:
def run_orca(name, insertion, pos, values, template):
result = []
for c, value in enumerate(values):
found = False
with open('%s_%s.inp' %(name, str(c)), 'w') as oinp:
newsertion = insertion[:]
newsertion[pos] = str(value)
newtemp = template % (newsertion[0], newsertion[1], newsertion[2])
oinp.write(newtemp)
p = subprocess.Popen("/srv/databases/orca/orca %s_%s.inp" %(name, str(c)),
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out=p.communicate()[0]
out = out.splitlines()
for outline in out:
if ('FINAL SINGLE' in outline):
found = True
result.append(float(outline.strip('\n').split()[-1]))
if not found:
result.append('error')
return result
In [41]:
def print_plot(name, insertion, pos, start, step, num, fitfunc, p0):
if (pos >= len(insertion)):
print("Wrong insertion index!\n")
return 9
x_o = [start+step*i for i in range(num+1)]
y_o = run_orca(name, insertion, pos, x_o, template)
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
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(math.floor(x_o[0]),math.ceil(x_o[-1]))
plt.show()
Оценим длину связи С-С
In [7]:
#function is f(x)=k(b-x)^2 + a
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2]
p0 = [1,1, -79] # Initial guess for the parameters
print_plot('CC_len', refvalues, 0, 1.32986, 0.02, 20, fitfunc, p0)
Проделаем аналогичное для валентного угла HCC
In [8]:
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2]
p0 = [1,1, -79]
print_plot('HCC_ang', refvalues, 1, 109.2, 0.2, 20, fitfunc, p0)
Проделаем то же самое для торсионного угла СС
In [3]:
template = '''!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 %s
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
In [14]:
tors = []
nrg = []
for i in range(31):
found = False
tors.append(-180 + i*12)
with open('CC_tors_%s.inp' % (str(i)), 'w') as oinp:
oinp.write(template % str(tors[i]))
p = subprocess.Popen("/srv/databases/orca/orca CC_tors_%s.inp" % (str(i)),
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out=p.communicate()[0]
out = out.splitlines()
for outline in out:
if ('FINAL SINGLE' in outline):
found = True
nrg.append(float(outline.strip('\n').split()[-1]))
if not found:
nrg.append('error')
print(zip(nrg, tors))
In [16]:
plt.scatter(tors,nrg)
plt.ylim([-79.200,-79.190])
plt.show()
Завсимость явно не квадратичная, а очень даже напоминает одну ояень известную тригонометрическую функцию. Попробуем зафитить с помощью косинуса:
In [27]:
fitfunc = lambda p, x: p[0]*(np.cos(p[1]*x)) + p[2]
errfunc = lambda p, x, y: fitfunc(p, x) - y
p0 = [1, 1, 1]
p1, success = optimize.leastsq(errfunc, p0[:], args=(tors, nrg))
print "Optimized params:", p1
#Plot it
plt.plot(tors, nrg, "ro", tors, fitfunc(p1,tors),"r-",c='blue',alpha=0.5)
plt.xlim(-190,190)
plt.show()
Получили 3 минимума, соотвествующие заторможенной конформации этана: один в 180 (он же -180) и два в -60, 60. Максимумы - заслоненная конформация - соотвественно на 120, -120 и 0. Зафитилось не очень хорошо.
*Увеличим шаг при расчете связи до 0.1
In [28]:
template = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 %s 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 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
In [30]:
lens = []
nrg = []
for i in range(20):
found = False
lens.append(0.52986 + i*0.1)
with open('CC_tors_new_%s.inp' % (str(i)), 'w') as oinp:
oinp.write(template % str(lens[i]))
p = subprocess.Popen("/srv/databases/orca/orca CC_tors_new_%s.inp" % (str(i)),
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out=p.communicate()[0]
out = out.splitlines()
for outline in out:
if ('FINAL SINGLE' in outline):
found = True
nrg.append(float(outline.strip('\n').split()[-1]))
if not found:
nrg.append('error')
print(zip(nrg, lens))
In [32]:
plt.scatter(lens,nrg)
plt.ylim([-80,-73])
plt.show()
Зависимость гиперболическая
Ссылка на архив со всеми созданными файлами.