HW5¶

In [1]:
%matplotlib inline
import subprocess
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
import math

Зададим общий вид матрицы, в которой будет возможно изменять три позиции:

* длину связи (length)
* угол связи (angle)
* торсионный угол (torsion)

И запишем их начальные значения для дальнейшего использования

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=1.52986
#angle=111.2
#torsion1=180
#torsion2=120
#torsion3=-120
In [3]:
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]
    for i in out.splitlines():
        if i.startswith("FINAL SINGLE POINT ENERGY"):
            return(float(i[len("FINAL SINGLE POINT ENERGY"):]))

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

    return out.splitlines()

Необходимые расчеты¶

1. Изменяем длину связи:¶

In [4]:
length_energies = []
lengths = np.arange(-10, 10+1)*0.02 + 1.52986
print lengths

[ 1.32986  1.34986  1.36986  1.38986  1.40986  1.42986  1.44986  1.46986
  1.48986  1.50986  1.52986  1.54986  1.56986  1.58986  1.60986  1.62986
  1.64986  1.66986  1.68986  1.70986  1.72986]
In [5]:
for length in lengths:
    length_energies.append(run_orca(inp.format(
        length = length,
        angle = 111.2,
        torsion1 = 180,
        torsion2 = 120,
        torsion3 = -120,
    )))

Энергия, соответствующая дефолтным значениям углов и различным значениям длины связи:

In [6]:
print length_energies

[-79.045965577984, -79.053603047534, -79.060094932944, -79.06554485476, -79.070047523843, -79.073689505282, -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]

2. Изменяем валентный угол:¶

In [7]:
angle_energies = []
angles = np.arange(-10, 10+1)*0.2 + 111.2
print angles

[ 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.
  113.2]
In [8]:
for angle in angles:
    angle_energies.append(run_orca(inp.format(
        length = 1.52986,
        angle = angle,
        torsion1 = 180,
        torsion2 = 120,
        torsion3 = -120,
    )))

Энергия, соответствующая различным значениям валентного угла, но дефолтным значениям длины и торсионных углов:

In [9]:
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]

3. Изменяем торсионный угол:¶

In [10]:
torsion_energies = []
torsions = np.arange(-180, 180+1, 12)
print torsions

[-180 -168 -156 -144 -132 -120 -108  -96  -84  -72  -60  -48  -36  -24  -12
    0   12   24   36   48   60   72   84   96  108  120  132  144  156  168
  180]
In [11]:
for torsion in torsions:
    torsion_energies.append(run_orca(inp.format(
        length = 1.52986,
        angle = 111.2,
        torsion1 = torsion,
        torsion2 = torsion - 60,
        torsion3 = torsion + 60,
    )))

Энергия, соответствующая различным знаечниям торсионного угла, при дефолтных значениях длины связи и валентного угла:

In [12]:
print torsion_energies

[-79.081531425441, -79.082613733243, -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]

Графики¶

1) Зависимость энергии от длины связи

In [13]:
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
print 'Function of approximation f(x)=k(b-x)^2 + a '

plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.title('Length vs Energie')
plt.xlabel('length')
plt.ylabel('energie')
plt.show()

Optimized params: [  0.63991836   1.55699293 -79.082194  ]
Function of approximation f(x)=k(b-x)^2 + a
Вывод по приближению: визуально аппроксимация неплохая¶

2) Зависимость энергии от валентного угла

In [14]:
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
print 'Function of approximation f(x)=k(b-x)^2 + a '

plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.title('Angle vs Energie')
plt.xlabel('angle')
plt.ylabel('energie')
plt.show()

Optimized params: [  3.31056886e-04   1.12270871e+02  -7.90819111e+01]
Function of approximation f(x)=k(b-x)^2 + a
Вывод по приближению: аппроксимация отличная¶

3) Зависимость энергии от торсионного угла

In [15]:
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
print 'Function of approximation f(x)=k(b-x)^2 + a '

plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.title('Torsion vs Energie')
plt.xlabel('torsion')
plt.ylabel('energie')
plt.show()

Optimized params: [  1.33284659e-08  -8.78147603e+01  -7.90827628e+01]
Function of approximation f(x)=k(b-x)^2 + a
Вывод по приближению: аппроксимация квадратичной функцией совершенно не подходит, вероятно подошла бы аппроксимация какой-нибудь периодичной функцией типа cos() или sin()¶
На заданном промежутке у функции три минимума¶

Можно попробовать провести аппроксимацию линейной комбинацией тригонометрических функций с константой:

$ f(x)=\alpha_0*\sin(\alpha_1x)+\alpha_2*\cos(\alpha_3x)+\alpha_4 $

In [16]:
x_o=torsions
y_o=torsion_energies

fitfunc = lambda p, x: p[0]*(np.sin(p[1]*x))+p[2]*(np.cos(p[3]*x))+p[4] # 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

plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.title('Torsion vs Energie')
plt.xlabel('torsion')
plt.ylabel('energie')
plt.show()

Optimized params: [ -1.50760519e-03   9.94841113e-01  -1.00591843e-03   9.94840464e-01
  -7.90825389e+01]
Вывод по приближению: визуально аппроксимация хорошая¶

Увеличим шаг до 0.1¶

In [17]:
new_energies = []
news = np.arange(-10, 10+1)*0.1 + 1.52986
news
Out[17]:
array([ 0.52986,  0.62986,  0.72986,  0.82986,  0.92986,  1.02986,
        1.12986,  1.22986,  1.32986,  1.42986,  1.52986,  1.62986,
        1.72986,  1.82986,  1.92986,  2.02986,  2.12986,  2.22986,
        2.32986,  2.42986,  2.52986])
In [18]:
for new in news:
    new_energies.append(run_orca(inp.format(
        length = new,
        angle = 111.2,
        torsion1 = 180,
        torsion2 = 120,
        torsion3 = -120,
    )))
In [19]:
print new_energies

[-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.013023521422, -78.993820247811, -78.974866462702, -78.956483580267, -78.938878296987, -78.922175247055]

Посмотрим на расположение точек

In [20]:
x_o=news
y_o=new_energies

plt.plot(x_o, y_o, "ro",c='blue',alpha=0.5)
plt.title('New length vs Energie')
plt.xlabel('new length')
plt.ylabel('energie')
plt.show()

Попробуем аппроксимацию гиперболой:

$ f(x)=\frac{\alpha_0}{\alpha_1 + x }+\alpha_2 $

In [21]:
x_o=news
y_o=new_energies

fitfunc = lambda p, x: p[0]*(1/(x+p[1]))+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('New length vs Energie')
plt.xlabel('new length')
plt.ylabel('energie')
plt.show()

Optimized params: [  0.76588462  -0.40958357 -79.60453662]

Попробуем усложнить аппроксимирующую функцию:

$ f(x)=\frac{\alpha_0}{\exp(\alpha_1 * x +\alpha_2)}+\alpha_3 $

In [22]:
x_o=news
y_o=new_energies

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

p0 = [1,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('New length vs Energie')
plt.xlabel('new length')
plt.ylabel('energie')
plt.show()

/home/shad/miniconda2/envs/hse/lib/python2.7/site-packages/ipykernel_launcher.py:4: RuntimeWarning: divide by zero encountered in divide
  after removing the cwd from sys.path.

Optimized params: [  8.55307655   5.24264199  -2.37456239 -79.04492105]
Вывод по приближению: визуально аппроксимация хорошая¶
Сравните полученные константы с данными из статьи о GAFF:¶

1Eh=627,5 ккал/моль 0.63991924 Eh=401.5493231

Константа из статьи: 311 ккал/моль

Несмотря на то, что константы одного порядка, все равно отличие сильное. Возможно это отличие обусловлено импользованием различных методов для подсчета. Приближение, полученное в данной работе хоть и неплохое, но далеко от идеала. Возможно orca не учитывает какие-то уникальные свойства молекулы