Вычисление параметров для молекулярной механики¶
import subprocess
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import optimize
from IPython.display import Image, display
import os
Параметры оптимизированной структуры:
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 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
Создадим функцию для получения структуры этана с заданными параметрами
def make_input(length, angle, torsion):
inp = f'''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 {length} 0 0
H 1 2 0 1.08439 {angle} 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 {torsion}
H 2 1 3 1.08439 111.200 {torsion-120}
H 2 1 3 1.08439 111.200 {torsion-240}
*
'''
return inp
Создадим функцию, которая будет сохранять ввод для программы Orca, запускать ее и сохранять результаты.
def run_orca(length, angle, torsion, experiment):
f_name = "orca_{}_{}_{}_{}.inp".format(experiment, length, angle, torsion,)
print(f_name)
inp = make_input(length, angle, torsion)
with open(f_name, 'w') as outfile:
outfile.write(inp)
with subprocess.Popen("/srv/databases/orca/orca {}".format(f_name),
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) as p:
out=p.communicate()[0].decode()
for l in out.splitlines():
if 'FINAL' in l:
return float(l.strip().split()[4])
return None
print(run_orca(1.52986, 111.200, 180, 'reference'))
1. Зависимость энергии от длины связи C-С¶
Перебираем связи C-C с длиной от 1.35 до 1.75 ангстрем с шагом 0.02.
Ссылка на файлы со всеми входными и выходными файлами: http://kodomo.fbb.msu.ru/~catherine.nesterenko/term8/HW6
x_o = np.arange(1.35,1.75,0.02)
y_o= [run_orca(i,111.200,180,'CC') for i in x_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,1.8)
plt.show()
Вывод: Длина связи в структуре, которую мы получили с помощью аппроксимации - 1.5499, что длиннее, чем 1.5299 для оптимизированной структуры.
2. Зависимость энергии от изменения валентного угла¶
Перебираем валентные углы от 109.2 до 113.2 градусов с шагом 0.2.
x_o = np.arange(109.2, 113.2,0.2)
y_o= [run_orca(1.52986,i,180,'angles') for i in x_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(109, 114)
plt.show()
Вывод: Получили угол в 111.95 градусов. Значение совпало со значением в оптимизированной структуре.
3. Зависимость энергии от торсионного угла C-C¶
Перебираем торсионные углы от -180 до 180 с шагом 12. Поменяли целевую функцию для аппроксимации.
x_o = np.arange(-180, 180, 12)
y_o= [run_orca(1.52986,111.200,i,'torsions') for i in x_o]
#function is f(x)=k*(1 + cos(a*x - b)) + c
fitfunc = lambda p, x: p[0]*(1 + np.cos(p[1]*x - p[2]))+p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [0.25,10,0, -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(-190, 190)
plt.show()
Вывод: Энергетические минимумы соответствуют заторможенной конформации, максимумы - заслоненной. Минимумы в -180 градусах и в +180 градусах можно считать за один, соответствует конформации.
4*. Улучшение аппроксимации длины C-С cвязи¶
Воспользуемся потенциалом Морзе в качестве аппроксимирующей функции.
x_o = np.arange(1.35,1.75,0.02)
y_o= [run_orca(i,111.200,180,'morse') for i in x_o]
#function is f(x)=k*(1 - exp(-a(x - b)))^2 + c
fitfunc = lambda p, x: p[0]*((1-np.exp(-p[1]*(x-p[2])))**2)+p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1,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,1.8)
plt.show()
Вывод: потенциал Морзе гораздо лучше описывает зависимость энергии от длины С-С связи. В оптимизированной структуре длина была 1.52986 ангстрем. Точка минимума, полученная при аппроксимации - 1.52981568. Совпадение вплоть до 4 знака после запятой.