Определение констант ковалентных взаимодействий для молекулярной механики.

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

inp = lambda length=1.52986, val=111.200, tors=180:'''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 {0} 0 0 
H 1 2 0 1.08439 {1} 0
H 1 2 3 1.08439 {1} 120
H 1 2 3 1.08439 {1} -120
H 2 1 3 1.08439 {1} {2}
H 2 1 5 1.08439 {1} {3}
H 2 1 5 1.08439 {1} {4}
*
'''.format(length, val, tors, tors-60, tors+60-360)

# запуск orca

def run_orca(inp):
    with open('inp', 'w') as outfile:
        outfile.write(inp)
    p = subprocess.Popen("orca inp", shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]

    for line in out.splitlines():
        if "FINAL SINGLE POINT ENERGY" in line:
            return float(line.split('FINAL SINGLE POINT ENERGY')[-1].strip())

bond_en = []
bond_lengths = np.linspace(length_default-step*9, length_default+step*10, num=20)
for bond_length in bond_lengths:
    bond_en.append(run_orca(inp(length=bond_length)))
    
val_en = []
vals = np.linspace(109.2, 113.2, num=20)
for val in vals:
    val_en.append(run_orca(inp(val=val)))
    
tors_en = []
tors = np.linspace(-180, 180, num=30)
for tor in tors:
    tors_en.append(run_orca(inp(tors=tor)))

def quadratic_aproximation(x_o, y_o, xlab, col):
    # апрокимируем квадратичной функцией 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
    print 'Function of approximation f(x)=k(b-x)^2 + a '

    plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c=col,alpha=0.5)
    plt.title('%s vs Energy' % xlab.title())
    plt.xlabel(xlab)
    plt.ylabel('energy')
    plt.show()

quadratic_aproximation(bond_lengths, bond_en, 'length', 'orange')

Optimized params: [  0.60466381   1.55589283 -79.08185347]
Function of approximation f(x)=k(b-x)^2 + a 

Зависимость энергии молекулы от длины С-С связи неплохо апрокисмируется квадратичной функцией.

quadratic_aproximation(vals, val_en, 'valence angle', 'red')

Optimized params: [  3.31057196e-04   1.12270856e+02  -7.90819111e+01]
Function of approximation f(x)=k(b-x)^2 + a 

Зависимость энергии молекулы от величины валентного угла HCC квадратична.

quadratic_aproximation(tors, tors_en, 'torsion angle', 'black')

Optimized params: [  1.35069094e-08  -8.60732198e+01  -7.90827614e+01]
Function of approximation f(x)=k(b-x)^2 + a 

Зависимость энергии молекулы от величины торсионного угла CC не является квадратичной.

Похоже, что зависимость периодическая.

# увеличим шаг изменения длины С-С связи до 0.1 А

step = 0.1
bond_en_1 = []
bond_lengths_1 = np.linspace(length_default-step*9, length_default+step*10, num=20)
for bond_length in bond_lengths_1:
    bond_en_1.append(run_orca(inp(length=bond_length)))

plt.plot(bond_lengths_1, bond_en_1, 'ro', c='orange',alpha=0.5)
plt.title('Length vs Energy')
plt.xlabel('length')
plt.ylabel('energy')
plt.show()

Апроксимируем степенной функцией:

x_o = bond_lengths_1
y_o = bond_en_1

fitfunc = lambda p, x: np.power(x, p[1])+p[0] 
errfunc = lambda p, x, y: fitfunc(p, x) - y

p0 = [-1,-79] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))

plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='orange',alpha=0.5)
plt.title('Length vs Energy')
plt.xlabel('length')
plt.ylabel('energy')
plt.show()