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

Цель задания: определение констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов.

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

Для того, чтобы преобразовать её в то, что прочитает psi4, посмотрим на пример с водой:
h2o = psi4.geometry("""
O
H 1 0.96
H 1 0.96 2 104.5
""")
Т.е.
O
H 1 длина_связи_O_H
H 1 длина_связи_O_H 2 угол_H_O_H

Таким образом для матрицы этана:

C
C 1 длина_связи_С_С
H 1 длина_связи_H_C 2 валентный_угол_H_C_C
H 1 длина_связи_H_C 2 H_C_C 3 угол_H_C_C_H
H 1 длина_связи_H_C 2 H_C_C 3 угол_H_C_C_H
H 2 длина_связи_H_C 1 H_C_C 3 торсионный_угол
H 2 длина_связи_H_C 1 H_C_C 5 угол_H_C_C_H
H 2 длина_связи_H_C 1 H_C_C 5 угол_H_C_C_H

C
C 1 1.52986
H 1 1.08439 2 111.200
H 1 1.08439 2 111.200 3 120
H 1 1.08439 2 111.200 3 -120
H 2 1.08439 1 111.200 3 0
H 2 1.08439 1 111.200 5 120
H 2 1.08439 1 111.200 5 -120

Необходимо:

  1. Получить 20 z-matrix c разными длиннами связи C-C с шагом 0.02 ангстрема для поиска варианта с минимальной энергией.

  2. Проделать аналогичные операции для валентного угла HCС, его значения должны изменяться от 109.2 до 113.2. Сравнить полученные константы с данными из статьи о GAFF.

  3. Проделать аналогичные операции для торсионного угла CC, его значения должны изменяться от -180 до 180 c шагом 12. Сохранить изображение и указать в отчёте количество минимумов функции.

In [1]:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
from scipy import optimize
from IPython.display import Image, display
In [2]:
import psi4

Ниже предложенный скрипт по запуску функций пакета Psi4, где

  1. psi4.geometry - представляет геометрию согласно вводу
  2. psi4.set_options maxiter - ограничивает макс число итераций
  3. psi4.energy - расчёт одноточечной энергии, в первом аргументенте - метод расчёта (SCF в данном случае) и базисный набор для представления волновых функций (в данном случае cc-pvtz - Correlation-Сonsistent Polarized Valence-only basis set Triple-Zeta)
In [5]:
def run_psi4(inp):
    m = psi4.geometry(inp)
    psi4.set_options({"maxiter": 200, "fail_on_maxiter" :  True})
    return psi4.energy('scf/cc-pvtz', molecule = m )

Расчёт минимум энергии для C-C связи¶

In [15]:
def calculate_energy(bond_length):
    geom = f'''
C 
C 1 {bond_length}
H 1 1.08439 2 111.200
H 1 1.08439 2 111.200 3 120
H 1 1.08439 2 111.200 3 -120
H 2 1.08439 1 111.200 3 0
H 2 1.08439 1 111.200 5 120
H 2 1.08439 1 111.200 5 -120
'''
    return run_psi4(geom)
In [16]:
print(calculate_energy(1))

-78.79588773541637
In [9]:
x_o=np.arange(1.3,1.7,0.02)
y_o = np.array(list(map(calculate_energy, x_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),"r-",c='blue',alpha=0.5)
plt.xlim(1,2)
plt.show()

Optimized params: [  0.75232916   1.55165175 -79.25655593]

C:\Users\днс\AppData\Local\Temp\ipykernel_6972\3011160597.py:13: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "ro" (-> color='r'). The keyword argument will take precedence.
  plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
C:\Users\днс\AppData\Local\Temp\ipykernel_6972\3011160597.py:13: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "r-" (-> color='r'). The keyword argument will take precedence.
  plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)

Оптимальная длина связи - 1.5517. В википедии - 1.5351, в статье GAFF - 1.526. Ну, немного больше.

Расчёт значения валентного угла¶

In [17]:
def calculate_energy(HCC_angle):
    geom = f'''
C 
C 1 1.52986
H 1 1.08439 2 {HCC_angle}
H 1 1.08439 2 111.200 3 120
H 1 1.08439 2 111.200 3 -120
H 2 1.08439 1 111.200 3 0
H 2 1.08439 1 111.200 5 120
H 2 1.08439 1 111.200 5 -120
'''
    return run_psi4(geom)
In [13]:
#Википедия говорит, что длина угла - 111.17, так что initial quess - 111
print(calculate_energy(111))

-78.85730475547705
In [22]:
x_o = np.arange(109.2,113.2,0.2)
y_o = np.array(list(map(calculate_energy, x_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, 111, -78] # 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),"r-",c='blue',alpha=0.5)
plt.xlim(108,114)
plt.show()

Optimized params: [ 3.93927516e-05  1.11938788e+02 -7.92550798e+01]

C:\Users\днс\AppData\Local\Temp\ipykernel_6972\2182522502.py:13: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "ro" (-> color='r'). The keyword argument will take precedence.
  plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
C:\Users\днс\AppData\Local\Temp\ipykernel_6972\2182522502.py:13: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "r-" (-> color='r'). The keyword argument will take precedence.
  plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)

Оптимальный валентный угол - 111.9388, на Википедии - 111.17, в статье GAFF значений валентных углов нет.

Расчёт значения торсионного угла¶

In [24]:
def calculate_energy_torsion(HCCH_angle):
    adjust = lambda x: (x + 180 + 360) % 360 - 180 if x != 180 else 180
    geom = f'''
C 
C 1 1.52986
H 1 1.08439 2 111.200
H 1 1.08439 2 111.200 3 120
H 1 1.08439 2 111.200 3 -120
H 2 1.08439 1 111.200 3 {HCCH_angle}
H 2 1.08439 1 111.200 3 {adjust(HCCH_angle + 120)}
H 2 1.08439 1 111.200 3 {adjust(HCCH_angle - 120)}
'''
    return run_psi4(geom)
x_torsion = np.arange(-180, 181, 12)
y_torsion = np.array(list(map(calculate_energy_torsion, x_torsion)))


#Plot it
plt.plot(x_torsion, y_torsion, "ro")
plt.show()

Четыре минимума (-180, -60, 60, 180), c периодичностью 120, что соответсвует заторможенной конформации.

In [25]:
from IPython import display
display.Image("https://cdn.masterorganicchemistry.com/wp-content/uploads/2020/02/15-energy-diagram-for-ethane-torsional-strain-versus-dihedral-angle.gif")
Out[25]:
<IPython.core.display.Image object>