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

import psi4

Немного подправим вид матрицы.

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 5 120
H 2 1.08439 1 111.200 5 -120
'''

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

run_psi4(inp)

Оптимизируем длину связи.

def calculate_energy(bond_length):
    geom = f'''
C 
C 1 {bond_length}
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 0
H 2 1.08439 1 111.200 3 120
H 2 1.08439 1 111.200 3 -120
'''
    return run_psi4(geom)

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

from IPython.display import clear_output

%%time
x_o = np.arange(1.3, 1.7, 0.02)
y_o = np.array(list(map(calculate_energy, x_o)))
clear_output()

#function is  f(x)=k(b-x)^2 + a
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] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print("Optimized params:", p1)

#Plot it
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.xlim(1,2)
plt.show()

*** tstop() called on ruslan-ThinkPad at Thu Mar 31 23:27:25 2022
Module time:
	user time   =       0.93 seconds =       0.02 minutes
	system time =       0.02 seconds =       0.00 minutes
	total time  =          1 seconds =       0.02 minutes
Total time:
	user time   =      23.39 seconds =       0.39 minutes
	system time =       0.72 seconds =       0.01 minutes
	total time  =         24 seconds =       0.40 minutes
Optimized params: [  0.75232986   1.55165169 -79.25655594]

/home/ruslan/miniconda3/envs/molsim/lib/python3.7/site-packages/ipykernel_launcher.py:14: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "ro" (-> color='r'). The keyword argument will take precedence.
  
/home/ruslan/miniconda3/envs/molsim/lib/python3.7/site-packages/ipykernel_launcher.py:14: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "r-" (-> color='r'). The keyword argument will take precedence.
  

CPU times: user 21.1 s, sys: 591 ms, total: 21.7 s
Wall time: 21.5 s

Оптимальная длина связи по солверу 1.5516 Å.

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

def calculate_energy_angle(angle):
    geom = f"""
C 
C 1 1.52986
H 1 1.08439 2 {angle}
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 0
H 2 1.08439 1 111.200 3 120
H 2 1.08439 1 111.200 3 -120
"""
    return run_psi4(geom)

x_angle = np.arange(109.2, 113.2, 0.2)
y_angle = np.array(list(map(calculate_energy_angle, x_angle)))
clear_output()

#function is  f(x)=k(b-x)^2 + a
p0 = [1, 111, -79] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_angle, y_angle))
print("Optimized params:", p1)

#Plot it
plt.plot(x_angle, y_angle, "ro", x_angle,fitfunc(p1,x_angle),"r-",c='blue',alpha=0.5)
plt.xlim(109, 114)
plt.show()

*** tstop() called on ruslan-ThinkPad at Thu Mar 31 23:27:43 2022
Module time:
	user time   =       0.87 seconds =       0.01 minutes
	system time =       0.03 seconds =       0.00 minutes
	total time  =          1 seconds =       0.02 minutes
Total time:
	user time   =      41.24 seconds =       0.69 minutes
	system time =       1.25 seconds =       0.02 minutes
	total time  =         42 seconds =       0.70 minutes
Optimized params: [ 3.93926846e-05  1.11938788e+02 -7.92550798e+01]

/home/ruslan/miniconda3/envs/molsim/lib/python3.7/site-packages/ipykernel_launcher.py:11: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "ro" (-> color='r'). The keyword argument will take precedence.
  # This is added back by InteractiveShellApp.init_path()
/home/ruslan/miniconda3/envs/molsim/lib/python3.7/site-packages/ipykernel_launcher.py:11: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "r-" (-> color='r'). The keyword argument will take precedence.
  # This is added back by InteractiveShellApp.init_path()

Оптимальным углом оказывается 111.94 °.

Далее исследуем торсионные углы.

def calculate_energy_torsion(torsion):
    adjust = lambda x: (x + 180 + 360) % 360 - 180 if x != 180 else 180
    geom = 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 {torsion}
H 2 1.08439 1 111.200 3 {adjust(torsion + 120)}
H 2 1.08439 1 111.200 3 {adjust(torsion - 120)}
'''
    return run_psi4(geom)

x_torsion = np.arange(-180, 181, 12)
y_torsion = np.array(list(map(calculate_energy_torsion, x_torsion)))
clear_output()

#Plot it
plt.plot(x_torsion, y_torsion, "ro")
plt.show()

*** tstop() called on ruslan-ThinkPad at Thu Mar 31 23:28:14 2022
Module time:
	user time   =       1.02 seconds =       0.02 minutes
	system time =       0.03 seconds =       0.00 minutes
	total time  =          1 seconds =       0.02 minutes
Total time:
	user time   =      71.19 seconds =       1.19 minutes
	system time =       2.17 seconds =       0.04 minutes
	total time  =         73 seconds =       1.22 minutes

Видим синусоиду с минимумами при значении параметра $\pm$60 ° и 180nbsp;°. Есть 3 минимума и максимума.

Сравним полученные оптимизированные значения ковалентных параметров со статьёй про GAFF. Непосредственно в статье приводится только значение длины связи, и оно составляет 1.526 — меньше моих 1.552. Если посмотреть на график, где происходила аппроксимация и оптимизация длины связи, то становится видно, что он не вполне хорошо приближается параболой и истинный минимум лежит как раз немного левее вершины параболы. С ковалентным углом проверить результат по статье не получится, так как там описывается только эвристика усреднения для вычисления угла между атомами трёх типов, если известны углы, когда первый и третий атомы одного типа. (Сами равновесные углы не приводятся).