import psi4
import numpy as np
psi4.core.set_output_file('output.dat')
from IPython.display import display,Image
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from math import pi
from math import e
import py3Dmol

В этом практикуме предлагается подобрать параметры для классической молекулярной динамики исходя из квантовомеханических расчетов: запараметризовать зависимость энергии ковалентного взаимодействия от межатомного расстояния,торсионных и валентных углов.

Молекула этана

Простейшая молекула, которая будет рассмотрена -- этан. В виде z-матрицы он выглядит так:

inp = '''
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 180
H 2 1.08439 1 111.200 6 120
H 2 1.08439 1 111.200 6 -120
'''

Вот наш герой:

m = psi4.geometry(inp)
view = py3Dmol.view()

view.addModel(m.save_string_xyz_file(),'xyz')
view.setStyle({'model': -1}, {'stick': {}})
view.zoomTo()
view.show()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

Расчет энергии молекулы этана:

def run_psi4(inp):
    
    m = psi4.geometry(g)
    psi4.set_options({"maxiter": 200, "fail_on_maxiter" :  True})
    enr = psi4.energy('scf/cc-pvtz', molecule = m )
    return enr

print(f'Энергия молекулы этана: {run_psi4(inp):.2f} Дж')

Энергия молекулы этана: -79.23 Дж

Зависимость энергии от длины связи С-С

Буду варьировать длину в районе исходной длины из z-матрицы (1.53 A): возьму 20 точек от 1.33 А до 1.73 А

r_arr = np.linspace(1.529-0.2,1.529+0.2,20)
r_arr

array([1.329     , 1.35005263, 1.37110526, 1.39215789, 1.41321053,
       1.43426316, 1.45531579, 1.47636842, 1.49742105, 1.51847368,
       1.53952632, 1.56057895, 1.58163158, 1.60268421, 1.62373684,
       1.64478947, 1.66584211, 1.68689474, 1.70794737, 1.729     ])

enrgs = np.empty((20,))
for i, r in enumerate(r_arr):
    g = template = f'''
                        C
                        C 1 {r}
                        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 180
                        H 2 1.08439 1 111.200 6 120
                        H 2 1.08439 1 111.200 6 -120
                        '''
    energy = run_psi4(g)
    enrgs[i] = energy

Посмотрим как получилось + построим аппроксимацию параболой:

def plot_val(enrgs, vals):
    x_o=vals
    y_o=enrgs

    # Аппроксимируем параболой
    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] # Инициализация случайными параметрами
    p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
    # p1 -- Оптимизированныее параметры
    plt.plot(x_o, fitfunc(p1,x_o),
             "r-",
             c='red',
             alpha=0.5,

             label = "Подгонка параболой")

    plt.plot(vals, enrgs, 
             "ro",
             c = "green",
             marker = 'o',
             label = "Сырые данные")

    plt.legend()
    return p1

p1 = plot_val(enrgs, r_arr)

Полученные параметры из аппроксимации параболой:

print(f"Cкорость роста энергии при удалении от оптимального расстояния : {p1[0]:.2f}\n\
Длина связи с минимальной энергией: {p1[1]:.2f} A\n\
Минимум энергии: {p1[2]:.2f} Дж")

Cкорость роста энергии при удалении от оптимального расстояния : 0.59
Длина связи с минимальной энергией: 1.55 A
Минимум энергии: -79.26 Дж

Зависимость энергии от величины валентного угла

Будем варьировать валентный угол CCH от 109.2 до 113.2 с шагом 0.2

enrgs_angl = np.empty((20,))
angls = np.arange(109.2, 113.2, 0.2)
for i, angle in enumerate(angls):
    g = template = f'''
                        C
                        C 1 1.52986
                        H 1 1.08439 2 {angle}
                        H 1 1.08439 2 {angle} 3 120
                        H 1 1.08439 2 {angle} 3 -120
                        H 2 1.08439 1 {angle} 3 180
                        H 2 1.08439 1 {angle} 6 120
                        H 2 1.08439 1 {angle} 6 -120
                        '''
    energy = run_psi4(g)
    enrgs_angl[i] = energy

p1 = plot_val(enrgs_angl, angls)

<matplotlib.legend.Legend at 0x7f857c5cd0d0>

Ну, идеально, что сказать. Параметры параболы:

p1

array([ 3.06271714e-04,  1.11143282e+02, -7.92600134e+01])

Зависимость энергии от величины торсионного угла

Теперь меняем торсионный угол от -180 до 180

enrgs_tors= np.empty((30,))
angls_tors = np.arange(-180, 180, 12)
for i, angle in enumerate(angls_tors):
    g = template = 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 {angle}
                        H 2 1.08439 1 111.200 6 120
                        H 2 1.08439 1 111.200 6 -120
                        '''

    energy = run_psi4(g)
    enrgs_tors[i] = energy

plt.plot(angls_tors, enrgs_tors, 
         "r-",
         c = "green",
         marker = 'o',
         label = "Сырые данные")

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

Попробуем аппроксимировать косинусом:

def plot_val2(enrgs, vals):
    x_o=vals
    y_o=enrgs

    # Аппроксимируем параболой
    fitfunc = lambda p, x: p[0]*(np.cos(p[1]*x*pi/180+p[2]))+p[3] # Target function
    errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

    p0 = [0,2.1, -10,0] # Инициализация случайными параметрами
    p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
    # p1 -- Оптимизированныее параметры
    plt.plot(x_o, fitfunc(p1,x_o),
             "r-",
             c='red',
             alpha=0.5,

             label = "Подгонка косинусом")

    plt.plot(vals, enrgs, 
             "ro",
             c = "green",
             marker = 'o',
             label = "Сырые данные")

    plt.legend()
    return p1

p1 = plot_val2(enrgs_tors, angls_tors)

Функция имеет 2 минимума, которые отвечают двум заторможенным конформациям. 3 максимума отвечают заслоненным конформациям. Парметры:

print(f"Амплитуда изменений энергии : {p1[0]:.3f}\n\
Угловая частота: {p1[1]:.2f}\n\
Начальная фаза: {p1[2]:.2f}\n\
Среднее значение энергии: {p1[2]:.2f}")

Амплитуда изменений энергии : -0.002
Угловая частота: 3.00
Начальная фаза: -9.42
Среднее значение энергии: -9.42