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

Дано: оптимизированная структура этана в виде Z-matrix (bond_length,bond angle, torsion angle)

'''!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 * '''

Воспользуемся функцией, которая запускает ORCA и возвращает рассчитанную энергию

import subprocess
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]
    out = out.splitlines()
    for x in out:
        if "FINAL" in x:
            return(float(x.split()[-1])*1.04*10**(-2))
        elif "Error:" in x:
            print inp
            print x
            print "~~~"

Рассчёт зависимости энергии от длины (C-C)-связи¶

Запишем шаблон, куда будем вставлять различные значения длины:

inp='''!HF RHF 6-31G
    * int 0 1
    C 0 0 0 0 0 0 
    C 1 0 0 %.5f 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
    *
    '''

Создадим списки со значениями длины, с помощью функции получим для каждого из них энергию в кДж/моль.

c_bond = 1.52986

bonds = [round(c_bond + x*0.02,5) for x in range(-10,11)]
energies = [run_orca(inp%b) for b in bonds]

Для построения графиков импортируем необходимые библиотеки

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

x_o=np.array(bonds)
y_o=np.array(energies)

print "Bond energy:"
print y_o

#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", x_o, fitfunc(p1,x_o),"b-",alpha=0.5)
plt.ylabel('Energy, kJ/mol')
plt.xlabel('Bond length, angstrem')
plt.ylim(-0.8225,-0.822)
plt.xlim(1.3,1.8)
plt.show()

Bond energy:
[-0.82207804 -0.82215747 -0.82222499 -0.82228167 -0.82232849 -0.82236637
 -0.82239612 -0.82241849 -0.82243417 -0.82244379 -0.82244793 -0.82244709
 -0.82244178 -0.82243241 -0.82241938 -0.82240307 -0.82238381 -0.82236188
 -0.82233758 -0.82231114 -0.82228281]
Optimized params: [ 0.00665515  1.55699291 -0.82245482]

В результате был получен график, отражающий зависимость значения энергии от длины С-С связи (рисунок выше). Помимо этого он был аппроксимирован параболой (синия линия).

Рассчёт зависимости энергии от значения валентного угла HCС¶

Шаблон для подставления значения валентного угла HCC:

inp_2='''!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 %.3f 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
    *
    '''

Создадим списки со значениями углов, получим для каждого энергию в кДж/моль.

angles = list(np.arange(109.2,113.3,0.2))
ang_energies = [run_orca(inp_2%a) for a in angles]

x_1=np.array(angles)
y_1=np.array(ang_energies)

print "Angle energy:"
print y_1

#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_1, y_1))
print "Optimized params:", p1
 
#Plot it
plt.plot(x_1, y_1, "ro", x_1, fitfunc(p1,x_1),"b-",alpha=0.5)
plt.ylabel('Energy, kJ/mol')
plt.xlabel('Bond angle, deg')
plt.ylim(-0.822445,-0.8224485)
plt.xlim(109,113.5)
plt.show()

Angle energy:
[-0.82244674 -0.82244702 -0.82244725 -0.82244746 -0.82244763 -0.82244776
 -0.82244787 -0.82244793 -0.82244796 -0.82244796 -0.82244793 -0.82244786
 -0.82244775 -0.82244761 -0.82244744 -0.82244723 -0.82244699 -0.82244672
 -0.82244641 -0.82244607 -0.82244569]
Optimized params: [  4.27092986e-07   1.10890551e+02  -8.22447966e-01]

В результате получаем график с зависимостью энергии от значения валентного угла HCC (рисунок выше). Синия линия - аппроксимация пораболой.

Рассчёт зависимости энергии от значения торсионного угла CC¶

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

inp_3='''!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 %d
    H 1 2 3 1.08439 111.200 %d
    H 2 1 3 1.08439 111.200 %d
    H 2 1 6 1.08439 111.200 %d
    H 2 1 6 1.08439 111.200 %d
    *
    '''

Значение угла предлагалось менять от -180° до 180° с шагом 12.

tors_angles = [180 + x*12 for x in range(-30,1)]
print len(tors_angles)
tors_energies = [run_orca(inp_3%(120,-120,t,120,-1*120)) for t in tors_angles]

31

x_2=np.array(tors_angles)
y_2=np.array(tors_energies)

print "Angle energy:"
print y_2

#function is  f(x)=cos(x * freq + phase) * A+ b

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

guess_freq = 1
guess_A = 3*np.std(y_2)/(2**0.5)
guess_phase = 0
guess_b = np.mean(y_2)

p0=[guess_freq, guess_A, guess_phase, guess_b] # Initial guess for the parameters

p1, success = optimize.leastsq(errfunc, p0[:], args=(x_2, y_2))
print "Optimized params:", p1
 
#Plot it
plt.plot(x_2, y_2, "ro", x_2, fitfunc(p1,x_2),"b-",alpha=0.5)
plt.ylabel('Energy, kJ/mol')
plt.xlabel('Bond angle, deg')
plt.ylim(-0.8236,-0.82366)
plt.xlim(-190,190)
plt.show()

Angle energy:
[-0.82365475 -0.82365023 -0.82363836 -0.82362363 -0.82361166 -0.82360707
 -0.82361166 -0.82362363 -0.82363836 -0.82365023 -0.82365475 -0.82365023
 -0.82363836 -0.82362363 -0.82361166 -0.82360707 -0.82361166 -0.82362363
 -0.82363836 -0.82365023 -0.82365475 -0.82365023 -0.82363836 -0.82362363
 -0.82361166 -0.82360707 -0.82361166 -0.82362363 -0.82363836 -0.82365023
 -0.82365475]
Optimized params: [  9.94834578e-01   2.38373316e-05  -3.30163066e-06  -8.23630956e-01]

Полученная зависимость энергии от величины торсионного угла получилась в отличае от предыдущих случаев - синусоидальной.

Наверх