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

In [1]:
import psi4
import numpy as np
psi4.core.set_output_file('output.dat')
In [2]:
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 [3]:
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 [4]:
import py3Dmol
view = py3Dmol.view()
view.addModel(psi4.geometry(inp).save_string_xyz_file(),'xyz')
Out[4]:
<py3Dmol.view at 0x7f992aba3b20>

Выглядит не очень (еле заметные линии, но пусть будет так)

In [5]:
run_psi4(inp)
Out[5]:
-79.26001266654337

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

Возьмем длину от 1.30 до 1.70, шаг 0.02 Ангстрем

In [11]:
CClengths = np.arange(1.30, 1.70, 0.02)
CClengths_energies = []
In [12]:
inp_CClength = '''
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 CClengths:
    CClengths_energies.append(run_psi4(inp_CClength.format(i)))
In [13]:
CClengths_energies
Out[13]:
[-79.21801576653898,
 -79.22644049209012,
 -79.23367565652505,
 -79.23982689782713,
 -79.24499065821256,
 -79.24925504732664,
 -79.25270061398953,
 -79.25540103771921,
 -79.2574237473626,
 -79.2588304762864,
 -79.25967776108736,
 -79.2600173904389,
 -79.25989681025317,
 -79.25935949020457,
 -79.25844525648515,
 -79.2571905944026,
 -79.25562892436501,
 -79.25379085425101,
 -79.25170441068752,
 -79.24939525100939]
In [14]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [19]:
x_o=CClengths
y_o=CClengths_energies

#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.25,1.7)
plt.show()

Optimized params: [  0.7068083    1.53860038 -79.26115339]

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

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

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

Здесь диапазон будет от 109.2 до 113.2. Шаг возьмем 0.2, тогда мы получим 20 разных значений

In [21]:
inp_HCCangle = '''
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
'''

HCCangles = np.arange(109.2, 113.2, 0.2)
HCCangles_energies = []
In [22]:
for i in HCCangles:
    HCCangles_energies.append(run_psi4(inp_HCCangle.format(i)))
In [23]:
HCCangles_energies
Out[23]:
[-79.2598705926487,
 -79.25989858162866,
 -79.25992349958007,
 -79.25994534962231,
 -79.25996413489553,
 -79.25997985846246,
 -79.25999252341924,
 -79.26000213274845,
 -79.2600086894439,
 -79.26001219656779,
 -79.26001266653942,
 -79.26001007384367,
 -79.26000444989225,
 -79.25999578811832,
 -79.2599840913957,
 -79.25996936264539,
 -79.25995160464004,
 -79.25993082032004,
 -79.25990701237967,
 -79.25988018372784]
In [25]:
x_o=HCCangles
y_o=HCCangles_energies

#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(109,113.6)
plt.show()

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

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

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

Значения будем менять от -180 до 180, шаг 12

In [27]:
inp_CCtorsion = '''
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
'''

CCtorsions = np.arange(-180, 181, 12)
CCtorsions_energies = []
In [28]:
for i in CCtorsions:
    CCtorsions_energies.append(run_psi4(inp_CCtorsion.format(i)))
In [29]:
CCtorsions_energies
Out[29]:
[-79.26001266654059,
 -79.25954606531678,
 -79.25831799033583,
 -79.25678696769599,
 -79.25553784348243,
 -79.25505826172113,
 -79.25553784347537,
 -79.25678696768469,
 -79.25831799033557,
 -79.25954606532011,
 -79.26001266653194,
 -79.25954606530043,
 -79.25831799032865,
 -79.25678696771084,
 -79.25553784347755,
 -79.25505826171099,
 -79.25553784348297,
 -79.25678696770294,
 -79.25831799033637,
 -79.25954606531977,
 -79.26001266652324,
 -79.25954606531015,
 -79.25831799034778,
 -79.25678696770709,
 -79.25553784347824,
 -79.25505826172416,
 -79.2555378434688,
 -79.25678696770235,
 -79.25831799034364,
 -79.2595460653219,
 -79.26001266654333]
In [40]:
plt.plot(CCtorsions, CCtorsions_energies)
plt.show()

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

In [35]:
x_o=CCtorsions
y_o=CCtorsions_energies

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, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.xlim(-200,200)
plt.show()

Optimized params: [ 2.47656780e-03  9.55040101e-01 -3.36858548e-06 -7.92600208e+01]

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

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

In [36]:
CClengths_new = np.arange(1, 3, 0.1)
CClengths_energies_new = []
In [37]:
for i in CClengths_new:
    CClengths_energies_new.append(run_psi4(inp_CClength.format(i)))
In [38]:
CClengths_energies_new
Out[38]:
[-78.83396995656011,
 -79.03621534656376,
 -79.15357633334051,
 -79.21801576656162,
 -79.24925504730011,
 -79.25967776108223,
 -79.25719059439633,
 -79.24688685833299,
 -79.23204753165604,
 -79.21477636767278,
 -79.19641550676525,
 -79.1778171025864,
 -79.15951735803988,
 -79.14184587002778,
 -79.12499449302214,
 -79.10906232032714,
 -79.09408657328227,
 -79.08006440488894,
 -79.06696829306388,
 -79.05475661660205]
In [43]:
plt.plot(CClengths_new, CClengths_energies_new)
plt.show()

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