In [ ]:
import subprocess
import time
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
Необходимо определить константы ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов.
In [2]:
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]
lines = out.splitlines()
energy = None
for line in lines:
line = str(line, 'utf-8')
if line.find("FINAL SINGLE POINT") > -1:
energy = float(line.split()[4])
return energy
Построим зависимость энергии молекулы от длины одной связи.
In [3]:
def iterative_run(name, values, inp_string):
energies = []
for i, x in enumerate(values):
inp = inp_string(x)
energies.append(run_orca(inp))
print("done")
return energies
In [4]:
cc_bond = [x/100 for x in range(120, 180, 2)]
inp_string_bond = lambda x : '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 {} 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
*
'''.format(x)
energies = iterative_run("bond", cc_bond, inp_string_bond)
In [5]:
x_ax=cc_bond
y_ax=energies
#аппроксимируем такой функцией -> f(x)=k(b-x)^2 + a
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2]
errfunc = lambda p, x, y: fitfunc(p, x) - y
p0 = [1,1, -79] # приблизительные значения параметра
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_ax, y_ax))
print("Optimized params:", p1)
#Отображение рассчитанных в Orca значений и минимум энергии
plt.plot(x_ax, y_ax, linewidth = 2, c = "blue",alpha=0.5)
i = np.argmin(y_ax)
plt.axvline(x=x_ax[i],linewidth = 3, c = "blue")
print("calculated energy minimum at {}".format(x_ax[i]))
#Аппроксимация и отметка минимума
plt.plot(x_ax, fitfunc(p1,x_ax), dashes=[30, 5, 10, 5], c='green')
plt.axvline(x= p1[1], c = "green")
# Минимум согласно Gaff
plt.axvline(x=1.526, linewidth = 3, c = "red")
plt.xticks([1.526, float("{0:.4f}".format(p1[1]))], rotation='vertical')
plt.show()
In [6]:
cc_bond_1 = [x/100 for x in range(50, 400, 10)]
energies_1 = iterative_run("bond", cc_bond_1, inp_string_bond)
In [7]:
x_ax=cc_bond_1
y_ax=energies_1
#Теперь в качестве аппроксимации используем сумму степенных функция
fitfunc1 = lambda t, x: np.multiply(t[0],pow(x,t[1])) + np.multiply(t[2], pow(x,t[3])) - t[4]
errfunc1 = lambda t, x, y: fitfunc1(t, x) - y
t0 = [1, -12, -10, -2, 79] # отправные точки для вычисления параметров
t1, success = optimize.leastsq(errfunc1, t0[:], args=(x_ax, y_ax))
print("Optimized params:", t1)
#Отобразим вычисленные значения и минимум энергии
plt.plot(x_ax, y_ax, linewidth = 2, c = "blue",alpha=0.5)
i = np.argmin(y_ax)
plt.axvline(x=x_ax[i], linewidth = 3, c = "blue")
print("minimum energy at {}".format(x_ax[i]))
#Отобразим оптимизированную функцию и ее минимум
plt.plot(x_ax, fitfunc1(t1,x_ax), dashes=[30, 5, 10, 5], c='green')
j = np.argmin(fitfunc1(t1,x_ax))
plt.axvline(x=x_ax[j], c = "green")
# длина связи сс в поле gaff
plt.axvline(x=1.526, linewidth = 3, c = "red")
plt.xticks([1.526, x_ax[j]], rotation = "vertical")
plt.show()
In [8]:
Fl_angles = [x/10 for x in range(1092, 1134, 2)]
inp_string = lambda x :'''!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 6 1.08439 111.200 120
H 2 1 6 1.08439 {} -120
*
'''.format(x)
energies_angles = iterative_run("angle", Fl_angles, inp_string)
In [9]:
x_ax=Fl_angles
y_ax=energies_angles
#Аппроксимируем -> f(x)=k(b-x)^2 + a
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2]
errfunc = lambda p, x, y: fitfunc(p, x) - y
p0 = [1,1, -79] # Отправная точка для параметров
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_ax, y_ax))
print("Optimized params:", p1)
#Отображаем на графике рассчитанных значений и минимума
plt.plot(x_ax, y_ax, linewidth = 2, c='blue',alpha=0.5)
i = np.argmin(y_ax)
plt.axvline(x=x_ax[i], linewidth = 3, c = "blue")
print("minimum energy at {}".format(x_ax[i]))
plt.plot( x_ax, fitfunc(p1,x_ax), dashes=[30, 5, 10, 5], c = "green")
plt.axvline(x=1.11195047e+02, c = "green")
#Отображаем величину угла в поле Amber
plt.axvline(x=109.5, c = "red")
plt.xticks([109.5, x_ax[i]], rotation = "vertical")
plt.show()
In [10]:
dih_angle = [x for x in range(-180, 181, 2)]
inp_string = lambda x :'''!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 {}
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''.format(x)
energies_dih = iterative_run("diherdal", dih_angle, inp_string)
In [11]:
import numpy as np
x_ax=dih_angle
y_ax=energies_dih
# Аппроксимируем функцией a*cos(b * x - c) - d, где а - амплитуда, b - периодичность, c - фаза, d - энергия, независимая от угла
fitfunc = lambda p, x: p[0] * np.cos(np.multiply(x, p[1]/2/np.pi) - p[2]) - p[3]
errfunc = lambda p, x, y: fitfunc(p, x) - y
p0 = [0.001, 0.33, 0, 79] # Затравки для вычисления параметров
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_ax, y_ax))
print("Optimized params:", p1)
#Отображаем на графике
plt.plot(x_ax, y_ax, linewidth = 2, c='blue',alpha=0.5)
plt.plot(x_ax, fitfunc(p1, x_ax), dashes=[30, 5, 10, 5], c = "green")
plt.show()
In [ ]: