Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов.
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
Оптимизированная структура этана в виде z-matrix:
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
# функция по запуску ORCA
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:
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()
Зависимость энергии молекулы от величины валентного угла, варьируем его от 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()
Зависимость энергии молекулы от величины торсионного угла, варьируем его значение от -180 до 180 с шагом 12
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()
Энергетические максимумы соответствуют заслоненной конформации, минимумы - заторможенной.