In [2]:
import psi4
import numpy as np
In [3]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize 
In [4]:
from IPython.display import Image
In [5]:
psi4.core.set_output_file('output.dat')

Ссылка на файл output.dat

In [6]:
import py3Dmol
In [7]:
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
In [8]:
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
'''
In [9]:
cc_len_real = 1.52986
hcc_angle_real = 111.2
cc_tors_real = -180
In [10]:
geom_ = geom_template.format(cc_len_real, hcc_angle_real, cc_tors_real)
e_real = run_psi4(geom_)
e_real
Out[10]:
-79.26001266654059
In [10]:
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()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

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

Проверим, что меняем действительно нужный параметр:

In [11]:
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()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

In [12]:
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
In [13]:
plt.plot(x_o, y_o, '-b.',  [cc_len_real], [e_real], 'r.', alpha=0.5)
plt.xlabel('Длина связи C-C')
plt.ylabel('Энергия')
plt.legend()
Out[13]:
<matplotlib.legend.Legend at 0x7ff189b079a0>

Попробуем зафитировать параболой:

In [14]:
# 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]
In [15]:
#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()

Фиттируется не очень. Если увеличить масштаб, видно, что это совсем не парабола:

In [16]:
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
In [17]:
plt.plot(x_o, y_o, '-b.',  [cc_len_real], [e_real], 'r.', alpha=0.5)
plt.xlabel('Длина связи C-C')
plt.ylabel('Энергия')
plt.legend()
Out[17]:
<matplotlib.legend.Legend at 0x7ff189a82760>

Она должна описываться ангармоническим потенциалом Морзе

In [18]:
Image(url='https://upload.wikimedia.org/wikipedia/commons/thumb/5/52/Morse-potential-ru.svg/1024px-Morse-potential-ru.svg.png', width=600)
Out[18]:
In [19]:
# 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 ]
In [20]:
#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()

Но фит всё ещё не очень хороший в левой части

Угол CCH¶

Проверим, что меняем действительно нужный параметр:

In [21]:
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()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

In [22]:
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_) 
In [23]:
plt.plot(x_o, y_o, '-b.',  [hcc_angle_real], [e_real], 'r.', alpha=0.5)
plt.xlabel('Угол C-C')
plt.ylabel('Энергия')
plt.legend()
Out[23]:
<matplotlib.legend.Legend at 0x7ff1898b6940>
In [24]:
# 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]
In [25]:
#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('Энергия')
Out[25]:
Text(0, 0.5, 'Энергия')

Торсион CC¶

Проверим, что меняем действительно нужный параметр:

In [26]:
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()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

In [11]:
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_) 
In [28]:
plt.plot(x_o, y_o, '-b.',  [cc_tors_real], [e_real], 'r.', alpha=0.5)
plt.xlabel('Угол CCH')
plt.ylabel('Энергия')
plt.legend()
Out[28]:
<matplotlib.legend.Legend at 0x7ff1898ffa60>
In [30]:
# 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]
In [31]:
#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('Энергия')
Out[31]:
Text(0, 0.5, 'Энергия')

Три минимума соответсвуют трём заторможенным конформациям