Вычисление параметров для молекулярной механики¶
Энергия молекулы в GAFF описывается суммой:
Найдем с помощью минимизации энергии из квантовых расчетов (orca) параметры r_eq, q_eq, n
length_start = 1.52986
HCC_start = 111.200
CC_start = 0
inp_template = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 {0} 0 0
H 1 2 0 1.08439 {1} 0
H 1 2 3 1.08439 {1} 120
H 1 2 3 1.08439 {1} -120
H 2 1 3 1.08439 {1} {2}
H 2 1 6 1.08439 {1} 120
H 2 1 6 1.08439 {1} -120
*
'''
# координаты атомов 7 и 8 заданы относительно 6,
# поэтому если его поворачивать, то 7 и 8 будут поворачиваться согласованно
def get_inp(length, HCC, CC):
return inp_template.format(length, HCC, CC)#, CC+120, CC-120, -CC+120, -CC-120)#CC+120, CC-120, CC+180, CC+120, CC-120)
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
import subprocess
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]
outlines = out.splitlines()
# extract energy: FINAL SINGLE POINT ENERGY'
# and return it as float
for l in outlines:
#print(l)
if 'FINAL SINGLE' in str(l):
return float(l.split()[-1])
run_orca(inp.format(length_start, HCC_start, CC_start))
Длина связи¶
Вычислим энергию для значений длины C-C связи от 1.3 до 1.7 ангстрем с шагом 0.02
x_o = np.arange(1.3, 1.7, 0.02)
x_o
y_o = np.array([run_orca(get_inp(i, HCC_start, CC_start)) for i in x_o])
y_o
x_o[np.argmin(y_o)]
Проведем квадратичную функцию через полученные пары длина связи - значение энергии.
fitfunc = lambda p, x: p[0] * (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
res = optimize.least_squares(errfunc, p0[:], args=(x_o, y_o), verbose=2, method='lm')
p1 = res.x
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)
plt.show()
Точка минимума энергии связи, приближенной кривой, равна параметру p[1]: ее значение равно 1.554. Видно, что парабола не очень хорошо приближает полученную функцию: истинный минимум энергии достигается при длине связи 1.54 ангстрема.
Валентный угол¶
Вычислим энергию для значений валентного угла H-C-C от 109.2 до 113.2 с шагом 0.2 и аналогично предыдущему пункту приблизим функцию квадратичной.
x_o1 = np.arange(109.2, 113.2, 0.2)
x_o1
y_o1 = np.array([run_orca(get_inp(length_start, i, CC_start)) for i in x_o1])
fitfunc = lambda p, x: p[0] * (p[1] - x) ** 2 + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1,1, 0] # Initial guess for the parameters
res = optimize.least_squares(errfunc, p0[:], args=(x_o1, y_o1), verbose=2, method='lm')
p1 = res.x
print("Optimized params:", p1)
#Plot it
plt.plot(x_o1, y_o1, "ro", x_o1,fitfunc(p1,x_o1),"r-",c='blue',alpha=0.5)
plt.xlim(109,114)
plt.show()
Получена точка минимума 111.7. В силовом поле GAFF2 значение этого угла считается равным 109.8 (C-sp3 : C-sp3: H связанный с алифатическим C без стягивающих электроны групп):
c3-c3-hc 46.8 109.80 SOURCE3_SOURCE5 179054 0.7972
Это значение ближе к углу тетраэдрического атома (около 109.5). Вероятно, отклонение от этого значения в наших данных вызвано сильной асимметрией молекулы (один из заместителей, C, крупнее, и от него отталкиваются атомы водорода).
Двугранный угол¶
Вычислим энергию для значений двугранного угла от -180 до 180 с шагом 6, и приблизим ее косинусом (параметры: амплитуда, частота, сдвиг фазы и постоянная энергия, не зависящая от угла).
x_o2 = np.arange(-180, 181, 6)
x_o2
y_o2 = np.array([run_orca(get_inp(length_start, HCC_start, i)) for i in x_o2])
y_o2
fitfunc = lambda p, x: p[0] * np.cos(x * p[1] * (np.pi) / 180 - p[2]) - p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1, 2.5, 0, 0] # Initial guess for the parameters
res = optimize.least_squares(errfunc, p0[:], args=(x_o2, y_o2), verbose=2, xtol=1e-15, method='lm')
p1 = res.x
print("Optimized params:", p1)
#Plot it
plt.plot(x_o2, y_o2, "ro", x_o2,fitfunc(p1,x_o2),"r-",c='blue',alpha=0.5)
plt.xlim(-180, 180)
plt.show()
У функции 3 максимума (n=3), соответствующих трем заслоненным конформациям этана.