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

import subprocess
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize 

%matplotlib inline

Важные функции

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
    )
    stdout = p.communicate()[0].decode()
    for s in stdout.splitlines():
        if "FINAL SINGLE POINT ENERGY" in s:
            return float(s.strip().split()[-1])

        
def fitfunc(p, x):
    """Target function"""
    return p[0]*pow(p[1]-x,2) + p[2]


def errfunc(p, x, y, fitfunc=fitfunc): 
    """Error function"""
    return fitfunc(p, x) - y


def draw(x, y, p0=(1,1, -79), fitfunc=fitfunc, xlims=None):
    p1, success = optimize.leastsq(lambda p, x, y: errfunc(p, x, y, fitfunc), p0[:], args=(x, y))
    print(success)
    print(f"Optimized params: {p1}")
    plt.plot(
        x, y, "ro",
        x, fitfunc(p1, x), "r-", 
        c="blue", alpha=0.5
    )
    if xlims is not None:
        plt.xlim(*xlims)
    plt.show()

template = """!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 {}
*
"""

def format_input(bond, angle, tors):
    return template.format(bond, angle, tors, tors-120, tors-240)

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

Ищем:

• длину связи C-C
• углы HCC
• торсионный угол связи C-C

C-C

sp = []
ls = np.arange(1.346, 1.74, 0.02)

for i in ls:
    i = round(i, 5)
    sp.append(run_orca(format_input(i, 111.200, 180)))

draw(ls, np.array(sp))

Optimized params: [  0.59409231   1.54979289 -79.19781436]

Лучшей длинной связи оказалась 1.5498 Å, чуть длиннее, чем исходное значение в 1.5299 Å.

Угол HCC

sp = []
ls = np.arange(109.2, 113.2, 0.2)

for i in ls:
    i = round(i, 5)
    sp.append(run_orca(format_input(1.52986, i, 180)))

draw(ls, np.array(sp))

Optimized params: [ 4.06308373e-05  1.11194916e+02 -7.91975723e+01]

Оптимальный угол равен 111.19, весьма приближен к исходному углу

Торсион

sp = []
ls = np.arange(-100, 100, 10)  # чуть меньше диапазон -- все равно функция переодичная с периодом в 120 градусов

for i in ls:
    sp.append(run_orca(format_input(1.52986, 111.200, i)))

draw(ls, np.array(sp), p0=(0.05, 10, -0.002), fitfunc=lambda p, x: p[0]*np.cos(p[1]*x)+p[2])

3
Optimized params: [ 2.29142799e-03  1.00007444e+01 -7.91952848e+01]

Аппроксимируем здесь косинусом, функция переодична, минимумы (незаслоненная конфигурация) при углах около 60, как и должно быть. Максимумы соответствуют заслоненной конфигурации.