Практикум 6

import subprocess
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import optimize
from IPython.display import Image, display

Определяеи функцию расчета энергии при помощи программы ORCA

def run_orca(inp):
    with open('orca.inp', 'w') as outfile:
        outfile.write(inp)
    with subprocess.Popen("orca orca.inp", shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) as p:
        out=p.communicate()[0].decode()
    for l in out.splitlines():
        if 'FINAL SINGLE POINT' in l:
            return float(l.strip().split()[4])
    return None

Предоставленная оптимизированная структура этана в виде z-matrix

def make_input(l,a,t):
    inp = f'''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 {l} 0 0 
H 1 2 0 1.08439 {a} 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 {t}
H 2 1 3 1.08439 111.200 {t-120}
H 2 1 3 1.08439 111.200 {t-240}
*
'''
    return inp

Теперь будем менять длину связи 1.2 до 1.5 с шагом 0.02.

x_o_b=np.arange(1.2,1.5,0.02)
y_o_b=[run_orca(make_input(l,111.200,180)) for l in x_o_b]

print(x_o_b)
print(y_o_b)

Построим зависимость энергии молекулы от длины одной связи

fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] # Целевая функция f(x)=k(b-x)^2 + a
errfunc = lambda p, x, y: fitfunc(p, x) - y # Функция ошибки, ее мы и минимизируем

p0 = [1,1, -79] # Наши начальные параметры
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o_b, y_o_b)) # Сама функция оптимизации, получаем из нее подобранные параметры
print("Optimized params:", p1)

plt.plot(x_o_b, y_o_b, "ro", x_o_b,fitfunc(p1,x_o_b),"r-",c='blue',alpha=0.5)
plt.xlim(1.2,1.8)
plt.show()

from IPython.display import Image
display(Image('/home/masha/Документы/Головин/прак6/im_1.png'))

Мы нашли минимум энергии и другие константы. Однако наша зависимость является не очень хорошим приближением.

Теперь будем изменять угол НСС от 109.2 до 113.2 с шагом 0.2

x_o_a=np.arange(109.2,113.2,0.2)
y_o_a=[run_orca(make_input(1.52986,a,180)) for a in x_o_a]

print(x_o_a)
print(y_o_a)

Построим зависимость энергии

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_a, y_o_a))
print("Optimized params:", p1)

plt.plot(x_o_a, y_o_a, "ro", x_o_a, fitfunc(p1,x_o_a), "r-", c='blue', alpha=0.5)
plt.xlim(109,114)
plt.show()

from IPython.display import Image
display(Image('/home/masha/Документы/Головин/прак6/im_2.png'))

Минимум по энергии тот же (-79.197), и при том же значении угла (111.2).

аналогичные операции для торсионного угла CC, его значения должны изменяться от -180 до 180 c шагом 12.

x_o_t=np.arange(-180,180,12)
y_o_t=[run_orca(make_input(1.52986,111.200,t)) for t in x_o_t]

print(x_o_t)
print(y_o_t)

fitfunc = lambda p, x: p[0]*(1 + np.cos(p[1]*x - p[2]))+p[3] # 
errfunc = lambda p, x, y: fitfunc(p, x) - y # Тут целевая функция задается как f(x)=k*(1 + cos(a*x - b)) + c

p0 = [0.25,10,0,-79]
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o_t, y_o_t))
print("Optimized params:", p1)

plt.plot(x_o_t, y_o_t, "ro", x_o_t, fitfunc(p1,x_o_t),"r-", c='blue', alpha=0.5)
plt.xlim(-200,200)
plt.show()

from IPython.display import Image
display(Image('/home/masha/Документы/Головин/прак6/im_3.png'))

Мы видим, что у нас минимумы энергии. Это связано с тем,что энергия не заслоненной конформации ниже, чем заслоненной. Минимумы расположены с периодом равным 120 градусам, а при сдвиге на 60 градусов мы видим максимум по энергии - заслоненную конформацию