In [27]:
import subprocess
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize 

%matplotlib inline
In [18]:
inp_default = '''!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 [19]:
inp = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 {length} 0 0 
H 1 2 0 1.08439 {angle} 0
H 1 2 3 1.08439 {angle} 120
H 1 2 3 1.08439 {angle} -120
H 2 1 3 1.08439 {angle} {torsion1}
H 2 1 5 1.08439 {angle} {torsion2}
H 2 1 5 1.08439 {angle} {torsion3}
*
'''

Исходные значения

In [37]:
default_length=1.52986
default_angle=111.2
default_torsion1=180
default_torsion2=120
default_torsion3=-120
difference=0.02 #шаг 0.02 ангстрема
dif_count=20 #количество разных длин
min_torsion=-180
max_torsion=180
steps=12

run_orca

In [10]:
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]
    energy_line = filter(lambda x: 'FINAL SINGLE POINT ENERGY' in x, out.splitlines())[0]
    return float(energy_line.split()[-1])

Функции для изменений

In [62]:
def append_length(lengths,default_angle,default_torsion1,default_torsion2,default_torsion3):
    length_energies = []
    for length in lengths:
        length_energies.append(run_orca(inp.format(
            length = length,
            angle = default_angle,
            torsion1 = default_torsion1,
            torsion2 = default_torsion2,
            torsion3 = default_torsion3,
        )))
    return length_energies
In [61]:
def append_angle(angles,default_length,default_torsion1,default_torsion2,default_torsion3):
    angle_energies = []
    for angle in angles:
        angle_energies.append(run_orca(inp.format(
            length = default_length,
            angle = angle,
            torsion1 = default_torsion1,
            torsion2 = default_torsion2,
            torsion3 = default_torsion3,
        )))
    return angle_energies
In [60]:
def append_torsion(torsions,default_length,default_angle):
    torsion_energies = []
    for torsion in torsions:
        torsion_energies.append(run_orca(inp.format(
            length = default_length,
            angle = default_angle,
            torsion1 = torsion,
            torsion2 = torsion - 60,
            torsion3 = torsion + 60,
        )))
    return torsion_energies

Изменение длины связи

In [50]:
lengths = np.arange(-dif_count/2, dif_count/2+1)*difference+default_length
length_energies=append_length(lengths,default_angle,default_torsion1,default_torsion2,default_torsion3)
In [51]:
print(length_energies)

[-79.045965590033, -79.053603047453, -79.060094932942, -79.06554485476, -79.070047523843, -79.073689505281, -79.076549916355, -79.078701061159, -79.080208995167, -79.081134076599, -79.081531431378, -79.081451404937, -79.080939963093, -79.080039062387, -79.078786990807, -79.077218675992, -79.075365973146, -79.073257919932, -79.070920984592, -79.068379281351, -79.065654773152]

Изменение угла валентности (angle)

In [52]:
angles = np.arange(-dif_count/2, dif_count/2+1)*difference*10 + default_angle
angle_energies=append_angle(angles,default_length,default_torsion1,default_torsion2,default_torsion3)
In [53]:
print(angle_energies)

[-79.078791138209, -79.079183400066, -79.079549472821, -79.079889311529, -79.080202891616, -79.080490187722, -79.080751172795, -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]

Изменение торсионного угла (torsion)

In [63]:
torsions = np.arange(min_torsion, max_torsion+1, steps)
torsion_energies=append_torsion(torsions,default_length,default_angle)
In [64]:
print(torsion_energies)

[-79.081531432144, -79.082613733241, -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.083542971544, -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]

Визуализация

length_energies(lengths)

In [69]:
x_o=lengths
y_o=length_energies

#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

plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.title('length_energies(lengths)')
plt.xlabel('lengths')
plt.ylabel('length_energies')
#plt.xlim(1,2)
plt.show()

Optimized params: [  0.63991838   1.55699292 -79.082194  ]

angles_energies(angles)

In [71]:
x_o=angles
y_o=angle_energies

#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

plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.title('angles_energies(angles)')
plt.xlabel('angles')
plt.ylabel('angles_energies')
#plt.xlim(1,2)
plt.show()

Optimized params: [  3.31056886e-04   1.12270871e+02  -7.90819111e+01]

torsion_energies(torsions)

In [72]:
x_o=torsions
y_o=torsion_energies

#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

plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.title('torsion_energies(torsions)')
plt.xlabel('torsions')
plt.ylabel('torsion_energies')
#plt.xlim(1,2)
plt.show()

Optimized params: [  1.33534816e-08  -8.75399665e+01  -7.90827626e+01]
Проблема. Попробуем другой вид функции (линейную комбинацию синусов и косинусов)
In [75]:
x_o=torsions
y_o=torsion_energies

fitfunc = lambda p, x: p[0]*np.sin(p[1]*x)+p[3]*np.cos(p[4]*x) + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [1,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
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.title('torsion_energies(torsions)')
plt.xlabel('torsions')
plt.ylabel('torsion_energies')
#plt.xlim(1,2)
plt.show()
plt.show()

Optimized params: [ -1.50760512e-03   9.94841117e-01  -7.90825389e+01  -1.00591822e-03
  -7.90110580e+01]

Наблюдаем три локальных минимума.

Увеличение шага до 0.1 ангстрема

In [76]:
difference=0.1
In [77]:
lengths_higher_difference = np.arange(-dif_count/2, dif_count/2+1)*difference+default_length
length_energies_higher_difference=append_length(lengths_higher_difference,default_angle,default_torsion1,default_torsion2,default_torsion3)
In [78]:
print (length_energies_higher_difference)

[-73.376057904506, -75.600038674083, -76.986925190758, -77.840020402069, -78.364137030277, -78.683902686411, -78.875495245128, -78.986233920144, -79.045965568718, -79.07368949118, -79.081531419723, -79.077218664924, -79.065654757297, -79.049930677904, -79.031982697254, -79.013023521421, -78.993820247811, -78.974866462701, -78.956483580267, -78.938878296987, -78.922175247055]
In [81]:
x_o=lengths_higher_difference
y_o=length_energies_higher_difference

fitfunc = lambda p, x: pow(x,p[1])+p[0] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [1, -79] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print "Optimized params:", p1

# Plot
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.title('length_energies(lengths)')
plt.xlabel('lengths')
plt.ylabel('length_energies')
#plt.xlim(1,2)
plt.show()

Optimized params: [-79.31281181  -2.79877937]
Функция $x^{-2.78}-79.31$. Не самая лучшая аппроксимация, но близко.