Задание 6

Проверяем, всё ли работает, и загружаем нужные нам модули.

!orca

This program requires the name of a parameterfile as argument
For example ORCA TEST.INP 

import subprocess
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import optimize

from IPython.display import Image, display

with subprocess.Popen('orca',shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) as p:
    print(p.communicate())

(b'This program requires the name of a parameterfile as argument\nFor example ORCA TEST.INP \n', b'')

Всё работает, так как экспортировали путь до того, как завели сервер.

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

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 line in out.splitlines():
        if 'FINAL SINGLE POINT' in line:
            return float(line.strip().split()[4])
    return None

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
*
'''

run_orca(inp)

-79.081531432268

Ура, работает!

Посмотрим на некоторые параметры молекулы этана в википедии.

display(Image('https://upload.wikimedia.org/wikipedia/commons/6/68/Ethane-staggered-CRC-MW-dimensions-2D.png', width=720))

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

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 связи

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

array([  1.5491822 ,   1.54783142, -79.08299495])

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()

Не очень хорошо фитируется, увы.

def morze_potential(length, d, a, l0, b):
    return d * (1 - np.exp(- a * (length - l0))) ** 2 + b

popt_m, pcov = optimize.curve_fit(morze_potential, bond_lengths, bond_length_energies)

popt_m

array([  0.19079672,   1.74329544,   1.53584299, -79.0815584 ])

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()

Прекрасно!

Почему возникает свободный член? У нас минимум функции не в 0, а меньше нуля, поэтому нам нужен свободный член для подгона.

from functools import partial

morze_for_optimization = partial(morze_potential, **dict(zip(['d', 'a', 'l0', 'b'], popt_m)))

optimization_results = optimize.minimize_scalar(morze_for_optimization, bounds=(1.3, 1.6))

optimization_results.x

1.535842940095419

В статье по разработке GAFF указано, что они выбрали длину связи C-C, равную 1.526. У нас, в случае с orca, слабо сходится. Зато очень похоже на данные википедии.

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

angle_values = np.arange(109.2, 113.21, 0.2)

angle_energies = np.array([run_orca(angle_template.format(value))
                           for value in angle_values])

def angle_potential(x, k, f0, a):
    return k / 2 * (x - f0) ** 2 + a

popt_a, pcov = optimize.curve_fit(angle_potential, angle_values, angle_energies)

popt_a

array([ 7.23517766e-05,  1.11037057e+02, -7.94918725e+01])

plt.plot(angle_values, angle_energies, 'o', label='orca')
plt.plot(angle_values, angle_potential(angle_values, *popt_a), label='fitted')
plt.legend()
plt.show()

Очень хорошо описали данные.

angle_for_optimization = partial(angle_potential, **dict(zip(['k', 'f0', 'a'], popt_a)))

optimization_results = optimize.minimize_scalar(angle_for_optimization, bounds=(110, 112))

optimization_results.x

111.03707096992738

У нас минимум в 111.037, а в статье ничего не нашли. Зато примерно похоже на данные из википедии.

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

torsion_values = np.arange(-180, 181, 12)

torsion_energies = np.array([run_orca(torsion_template.format(value))
                             for value in torsion_values])

plt.plot(torsion_values, torsion_energies, 'o')
plt.xticks(np.arange(-180, 181, 60))
plt.show()

def torsion_potential(x, v, n, g, a):
    return v / 2 * (1 + np.cos(np.radians(n * np.pi * x - g))) + a

popt_t, pcov = optimize.curve_fit(torsion_potential, torsion_values, torsion_energies)

popt_t

array([ 5.19887996e-03,  9.55065908e-01, -1.46395981e-02, -7.95705882e+01])

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. Неудивительно, ведь это заторможенная конформация этана.

plt.plot(np.arange(-180, 180, 1), torsion_potential(np.arange(-180, 180, 1), *popt_t))

[<matplotlib.lines.Line2D at 0x7fd468301438>]

Да, получилось правильно.

torsion_for_optimization = partial(torsion_potential, **dict(zip(['v', 'n', 'g', 'a'], popt_t)))

optimization_results = optimize.minimize_scalar(torsion_for_optimization, bounds=(0, 120))

optimization_results.x

59.98655403993325

Ну, почти 60!

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

long_bond_lengths = np.arange(1, 2.6, 0.1)

long_bond_length_energies = np.array([run_orca(bond_length_template.format(value))
                                      for value in long_bond_lengths])

plt.plot(long_bond_lengths, long_bond_length_energies)

[<matplotlib.lines.Line2D at 0x7fd467e1d940>]

popt_long, pcov = optimize.curve_fit(morze_potential, long_bond_lengths, long_bond_length_energies)

popt_long

array([  0.2399479 ,   1.65276024,   1.53227024, -79.08500183])

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()

Ну конечно, Морзе справился!

morze_long_for_optimization = partial(morze_potential, **dict(zip(['d', 'a', 'l0', 'b'], popt_long)))

optimization_results = optimize.minimize_scalar(morze_long_for_optimization, bounds=(1.4, 1.6))

optimization_results.x

1.5322701588167116

Ну, чуть дальше получилось, чем при небольшом шаге, но всё ещё не совпадает с длиной связи для GAFF.

Обсуждение

Разница со статьей получилась, видимо, потому, что авторы статьи исходили из других предположений для ab initio, нежели мы, используя другие модели для потенциалов.