Вычисление параметров для молекулярной механики¶
0. Импорты и константы¶
from textwrap import dedent
from fabric import Connection
import keyring
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
import seaborn
from jupyterthemes import jtplot
jtplot.style(figsize=(13, 11), fscale=1.5)
BOND, ANGLE, DIHEDRAL, LONGBOND = 0, 1, 2, 3
1. Материалы и методы¶
1) Генераторы входных файлов для орки¶
def inp_bond(bond_len):
return dedent(f'''
!MP2 6-31G*
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 {bond_len:.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 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*''').strip()
def inp_ang(angle):
return dedent(f'''
!MP2 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 {angle:.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 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*''').strip()
def inp_dih(dih):
return dedent(f'''
!MP2 6-311G**
* 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 {dih:d}
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*''').strip()
2) Запускатор орки¶
def run_orca(conn, inp):
with conn.cd('~/tmp/term8'):
conn.run(f'echo "{inp}" > orca.inp', hide='stdout')
energy = conn.run('/srv/databases/orca/orca orca.inp | grep -E "FINAL .*? ENERGY"',
hide='stdout').stdout.split()[-1]
return float(energy)
3) Удобный интерфейс для расчета энергий (и их сохранения)¶
def calc_energies(mode, values):
get_inp = {BOND: inp_bond,
ANGLE: inp_ang,
DIHEDRAL: inp_dih,
LONGBOND: inp_bond
}[mode]
energies = np.zeros(len(values))
with Connection('kodomo.fbb.msu.ru', user='aidar',
connect_kwargs={'password': keyring.get_password('kodomo', 'aidar')}) as kodomo:
for i, value in enumerate(values):
energy = run_orca(kodomo, get_inp(value))
energies[i] = energy
energies *= 627.5 # Convert Hartree to kcal/mol
file_name = {BOND: 'bond.txt',
ANGLE: 'angle.txt',
DIHEDRAL: 'dih.txt',
LONGBOND: 'longbond.txt'
}[mode]
with open(file_name, 'w') as f:
f.write('value\tenergy(kcal/mol)\n')
for v, e in zip(values, energies):
f.write(f'{v}\t{e:.3f}\n')
return energies
4) Еще более удобный интерфейс для оптимизации функций и построения графика¶
def fit_plt(mode, values):
energies = calc_energies(mode, values)
fit_func = {BOND: lambda x, k, r, a: k * (x - r)**2 + a,
ANGLE: lambda x, k, r, a: k * (x - r)**2 + a,
DIHEDRAL: lambda x, k, r, a: k * (1 + np.cos(np.radians(3*x - r))) + a,
LONGBOND: lambda x, k, r, a: k * (np.exp(-2 * 1.836 * (x - r))
- 2 * np.exp(-1.836 * (x - r))) + a
}[mode]
(k, r, a), _ = optimize.curve_fit(fit_func, values, energies)
jtplot.style(figsize=(14, 12), fscale=1.4)
plt.plot(values, energies, label='orca', marker='o')
newvalues = np.linspace(min(values), max(values), 100)
func_label = {BOND: f'y = {k:.3f} * (x - {r:.3f})^2 - {-a:.3f}',
ANGLE: f'y = {k:.3f} * (x - {r:.3f})^2 - {-a:.3f}',
DIHEDRAL: f'y = {k:.3f} * (1 + cos(3*x - {r:.3f})) - {-a:.3f}',
LONGBOND: f'y = {k:.3f} * (e^(-2 * 1.836 * (x - {r:.3f}))'
f' - 2 * e^(-1.836 * (x - {r:.3f}))) - {-a:.3f}'
}[mode]
plt.plot(newvalues, fit_func(newvalues, k, r, a), label=func_label)
xlabel = {BOND: 'Bond length (Angstroms)',
ANGLE: 'Angle (degrees)',
DIHEDRAL: 'Dihedral angle (degrees)',
LONGBOND: 'Bond length (Angstroms)'
}[mode]
plt.xlabel(xlabel)
plt.ylabel('Energy (kcal/mol)')
plt.legend(fontsize=16)
plt.ticklabel_format(useOffset=False)
2. Результаты¶
1) Параметры для связи C-C¶
bond_lens = np.arange(1.52986 - 0.2, 1.52986 + 0.21, 0.02)
fit_plt(BOND, bond_lens)
Файл с энергиями
Как видно из графика, модель пружины плохо описывает энергию связи. У нас получились параметры 1.55 и 373.8 для оптимальной длины и жесткости связи соответственно. В статье указаны 1.526 и 2086.
2) Параметры для угла H-C-C¶
angles = np.arange(109.2, 113.3, 0.2)
fit_plt(ANGLE, angles)
Файл с энергиями
А вот для угла модель пружины работает на удивление хорошо. У нас получились параметры 111 и 0.023 для оптимального угла и жесткости угла соответственно. В статье указана только жесткость -- 0.015.
3) Параметры для торсионного угла H-C-C-H¶
dihedrals = np.arange(-180, 181, 12)
fit_plt(DIHEDRAL, dihedrals)
Файл с энергиями
Для торсионного угла косинус великолепно описывает энергию. У нас получились параметры 0 и 1.63 для максимума энергии и жесткости угла соответственно. В статье указаны 0 и 1.2.
4) Улучшенные параметры для связи C-C¶
longbond_lens = np.arange(1.52986 - 0.7, 1.52986 + 1, 0.05)
fit_plt(LONGBOND, longbond_lens)
Файл с энергиями
Для аппроксимации энергии связи очень хорошо подошел потенциал Морзе. У нас получились параметры 1.5 и 121 для оптимальной длины и "жесткости" связи.
3. Обсуждения¶
Наши параметры в большинстве случаев были очень похожи на статейные (за исключением жесткости связи, но я вообще не представляю почему там так получилось). А небольшие различия можно объяснить следующими факторами:
- очевидно, что мы брали разные точки (и разные интервалы) для фиттинга
- в статье использовали другие базисы и методы расчета (MP2/6-31G* для связей и угла; MP4/6-311G** для торсионного угла)
- в статье делалил непонятную поправку к значениям, полученным ab initio
- для получения финальных констант они так же использовали данные из поля AMBER и из структур кристаллов.
Также стоит отметить, что пружина плохо описывает ковалентную связь, а потенциал Морзе хорошо, но видимо их хорошо моделировать не надо (нам же скорее нековалентные взаимодействия важны), поэтому и пружинка сойдет (давая сильный прирост к скорости расчета).