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

In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import optimize
In [14]:
ORCA_PATH = "./orca/orca"  # На kodomo не грузится numpy (???) + мало памяти.
In [15]:
import subprocess
In [16]:
from IPython.display import Image, display

напишем функцию для запуска orca:

In [18]:
def run_orca(inp):
    with open('orca.inp', 'w') as outfile:
        outfile.write(inp)
    with subprocess.Popen(f"{ORCA_PATH} 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 np.nan

проверяем

In [19]:
input_orca = '''!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 [8]:
run_orca(input_orca)
Out[8]:
-79.081531418112

Это успех!

Сделаем шаблоны файлов

In [20]:
bond = '''!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 = '''!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
*'''
tors = '''!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
*'''

Связи

In [7]:
bonds = np.linspace(1.30986, 1.72986, 22)
bonds
Out[7]:
array([1.30986, 1.32986, 1.34986, 1.36986, 1.38986, 1.40986, 1.42986,
       1.44986, 1.46986, 1.48986, 1.50986, 1.52986, 1.54986, 1.56986,
       1.58986, 1.60986, 1.62986, 1.64986, 1.66986, 1.68986, 1.70986,
       1.72986])
In [22]:
run_orca_bonds = np.vectorize(lambda x: run_orca(bond.format(x)), otypes=[np.float])
In [5]:
bonds_e = run_orca_bonds(bonds)
bonds_e
Out[5]:
array([-79.03706923, -79.04596559, -79.05360305, -79.06009493,
       -79.06554486, -79.07004752, -79.07368951, -79.07654992,
       -79.07870106, -79.080209  , -79.08113408, -79.08153143,
       -79.08145141, -79.08093996, -79.08003906, -79.07878699,
       -79.07721868, -79.07536597, -79.07325792, -79.07092099,
       -79.06837928, -79.06565477])

Отрисуем результаты и зафитим функцию по ним.

In [28]:
#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=(bonds, bonds_e))
print("Optimized params:", p1)

#Plot it
plt.plot(bonds, bonds_e, "x", c='green', alpha=0.5)
plt.plot(bonds, fitfunc(p1, bonds), "r-", c='blue')
#plt.xlim(1,2)
plt.show()

Optimized params: [  0.67723276   1.55759168 -79.08260194]
Давайте сравним со статьей GAFF - там указана длина связи 1.526. У нас - 1.557. Можно списать на неточный фитинг, по точкам получается около 1.553

Углы

In [12]:
angles = np.linspace(109.2, 113.2, 21)
angles
Out[12]:
array([109.2, 109.4, 109.6, 109.8, 110. , 110.2, 110.4, 110.6, 110.8,
       111. , 111.2, 111.4, 111.6, 111.8, 112. , 112.2, 112.4, 112.6,
       112.8, 113. , 113.2])
In [44]:
run_orca_angles = np.vectorize(lambda x: run_orca(angle.format(x)), otypes=[np.float])
In [8]:
angles_e = run_orca_angles(angles)
angles_e
Out[8]:
array([-79.08141768, -79.08144397, -79.08146694, -79.08148658,
       -79.0815029 , -79.08151591, -79.08152561, -79.08153201,
       -79.08153511, -79.08153492, -79.08153143, -79.08152466,
       -79.08151461, -79.08150128, -79.08148468, -79.0814648 ,
       -79.08144167, -79.08141527, -79.08138562, -79.08135271,
       -79.08131656])
In [10]:
def angle_potential(x, k, f0, a):
    return k / 2 * (x - f0) ** 2 + a
In [29]:
#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, -8] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(angles, angles_e))
print("Optimized params:", p1)

#Plot it
plt.plot(angles, angles_e, "x", c='green', alpha=0.5)
plt.plot(angles, fitfunc(p1, angles), "r-", c='blue')
#plt.xlim(1,2)
plt.show()

Optimized params: [ 4.10776520e-05  1.10890769e+02 -7.90815354e+01]

Данных в сатье нет. Фитинг вполне хорош.

Торсионные углы

In [4]:
torss = np.linspace(-180, 180, 31).astype(int)
torss
Out[4]:
array([-180, -168, -156, -144, -132, -120, -108,  -96,  -84,  -72,  -60,
        -48,  -36,  -24,  -12,    0,   12,   24,   36,   48,   60,   72,
         84,   96,  108,  120,  132,  144,  156,  168,  180])
In [3]:
torss_e = run_orca_torss(torss)
torss_e
Out[3]:
array([-79.1975724 , -79.19713771, -79.19599652, -79.19457972,
       -79.19342857, -79.19298772, -79.19342857, -79.19457972,
       -79.19599652, -79.19713771, -79.1975724 , -79.19713771,
       -79.19599652, -79.19457972, -79.19342857, -79.19298772,
       -79.19342857, -79.19457972, -79.19599652, -79.19713771,
       -79.1975724 , -79.19713771, -79.19599652, -79.19457972,
       -79.19342857, -79.19298772, -79.19342857, -79.19457972,
       -79.19599652, -79.19713771, -79.1975724 ])
In [5]:
fitfunc = lambda p, x: p[0] * (1 + np.cos(p[1] * np.pi * np.radians(x) - p[2])) + p[3]  # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y  # Error function

p0 = [0, 1, -1, -1] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(torss, torss_e))
print("Optimized params:", p1)

#Plot it
plt.plot(torss, torss_e, "x", c='darkred', alpha=0.5)
good_torss = np.linspace(-180, 180, 361)
plt.plot(good_torss, fitfunc(p1, good_torss), "r-", c='red')
#plt.xlim(1,2)
plt.show()

Optimized params: [ 2.29204658e-03  9.54986114e-01 -2.45630292e-06 -7.91975763e+01]

Получается, что значения минимумов соответствуют -180, -60, 60, 180 - т.е. анти (заторможенным)-конформациям этана. Максимумы - заслонённая, наименее выгодная конформация.