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

Алгоритм для Orca

def run_orca(inp):
    with open('orca.inp', 'w') as outfile:
        outfile.write(inp)
    with subprocess.Popen("orca orca.inp", shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) as p:
        out=p.communicate()[0].decode()
    for l in out.splitlines():
        if 'FINAL SINGLE POINT' in l:
            return float(l.strip().split()[4])
    return None

Оптимизированная структура этана в виде z-matrix

inp = '''!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
*
'''

def help_func(name, values, inp_str):
    energies = []
    times = []
    length = len(values) 
    for i, x in enumerate(values):
        print("{} {:>7.2f}".format(name, x), end = "  ")
        if i != 0:
            t = np.mean(times) * (length - i)
            print("remaining time {:.0f}:{:02.0f}".format(t//60, t % 60), end = " ")
        print("", end = "\r")
        inp = inp_str(x)
        start_time = time.time()
        energies.append(run_orca(inp))
        times.append(time.time() - start_time)
    print("\ndone")
    return energies

matrix = [x/100 for x in range(120, 180, 2)]
inp_str_bond = lambda x : '''!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 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 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''.format(x)
energy = help_func("bond", matrix, inp_str_bond)

x_o=matrix
y_o=energy

#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)

plt.plot(x_o, y_o, "ro", x_o, fitfunc(p1,x_o),"r-", c = "blue",alpha=0.5)

Проделайте аналогичные операции для валентного угла HCС, его значения должны изменяться от 109.2 до 113.2.

angles = [x/10 for x in range(109.2, 113.4, 0.2)]
inp_str = lambda x :'''!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 6 1.08439 111.200 120
H 2 1 6 1.08439 {} -120
*
'''.format(x)

energy = help_func("angle", angles, inp_str)

x_o=angles
y_o=energy

#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)

plt.plot(x_o, y_o, "ro", c='blue',alpha=0.5)

Проделайте аналогичные операции для торсионного угла CC, его значения должны изменяться от -180 до 180 c шагом 12

angles = [x/10 for x in range(-180, 180, 12)]
inp_str = lambda x :'''!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 6 1.08439 111.200 120
H 2 1 6 1.08439 {} -120
*
'''.format(x)

energy = help_func("angle", angles, inp_str)

x_o=angles
y_o=energy

fitfunc = lambda p, x: p[0] * np.cos(np.multiply(x, p[1]/2/np.pi) - p[2]) - p[3]
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [0.001, 0.33, 0,  79] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print("Optimized params:", p1)

plt.plot(x_o, y_o, "ro",c='gray',alpha=0.5)