Вычисление параметров для молекулярной механики¶
Данная работа заключалась в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов. Это производилось на примере молекулы этана с использованием ORCA.
Имеем оптимизированную структура этана в виде z-matrix (bond_length,bond angle, torsion 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 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
*
Далее с помощью ORCA были расчитаны энергии при разных значениях длины связи (C-C). Были подставлены 20 разных длинн связи (C-C) с шагом 0.02 ангстрема от оптимального значения.
import subprocess, os
def run_orca(inp):
with open('orca.inp', 'w') as outfile:
outfile.write(inp)
my_env = os.environ.copy()
p = subprocess.Popen("/srv/databases/orca/orca orca.inp",
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=my_env)
p.wait()
out=p.communicate()[0]
lines=out.splitlines()
for line in lines:
if "FINAL SINGLE POINT ENERGY" in line:
line=line.split()[4]
return float(line)
ccbond=1.52986
bonds=[]
b_energies=[]
for n in range(-10,10):
newbond=ccbond+n*0.02
bonds.append(newbond)
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 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
b_energy=run_orca(inp %newbond)
b_energies.append(b_energy)
print "bond length:"
print bonds
print "bond energy:"
print b_energies
С помощью matplotlib построена зависимость энергии молекулы от длинны С-С связи в молекуле этана.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
x_o=bonds
y_o=b_energies
#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 parameters:", p1
#Plot it
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.ylabel('Energy, eV')
plt.xlabel('C-C bond length, angstrem')
plt.xlim(1.3,1.75)
plt.show()
График зависимости в приближении имеет вид параболы. Минимуму энергии соответствует длина С-С связи равная 1.52986 Å, именно такое значение эта связи имеет в оптимизированной структуре этана. Также можно заметить, что более маленькие значения длины связи имеют большее влияние на энергию, чем симмитричные значения связей в большую сторону от значения 1.52986 Å.
Аналогичные операции были проведены для валентного угла HCС, его значения изменятся от 109.2° до 113.2°.
hccangle=111.200
angles=[]
a_energies=[]
for n in range(-10,10):
newangle=hccangle+n*0.2
angles.append(newangle)
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
*
'''
a_energy=run_orca(inp %newangle)
a_energies.append(a_energy)
print "hcc angle:"
print angles
print "angle energy:"
print a_energies
x_o=angles
y_o=a_energies
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_o, y_o))
print "Optimized parameters:", p1
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.ylabel('Energy, eV')
plt.xlabel('Angle, degree')
plt.xlim(109,113.4)
plt.show()
График зависимости энергии молекулы этана от значения HCC угла аналогично с предыдущим случаем приближается параболой. Минимум энергии соответствует значению HCC угла 111.2°, именно такое значение этот угол имеет в оптимизированной структуре этана.
Аналогичные операции были проделаны для торсионного угла CC, его значения изменятся от -180° до 180° c шагом 12.
tors_angle=180
tors_angles=[]
t_energies=[]
for n in range(-30,1):
newtors=tors_angle+n*12
tors_angles.append(newtors)
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 %f
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
t_energy=run_orca(inp %newtors)
t_energies.append(t_energy)
print "torsion angle:"
print tors_angles
print "torsion angle energy:"
print t_energies
x_o=tors_angles
y_o=t_energies
plt.plot(x_o, y_o, "ro", c='blue',alpha=0.5)
plt.ylabel('Energy, eV')
plt.xlabel('Torsion angle, degree')
plt.xlim(-180,180)
plt.show()
plt.savefig("energy_torsangle")
Этот график зависимости энергии от значений торсионного угла в молекуле этана имеет синусоидальную форму. Он имеет 4 минимума: при значениях -180, -60, 60 и 180. Это соответствует торсионным углам H-C-C-H с различными парами водородов в оптимизированной молекуле этана.