In [ ]:
Занятие 6. Вычисление параметров для молекулярной механики
In [1]:
#Definig orca function
def run_orca(inp):
    with open('orca.inp', 'w') as outfile:
        outfile.write(inp)
    p = subprocess.Popen("/srv/databases/orca/orca orca.inp", 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]
    outlines=out.splitlines()

    # extract energy: FINAL SINGLE POINT ENERGY'
    # and return it as float
    for i in outlines:
        if "FINAL SINGLE POINT ENERGY" in i:
            return float(i.split()[4])
In [ ]:
Теперь рассчитаем энергии для 20 разных длинн связи с шагом 0.02 ангстрема для расчёта энергии в ORCA с разными значениями по длине одной из связей (C-C).
In [2]:
import subprocess, os
bond=1.52986
x_axis=[]
y_axis=[]

print "bond\tenergy"
for i in range(-10,10):
    bond = 1.52986 + i*0.02
    x_axis.append(bond) 
    inp = '''!HF RHF 6-31G
    * int 0 1
    C 0 0 0 0 0 0 
    C 1 0 0 %f 0 0 
    H 1 2 0 1.08439 111.200 0
    H 1 2 3 1.08439 111.200 120
    H 1 2 3 1.08439 111.200 -120
    H 2 1 3 1.08439 111.200 180
    H 2 1 6 1.08439 111.200 120
    H 2 1 6 1.08439 111.200 -120
    *
    ''' % bond
    energy = run_orca(inp)
    print str(bond)+"\t"+str(energy)
    y_axis.append(energy)
bond	energy
1.32986	-79.1646310741
1.34986	-79.1718840752
1.36986	-79.1780232087
1.38986	-79.1831502051
1.40986	-79.187357994
1.42986	-79.1907314465
1.44986	-79.1933480709
1.46986	-79.1952786405
1.48986	-79.1965877554
1.50986	-79.1973343976
1.52986	-79.1975723852
1.54986	-79.1973508238
1.56986	-79.1967145065
1.58986	-79.1957042783
1.60986	-79.1943573951
1.62986	-79.1927077537
1.64986	-79.1907862852
1.66986	-79.1886211376
1.68986	-79.186237935
1.70986	-79.1836599929
In [ ]:
Построим график и линию тренда
In [3]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

x_o=np.array(x_axis)
y_o=np.array(y_axis)
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, "bo", x_o,fitfunc(p1,x_o),"g-",alpha=0.5)
plt.xlim(1.3,1.73)
plt.ylabel('Energy, eV')
plt.xlabel('C-C bond length, angstrem')
plt.show()
Optimized params: [  0.64460622   1.54848049 -79.19820325]
In [ ]:
Сделаем то же самое, но для валентного угла HCС. Его значения должны изменяться от 109.2 до 113.2
In [8]:
x_axis=[]
y_axis=[]

print "angle\tenergy"
for i in range(-10,10):
    angle = 111.2 + i*0.2
    x_axis.append(angle) 
    inp = '''!HF RHF 6-31G
    * int 0 1
    C 0 0 0 0 0 0 
    C 1 0 0 1.52986 0 0 
    H 1 2 0 1.08439 %f 0
    H 1 2 3 1.08439 111.200 120
    H 1 2 3 1.08439 111.200 -120
    H 2 1 3 1.08439 111.200 180
    H 2 1 6 1.08439 111.200 120
    H 2 1 6 1.08439 111.200 -120
    *
    ''' % angle
    energy = run_orca(inp)
    print str(angle)+"\t"+str(energy)
    y_axis.append(energy)
angle	energy
109.2	-79.1974104866
109.4	-79.1974413122
109.6	-79.197468918
109.8	-79.1974932689
110.0	-79.1975143495
110.2	-79.1975321602
110.4	-79.1975467038
110.6	-79.1975579838
110.8	-79.1975660039
111.0	-79.1975707674
111.2	-79.1975722779
111.4	-79.1975705389
111.6	-79.1975655538
111.8	-79.1975573261
112.0	-79.1975458591
112.2	-79.1975311562
112.4	-79.1975132207
112.6	-79.1974920559
112.8	-79.1974676651
113.0	-79.1974400515
In [9]:
x_o=np.array(x_axis)
y_o=np.array(y_axis)
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, "bo", x_o,fitfunc(p1,x_o),"g-",alpha=0.5)
plt.xlim(108,115)
plt.ylabel('Energy, eV')
plt.xlabel('Angle')
plt.show()
Optimized params: [  4.06308675e-05   1.11194916e+02  -7.91975723e+01]
In [ ]:
И тепеь то же самое, но для торсионного угла CC, его значения должны изменяться от -180 до 180 c шагом 12.
In [15]:
x_axis=[]
y_axis=[]

print "dihedral\tenergy"
for i in range(-15,16):
    dihedral = 0 + i*12
    x_axis.append(dihedral) 
    inp = '''!HF RHF 6-31G
    * int 0 1
    C 0 0 0 0 0 0 
    C 1 0 0 1.52986 0 0 
    H 1 2 0 1.08439 111.200 0
    H 1 2 3 1.08439 111.200 120
    H 1 2 3 1.08439 111.200 -120
    H 2 1 3 1.08439 111.200 %f
    H 2 1 6 1.08439 111.200 120
    H 2 1 6 1.08439 111.200 -120
    *
    ''' % dihedral
    energy = run_orca(inp)
    print str(dihedral)+"\t"+str(energy)
    y_axis.append(energy)
dihedral	energy
-180	-79.1975724034
-168	-79.1971376991
-156	-79.1959965114
-144	-79.194579711
-132	-79.1934285535
-120	-79.1929877026
-108	-79.1934285529
-96	-79.1945797095
-84	-79.1959965099
-72	-79.1971376989
-60	-79.1975723928
-48	-79.1971376991
-36	-79.1959965114
-24	-79.194579711
-12	-79.1934285535
0	-79.1929877026
12	-79.1934285529
24	-79.1945797095
36	-79.1959965099
48	-79.1971376989
60	-79.1975723928
72	-79.1971376991
84	-79.1959965114
96	-79.194579711
108	-79.1934285535
120	-79.1929877026
132	-79.1934285529
144	-79.1945797095
156	-79.1959965099
168	-79.1971376989
180	-79.1975723928
In [21]:
x_o=np.array(x_axis)
y_o=np.array(y_axis)
#I was forced to think about a new function, the old one looked like ... a power, but here is definetly a periodic one
#It was strange that sin didn't work, but cos did.
fitfunc = lambda p, x: p[0]*np.cos(p[1]*x) + 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, "bo", x_o,fitfunc(p1,x_o),"g-",alpha=0.5)
plt.xlim(-200,200)
plt.ylabel('Energy, eV')
plt.xlabel('C-C bond length, angstrem')
plt.show()
Optimized params: [  2.29205061e-03   9.94834578e-01  -7.91952842e+01]
In [ ]:
Теперь увеличим шаг до 0.1 при расчёте энергии для длин связи.
In [22]:
print "bond\tenergy"
for i in range(-20,20):
    bond = 1.52986 + i*0.01
    x_axis.append(bond) 
    inp = '''!HF RHF 6-31G
    * int 0 1
    C 0 0 0 0 0 0 
    C 1 0 0 %f 0 0 
    H 1 2 0 1.08439 111.200 0
    H 1 2 3 1.08439 111.200 120
    H 1 2 3 1.08439 111.200 -120
    H 2 1 3 1.08439 111.200 180
    H 2 1 6 1.08439 111.200 120
    H 2 1 6 1.08439 111.200 -120
    *
    ''' % bond
    energy = run_orca(inp)
    print str(bond)+"\t"+str(energy)
    y_axis.append(energy)
bond	energy
1.32986	-79.1646310741
1.33986	-79.1684035451
1.34986	-79.1718840913
1.35986	-79.1750863188
1.36986	-79.1780232254
1.37986	-79.1807072209
1.38986	-79.1831502084
1.39986	-79.1853635195
1.40986	-79.1873579996
1.41986	-79.1891440095
1.42986	-79.1907314538
1.43986	-79.1921297963
1.44986	-79.1933480798
1.45986	-79.1943949446
1.46986	-79.1952786479
1.47986	-79.1960070756
1.48986	-79.1965877649
1.49986	-79.1970279164
1.50986	-79.1973344072
1.51986	-79.1975138082
1.52986	-79.1975723951
1.53986	-79.197516162
1.54986	-79.1973508341
1.55986	-79.197081879
1.56986	-79.1967145161
1.57986	-79.1962537335
1.58986	-79.195704289
1.59986	-79.1950707285
1.60986	-79.194357389
1.61986	-79.1935684116
1.62986	-79.1927077488
1.63986	-79.1917791734
1.64986	-79.190786282
1.65986	-79.1897325104
1.66986	-79.188621136
1.67986	-79.1874552828
1.68986	-79.1862379334
1.69986	-79.1849719324
1.70986	-79.1836599913
1.71986	-79.1823046971
In [27]:
x_o=np.array(x_axis)
y_o=np.array(y_axis)
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, "bo", x_o,fitfunc(p1,x_o),"g-",alpha=0.5)
plt.xlim(1.3,1.73)
plt.ylabel('Energy, eV')
plt.xlabel('C-C bond length, angstrem')
plt.show()
x_o_1=np.array(x_axis)
y_o_1=np.array(y_axis)
Optimized params: [ -2.34105073e-07   4.25627517e-01  -7.91907893e+01]
In [ ]:
При построении линии тренда выходит что-то странное. Я построила в Excel то же самое и там точки прекрасно аппроксимируется линией тренда y = 6E-05x2 - 0,0029x - 79,165. Но вдруг получится аппроксимировать потенциалом Морзе? Потенциал Морзе  функция потенциальной энергии элетростатического поля, предложенная американским физиком Филипом Морзе, как приближение для энергии двухатомной молекулы.  
In [26]:
x_o=np.array(x_axis)
y_o=np.array(y_axis)

#function is  f(x)=(1-e^(-k(x-a)))^2
fitfunc = lambda p, x: p[0]*pow(1-np.e**(-p[1]*(x-p[2])),2) # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [-79,45,0.5] # 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.title('Energy(bond)')
plt.show()
/srv/databases/orca/miniconda2/lib/python2.7/site-packages/ipykernel/__main__.py:5: RuntimeWarning: overflow encountered in power
/srv/databases/orca/miniconda2/lib/python2.7/site-packages/ipykernel/__main__.py:5: RuntimeWarning: overflow encountered in square
Optimized params: [-79.   45.    0.5]
In [29]:
print "bond\tenergy"
for i in range(-10,10):
    bond = 1.52986 + i*0.01
    x_axis.append(bond) 
    inp = '''!HF RHF 6-31G
    * int 0 1
    C 0 0 0 0 0 0 
    C 1 0 0 %f 0 0 
    H 1 2 0 1.08439 111.200 0
    H 1 2 3 1.08439 111.200 120
    H 1 2 3 1.08439 111.200 -120
    H 2 1 3 1.08439 111.200 180
    H 2 1 6 1.08439 111.200 120
    H 2 1 6 1.08439 111.200 -120
    *
    ''' % bond
    energy = run_orca(inp)
    print str(bond)+"\t"+str(energy)
    y_axis.append(energy)
bond	energy
1.42986	-79.1907314552
1.43986	-79.1921297964
1.44986	-79.1933480798
1.45986	-79.1943949446
1.46986	-79.1952786479
1.47986	-79.1960070756
1.48986	-79.1965877649
1.49986	-79.1970279164
1.50986	-79.1973344072
1.51986	-79.1975138082
1.52986	-79.1975723951
1.53986	-79.197516162
1.54986	-79.1973508341
1.55986	-79.197081879
1.56986	-79.1967145161
1.57986	-79.1962537335
1.58986	-79.195704289
1.59986	-79.1950707285
1.60986	-79.194357389
1.61986	-79.1935684116
In [30]:
x_o=np.array(x_axis)
y_o=np.array(y_axis)
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,0.029, -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, "bo", x_o,fitfunc(p1,x_o),"g-",alpha=0.5)
plt.xlim(1.3,1.73)
plt.ylabel('Energy, eV')
plt.xlabel('C-C bond length, angstrem')
plt.show()
Optimized params: [ -2.18102970e-07   3.03609901e-01  -7.91911209e+01]
In [ ]: