Вычисление ковалентных параметров для молекулярной механики¶
Нам нужно определить константы ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов.
In [1]:
# импортируем библиотеки
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
import pandas as pd
import subprocess
In [2]:
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
*
'''
In [3]:
# модифицируем функцию по запуску ORCA
def run_orca(inp):
with open('/home/pavel/Desktop/8_семестр/Models_polymers/orca.inp', 'w') as outfile:
outfile.write(inp)
with subprocess.Popen("/home/pavel/orca/orca /home/pavel/Desktop/8_семестр/Models_polymers/orca.inp", shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) as p:
out=p.communicate()[0].decode()
for line in out.splitlines():
if 'FINAL SINGLE POINT' in line:
return float(line.strip().split()[4])
return None
In [4]:
run_orca(inp)
Out[4]:
Построим зависимость энергии молекулы от длины одной связи. Нужно брать разные значения по длине одной связи С-С и рассчитать 20 разных длинн связи с шагом 0.02 ангстрема.
In [5]:
# сделаем сразу для валентного угла HCС, торсионного угла CC, длины связи СС
inp1 = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 {:.5f} 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 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
inp2 = '''!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 {:.3f} 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 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
inp3 = '''!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 {:.1f}
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
In [10]:
def fitparam(x_o,y_o,p0):
#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
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print ("Optimized params:", p1)
plt.plot(x_o, y_o, '.')
plt.plot(x_o,fitfunc(p1,x_o))
plt.legend()
plt.show()
def fitparamsin(x_o,y_o,p0):
#function is f(x)=k(b-x)^2 + a
fitfunc = lambda p, x: np.sin(x*p[0]+p[1])*p[2]+p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
res = optimize.least_squares(errfunc, p0[:], args=(x_o, y_o), verbose=2, xtol=1e-15, method='lm')
p1 = res.x
print("Optimized params:", p1)
plt.plot(x_o, y_o, '.')
plt.plot(x_o,fitfunc(p1,x_o))
plt.legend()
plt.show()
In [7]:
bond_lengths = np.arange(1.3, 1.7, 0.02)
bond_length_energies = np.array([run_orca(inp1.format(value)) for value in bond_lengths])
fitparam(bond_lengths,bond_length_energies, [1,1.5, -79])
Получили близкое к статье значение - 1.53 вместо 1.526. Возможно, они использовали другие модели для потенциала взаимодействия.
In [8]:
angle_values = np.arange(109.2, 113.4, 0.2)
angle_energies = np.array([run_orca(inp2.format(value)) for value in angle_values])
fitparam(angle_values,angle_energies,[1,1, 0])
In [11]:
torsion_values = np.arange(-180, 192, 12)
torsion_energies = np.array([run_orca(inp3.format(value)) for value in torsion_values])
fitparamsin(torsion_values,torsion_energies,[1, 2.5, 0, 0])
Видно, что есть четыре минимума -180, 180, -60, 60. Если учесть периодичность, то минимумов три. Они должны соответствовать заторможенной конформации этана. Соответственно максимумы - заслоненной.