В этом таске мы вертим, растягиваем и гнём молекулу этана, и смотрим, как у неё от этого меняется энергия.
Издеваться над молекулой мы будем тремя способами:

  • Варьируя длину связи C—C;
  • Варьируя угол связи HCC (общий для всех водородов);
  • Вращая две тройки водородов относительно друг друга.

Энергия изменённой молекулы будет считаться с помощью ORCA.
Значения энергии мы сохраним и нарисуем как функции вышеозначенных переменных с помощью matplotlib.

Начнём сначала: мы будем генерировать входные файлы для ORCA по шаблону. Шаблон и стандартные значения — вот:

In [2]:
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}
*
'''
length_base = 1.52986
angle_base = 111.2
torsion1_base = 180
# torsion2_base = torsion1_base - 60
# torsion3_base = torsion1_base + 60

Функция полного цикла работы с ORCA:

In [20]:
def run_orca(inp):
    import subprocess
    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]
    for s in out.splitlines():
        if s.startswith("FINAL SINGLE POINT ENERGY"):
            return(float(s[len("FINAL SINGLE POINT ENERGY"):]))

    # extract energy: FINAL SINGLE POINT ENERGY'
    # and return it as float

    return out.splitlines()

Циклы перебора значений.

  • Длина связи CC варьируется от -0.2Å до +0.2Å к стандартной, с шагом 0.02;

  • Угол CCH варьируется от -2° до +2° к стандартному, с шагом 0.2;

  • Торсионный угол torsion1 варьируется от -180° до 180° с шагом 12°, torsion2 и torsion3 подкручиваются вместе с ним.
In [36]:
import numpy as np
In [33]:
length_energies = []
lengths = np.arange(-10, 10+1)*0.02 + length_base
for length in lengths:
    length_energies.append(run_orca(inp.format(
        length = length,
        angle = angle_base,
        torsion1 = torsion1_base,
        torsion2 = torsion1_base - 60,
        torsion3 = torsion1_base + 60,
    )))
print("Lengths done")
angle_energies = []
angles = np.arange(-10, 10+1)*0.2 + angle_base
for angle in angles:
    angle_energies.append(run_orca(inp.format(
        length = length_base,
        angle = angle,
        torsion1 = torsion1_base,
        torsion2 = torsion1_base - 60,
        torsion3 = torsion1_base + 60,
    )))
print("Angles done")
torsion_energies = []
torsions = np.arange(-180, 180+1, 12)
for torsion in torsions:
    torsion_energies.append(run_orca(inp.format(
        length = length_base,
        angle = angle_base,
        torsion1 = torsion,
        torsion2 = torsion - 60,
        torsion3 = torsion + 60,
    )))
print("Torsions done")
print("All done")

Lengths done
Angles done
Torsions done
All done

Импортируем matplotlib (и scipy.optimize).

In [23]:
import matplotlib.pyplot as plt
%matplotlib inline
from scipy import optimize

Рисуем графики зависимости энергии от.
Интереса ради аппроксимируем эти графики параболами.

In [43]:
variables = [lengths,angles,torsions]
energies = [length_energies, angle_energies, torsion_energies]
for x_o, y_o in zip(variables, 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("Formula: k(x-b)^2 + a")
    print "Optimized params:", p1
    print "R-squared:", 1-np.mean(errfunc(p1, x_o, y_o)**2)/np.var(y_o)


    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()

Formula: k(x-b)^2 + a
Optimized params: [  0.63991807   1.55699293 -79.082194  ]
R-squared: 0.976335735993

Formula: k(x-b)^2 + a
Optimized params: [  3.31056886e-04   1.12270871e+02  -7.90819111e+01]
R-squared: 0.999998813909

Formula: k(x-b)^2 + a
Optimized params: [  1.33569535e-08  -8.76360807e+01  -7.90827629e+01]
R-squared: 0.0505754744814

О качестве аппроксимаций:

  • Зависимость от длины связи приближается параболой неплохо, но не идеально: невооружённым взглядом видно, что подлежащая функция несимметрична.
  • Зависимость от угла HCC приближается идеально, ничего менять не надо.
  • Зависимость от торсионного угла, как и ожидалось, имеет период 120°, а посему параболой приближается из рук вон плохо. Взять хотя бы три минимума, которые встречаются от -180° до +180°.

Посчитаем зависимость от длины связи на большем диапазоне (с шагом 0.1) и приблизим более хитрыми потенциалами.

In [73]:
length_energies = []
lengths = np.arange(-10, 10+1)*0.1 + length_base
for length in lengths:
    length_energies.append(run_orca(inp.format(
        length = length,
        angle = angle_base,
        torsion1 = torsion1_base,
        torsion2 = torsion1_base - 60,
        torsion3 = torsion1_base + 60,
    )))
print("Lengths done")
print("All done")

Lengths done
All done

Пробуем потенциал Морзе.

In [102]:
x_o = lengths
y_o = length_energies
fitfunc = lambda p, x: p[0]+p[1]*(1-np.exp(-p[2]*(x-p[3])))**2 # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [-79, 1, 2, 1.5] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print("Formula: E_0+D_e*(1-e^(-a(l-l_0)))")
print "Optimized params:", p1
print "R-squared:", 1-np.mean(errfunc(p1, x_o, y_o)**2)/np.var(y_o)


plt.plot(x_o, y_o, "ro", 
         x_o,fitfunc(p1,x_o),"r-",
         c='blue',alpha=0.5)
#plt.plot(x_o,fitfunc(p0,x_o),"r-",c='red',alpha=0.5)
#plt.xlim(1,2)
plt.show()

Formula: E_0+D_e*(1-e^(-a(l-l_0)))
Optimized params: [-79.06965375   0.15032109   2.05354242   1.48749283]
R-squared: 0.999872305235

Хорошая подгонка! Радиус не совсем точен (1.49 потив истинных 1.53), но это простительно: точность ограничена разрешением в 0.1Å.