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

Молекулярное моделирование биополимеров, 2017. Факультет биоинженерии и биоинформатики, МГУ им. М. В. Ломоносова

In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import optimize 
import subprocess
import numpy as np
In [88]:
inp = '''!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 %.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 %.1f
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
In [4]:
ref_bond = 1.52986
ref_angle = 111.2
ref_dihed = 180
In [66]:
def run_orca(inp):
    with open('orca.inp', 'w') as outfile:
        outfile.write(inp)
    p = subprocess.Popen("/srv/databases/orca/orca orca.inp", 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]
    O = out.splitlines()
    result = None
    for line in O:
        if 'FINAL SINGLE POINT ENERGY' in line:
            result = line

    # extract energy: FINAL SINGLE POINT ENERGY'
    # and return it as float
    
    return result
In [55]:
def fitparam(my_x,my_y,xlim=(0,3)):
    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=(my_x,my_y))
    print "Optimized params:", p1

    #Plot it
    plt.plot(my_x,my_y, "ro", my_x ,fitfunc(p1,my_x),"r-",c='blue',alpha=0.5)
    plt.xlim(xlim)
    plt.show()

Связь

In [7]:
res = []
ref = 1.52986
delta = 0.04
points = 70
if ref-(points/2*delta) < 0:
    print("Negative bond length! You are idiot!")
L = np.arange(ref-delta*points/2,ref+delta*points/2,delta)
In [8]:
for i in L:
    inpstr = inp % (i, ref_angle, ref_dihed)
    res.append(run_orca(inpstr))
In [9]:
res = [float(x.split()[-1]) for x in res]
In [10]:
plt.scatter(L[3:],res[3:])
plt.show()

Похоже на потенциал Морзе, если присмотреться поближе

In [11]:
plt.scatter(L[15:],res[15:])
plt.show()
In [12]:
fitparam(L[15:-5],res[15:-5])

Optimized params: [  0.99731425   1.92018577 -79.20221602]

Вообще не туда

In [13]:
fitparam(L[18:-20],res[18:-20])

Optimized params: [  1.45881555   1.64959812 -79.1545218 ]

Уже поближе

In [14]:
fitparam(L[28:-25],res[28:-25])

Optimized params: [  0.58547809   1.59265437 -79.08238855]

Что-то похожее

Вспомним теперь о потенциале Морзе и попробуем его параметризовать. Начальная догадка для p[0] - глубина ямы, p[2] - длина связи. Хотя при фитировании это оказалась вовсе не длина связи

In [40]:
def fitmorse(my_x,my_y):
    my_x = np.array(my_x)
    my_y = np.array(my_y)
    
    
    # f(x) = D*(1-e^(-a(x - x0)))^2
    fitfunc = lambda p, x: p[0]*((1 - np.exp(-p[1]*(x-p[2])))**2) # Target function
    errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

    p0 = [-79, 1, -1.5] # Initial guess for the parameters
    p1, success = optimize.leastsq(errfunc, p0[:], args=(my_x,my_y))
    print "Optimized params:", p1

    #Plot it
    plt.plot(my_x,my_y, "ro", my_x ,fitfunc(p1,my_x),"r-",c='blue',alpha=0.5)
    plt.xlim(0,3)
    plt.show()
In [41]:
fitmorse(L[3:],res[3:])

Optimized params: [ -7.89505081e+01   5.87319306e+00  -5.26663458e-02]

Угол

In [48]:
def angle():
    res = []
    ref = ref_angle
    delta = 0.2
    points = 20
    if ref-(points/2*delta) < 0:
        print("Negative bond length! You are idiot!")
    L = np.arange(ref-delta*points/2,ref+delta*points/2,delta)
    for i in L:
        inpstr = inp % (ref_bond, i, ref_dihed)
        res.append(run_orca(inpstr))

    res = [float(x.split()[-1]) for x in res]
    return L,res
In [49]:
lA, resA = angle()
In [53]:
plt.scatter(lA[:],resA[:])
plt.ylim([-79.082,-79.081])
plt.show()

Отличия в энергии незначительны. Опишем параболически

In [56]:
fitparam(lA,resA,(109,113))

Optimized params: [  4.10917975e-05   1.10890455e+02  -7.90815353e+01]

Двугранный угол

In [96]:
def torsion():
    res = []
    L = np.arange(-180,181,6)
    for i in L:
        inpstr = inp % (ref_bond, ref_angle, i)
        res.append(run_orca(inpstr))

    return L,res
In [97]:
lT, resT = torsion()
In [98]:
resT = [float(x.split()[-1]) for c,x in enumerate(resT)]
In [99]:
plt.scatter(lT,resT)
Out[99]:
<matplotlib.collections.PathCollection at 0x7f65e8cb3510>

Будем приближать косинусом

In [102]:
def fitcosine(my_x,my_y):
    my_x = np.array(my_x)
    my_y = np.array(my_y)
    
    
    # f(x) = D*(1-e^(-a(x - x0)))^2
    fitfunc = lambda p, x: p[0]*(np.cos(p[1]*x)) + p[2] # Target function
    errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

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

    #Plot it
    plt.plot(my_x,my_y, "ro", my_x ,fitfunc(p1,my_x),"r-",c='blue',alpha=0.5)
    plt.xlim(-180,180)
    plt.show()
In [103]:
fitcosine(lT, resT)

Optimized params: [  2.29217272e-03   9.94834910e-01  -7.91952843e+01]

Получилось