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

import psi4
import numpy as np
import py3Dmol
from scipy import optimize
from matplotlib import pyplot as plt

psi4.core.set_output_file('output.dat')

Представим структуру молекулы этана в виде z-матрицы, в которой удобно изменять значения длины C-C связи, плоского угла H-C-C и торсиона вокруг C-C:

ethane_template = '''
C
C 1 {0}
H 1 1.08439 2 {1}
H 1 1.08439 2 111.2 3 120
H 1 1.08439 2 111.2 3 -120
H 2 1.08439 1 111.2 3 {2}
H 2 1.08439 1 111.2 6 120
H 2 1.08439 1 111.2 6 -120
'''

Визуализируем молекулу:

mol = psi4.geometry(ethane_template.format(1.52986, 111.2, 180))
view = py3Dmol.view()
view.addModel(mol.save_string_xyz_file(), 'xyz')
view.setStyle({'stick': {}})
view.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

def scan_ethane(cc_bond=[1.52986], hcc_angle=[111.2], hcch_torsion=[180]): # единая функция для сканирования длины связи, угла и торсиона
    scan_length = max(len(cc_bond), len(hcc_angle), len(hcch_torsion))
    bonds = cc_bond * scan_length if len(cc_bond) == 1 else cc_bond
    angles = hcc_angle * scan_length if len(hcc_angle) == 1 else hcc_angle
    torsions = hcch_torsion * scan_length if len(hcch_torsion) == 1 else hcch_torsion
    energies = []
    for bond, angle, torsion in zip(bonds, angles, torsions):
        mol = psi4.geometry(ethane_template.format(bond, angle, torsion))
        psi4.set_options({'maxiter': 200, 'fail_on_maxiter' :  True})
        energies.append(psi4.energy('scf/cc-pvtz', molecule=mol))
    return energies

Расчет параметров связи C-C

bond_range = np.arange(1.31, 1.71, 0.02)
bond_energy = scan_ethane(cc_bond=bond_range)

Аппроксимируем зависимость энергии от длины связи гармоническим потенциалом:

harmonic = lambda x, k, x_0, c: k / 2 * (x - x_0) ** 2 + c
bond_param = optimize.curve_fit(harmonic, bond_range, bond_energy, method='lm')[0]
print('k = {:.4f} hartree / Ų, r₀ = {:.3f} Å'.format(*bond_param))

k = 1.3467 hartree / Ų, r₀ = 1.541 Å

plt.plot(bond_range, bond_energy, 'ko',
         bond_range, harmonic(bond_range, *bond_param), 'k',
         alpha=0.5)
plt.xlabel('C-C bond length, Å')
plt.ylabel('Energy, hartree')
plt.legend(['QM energy', 'Fit'])
plt.show()

Из графика видно, что на данном масштабе энергия уже не очень точно описывается параболой и заметна асимметрия относительно минимума. Чтобы учесть это, применим потенциал Морзе:

bond_range_1 = np.arange(1, 2.5, 0.05)
bond_energy_1 = scan_ethane(cc_bond=bond_range_1)

morse = lambda x, d, a, x_0, c: d * (1 - np.exp(-a * (x - x_0))) ** 2 + c
bond_param_1 = optimize.curve_fit(morse, bond_range_1, bond_energy_1, method='lm', p0=[0.1, 1, 1.5, -80])[0]
print('D = {:.4f}, a = {:.4f}, r₀ = {:.3f} Å'.format(*bond_param_1))

D = 0.2300, a = 1.6479, r₀ = 1.520 Å

plt.plot(bond_range_1, bond_energy_1, 'ko',
         np.arange(1, 2.5, 0.01), morse(np.arange(1, 2.5, 0.01), *bond_param_1), 'k',
         np.arange(1, 2.5, 0.01), harmonic(np.arange(1, 2.5, 0.01), *bond_param), 'k--',
         alpha=0.5)
plt.xlabel('C-C bond length, Å')
plt.ylabel('Energy, hartree')
plt.legend(['QM energy', 'Morse fit', 'Harmonic fit'])
plt.show()

Действительно, в широком диапазоне длины связи потенциал Морзе работает лучше, чем гармонический.

Расчет параметров угла H-C-C

angle_range = np.arange(109.2, 113.2, 0.2)
angle_energy = scan_ethane(hcc_angle=angle_range)

Для угла также используем гармонический потенциал:

angle_param = optimize.curve_fit(harmonic, angle_range, angle_energy, method='lm')[0]
print('k = {:.9f} hartree / deg², r₀ = {:.3f}°'.format(*angle_param))

k = 0.000076133 hartree / deg², r₀ = 111.132°

plt.plot(angle_range, angle_energy, 'ko',
         angle_range, harmonic(angle_range, *angle_param), 'k',
         alpha=0.5)
plt.xlabel('H-C-C angle, °')
plt.ylabel('Energy, hartree')
plt.legend(['QM energy', 'Fit'])
plt.show()

Расчет параметров торсионного угла H-C-C-H

torsion_range = np.arange(-180, 180, 10)
torsion_energy = scan_ethane(hcch_torsion=torsion_range)

Энергию торсионного угла будем описывать синусоидой:

periodic = lambda x, k, n, x_0, c: k * (1 + np.cos(n * x - x_0)) + c
torsion_param = optimize.curve_fit(periodic, torsion_range, torsion_energy, method='lm', p0=[0.001, 0.05, 0, -80])[0]
print('k = {:.6f}, n = {:.6f}, phi₀ = {:.3f}°'.format(*torsion_param))

k = 0.002477, n = 0.052366, phi₀ = -0.000°

plt.plot(torsion_range, torsion_energy, 'ko',
         np.arange(-180, 180, 1), periodic(np.arange(-180, 180, 1), *torsion_param), 'k',
         alpha=0.5)
plt.xlabel('H-C-C-H torsion, °')
plt.ylabel('Energy, hartree')
plt.legend(['QM energy', 'Fit'])
plt.xticks(range(-180, 240, 60))
plt.show()

Зависимость потенциала от значения торсионного угла имеет 3 максимума (-120°, 0°, 120° — заслоненные конформации) и 3 минимума (-60°, 60°, 180° — заторможенные конформации).