Вычисление параметров для молекулярной механики¶
Все результаты вычислений программой Orca доступны по ссылке: https://kodomo.fbb.msu.ru/~anton.vlasov/term8/prac6/.
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
import sys
from tqdm import tqdm
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 get_inp(l,val,t):
inp = f'''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 {l} 0 0
H 1 2 0 1.08439 {val} 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 {t}
H 2 1 3 1.08439 111.200 {t-120}
H 2 1 3 1.08439 111.200 {t-240}
*
'''
return inp
def run_orca(inp, subfolder):
with open(os.path.join(subfolder, 'orca.inp'), 'w') as outfile:
outfile.write(inp)
with subprocess.Popen("/srv/databases/orca/orca orca.inp", cwd=subfolder,
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) as p:
out=p.communicate()[0].decode()
for l in out.splitlines():
if 'FINAL SINGLE POINT' in l:
return float(l.strip().split()[4])
return None
def get_energy(l, val, t, subfolder):
path = os.path.join(subfolder, f'{l}_{val}_{t}')
if not os.path.exists(path):
os.makedirs(path)
return run_orca(get_inp(l, val, t), path)
print(get_energy(1.52986, 111.200, 180, 'test'))
Изменение длины связи C-C¶
Варьируем длину связи C-C от 1.35 до 1.75 с шагом 0.02 ангстрема.
x_o = np.arange(1.35,1.75,0.02)
y_o= [get_energy(i,111.200,180,'lengths') for i in tqdm(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.52986, а длина связи, полученная аппроксимацией - 1.5499,что немного длиннее.
Изменение валентного угла¶
Варьируем изменение валентного угла с от 109.2 до 113.2 с шагом 0.2 градуса.
x_o = np.arange(109.2, 113.2,0.2)
y_o= [get_energy(1.52986,i,180,'val_angles') for i in tqdm(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.950 градусов, что практически соответствует значению в оптимизированной структуре, которое составляет 111.200 градусов.
Изменение торсионного угла C-C¶
Варьируем значение торсионного угла с шагом 12 от -180 до 180
x_o = np.arange(-180, 180, 12)
y_o= [get_energy(1.52986,111.200,i,'torsion_angles') for i in tqdm(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()
Наблюдаем 4 минимума у функции, но на самом деле минимумы в -180 градусах и в +180 градусах можно считать за один.
Энергетические максимумы и минимумы соответствуют заслоненной и заторможенной конформациям. Начальный угол, равный 180 градусам, соответствует заторможенной конформации.