Задание 2: Вычисление параметров при изменении значения одного валентного угла H-C-C в диапозоне 109.2-113.2° с шагом 0.2°.
In [2]:
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
*
'''
In [3]:
j=109.2
name="molecula_tri"
for i in range(20):
    out=open(name+str(i),"w")
    out.write(inp %(j))
    out.close()
    j=j+0.2
In [5]:
import subprocess
def energy(fineexe):
    p = subprocess.Popen("/srv/databases/orca/orca "+fineexe, 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]
    out.splitlines()
    for li in out.splitlines():
        data=li.split("FINAL SINGLE POINT ENERGY")
        if len(data)==2:
            return(data[1].strip())
In [6]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [21]:
x_o=np.arange(109.2, 113.2, 0.2)
y_o=[]

for i in range(20):
    name = "molecula_tri"+str(i)
    it = energy(name)
    y_o.append(float(it))

y_o = np.array(y_o)
print x_o    
print y_o

[ 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. ]
[-79.19741049 -79.1974414  -79.19746904 -79.19749339 -79.19751448
 -79.19753229 -79.19754683 -79.19755811 -79.19756613 -79.19757089
 -79.1975724  -79.19757067 -79.19756568 -79.19755745 -79.19754599
 -79.19753128 -79.19751335 -79.19749218 -79.19746779 -79.19744018]
In [22]:
#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.5)

Optimized params: [  4.06431971e-05   1.11195037e+02  -7.91975724e+01]
Out[22]:
(109, 113.5)

Был получен график зависимости энергии соединения от величены угла. Min - 111.57°