Практикум 6: Вычисление ковалентных параметров для молекулярной механики

#подгружаем библиотеки
import subprocess
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import optimize
from IPython.display import Image, display
import os
import sys
from tqdm import tqdm

Исходная оптимизированная структура этана в виде 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
*
'''

Функция, возвращающая матрицу (с переменными длины связи СС, угла НСС и торсионного угла СС).

def mk_inp(l,v,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 {v} 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

Прописываем функцию для расчета энергии молекулы исходя из ее z-матрицы с помощью ORCA.

def run_orca(inp, subfolder):
    with open(os.path.join(subfolder, 'orca.inp'), 'w') as outfile:
        outfile.write(inp)
    with subprocess.Popen("/srv/databases/orca/orca orca.inp", cwd=subfolder,
                          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

Проверка работы программы до варьирования параметров:

print(run_orca(1.52986, 111.200, 180))

-79.197572404432

Сначала попробуем изменить длину связи СС. Варьируем длину связи от 1.35 до 1.75 с шагом 0.02 ангстрем.

x_o = np.arange(1.35,1.75,0.02)
y_o= [run_orca(i,111.200,180) for i in tqdm(x_o)]

#Целевая функция и функция ошибки:
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 params:", p1)

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

100%|██████████| 20/20 [00:39<00:00,  1.98s/it]

Optimized params: [  0.58214253   1.54992786 -79.19771957]

В итоге значение энергии после оптимизации: -79.19771957. Длина связи немного увеличивается по сравнению с параметрами в оптимизированной структуре. Далее рассмотрим изменение валентного кгла (от 109.2 до 113.2 градосов с шагом 0.2 градуса):

x_o = np.arange(109.2, 113.2,0.2)
y_o= [run_orca(1.52986,i,180) for i in tqdm(x_o)]

#Целевая функция и функция ошибки
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 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, 114)
plt.show()

100%|██████████| 20/20 [00:34<00:00,  1.71s/it]

Optimized params: [ 4.06431971e-05  1.11195037e+02 -7.91975724e+01]

Аппроксимация прошла очень хорошо, из графика видно, что кривые прктически идеально совпадают. Минимум энергии остается таким же. Торсионный угол:

x_o = np.arange(-180, 180, 12)
y_o= [run_orca(1.52986,111.200,i) for i in tqdm(x_o)]

#Целевая функция и функция ошибки
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 

p0 = [0.25,10,0, -79] 
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(-190, 190)
plt.show()

100%|██████████| 30/30 [01:00<00:00,  2.03s/it]

Optimized params: [ 2.29231443e-03  1.00007394e+01 -6.15060266e-07 -7.91975767e+01]

На графике минимумы по энергии соответствуют заторможенной конформации, а максимумы - заслоненной.