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

In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import optimize
In [3]:
ORCA_PATH = "./orca/orca"  # На kodomo не грузится numpy (???) + мало памяти.
In [4]:
import subprocess
In [5]:
from IPython.display import Image, display
In [6]:
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 [7]:
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 [12]:
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
*'''

Посчитаем минимум для C-C связи

In [36]:
bonds = np.linspace(1.3, 1.7, 21)
bonds
Out[36]:
array([1.3 , 1.32, 1.34, 1.36, 1.38, 1.4 , 1.42, 1.44, 1.46, 1.48, 1.5 ,
       1.52, 1.54, 1.56, 1.58, 1.6 , 1.62, 1.64, 1.66, 1.68, 1.7 ])
In [37]:
run_orca_bonds = np.vectorize(lambda x: run_orca(bond.format(x)), otypes=[np.float])
In [38]:
bonds_e = run_orca_bonds(bonds)
In [39]:
bonds_e
Out[39]:
array([-79.03218253, -79.04174446, -79.04998776, -79.05703088,
       -79.06298212, -79.06794047, -79.07199644, -79.07523272,
       -79.0777249 , -79.07954205, -79.08074728, -79.08139822,
       -79.08154754, -79.0812433 , -79.08052939, -79.07944588,
       -79.07802931, -79.07631301, -79.07432738, -79.07210013,
       -79.06965651])

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

In [40]:
#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='darkred', alpha=0.5)
plt.plot(bonds, fitfunc(p1, bonds), "r-", c='red')
#plt.xlim(1,2)
plt.show()

Optimized params: [  0.74214241   1.55193405 -79.08296748]

В статье GAFF указана длина связи 1.526. У нас - 1.552, но с учетом кривого фитирования, скорее всего, должно получаться ~= 1.54

Аналогичную операцию повторим для углов

In [43]:
angles = np.linspace(109.2, 113.2, 21)
angles
Out[43]:
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 [45]:
angles_e = run_orca_angles(angles)
In [ ]:
def angle_potential(x, k, f0, a):
    return k / 2 * (x - f0) ** 2 + a
In [47]:
#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='darkred', alpha=0.5)
plt.plot(angles, fitfunc(p1, angles), "r-", c='red')
#plt.xlim(1,2)
plt.show()

Optimized params: [ 3.61815365e-05  1.11037155e+02 -7.94918725e+01]

Зафитировалось здесь лучше. Получили минимум в 111.0. В статье не указан.

In [60]:
torss = np.linspace(-180, 180, 31).astype(int)
torss
Out[60]:
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 [61]:
run_orca_torss = np.vectorize(lambda x: run_orca(tors.format(x)), otypes=[np.float])
In [63]:
torss_e = run_orca_torss(torss)
In [64]:
plt.plot(torss, torss_e, "x", c='darkred', alpha=0.5)
Out[64]:
[<matplotlib.lines.Line2D at 0x7f82e255de10>]

Кажется, можно приблизить синусоидой.

In [83]:
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: [ 5.19919783e-03  9.55067198e-01  1.42276031e-05 -7.95705884e+01]

Ура, получилось! Следовательно, значения минимумов соответствуют -180, -60, 60, 180 - т.е. анти-конформациям этана. Довольно неплохо.

И проделаем для длины связи с бОльшим шагом и на большем диапазоне.

In [84]:
bonds_l = np.linspace(1.0, 3.0, 21)
bonds_l
Out[84]:
array([1. , 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2. , 2.1, 2.2,
       2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3. ])
In [85]:
bonds_le = run_orca_bonds(bonds_l)
In [86]:
bonds_le
Out[86]:
array([-78.60460671, -78.82848398, -78.95960734, -79.0321825 ,
       -79.06794046, -79.08074726, -79.07944587, -79.06965649,
       -79.05492949, -79.03749227, -79.01873861, -78.99954659,
       -78.98047826, -78.96189884, -78.94404506, -78.92706354,
       -78.91103378, -78.89598374, -78.8819023 , -78.86875013,
       -78.85646952])
In [87]:
plt.plot(bonds_l, bonds_le, "x", c='darkred', alpha=0.5)
Out[87]:
[<matplotlib.lines.Line2D at 0x7f82e17d77d0>]

Я бы сказал, похоже на a/x^n + kx + b

In [89]:
#function is  f(x)=k(b-x)^2 + a
fitfunc = lambda p, x: p[0] / np.power(x, p[1]) + p[2] * x + p[3]  # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y  # Error function

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

#Plot it
plt.plot(bonds_l, bonds_le, "x", c='darkred', alpha=0.5)
plt.plot(bonds_l, fitfunc(p1, bonds_l), "r-", c='red')
#plt.xlim(1,2)
plt.show()

Optimized params: [  0.63099626   5.71099485   0.19319248 -79.4215608 ]

Может, и похоже. С 4 степенями свободы-то трудно сделать, чтобы плохо вышло.

In [90]:
from scipy.optimize import fmin
In [92]:
fmin(lambda x: fitfunc(p1, x), 1.5)

Optimization terminated successfully.
         Current function value: -79.070471
         Iterations: 11
         Function evaluations: 22
Out[92]:
array([1.54650879])

С этой не очень точной оптимизацией получилось примерно 1.55. Примерно так же, как в прошлый раз.