In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import optimize
import subprocess
import os
import subprocess
In [2]:
def run_orca(inp, filename):
    with open(filename, 'w') as outfile:
        outfile.write(inp)
    p = subprocess.Popen("/srv/databases/orca/orca %s"%(filename), 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]
    lines = out.splitlines()
    e = 0
    for line in lines:
        #print(line[1:])
        if "FINAL SINGLE POINT" in str(line):
        #if line.find("FINAL SINGLE POINT") > -1:
            FSP = line.strip().split()
            e = float(FSP[4])
    return e
In [3]:
def change_by_step(step, start, end, const, flag):
    values = []
    energ = []
    for n in range(start, end):
        value = const + step*n
        values.append(value)
        if flag == 0:
            inp = '''!HF RHF 6-31G
            * int 0 1
            C 0 0 0 0 0 0 
            C 1 0 0 %s 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
            *
            '''%(str(value))
            energ.append(run_orca(inp, "orca_bond_%s.inp" % (str(n))))
        if flag == 1:
            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 %s 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
            *
            '''%(str(value))
            energ.append(run_orca(inp, "orca_hcc_%s.inp" % (str(n))))
    return values, energ
In [5]:
def show_graph(value_1, value_2, minimum, maximum, y_name, x_name, title_name):
    #Построение графика
    x_o=np.array(value_1)
    y_o=np.array(value_2)

    fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] 
    errfunc = lambda p, x, y: fitfunc(p, x) - y 
    p0 = [1,1, -79] 
    if x_name == 'Торсионный угол C-C (gegrees)':
        p0 = [1,1, -79, 0]
    p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
    print("Оптимизированные параметры:", p1)

    plt.plot(x_o, y_o, "ro", x_o, fitfunc(p1,x_o), "r-", c="teal", alpha=0.9)
    plt.xlim(minimum, maximum) 

    plt.axvline(x = p1[1], c = "slategrey", alpha=0.45)
    plt.axhline(y = p1[2], c = "slategrey", alpha=0.45)
    plt.plot(p1[1], p1[2],'ro', c = "slategrey", alpha=0.85)

    plt.ylabel(y_name, fontsize = 12)
    plt.xlabel(x_name, fontsize = 12)
    plt.title(title_name, pad = 25, fontsize = 14)

    plt.show()
In [6]:
lengths, energies = change_by_step(0.02, -11, 11, 1.52986, 0)
In [7]:
show_graph(lengths, energies, 1.3,1.75, 
           'Energy (kcal/mol)',
           'Length С-С (A)', 
           'Energy-C-C lenght dependance')

Оптимизированные параметры: [ 8.70040630e-36  9.99996332e-01 -3.13680513e-36]
In [8]:
lengths_new, energies_new = change_by_step(0.1, -11, 11, 1.52986, 0)
In [11]:
x_o=np.array(lengths_new)
y_o=np.array(energies_new)

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

p0 = [0.5,-78, -3, -6]
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print("Оптимизированные параметры:", p1)


#Plot it
plt.plot(x_o, y_o, "ro", c="teal", alpha=0.9)
plt.plot(x_o, fitfunc(p1, x_o), "r-", c="slategrey", alpha=0.45)
plt.xlim(0.4,2.6)
plt.ylabel('Energy (kcal/mol)', fontsize = 12)
plt.xlabel('Length С-С (A)', fontsize = 12)
plt.title('Dependance (step = 0.1)', pad = 25, fontsize = 14)
plt.show()

Оптимизированные параметры: [ 6.03369168e-29  4.40267548e-30 -3.00000631e+00 -5.99999837e+00]
In [12]:
angles, an_energies = change_by_step(0.2, 0, 21, 109.2, 1)
In [13]:
show_graph(angles, an_energies, 109,113.4, 
           'Energy (kcal/mol)',
           'Angle H-C-C (degrees)', 
           'Dependance')

Оптимизированные параметры: [ 6.43345150e-168  9.80256094e-001 -7.79120286e-164]
In [14]:
tor_angles = []
tor_energies = []
step = 12
for n in range(31):
    tor_angle = -180 + step*n
    tor_angles.append(tor_angle)
    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 %s
    H 2 1 6 1.08439 111.200 120
    H 2 1 6 1.08439 111.200 -120
    *
    '''%(str(tor_angle))
    with open('orca_cc_%s.inp' % (str(n)), 'w') as outfile:
        outfile.write(inp)
    p = subprocess.Popen("/srv/databases/orca/orca orca_cc_%s.inp" % (str(n)), 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]
    outlines = out.splitlines()
    for outline in outlines:
        if 'FINAL SINGLE POINT' in str(outline):
            en = (str(outline).strip('\n').split()[-1])[:-1]
            tor_energies.append(float(en))
In [17]:
x_o=np.array(tor_angles)
y_o=np.array(tor_energies)

fitfunc = lambda p, x: np.sin(x*p[0]+p[1])*p[2]+p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [1,1, -79, 0] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print("Оптимизированные параметры:", p1)


#Plot it
plt.plot(x_o, y_o, "ro", c="teal", alpha=0.9)
plt.plot(x_o, fitfunc(p1, x_o), "r-", c="slategrey", alpha=0.45)
plt.xlim(-185,185)
plt.ylabel('Energy (kcal/mol)', fontsize = 12)
plt.xlabel('Angle C-C (degrees)', fontsize = 12)
plt.title('Dependance', pad = 25, fontsize = 14)
plt.show()

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-17-29b7741e080f> in <module>
      6 
      7 p0 = [1,1, -79, 0] # Initial guess for the parameters
----> 8 p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
      9 print("Оптимизированные параметры:", p1)
     10 

/opt/anaconda3/lib/python3.7/site-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    381     if not isinstance(args, tuple):
    382         args = (args,)
--> 383     shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
    384     m = shape[0]
    385 

/opt/anaconda3/lib/python3.7/site-packages/scipy/optimize/minpack.py in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     24 def _check_func(checker, argname, thefunc, x0, args, numinputs,
     25                 output_shape=None):
---> 26     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     27     if (output_shape is not None) and (shape(res) != output_shape):
     28         if (output_shape[0] != 1):

<ipython-input-17-29b7741e080f> in <lambda>(p, x, y)
      3 
      4 fitfunc = lambda p, x: np.sin(x*p[0]+p[1])*p[2]+p[3] # Target function
----> 5 errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
      6 
      7 p0 = [1,1, -79, 0] # Initial guess for the parameters

ValueError: operands could not be broadcast together with shapes (31,) (0,)
In [18]:
In [19]:
In [ ]:
In [ ]: