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]
Т.е. мы аппроксимировали зависимость с помощью параболы. Вышло, в целом, не так уж хорошо, но и не сильно плохо.
В статье про 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]
Здесь зависимость энергии от угла ложится на параболу практически идеально
Зависимость энергии от торсионного угла 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()
Попробуем аппроксимировать на синусоиду.
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]
На графике 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()
Получилась совсем не парабола. Зато такая зависимость могла бы аппроксимироваться с помощью функции Морзе.