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

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from IPython.display import Image, display
import subprocess
%matplotlib inline

Задание функций для работы

def make_val(inp, bond, angle, tors):
    return inp.format(bond, angle, tors, tors-120, tors-240)

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)
    #print(p.communicate())
    out=p.communicate()[0]
    for i in range(len(out.splitlines())):
        s=out.splitlines()[i]
        s=a.decode()
        if 'FINAL SINGLE POINT ENERGY' in s:
            return float(s.strip().split()[-1])

def draw_graph(X, Y, xmin, xmax):
    # fake x array , replace with real one
    x_o=X
    # fake y array, replace with energies
    y_o=Y

    #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(xmin, xmax)
    plt.show()

Вычисление параметров

Имеет оптимизированную молекулу этана, у которой пропущены значения для:
- длины связи С-С,
- валентного угла HCC,
- торсионного угла СС

inp = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 {} 0 0 
H 1 2 0 1.08439 {} 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 {}
H 2 1 3 1.08439 111.200 {}
H 2 1 3 1.08439 111.200 {}
*
'''

Варьирование длины связи С-С

sp = []
j = 0

for i in np.arange(1.346, 1.74, 0.02):
    i = round(i, 5)
    sp.append(run_orca(make_val(inp, i, 111.200, 180)))

draw_graph(np.arange(1.346, 1.74, 0.02), np.array(sp), 1.3, 1.8)

Optimized params: [  0.60674062   1.5531397  -78.99305146]

Мы варьировали значение длины связи С-С от 1.346 до 1.740 с шагом 0.02 Å (по итогу - 20 вариантов длины связи).

Полученное значение минимума по энергии соответствует длине связи, равной 1.5531397 Å. Исходное значение в оптимизированной структуре - 1.52986 Å. Наша связь получилась длинноватой, не самая удачная оптимизация.

Варьирование угла НСС

sp = []
j = 0

for i in np.arange(109.2,113.2,0.2):
    i = round(i, 5)
    sp.append(run_orca(make_val(inp, 1.52986, i, 180)))

draw_graph(np.arange(109.2,113.2,0.2), np.array(sp), 109, 113)

Optimized params: [ 4.09206083e-05  1.09934889e+02 -7.89927841e+01]

Мы варьировали значение угла HCC от 109.2 до 113.2 с шагом 0.2.

Полученное значение минимума по энергии соответствует углу, равному 109.93, несколько ниже исходного оптимального значения.

Варьирование торсиона СС

sp = []
j = 0

for i in np.arange(-180,180,12):
    sp.append(run_orca(make_val(inp, 1.52986, 111.200, 180)))

fitfunc = lambda p, x: p[0]*(1 + np.cos(p[1]*x - p[2]))+p[3] 
errfunc = lambda p, x, y: fitfunc(p, x) - y 

p0 = [0.25,10,0,-79]
p1, success = optimize.leastsq(errfunc, p0[:], args=(np.arange(-180,180,12), np.array(sp)))
print("Optimized params:", p1)

plt.plot(np.arange(-180,180,12), np.array(sp), "ro", np.arange(-180,180,12), fitfunc(p1,np.arange(-180,180,12)),"r-", c='blue', alpha=0.5)
plt.xlim(-200,200)
plt.show()

Optimized params: [ 2.29231461e-03  1.00007394e+01 -6.15056377e-07 -7.91975767e+01]

Как ни странно, мы видим периодическую функцию, которую аппроксимировали с помощью косинуса. Энергетическим минимумам соответствуют заторможенной конформации, а максимумы - заслоненной. Период между максимумами/минимумами равен 120°.