Занятие 6. Вычисление параметров для молекулярной механики.

In [31]:
%matplotlib inline
import subprocess
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [19]:
structure = '''!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 [20]:
def run_orca(inp):
    with open('orca.inp', 'w') as outfile:
        outfile.write(inp)
    p = subprocess.Popen("/srv/databases/orca/orca orca.inp", 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out = p.communicate()[0]
    for i in out.splitlines():
        if i.startswith("FINAL"):
            return float(i.split(' ')[-1])
In [21]:
# эта функция будет строить график с аппроксимацией параболчиеской функцией
def plot_with_app(x_axis, y_axis):
    # добавляем аппроксимацию полученной зависимости функцией 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_axis, y_axis))
    plt.figure(figsize=(10, 6))
    plt.plot(x_axis, y_axis, "ro", x_axis, fitfunc(p1, x_axis), "r-", c='green', alpha=0.5)
    print "Optimized params:", p1

Зависимость энергии молекулы от длины одной C-C связи.

In [22]:
# генерим входные файлы для ORCA с новой длиной C-C связи
# длина связи изменяется от исходного значения 1.52986 с шагом 0.02

mid_len = 1.52986
incr_len = 0.02
start_len = mid_len - incr_len * 10
stop_len = mid_len + incr_len * 10
en_len = []  # energies
j = start_len

for i in range(21):
    struct_len = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 %s 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
*
''' % str(j)
    en_len.append(run_orca(struct_len))
    j += incr_len
In [42]:
x_len = np.arange(start_len, stop_len + incr_len, incr_len)
y_len = np.array(en_len)
In [56]:
plot_with_app(x_len, y_len)

Optimized params: [  0.63991931   1.55699288 -79.08219401]

Зависимость энергии молекулы от валентного угла HCC.

In [76]:
# генерим входные файлы для ORCA с новым значением валентного угла HCC
# угол изменяется от исходного значения 111.2 с шагом 0.2

mid_ang = 111.2
incr_ang = 0.2
start_ang = mid_ang - incr_ang * 10
stop_ang = mid_ang + incr_ang * 10
en_ang = []
k = start_ang

for i in range(21):
    struct_ang = '''!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 %s 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
*
''' % str(k)
    en_ang.append(run_orca(struct_ang))
    k += incr_ang
In [62]:
x_ang = np.arange(start_ang, stop_ang + 0.02, incr_ang)
y_ang = np.array(en_ang)
In [64]:
plot_with_app(x_ang, y_ang)

Optimized params: [  4.10665715e-05   1.10890550e+02  -7.90815352e+01]

Зависимость энергии молекулы от торсионного угла CC.¶

In [6]:
# генерим входные файлы для ORCA с новым значением торсионного угла CC
# угол изменяется от исходного значения -180 до 180 с шагом 12

mid_tang = 0
incr_tang = 12
start_tang = mid_tang - incr_tang * 15
stop_tang = mid_tang + incr_tang * 15
en_tang = []
l = start_tang

for i in range(31):
    struct_tang = '''!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 %s
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
''' % (str(l))
    en_tang.append(run_orca(struct_tang))
    l += incr_tang
In [10]:
x_tang = np.arange(start_tang, stop_tang + 0.02, incr_tang)
y_tang = np.array(en_tang)
In [16]:
plt.figure(figsize=(10, 6))
plt.scatter(x_tang, y_tang, c='green')
Out[16]:
<matplotlib.collections.PathCollection at 0x7f2596dfcc50>

На графике формально видно 4 минимума функции, но так как функция явно периодическая, минимумы на 180 и на -180 - это на самом деле один и тот же минимум, поэтому всего минимумов 3. Эти минимумы соответствуют заторможенным конформациям этана (-180, -60, 60). Максимумы функции (-120, 0, 120) соответствуют заслоненным конформациям.

Зависимость энергии молекулы от длины одной C-C связи (с шагом 0.1).

In [32]:
# генерим входные файлы для ORCA с новой длиной C-C связи
# длина связи изменяется от исходного значения 1.52986 с шагом 0.1

incr_llen = 0.1
start_llen = mid_len - incr_llen * 10
stop_llen = mid_len + incr_llen * 10
en_llen = []  # energies
m = start_llen

for i in range(21):
    struct_llen = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 %s 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
*
''' % str(m)
    en_llen.append(run_orca(struct_llen))
    m += incr_llen
In [27]:
x_llen = np.arange(start_llen, stop_llen + 0.02, incr_llen)
y_llen = np.array(en_llen)
In [28]:
plt.figure(figsize=(10, 6))
plt.scatter(x_llen, y_llen, c='green')
Out[28]:
<matplotlib.collections.PathCollection at 0x7f258fc9fc10>

Зависимость похожа на гиперболическую, но в конце не стремится к нулю, а отклоняется от Ox.