Задание 1: Вычисление параметров при изменении значения длины связи C-C в диапазоне 1.33-1.71Å с шагом 0.2Å.
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 [2]:
print inp
data = inp.split("1.52986")
j=1.52986-0.2
print j
name="molecula"
for i in range(20):
    out=open(name+str(i),"w")
    out.write(data[0]+str(j)+data[1])
    out.close()
    j=j+0.02

!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
*

1.32986
In [9]:
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:
            #print data[1].strip()
            return(data[1].strip())
In [10]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [11]:
x_o=np.arange(1.33, 1.73, 0.02)
y_o=[]

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

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

[ 1.33  1.35  1.37  1.39  1.41  1.43  1.45  1.47  1.49  1.51  1.53  1.55
  1.57  1.59  1.61  1.63  1.65  1.67  1.69  1.71]
[-79.04596559 -79.05360305 -79.06009493 -79.06554486 -79.07004752
 -79.07368951 -79.07654992 -79.07870106 -79.080209   -79.08113408
 -79.08153143 -79.08145141 -79.08093996 -79.08003906 -79.07878699
 -79.07721868 -79.07536598 -79.07325792 -79.07092099 -79.06837928]
In [12]:
#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)

Optimized params: [  0.66824655   1.55403272 -79.08233322]
Out[12]:
(1, 2)

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