Задание 3: Вычисление параметров при изменении значения одного торсионного угла (от -180° до 180° с шагом 12°)
In [3]:
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 %i
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
In [16]:
for i in range(-180,181, 12):
    #l = 111.2 + i*0.02
    name = 'mol'+str(i)+'_tor.inp'
    out_file = open (name, 'w')
    out_file.write(inp %(i))
    out_file.close()
In [17]:
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 [18]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [21]:
x_o=np.arange(-180,181, 12)
y=[]
for i in range(-180,181, 12):
    name = 'mol'+str(i)+'_tor.inp'
    en = energy(name)
    y.append(float(en))
y_o = np.array(y)
In [24]:
#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",c='blue',alpha=0.5)
plt.xlim(-170,190)
Out[24]:
(-170, 190)

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