Домашняя работа №5¶
Загрузим необходимые пакеты и создадим z-матрицу этана, которая в дальнейшем будет использоваться как основа для расчета энергии других конформаций этана.
import subprocess
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
base_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
*
'''
Реализуем функции, осуществляющие расчет энергии молекулы по ее z-матрице и генерацию z-матриц с различными длинами C-C связи, различными значениями валентного угла HCC и различными значениями торсионного угла HCCH.
def get_energy(inp):
with open('orca.inp', 'w') as outfile:
outfile.write(inp)
p = subprocess.Popen("orca/orca orca.inp",
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out=p.communicate()[0]
energy_line = filter(lambda x: 'FINAL SINGLE POINT ENERGY' in x, out.splitlines())[0]
return float(energy_line.split()[-1])
def gen_bond_inp(base_inp, bond_length):
base_inp = base_inp.splitlines()
base_inp[3] = base_inp[3].split()
base_inp[3][4] = str(bond_length)
base_inp[3] = ' '.join(base_inp[3])
return '\n'.join(base_inp)
def gen_angle_inp(base_inp, val_angle):
base_inp = base_inp.splitlines()
base_inp[4] = base_inp[4].split()
base_inp[4][5] = str(val_angle)
base_inp[4] = ' '.join(base_inp[4])
return '\n'.join(base_inp)
def gen_tors_inp(base_inp, tors_angle):
base_inp = base_inp.splitlines()
base_inp[7] = base_inp[7].split()
base_inp[7][6] = str(tors_angle)
base_inp[7] = ' '.join(base_inp[7])
return '\n'.join(base_inp)
def plot_energy(points, energies, fitfunc, initial_guess):
x_o=points
y_o=energies
errfunc = lambda p, x, y: fitfunc(p, x) - y
p1, success = optimize.leastsq(errfunc, initial_guess[:], args=(x_o, y_o))
print "Optimized params:", p1
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.show()
Переберем 20 значений длины связи CC по сетке с шагом 0.02 начиная со значения 1.32986. И расчитаем значения энергии этана для каждой длины связи.
start_bond_length = 1.32986
step = 0.02
step_num = 20
bond_length_list = np.linspace(start_bond_length, start_bond_length + step * (step_num - 1), step_num)
energies = [get_energy(gen_bond_inp(base_inp, bond_length)) for bond_length in bond_length_list]
Изобразим получившиеся значения энергии на графике. А также подберем методом наименьших квадратов значения параметров ф-ции $E(r) = k * (r - r_{eq})^2 + C$, при которых она наилучшим образом аппроксимирует полученные точки.
%matplotlib inline
plot_energy(bond_length_list, energies, lambda p, x: p[0]*pow(p[1]-x,2) + p[2], [1., 1., -79.])
Получили, что оптимальными значениями параметров являются $k = 0.64, r_{eq} = 1.55, C = -79.20$.
Теперь расчитаем 50 значений валентного угла HCC от 109.2 до 113.2 и вычислим значения энергий соответствующих конформаций этана.
start_angle = 109.2
stop_angle = 113.2
step_num = 50
val_angle_list = np.linspace(start_angle, stop_angle, step_num)
energies = [get_energy(gen_angle_inp(base_inp, val_angle)) for val_angle in val_angle_list]
Изобразим получившиеся значения энергии на графике. А также подберем методом наименьших квадратов значения параметров ф-ции $E(\theta) = k * (\theta - \theta_{eq})^2 + C$, при которых она наилучшим образом аппроксимирует полученные точки.
%matplotlib inline
plot_energy(val_angle_list, energies, lambda p, x: p[0]*pow(p[1]-x,2) + p[2], [1., 1., -79.])
Получили, что оптимальными значениями параметров являются $k = 4.06e-05, \theta_{eq} = 1.11e+02, C = -79.20$.
Теперь расчитаем значения торсионного угла HCCH от -180 до 180 с шагом 12 и вычислим значения энергий соответствующих конформаций этана.
start_angle = -180.0
stop_angle = 180.0
step = 12.0
tors_angle_list = np.linspace(start_angle, stop_angle, (stop_angle - start_angle) / step + 1)
energies = [get_energy(gen_tors_inp(base_inp, tors_angle)) for tors_angle in tors_angle_list]
Изобразим получившиеся значения энергии на графике. А также подберем методом наименьших квадратов значения параметров ф-ции $E(\phi) = \frac{k}{2} * (1 + cos(n*\phi - \gamma)) + C$, при которых она наилучшим образом аппроксимирует полученные точки.
%matplotlib inline
plot_energy(tors_angle_list, energies, lambda p, x: p[0] / 2 * (1 + np.cos(p[1] * x - p[2])) + p[3], [1., 1., 1., -79.])
Получили, что оптимальными значениями параметров являются $k = 4.58e-03, n = 9.95e-01, \gamma = -1.65e-06, C = -79.20$.
Теперь снова вычислим энергии конформаций этана для 200 различных длин связей CC с шагом 0.02 начиная с 1.
start_bond_length = 1
step = 0.02
step_num = 200
bond_length_list = np.linspace(start_bond_length, start_bond_length + step * (step_num - 1), step_num)
energies = [get_energy(gen_bond_inp(base_inp, bond_length)) for bond_length in bond_length_list]
Попробуем аппроксимировать полученные точки потенциалом Морса: $V(r) = D_e (1-e^{-a(r-r_e)})^2 + C$
%matplotlib inline
plot_energy(bond_length_list, energies, lambda p, x: p[0]*pow(1 - np.exp(p[1] * (p[2] - x)),2) + p[3], [1., 1., 1., -79.])
