%matplotlib inline
import subprocess
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
import warnings
warnings.filterwarnings('ignore')
Оптимизированная структура этана в виде z-matrix с параметрами:¶
- cc - длина связи $C-C$
- ch - длина связи $C-H$
- ah - валентный угол $HCC$
- t1, t2, t3 - торсионные углы
Функция генерации файла orca.inp для расчёта энергии в ORCA:
def create_inp(cc = 1.52986, ch = 1.08439, ah = 111.200, t1 = 180, t2 = 120, t3 = -120):
return '''!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} {2} 0
H 1 2 3 {1} {2} 120
H 1 2 3 {1} {2} -120
H 2 1 3 {1} {2} {3}
H 2 1 5 {1} {2} {4}
H 2 1 5 {1} {2} {5}
*
'''.format(cc, ch, ah, t1, t2, t3)
Функция по запуску ORCA:
def run_orca(inp, inp_name):
with open(inp_name + ".inp", 'w') as outfile:
outfile.write(inp)
p = subprocess.Popen("/home/shad/progs/bin/orca " + inp_name + ".inp",
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out = p.communicate()[0]
# extract energy: FINAL SINGLE POINT ENERGY'
# and return it as float
key_line = b"FINAL SINGLE POINT ENERGY"
energy = [float(line[len(key_line):]) for line in out.splitlines() if line.startswith(key_line)]
return energy[0]
Меняем параметр длины связи $C-C$ и смотрим, какая энергия получается при запуске ORCA:
cc_lengths = np.arange(-10, 10+1)*0.02 + 1.52986
cc_energies = [run_orca(create_inp(cc), "cc_lenght_" + str(cc)) for cc in cc_lengths]
print(cc_lengths)
print(cc_energies)
Функция, отображающая зависимость энергии от изменяемого параметра:
def graph(x_param, y_energy, x_str = 'Length CC', y_str = 'Energy', \
fitfunc_str = 'f(x) = k(b-x)^2 + a',
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2], \
p0 = [1, 1, -79]):
# Error function
errfunc = lambda p, x, y: fitfunc(p, x) - y
# Optimize
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_param, y_energy))
# Plot it
print("Finction:", fitfunc_str)
print("Optimized params:", p1)
err = errfunc(p1, x_param, y_energy)
print("MSE:", np.sum(err**2) / len(err))
plt.plot(x_param, y_energy, "ro", x_param, fitfunc(p1, x_param), "r-", c='blue', alpha=0.5)
plt.xlabel(x_str)
plt.ylabel(y_str)
plt.xlim(min(x_param), max(x_param))
plt.show()
Построим зависимость энергии молекулы от длины $C-C$ связи:
graph(cc_lengths, cc_energies)
Получилась довольно хорошая аппроксимация.
В нашем случае оптимльная длина связи $r_{eq} = 1.56$, в статье $r_{eq} = 1.526$.
Проделаем аналогичные операции для валентного угла $HCC$:
hcc_angles = np.arange(-10, 10+1)*0.2 + 111.2
hcc_energies = [run_orca(create_inp(ah = hcc), "hcc_angle_" + str(hcc)) for hcc in hcc_angles]
print(hcc_angles)
print(hcc_energies)
graph(hcc_angles, hcc_energies, 'Angle HCC', 'Energy')
Аппроксимация идеальна, очень маленький MSE.
Проделаем аналогичные операции для торсионного угла:
t_angles = np.arange(-180, 180, 12)
t_energies = [run_orca(create_inp(t1 = t, t2 = t - 60, t3 = t + 60), "t_angles_" + str(t)) for t in t_angles]
print(t_angles)
print(t_energies)
graph(t_angles, t_energies, 'Torsion', 'Energy')
Видно, что аппроксимация очень плохая.
Также, стоит заметить, что на заданном промежутке функция имеет три максимума и три минимума, что наводит на мысль о периодичности.
Попробуем аппроксимировать какой-нибудь периодической функцией: $f(x) = a \sin (bx + c) + d$
graph(t_angles, t_energies, 'Torsion', 'Energy', 'f(x) = a*sin(bx + c) + d', \
lambda p, x: p[0]*np.sin(p[1]*x + p[2]) + p[3], \
p0 = [1, 1, 1, -79])
В данном случае получилось довольно хорошее приближение.
Теперь увеличим шаг до $0.1$ при расчете длины связи:
new_cc_lengths = np.arange(-10, 10+1)*0.1 + 1.52986
new_cc_energies = [run_orca(create_inp(cc), "cc_lenght_" + str(cc)) for cc in new_cc_lengths]
print(new_cc_lengths)
print(new_cc_energies)
graph(new_cc_lengths, new_cc_energies)
Аппроксимация не удалась.
Попробуем описать экспоненциальной функцией: $f(x) = a e^{−bx+c}+d$
graph(new_cc_lengths, new_cc_energies, 'Length', 'Energy', 'f(x) = a e^{−bx+c}+d', \
lambda p, x: p[0]*1/(np.exp(p[1]*x+p[2]))+p[3],
p0 = [1,1,1, -79])
Теперь аппроксимация получилась довольно хорошей.