Вычисление параметров для молекулярной механики¶
Суть этой работы состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов. Нам предоставлена оптимизированная структура этана в виде 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
*
'''
Наша цель состоит в том, что бы рассчитать 20 разных длинн связи с шагом 0.02 ангстрема для расчёта энергии в ORCA с разными значениями по длине одной из связей (C-C). Создаем функцию для запуска ORCA:
import subprocess, os
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()
for line in lines:
if "FINAL SINGLE POINT ENERGY" in line:
energy=line.split()[4]
return float(energy)
Согласно данному нам файлу длина С-С связи 1.52986 ангстрема. Расчитаем энергию для длины связи от 1.32986 до 1.72986.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
energies=[]
bonds=[]
bond0=1.52986
for i in range(-10,11):
bond=bond0+i*0.02
bonds.append(bond)
#создаем z-матрицу
inp = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 %6.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 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
energies.append(run_orca(inp %bond))
print bonds
print energies
Строим графики зависимости энергии от длины С-С связи:
x_o=np.array(bonds)
y_o=np.array(energies)
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.3,1.75)
plt.ylabel('Energy')
plt.xlabel('C-C bond length, angstrem')
plt.show()
Зависимость энергии от длины С-С связи может быть апроксимирована параболой с минимум, соответствующим длине связи в оптимизированной структуре этана - 1.55699284 ангстрем.
Проделаем аналогичные операции для валентного угла HCС, его значения должны изменяться от 109.2° до 113.2° с шагом 0.2°.
energies_angle=[]
angles=[]
angle0=111.2
for i in range(-10,11):
angle=angle0+i*0.2
angles.append(angle)
#создаем z-матрицу
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 %6.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
*
'''
energies_angle.append(run_orca(inp %angle))
print angles
print energies_angle
x_o=np.array(angles)
y_o=np.array(energies_angle)
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,113.4)
plt.ylabel('Energy')
plt.xlabel('Angle HCC')
plt.show()
Зависимость энергии от величины угла НСС получилось идеально апроксимировать параболой с минимум, соответствующим углу в оптимизированной структуре - 1.11195140e+02°.
Проделываем аналогичные операции для торсионного угла HCСH, его значения должны изменяться от -180° до 180° с шагом 12°.
energies_tors_angle=[]
tors_angles=[]
tors_angle0=0
for i in range(-15,16):
tors_angle=tors_angle0+i*12
tors_angles.append(tors_angle)
#создаем z-матрицу
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 %i
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
energies_tors_angle.append(run_orca(inp %tors_angle))
print tors_angles
print energies_tors_angle
x_o=np.array(tors_angles)
y_o=np.array(energies_tors_angle)
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
p0 = [1,1, -79, 0] # 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(-185,185)
plt.ylabel('Energy')
plt.xlabel('Torsion angle HCCH')
plt.show()
Зависимость энергии от торсионного угла НССН апроксимируется функцией синуса. Функция имеет 4 минимума: -180°, -60°, 60°, 180°, которые соответствуют различным парам водорода в структуре этана.
Увеличим шаг до 0.1 ангстрема при расчёте связи, длина связи будет варьировать в пределах 0.52986 до 2.52986. Построим зависимость.
energies_2=[]
bonds_2=[]
bond0=1.52986
for i in range(-10,11):
bond=bond0+i*0.1
bonds_2.append(bond)
#создаем z-матрицу
inp = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 %6.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 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
energies_2.append(run_orca(inp %bond))
print bonds_2
print energies_2
x_o=np.array(bonds_2)
y_o=np.array(energies_2)
fitfunc = lambda p, x: p[0]*(pow(x,p[2])-pow(x,p[3])) + p[1] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [0.5,-78, -3, -6] # 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(0.5,2.55)
plt.ylabel('Energy')
plt.xlabel('C-C bond length, angstrem')
plt.show()
Зависимость энергии от длины С-С связи с шагом в 0,1 ангстрем можно апроксимировать функцией x^(-1.83377971)-x^(-1.78607697).