6. Вычисление параметров для молекулярной механики¶

In [2]:
import psi4
import numpy as np
psi4.core.set_output_file('output.dat')
In [4]:
def run_psi4(inp):
    m = psi4.geometry(inp)
    psi4.set_options({"maxiter": 200, "fail_on_maxiter" :  True})
    x = psi4.energy('scf/cc-pvtz', molecule = m )
    return x
In [5]:
inp = '''
C
C 1 1.52986
H 1 1.08439 2 111.200
H 1 1.08439 2 111.200 3 120
H 1 1.08439 2 111.200 3 -120
H 2 1.08439 1 111.200 3 180
H 2 1.08439 1 111.200 6 120
H 2 1.08439 1 111.200 6 -120
'''
In [9]:
import py3Dmol
view = py3Dmol.view()
view.addModel(psi4.geometry(inp).save_string_xyz_file(),'xyz')
view.setStyle({'stick':{'colorscheme':'cyanCarbon'}})

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Out[9]:
<py3Dmol.view at 0x13129ec1ed0>
In [10]:
run_psi4(inp)
Out[10]:
-79.26001266653356

Зависимость энергии от длины связи С-С¶

In [11]:
CCl = np.arange(1.28, 1.72, 0.02)
CCl_e = []

inp_CC = '''
C
C 1 {}
H 1 1.08439 2 111.200
H 1 1.08439 2 111.200 3 120
H 1 1.08439 2 111.200 3 -120
H 2 1.08439 1 111.200 3 180
H 2 1.08439 1 111.200 6 120
H 2 1.08439 1 111.200 6 -120
'''

for i in CCl:
    CCl_e.append(run_psi4(inp_CC.format(i)))
In [12]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize 
In [15]:
x_o=CCl
y_o=CCl_e

#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, "o", x_o,fitfunc(p1,x_o),"-",c='blue',alpha=0.5)
plt.xlim(1.25,1.7)
plt.show()
Optimized params: [  0.71579156   1.5420268  -79.26149991]
No description has been provided for this image

Т.е. мы аппроксимировали зависимость с помощью параболы. Вышло, в целом, не так уж хорошо, но и не сильно плохо.

В статье про GAFF указана длина связи 1.526, у нас же вышло 1.542. Возможно, тут дело в не самой удачной аппроксимации...

Зависимость энергии от валентного угла HCC¶

In [20]:
inp_HCC = '''
C
C 1 1.52986
H 1 1.08439 2 {}
H 1 1.08439 2 111.200 3 120
H 1 1.08439 2 111.200 3 -120
H 2 1.08439 1 111.200 3 180
H 2 1.08439 1 111.200 6 120
H 2 1.08439 1 111.200 6 -120
'''
109.5
HCCa = np.arange(109.2, 113.2, 0.2)
HCCa_e = []

for i in HCCa:
    HCCa_e.append(run_psi4(inp_HCC.format(i)))
In [21]:
x_o=HCCa
y_o=HCCa_e

#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, "o", x_o,fitfunc(p1,x_o),"-",c='blue',alpha=0.5)
plt.xlim(109,113.6)
plt.show()
Optimized params: [ 3.80666172e-05  1.11132163e+02 -7.92600128e+01]
No description has been provided for this image

Здесь зависимость энергии от угла ложится на параболу практически идеально

Зависимость энергии от торсионного угла CC¶

In [19]:
inp_CCt = '''
C
C 1 1.52986
H 1 1.08439 2 111.200
H 1 1.08439 2 111.200 3 120
H 1 1.08439 2 111.200 3 -120
H 2 1.08439 1 111.200 3 {}
H 2 1.08439 1 111.200 6 120
H 2 1.08439 1 111.200 6 -120
'''

CCt = np.arange(-180, 181, 12)
CCt_e = []

for i in CCt:
    CCt_e.append(run_psi4(inp_CCt.format(i)))

Помотрим на то, как выглядит зависимость энергии от торсиона.

In [22]:
plt.plot(CCt, CCt_e)
plt.show()
No description has been provided for this image

Попробуем аппроксимировать на синусоиду.

In [24]:
x_o=CCt
y_o=CCt_e

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

p0 = [0, 1, -1, -1] # 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, "o", x_o,fitfunc(p1,x_o),"-",c='blue',alpha=0.5)
plt.xlim(-200,200)
plt.show()
Optimized params: [ 2.47654449e-03  9.55039759e-01  9.68791789e-06 -7.92600209e+01]
No description has been provided for this image

На графике 4 минимума: в -180, -60, 60, 180 - заторможенные конформации.

Опять зависимость энергии от длины связи СС, но с шагом 0.1 ангстрем¶

In [25]:
CCl_new = np.arange(1, 3, 0.1)
CCl_e_new = []

for i in CCl_new:
    CCl_e_new.append(run_psi4(inp_CC.format(i)))
In [27]:
plt.plot(CCl_new, CCl_e_new)
plt.show()
No description has been provided for this image

Получилась совсем не парабола. Зато такая зависимость могла бы аппроксимироваться с помощью функции Морзе.