Загрузка пакетов
In [5]:
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import optimize
import subprocess
import numpy as np
import math
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit import RDConfig
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from IPython.display import display,Image
Задача -> Определение констант ковалентных взаимодействий для молекулярной механики этана на основе квантово-химических расчётов
In [7]:
etan=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(etan)
display(etan)
In [ ]:
# Имеем на входе оптимизированную структуру этана:
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
*
'''
# Из кода на "птичьем" языке видим, что расчет оптимальной структуры был сделан так,
# что для внутренних электронов использовались 6 гауссианов с фиксированными коэффициентами,
# для "внешних", валентных, электронов, участвующих в ковалентных взаимодействиях, коэффициенты в 3 гауссианах фиксированы,
# в 1 гауссиане - вариабельные и оптимизировались.
# В Z-matrix:
# 5 колонка - длина связи
# 6 колонка - валентный угол
# 7 колонка - торсионный угол
Функция для генерации файлов для ORCA (.inp):
In [36]:
def inp_len(length):
inp = '''!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 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''.format(length)
return inp
Функция по запуску ORCA:
In [10]:
def run_orca(inp):
with open('orca.inp', 'w') as outfile:
outfile.write(inp)
p = subprocess.Popen("/home/shad/progs/bin/orca orca.inp",
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out=p.communicate()[0]
# extract energy: FINAL SINGLE POINT ENERGY'
for i in out.splitlines():
if i.startswith("FINAL SINGLE POINT ENERGY"):
return(float(i[len("FINAL SINGLE POINT ENERGY"):]))
# and return it as float
return out.splitlines()
Зависимость энергии молекулы от длины одной связи:
In [41]:
def energy_dependency(x,y,title,lim1,lim2):
# fake x array , replace with real one
x_o=x
# fake y array, replace with energies
y_o=y
#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(lim1,lim2)
plt.title(title)
plt.show()
Формируем массив 20 различных длинн связи с шагом в 0.02 ангстрем
In [34]:
link_lengths = np.arange(-10,11)*0.02 + 1.52986
print (link_lengths)
In [37]:
energy_len = list()
for i in link_lengths:
energy_len.append(run_orca(inp_len(i)))
print(energy_len)
Рисуем график зависимости энергии молекулы от длины связи C-C
In [42]:
energy_dependency(link_lengths,energy_len, 'Link length and energy',1.3,1.8)
Формируем массив 20 различных валентных углов в пределах от 109.2 до 113.2
In [39]:
def inp_angle(angle):
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 {0} 0
H 1 2 3 1.08439 {0} 120
H 1 2 3 1.08439 {0} -120
H 2 1 3 1.08439 {0} 180
H 2 1 5 1.08439 {0} 120
H 2 1 5 1.08439 {0} -120
*
'''.format(angle)
return inp
In [26]:
angles = np.arange(109.2,113.2,(113.2-109.2)/20)
print(angles)
In [27]:
len(angles)
Out[27]:
In [29]:
energy_angles = list()
for i in angles:
energy_angles.append(run_orca(inp_angle(i)))
print(energy_angles)
Рисуем график зависимости энергии молекулы от валентного угла
In [43]:
energy_dependency(angles,energy_angles, 'valence angles and energy',109,113.4)
Отличная аппроксимация, как и ожидалось (из лекции).
Формируем массив различных торсионных углов в пределах от -180 до 180 с шагом 12
In [65]:
def inp_torsion(torsion):
torsion1 = torsion - 60.
torsion2 = torsion + 60.
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 {0}
H 2 1 5 1.08439 111.200 {1}
H 2 1 5 1.08439 111.200 {2}
*
'''.format(torsion, torsion1, torsion2)
return inp
In [54]:
torsions = np.arange(-180,192,12)
print(torsions)
In [55]:
len(torsions)
Out[55]:
In [66]:
energy_torsions = []
for i in torsions:
energy_torsions.append(run_orca(inp_torsion(i)))
print(energy_torsions)
Рисуем график зависимости энергии молекулы от валентного угла
In [67]:
energy_dependency(torsions,energy_torsions, 'torsions and energy',-180,180)
Зафитим функцию через Фурье
In [69]:
lim1 = -180
lim2 = 180
title = 'torsions and energy'
# fake x array , replace with real one
x_o=torsions
# fake y array, replace with energies
y_o=energy_torsions
#function is f(x)=k(b-x)^2 + a
fitfunc = lambda p, x: p[0]*np.sin(p[1]*x)+p[3]*np.cos(p[4]*x) + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1,1,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(lim1,lim2)
plt.title(title)
plt.show()
Увеличиваем шаг до 0.1 ангстрема
In [72]:
link_lengths_2 = np.arange(-10,10)*0.1 + 1.52986
print (link_lengths_2)
In [73]:
len(link_lengths_2)
Out[73]:
In [74]:
energy_len_2 = list()
for i in link_lengths_2:
energy_len_2.append(run_orca(inp_len(i)))
print(energy_len_2)
In [75]:
energy_dependency(link_lengths_2,energy_len_2, 'New link length and energy',0.52,2.43)
Зафитим гиперболой:
In [76]:
lim1 = 0.52
lim2 = 2.43
title = 'New link length and energy'
# fake x array , replace with real one
x_o=link_lengths_2
# fake y array, replace with energies
y_o=energy_len_2
fitfunc = lambda p, x: pow(1/x,p[1])+p[0] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [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(lim1,lim2)
plt.title(title)
plt.show()
