Рудаков Кирилл ДЗ 5

In [4]:
import subprocess
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
%matplotlib inline
In [1]:
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 180
H 2 1 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
In [3]:
def run_orca(inp):
    with open('orca.inp', 'w') as outfile:
        outfile.write(inp)
    p = subprocess.Popen("/home/shad/progs/bin/orca orca.inp", 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]
    out.splitlines()

    # extract energy: FINAL SINGLE POINT ENERGY'
    for i in out.splitlines():
        if i.startswith("FINAL SINGLE POINT ENERGY"):
            return(float(i[len("FINAL SINGLE POINT ENERGY"):]))
    # and return it as float

    return out.splitlines()
In [5]:
run_orca(inp)
Out[5]:
-79.081531418111
In [6]:
# fake x array , replace with real one
x_o=np.arange(1,2,0.1)
# fake y array, replace with energies
y_o=-70+2*(x_o -1.58)**2 

#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: [  2.     1.58 -70.  ]
In [9]:
links = np.arange(-10,11)*0.02 + 1.52986
links
Out[9]:
array([ 1.32986,  1.34986,  1.36986,  1.38986,  1.40986,  1.42986,
        1.44986,  1.46986,  1.48986,  1.50986,  1.52986,  1.54986,
        1.56986,  1.58986,  1.60986,  1.62986,  1.64986,  1.66986,
        1.68986,  1.70986,  1.72986])
In [10]:
gen_inp = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 {} 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 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
In [14]:
orcas = [run_orca(gen_inp.format(x)) for x in links]

Получившиеся значения энергии

In [15]:
# fake x array , replace with real one
x_o=links
# fake y array, replace with energies
y_o=orcas

#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.63991887   1.5569929  -79.08219401]

Аппроксимация приблизительно схожа с колебаниями

Аналогичные операции для валентного угла HCС

In [22]:
angles = np.arange(-10,11)*0.2+111.2
angles
Out[22]:
array([ 109.2,  109.4,  109.6,  109.8,  110. ,  110.2,  110.4,  110.6,
        110.8,  111. ,  111.2,  111.4,  111.6,  111.8,  112. ,  112.2,
        112.4,  112.6,  112.8,  113. ,  113.2])
In [21]:
gen_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 {} 0
H 1 2 3 1.08439 {} 120
H 1 2 3 1.08439 {} -120
H 2 1 3 1.08439 {} 180
H 2 1 5 1.08439 {} 120
H 2 1 5 1.08439 {} -120
*
'''
In [24]:
orcas = [run_orca(gen_inp.format(x)) for x in angles]
orcas
Out[24]:
[-79.078791124966,
 -79.079183399887,
 -79.079549472852,
 -79.07988931153,
 -79.080202891616,
 -79.080490187722,
 -79.080751172794,
 -79.080985818069,
 -79.08119409307,
 -79.081375965653,
 -79.081531402161,
 -79.081660366819,
 -79.081762822554,
 -79.081838730525,
 -79.081888050195,
 -79.081910739351,
 -79.081906754045,
 -79.081876048815,
 -79.081818576377,
 -79.081734287686,
 -79.081623132149]
In [28]:
# fake x array , replace with real one
x_o=angles
# fake y array, replace with energies
y_o=orcas

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

Optimized params: [  3.31057709e-04   1.12270870e+02  -7.90819111e+01]

Как видно из графика, аппроксимация отлично описывает колебания

In [31]:
tors = np.linspace(-180, 180, (180 + 180) / 12 + 1)
tors
Out[31]:
array([-180., -168., -156., -144., -132., -120., -108.,  -96.,  -84.,
        -72.,  -60.,  -48.,  -36.,  -24.,  -12.,    0.,   12.,   24.,
         36.,   48.,   60.,   72.,   84.,   96.,  108.,  120.,  132.,
        144.,  156.,  168.,  180.])
In [36]:
gen_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 {}
H 2 1 5 1.08439 111.200 {}
H 2 1 5 1.08439 111.200 {}
*
'''
In [37]:
orcas = [run_orca(gen_inp.format(x,x-60,x+60)) for x in tors]
orcas
Out[37]:
[-79.081531418111,
 -79.082613733222,
 -79.083665197014,
 -79.084283343101,
 -79.084235490395,
 -79.083542971542,
 -79.082468786938,
 -79.081419140127,
 -79.080793992197,
 -79.080835703972,
 -79.081531425731,
 -79.082613733261,
 -79.083665196902,
 -79.084283343103,
 -79.084235490374,
 -79.083542971543,
 -79.082468786883,
 -79.081419140194,
 -79.080793992237,
 -79.080835704057,
 -79.081531425784,
 -79.082613733277,
 -79.083665197014,
 -79.084283343103,
 -79.084235490391,
 -79.083542971545,
 -79.082468786931,
 -79.081419140068,
 -79.080793992183,
 -79.080835703989,
 -79.081531425729]
In [38]:
# fake x array , replace with real one
x_o=tors
# fake y array, replace with energies
y_o=orcas

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

Optimized params: [  1.33350435e-08  -8.78365661e+01  -7.90827630e+01]

У функции 3 минимума, аппроксимация не описывает колебания ( функция похожа на cos/sin)

In [42]:
# fake x array , replace with real one
x_o=tors
# fake y array, replace with energies
y_o=orcas

#function is  f(x)=k/2*sin(b*x+c) + d
fitfunc = lambda p, x: p[0]/2*np.sin(p[1]*x+p[2]) + p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [1,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(-182,182)
plt.show()

Optimized params: [ -3.62478006e-03   9.94840947e-01   5.88393389e-01  -7.90825389e+01]

Увеличим шаг до 0.1 ангстрема при расчёте связей (C-C).

In [45]:
links = np.arange(-10,11)*0.1 + 1.52986
links
Out[45]:
array([ 0.52986,  0.62986,  0.72986,  0.82986,  0.92986,  1.02986,
        1.12986,  1.22986,  1.32986,  1.42986,  1.52986,  1.62986,
        1.72986,  1.82986,  1.92986,  2.02986,  2.12986,  2.22986,
        2.32986,  2.42986,  2.52986])
In [46]:
gen_inp = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 {} 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 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
In [47]:
orcas = [run_orca(gen_inp.format(x)) for x in links]
orcas
Out[47]:
[-73.376057904506,
 -75.600038674084,
 -76.986925190758,
 -77.840020402069,
 -78.364137030278,
 -78.683902686411,
 -78.875495245128,
 -78.986233920144,
 -79.045965568718,
 -79.07368949118,
 -79.081531419723,
 -79.077218664924,
 -79.065654757297,
 -79.049930677904,
 -79.031982697254,
 -79.013023521422,
 -78.993820247811,
 -78.974866462702,
 -78.956483580267,
 -78.938878296987,
 -78.922175247055]
In [49]:
# fake x array , replace with real one
x_o=links
# fake y array, replace with energies
y_o=orcas

#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(0,3)
plt.show()

Optimized params: [  2.64135253   1.81379259 -79.55943724]

Функция похожа на экспоненциальную

In [58]:
# fake x array , replace with real one
x_o=links
# fake y array, replace with energies
y_o=orcas

#function is  f(x)=k*exp{a*x+b}+c
fitfunc = lambda p, x: p[0]*np.exp(p[1]*x + p[2]) + p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [1,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(0,3)
plt.show()

Optimized params: [  7.38100677  -5.24263473   2.52193897 -79.04492159]
In [ ]:
!jupyter nbconvert --to html rudakov_kirill_hw5.ipynb