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

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

Задание: https://kodomo.fbb.msu.ru/wiki/2016/8/MolSim/task6

Скачать ORCA локально не получилось, поэтому пробуем запустить jupyter notebook на kodomo

ssh -L 9200:localhost:8888 potanina.darya@kodomo.fbb.msu.ru
jupyter notebook --ip=0.0.0.0 --no-browser
// The Jupyter Notebook is running at....
// Открыть в браузере http://localhost:9200/, ввести токен в качестве пароля.
// ура! jupyter notebook

Ставим необходимые библиотеки:

jupyter kernelspec list
python3.8 -m pip install ipykernel --user
python3.8 -m ipykernel install --user

pip3.8 install matplotlib --user
pip3.8 install scipy --user
pip3.8 install numpy --user
In [1]:
!which python

/usr/bin/python
In [2]:
!/srv/databases/orca/orca

This program requires the name of a parameterfile as argument
For example ORCA TEST.INP
In [10]:
import subprocess
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import optimize
import matplotlib as mpl

from IPython.display import Image, display
In [4]:
def run_orca(inp):
    with open('orca.inp', 'w') as outfile:
            outfile.write(inp)
    p = subprocess.Popen("/srv/databases/orca/orca orca.inp", shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out = p.communicate()[0]
    lines = out.splitlines()
    energy = None
    for line in lines:
        line = str(line, 'utf-8')
        if line.find("FINAL SINGLE POINT") > -1:
            energy = float(line.split()[4])
    return energy
In [5]:
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
*
'''
In [6]:
run_orca(inp)
Out[6]:
-79.081531429795

Подготовим шаблоны файлов

In [7]:
bond_length_template = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 {:.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 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
angle_template = '''!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 {:.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
*'''
torsion_template = '''!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 {:d}
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*'''

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

In [8]:
bond_lengths = np.arange(1.3, 1.7, 0.02)
bond_length_energies = np.array([run_orca(bond_length_template.format(value))
                                 for value in bond_lengths])

def length_energy(length, k, l0, a):
    return k / 2 * (length - l0) ** 2 + a


popt_l, pcov = optimize.curve_fit(length_energy, bond_lengths, bond_length_energies)


popt_l
Out[8]:
array([  1.54918147,   1.54783146, -79.08299495])
In [12]:
mpl.style.use('seaborn')
plt.plot(bond_lengths, bond_length_energies, 'o', label='orca')
plt.plot(bond_lengths, length_energy(bond_lengths, *popt_l), label='fitted')
plt.legend()
plt.show()

Не идеально, попробуем улучшить.

In [13]:
def morze_potential(length, d, a, l0, b):
    return d * (1 - np.exp(- a * (length - l0))) ** 2 + b
In [14]:
popt_m, pcov = optimize.curve_fit(morze_potential, bond_lengths, bond_length_energies)
In [15]:
popt_m
Out[15]:
array([  0.19079564,   1.74329991,   1.53584295, -79.0815584 ])
In [16]:
plt.plot(bond_lengths, bond_length_energies, 'o', label='orca')
plt.plot(bond_lengths, morze_potential(bond_lengths, *popt_m), label='fitted')
plt.legend()
plt.show()

Отлично. Со свободным членом фитится сильно лучше.

In [17]:
from functools import partial
In [18]:
morze_for_optimization = partial(morze_potential, **dict(zip(['d', 'a', 'l0', 'b'], popt_m)))
In [19]:
optimization_results = optimize.minimize_scalar(morze_for_optimization, bounds=(1.3, 1.6))
In [20]:
optimization_results.x
Out[20]:
1.5358429073226942

По рассчетам Orca связь получилась длиной около 1.54 ангстрем. В статье о разработке поля GAFF указано, что они выбрали длину C-C связи равной 1.526. Не знаю, насколько допустима такая точность, если они правы6 то наша ошибка может оказаться критичной в моделировании, но в целом это неплохой результат.

Расчёт минимума энергии для валентного угла HCC

In [33]:
angle_values = np.arange(109.2, 113.21, 0.2)
In [45]:
angle_energies = np.array([run_orca(angle_template.format(value))
                           for value in angle_values])
In [46]:
def angle_potential(x, k, f0, a):
    return k / 2 * (x - f0) ** 2 + a
In [ ]:
popt_a, pcov = optimize.curve_fit(angle_potential, angle_values, angle_energies)
In [ ]:
plt.plot(angle_values, angle_potential(angle_values, *popt_a), label='fitted')
plt.legend()
plt.show()
In [ ]:
angle_for_optimization = partial(angle_potential, **dict(zip(['k', 'f0', 'a'], popt_a)))
In [ ]:
optimization_results = optimize.minimize_scalar(angle_for_optimization, bounds=(110, 112))
In [ ]:
optimization_results.x

Вот. У меня это не заработало, я просто сохраню код. Ошибка была в том, что функция optimize.curve_fit не принимает на вход значения с NaN, но таких вроже не должно возникать, и во всяком случае у меня не получилось это поправить. Так что я не знаю, какой у меня минимум энергии для валентного угла HCC.

Рассчёт минимума энергии для торсионного угла CC

In [50]:
torsion_values = np.arange(-180, 181, 12)
In [51]:
torsion_energies = np.array([run_orca(torsion_template.format(value))
                             for value in torsion_values])
In [52]:
plt.plot(torsion_values, torsion_energies, 'o')
plt.xticks(np.arange(-180, 181, 60))
plt.show()
In [53]:
def torsion_potential(x, v, n, g, a):
    return v / 2 * (1 + np.cos(np.radians(n * np.pi * x - g))) + a
In [54]:
popt_t, pcov = optimize.curve_fit(torsion_potential, torsion_values, torsion_energies)
In [55]:
popt_t
Out[55]:
array([ 5.19916891e-03,  9.55075462e-01, -8.58960511e-03, -7.95705884e+01])
In [56]:
plt.plot(torsion_values, torsion_energies, 'o', label='orca')
plt.plot(torsion_values, torsion_potential(torsion_values, *popt_t), label='fitted')
plt.xticks(np.arange(-180, 181, 60))
plt.legend()
plt.show()

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

In [57]:
plt.plot(np.arange(-180, 180, 1), torsion_potential(np.arange(-180, 180, 1), *popt_t))
Out[57]:
[<matplotlib.lines.Line2D at 0x7fab87577250>]
In [58]:
torsion_for_optimization = partial(torsion_potential, **dict(zip(['v', 'n', 'g', 'a'], popt_t)))
In [59]:
optimization_results = optimize.minimize_scalar(torsion_for_optimization, bounds=(0, 120))
In [60]:
optimization_results.x
Out[60]:
59.98797623551164

Ага, отлично, получили значение очень близкое к 60.

Рассчёт минимума энергии для длины связи С-С с бОльшим шагом

In [63]:
long_bond_lengths = np.arange(1, 2.6, 0.1)
In [64]:
long_bond_length_energies = np.array([run_orca(bond_length_template.format(value))
                                      for value in long_bond_lengths])
In [65]:
plt.plot(long_bond_lengths, long_bond_length_energies)
Out[65]:
[<matplotlib.lines.Line2D at 0x7fab876091f0>]
In [66]:
popt_long, pcov = optimize.curve_fit(morze_potential, long_bond_lengths, long_bond_length_energies)
In [67]:
popt_long
Out[67]:
array([  0.23994791,   1.65276025,   1.53227023, -79.08500183])
In [68]:
plt.plot(long_bond_lengths, long_bond_length_energies, 'o', label='orca')
plt.plot(long_bond_lengths, morze_potential(long_bond_lengths, *popt_long), label='fitted')
plt.legend()
plt.show()

О, оно отлично фитится. Посмотрим, какое значение получили в этот раз.

In [69]:
morze_long_for_optimization = partial(morze_potential, **dict(zip(['d', 'a', 'l0', 'b'], popt_long)))
In [70]:
optimization_results = optimize.minimize_scalar(morze_long_for_optimization, bounds=(1.4, 1.6))
In [71]:
optimization_results.x
Out[71]:
1.5322701552817175

Хм, ну это опять не 1.526, но ближе, чем в прошлый раз.

Итог: наши полученные значения похожи на те, что получены в статье, но все-таки отличаются в 1-2 знаках, что может оказаться критичным при моделировании. Не знаю, как определить, чьи параметры лучше, возможно, они использовали другие исходные данные для рассчета ab initio.

In [ ]: