Вычисление параметров для молекулярной механики
Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов.
Рассчёт зависимости энергии от длины (C-C)-связи
In [ ]:
##исходные данные
'''!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
*
'''
In [1]:
##функция, которая запускает ORCA и возвращает рассчитанную энергию
import subprocess
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]
out = out.splitlines()
for x in out:
if "FINAL" in x:
return(float(x.split()[-1])*1.04*10**(-2))
elif "Error:" in x:
print inp
print x
print "~~~"
In [2]:
##шаблон Z-matrix
inp='''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 %.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
*
'''
In [3]:
##списки со значениями длины, с помощью функции можно получить для каждого из них энергию в кДж/моль
c_bond = 1.52986
bonds = [round(c_bond + x*0.02,5) for x in range(-10,11)]
energies = [run_orca(inp%b) for b in bonds]
In [4]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [8]:
x_o=np.array(bonds)
y_o=np.array(energies)
print "Bond energy:"
print y_o
#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),"b-",alpha=0.5)
plt.ylabel('Energy, kJ/mol')
plt.xlabel('Bond length, angstrem')
plt.ylim(-0.8225,-0.822)
plt.xlim(1.3,1.8)
plt.show()
plt.savefig("pr6")
##Получен график, отражающий зависимость значения энергии от длины С-С связи и аппроксимированный параболой (синия линия).
Рассчёт зависимости энергии от значения валентного угла HCС
In [10]:
##данные
inp_2='''!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 %.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 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
In [11]:
##списки со значениями углов, для каждого энергию в кДж/моль.
angles = list(np.arange(109.2,113.3,0.2))
ang_energies = [run_orca(inp_2%a) for a in angles]
In [12]:
x_1=np.array(angles)
y_1=np.array(ang_energies)
print "Angle energy:"
print y_1
#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_1, y_1))
print "Optimized params:", p1
#Plot it
plt.plot(x_1, y_1, "ro", x_1, fitfunc(p1,x_1),"b-",alpha=0.5)
plt.ylabel('Energy, kJ/mol')
plt.xlabel('Bond angle, deg')
plt.ylim(-0.822445,-0.8224485)
plt.xlim(109,113.5)
plt.show()
##Получили график с зависимостью энергии от значения валентного угла HCC и аппроксимированный параболой (синия линия).
Рассчёт зависимости энергии от значения торсионного угла CC
In [18]:
inp_3='''!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 %d
H 1 2 3 1.08439 111.200 %d
H 2 1 3 1.08439 111.200 %d
H 2 1 6 1.08439 111.200 %d
H 2 1 6 1.08439 111.200 %d
*
'''
In [19]:
tors_angles = [180 + x*12 for x in range(-30,1)]
print len(tors_angles)
tors_energies = [run_orca(inp_3%(120,-120,t,120,-1*120)) for t in tors_angles]
In [20]:
x_2=np.array(tors_angles)
y_2=np.array(tors_energies)
print "Angle energy:"
print y_2
#function is f(x)=cos(x * freq + phase) * A+ b
fitfunc = lambda p, x: np.cos(x * p[0] + p[2]) * p[1] + p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
guess_freq = 1
guess_A = 3*np.std(y_2)/(2**0.5)
guess_phase = 0
guess_b = np.mean(y_2)
p0=[guess_freq, guess_A, guess_phase, guess_b] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_2, y_2))
print "Optimized params:", p1
#Plot it
plt.plot(x_2, y_2, "ro", x_2, fitfunc(p1,x_2),"b-",alpha=0.5)
plt.ylabel('Energy, kJ/mol')
plt.xlabel('Bond angle, deg')
plt.ylim(-0.8236,-0.82366)
plt.xlim(-190,190)
plt.show()
##Получена зависимость энергии от величины торсионного угла - синусоидальная.
In [ ]: