In [1]:
#Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов.
In [7]:
import psi4
import numpy as np
psi4.core.set_output_file('output.dat')
In [50]:
#z-matrix
inp = '''
0 1
C
C 1 1.52986
H 1 1.08439 2 111.200
H 1 1.08439 2 111.200 3 120.0
H 1 1.08439 2 111.200 3 -120.0
H 2 1.08439 1 111.200 3 180.0
H 2 1.08439 1 111.200 6 120.0
H 2 1.08439 1 111.200 6 -120.0
'''
In [51]:
def run_psi4(inp):
m = psi4.geometry(inp)
psi4.set_options({"maxiter": 200, "fail_on_maxiter" : True})
e = psi4.energy('scf/cc-pvtz', molecule = m)
return e, m
In [52]:
e, m = run_psi4(inp)
In [54]:
import py3Dmol
view = py3Dmol.view()
view.addModel(m.save_string_xyz_file(), 'xyz')
view.setStyle({'stick': {}, 'sphere': {'radius': 0.35}})
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
Out[54]:
<py3Dmol.view at 0x7f1682157e50>
In [55]:
# Построим зависимость энергии молекулы от длины связи
In [56]:
import pandas as pd
In [76]:
d = [1.34 + i * 0.02 for i in range(20)]
df = pd.DataFrame({'distances': d})
def diff_dist(d):
inp = f'''
0 1
C
C 1 {d}
H 1 1.084 2 111.2
H 1 1.084 2 111.2 3 120
H 1 1.084 2 111.2 3 -120
H 2 1.084 1 111.2 3 180
H 2 1.084 1 111.2 6 120
H 2 1.084 1 111.2 6 -120
'''
psi4.core.clean()
e, m = run_psi4(inp)
return e
df['energy'] = df['distances'].apply(diff_dist)
print(df.head())
distances energy 0 1.34 -79.233665 1 1.36 -79.239818 2 1.38 -79.244983 3 1.40 -79.249248 4 1.42 -79.252695
In [58]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [59]:
x_o=df['distances'].values
y_o=df['energy'].values
#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)
#Plot it
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.xlim(1,2)
plt.show()
Optimized params: [ 0.58152001 1.54385331 -79.26026046]
In [ ]:
# Проделаем аналогичные операции для валентного угла HCC
In [42]:
angle = np.linspace(109.2, 113.2, 20)
df2 = pd.DataFrame({'angle': angle})
def diff_angles(angle):
inp = f'''
0 1
C
C 1 1.52986
H 1 1.084 2 {angle}
H 1 1.084 2 {angle} 3 120
H 1 1.084 2 {angle} 3 -120
H 2 1.084 1 {angle} 3 180
H 2 1.084 1 {angle} 6 120
H 2 1.084 1 {angle} 6 -120
'''
psi4.core.clean()
e, m = run_psi4(inp)
return e
df2['energy'] = df2['angle'].apply(diff_angles)
In [43]:
x_o=df2['angle'].values
y_o=df2['energy'].values
#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)
#Plot it
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.show()
Optimized params: [ 3.06678287e-04 1.11140298e+02 -7.92600143e+01]
In [ ]:
#Проделайте аналогичные операции для торсионного угла CC, его значения должны изменяться от -180 до 180 c шагом 12.
In [46]:
tors = [-180 + i * 12 for i in range(31)]
df3 = pd.DataFrame({'tors': tors})
def diff_tors(tors):
inp = f'''
0 1
C
C 1 1.52986
H 1 1.084 2 111.2
H 1 1.084 2 111.2 3 120
H 1 1.084 2 111.2 3 -120
H 2 1.084 1 111.2 3 {tors}
H 2 1.084 1 111.2 3 {tors + 120}
H 2 1.084 1 111.2 3 {tors - 120}
'''
psi4.core.clean()
e, m = run_psi4(inp)
return e
df3['energy'] = df3['tors'].apply(diff_tors)
In [86]:
x_o=df3['tors'].values
y_o=df3['energy'].values
plt.plot(x_o, y_o, c='blue', alpha=0.5)
plt.scatter(x_o, y_o, c='blue', alpha=0.5)
plt.show()
In [ ]:
#Увеличьте шаг до 0.1 ангстрема при расчёте связи.
#Постройте зависимость.
#Укажите какой функцией можно было бы аппроксимировать наблюдаемую зависимость.
In [66]:
d2 = [1.34 + i * 0.1 for i in range(20)]
df4 = pd.DataFrame({'distances': d2})
def diff_dist(d):
inp = f'''
0 1
C
C 1 {d}
H 1 1.084 2 111.2
H 1 1.084 2 111.2 3 120
H 1 1.084 2 111.2 3 -120
H 2 1.084 1 111.2 3 180
H 2 1.084 1 111.2 6 120
H 2 1.084 1 111.2 6 -120
'''
psi4.core.clean()
e, m = run_psi4(inp)
return e
df4['energy'] = df4['distances'].apply(diff_dist)
In [68]:
x_o=df4['distances'].values
y_o=df4['energy'].values
#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)
#Plot it
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.show()
Optimized params: [ 1.27749409e-02 -3.19923151e+00 -7.95345584e+01]
In [73]:
# Из-за большого расстояния связь C-C разрывается
# Парабола больше не подходит для аппроксимации
# Можно использовать потенциал Морзе
In [75]:
morse_func = lambda p, x: p[0] * (1 - np.exp(-p[1] * (x - p[2])))**2 + p[3]
err_morse = lambda p, x, y: morse_func(p, x) - y
p0_m = [0.1, 2.0, 1.53, y_o.min()]
p_m, success = optimize.leastsq(err_morse, p0_m, args=(x_o, y_o))
plt.scatter(x_o, y_o, color='blue', alpha = 0.5)
x_fine = np.linspace(x_o.min(), x_o.max(), 100)
plt.plot(x_fine, morse_func(p_m, x_fine), 'b-', alpha = 0.5)
plt.show()
In [ ]: