
import psi4
import numpy as np


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


from IPython.display import Image


psi4.core.set_output_file('output.dat')


import py3Dmol


def run_psi4(g):
    m = psi4.geometry(g)  # 
    psi4.set_options({"maxiter": 200, "fail_on_maxiter" :  True})
    e = psi4.energy('scf/cc-pvtz', molecule = m )

    return e


geom_template = '''
C
C 1 {0}
H 1 1.08439 2 111.2
H 1 1.08439 2 111.2 3 120
H 1 1.08439 2 111.2 3 -120
H 2 1.08439 1 111.2 3 {2}
H 2 1.08439 1 111.2 6 120
H 2 1.08439 1 {1} 6 -120
'''


cc_len_real = 1.52986
hcc_angle_real = 111.2
cc_tors_real = -180


geom_ = geom_template.format(cc_len_real, hcc_angle_real, cc_tors_real)
e_real = run_psi4(geom_)
e_real

-79.26001266654059


view = py3Dmol.view()
g = geom_template.format(cc_len_real, hcc_angle_real, cc_tors_real)
m = psi4.geometry(g)
view.addModel(m.save_string_xyz_file(),'xyz')
view.zoomTo()
view.setStyle({'model': -1}, {"stick": {}})

view.show()


view = py3Dmol.view()
g = geom_template.format(3, hcc_angle_real, cc_tors_real)
m = psi4.geometry(g)
view.addModel(m.save_string_xyz_file(),'xyz')
view.zoomTo()
view.setStyle({'model': -1}, {"stick": {}})

view.show()


x_o = np.arange(1.52986 - 0.02*10, 1.52986 + 0.02*10, 0.02)
y_o = np.zeros(x_o.shape)

for i, x in enumerate(x_o):
    geom_ = geom_template.format(x, hcc_angle_real, cc_tors_real)
    e = run_psi4(geom_) 
    y_o[i] = e


plt.plot(x_o, y_o, '-b.',  [cc_len_real], [e_real], 'r.', alpha=0.5)
plt.xlabel('Длина связи C-C')
plt.ylabel('Энергия')
plt.legend()

<matplotlib.legend.Legend at 0x7ff189b079a0>


# 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, cc_len_real, -79] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print("Optimized params:", p1)

Optimized params: [  0.61117515   1.54313697 -79.26049059]


#Plot it
plt.plot(x_o, y_o, "ro",    x_o, fitfunc(p1,x_o), "-r",    c='blue', alpha=0.5)
plt.xlabel('Длина связи C-C')
plt.ylabel('Энергия')
plt.show()


x_o = np.arange(1.52986 - 0.1*10, 1.52986 + 0.1*10, 0.02)
y_o = np.zeros(x_o.shape)

for i, x in enumerate(x_o):
    geom_ = geom_template.format(x, hcc_angle_real, cc_tors_real)
    e = run_psi4(geom_) 
    y_o[i] = e


plt.plot(x_o, y_o, '-b.',  [cc_len_real], [e_real], 'r.', alpha=0.5)
plt.xlabel('Длина связи C-C')
plt.ylabel('Энергия')
plt.legend()

<matplotlib.legend.Legend at 0x7ff189a82760>


Image(url='https://upload.wikimedia.org/wikipedia/commons/thumb/5/52/Morse-potential-ru.svg/1024px-Morse-potential-ru.svg.png', width=600)


# Morze potential f(x) = c * ( 1 - exp(-a(x-b)) )^2

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

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

Optimized params: [ -0.08766934   5.43734272 -79.2249846 ]


#Plot it
plt.plot(x_o, y_o, "b-", label='real', alpha=0.5)
plt.plot(x_o, fitfunc(p1,x_o), "r-", alpha=0.5, label='fit')
plt.xlabel('Длина связи C-C')
plt.ylabel('Энергия')
plt.legend()
plt.show()


view = py3Dmol.view()
g = geom_template.format(cc_len_real, 10, cc_tors_real)
m = psi4.geometry(g)
view.addModel(m.save_string_xyz_file(),'xyz')
view.zoomTo()
view.setStyle({'model': -1}, {"stick": {}})
view.show()


x_o = np.arange(109.2, 113.2, 0.2)
y_o = np.zeros(x_o.shape)

for i, x in enumerate(x_o):
    geom_ = geom_template.format(cc_len_real, x, cc_tors_real)
    y_o[i] = run_psi4(geom_)


plt.plot(x_o, y_o, '-b.',  [hcc_angle_real], [e_real], 'r.', alpha=0.5)
plt.xlabel('Угол C-C')
plt.ylabel('Энергия')
plt.legend()

<matplotlib.legend.Legend at 0x7ff1898b6940>


# 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 = [0.1, 111, 0] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))

print("Optimized params:", p1)

Optimized params: [ 3.80665270e-05  1.11132163e+02 -7.92600128e+01]


#Plot it
plt.plot(x_o, y_o, "ro",    x_o, fitfunc(p1, x_o), "r-", c='blue', alpha=0.5)
plt.xlabel('Длина связи C-C')
plt.ylabel('Энергия')

Text(0, 0.5, 'Энергия')


view = py3Dmol.view()
g = geom_template.format(cc_len_real, hcc_angle_real, 0)
m = psi4.geometry(g)
view.addModel(m.save_string_xyz_file(),'xyz')
view.zoomTo()
view.setStyle({'model': -1}, {"stick": {}})

view.show()


x_o = np.arange(-180, 180, 12)
y_o = np.zeros(x_o.shape)

for i, x in enumerate(x_o):
    geom_ = geom_template.format(cc_len_real, hcc_angle_real, x)
    y_o[i] = run_psi4(geom_)


plt.plot(x_o, y_o, '-b.',  [cc_tors_real], [e_real], 'r.', alpha=0.5)
plt.xlabel('Угол CCH')
plt.ylabel('Энергия')
plt.legend()

<matplotlib.legend.Legend at 0x7ff1898ffa60>


# function is  f(x) = a * cos(bx + c) + d

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

p0 = [1, 3.14/180*3, 0, 0]   # Initial guess for the parameters

p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))

print("Optimized params:", p1)

Optimized params: [ 2.47709772e-03  5.23648982e-02  1.21428853e-09 -7.92575446e+01]


#Plot it
plt.plot(x_o, y_o, "ro", c='blue', alpha=0.5)
plt.plot(x_o, fitfunc(p1, x_o), "r-", c='blue', alpha=0.5)
plt.xlabel('Торсион')
plt.ylabel('Энергия')

Text(0, 0.5, 'Энергия')

Длина связи С-С¶

Угол CCH¶

Торсион CC¶