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

В этом задании нужно было научиться определять константы ковалентных взаимодействий на примере этана. Молекула этана задана в виде Z-matrix.

Для начала мы будем менять значения длины связи C-C в диапазоне 1.33-1.71Å с шагом 0.2Å. Для этого создадим 20 файлов с соответствующими Z-matrix.

In [ ]:
text = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 %f 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 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
for i in range(-10,10):
    name = 'A'+str(i)+'_CC.inp'
    l = 1.53 + i*0.02
    out_file = open (name, 'w')
    out_file.write(text %(l))
    out_file.close()

Напишем функцию для запуска ORCA и получения значения энергии от входной Z-matrix.

In [ ]:
import subprocess
def run_orca(inp):
    p = subprocess.Popen("/srv/databases/orca/orca %s" %(inp), 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]
    out.splitlines()
    
    #print out
    
    for line in out.splitlines():
        if 'FINAL SINGLE POINT ENERGY ' in line:
            en = line.split(' ')[-1]
            return en

    # extract energy: FINAL SINGLE POINT ENERGY'
    # and return it 
    
    #return ....
In [ ]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [ ]:
x_o=np.arange(1.33, 1.73, 0.02)
y=[]
for i in range(-10,10):
    name = 'A'+str(i)+'_CC.inp'
    en = run_orca(name)
    y.append(float(en))
y_o = np.array(y)
In [ ]:
#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),"r-",c='blue',alpha=0.5)
plt.xlim(1,2)
#plt.show()
plt.savefig("pr6_gr1.png")


Рис 1. График зависимости энергии соединения от длины связи C-C.

Мы получили график зависимости энергии соединения от длины связи. Минимум на графике примерно соответствует значению 1.53Å, чего и следовало ожидать.

Далее фиксируем все, кроме значения одного валентного угла H-C-C. Будем менять его в диапазоне 109.2-113.2° с шагом 0.2°. Для этого создадим 21 файл с Z-matrix.

In [ ]:
text = '''!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 %f 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 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
for i in range(-10,11):
    name = 'B'+str(i)+'_CC.inp'
    l = 111.2 + i*0.2
    out_file = open (name, 'w')
    out_file.write(text %(l))
    out_file.close()
In [ ]:
x_o=np.arange(109.2, 113.3, 0.2)
y=[]
for i in range(-10,11):
    name = 'B'+str(i)+'_CC.inp'
    en = run_orca(name)
    y.append(float(en))
y_o = np.array(y)
In [ ]:
#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),"r-",c='blue',alpha=0.5)
plt.xlim(108,114)
#plt.show()
plt.savefig("pr6_gr2")


Рис 2. График зависимости энергии соединения от величины валентного угла H-C-C.

Получился график зависимости энергии соединения от величины валентного угла. Опять же минимум энергии соответствует наиболее оптимальной структуре. В данном случае это 111.2°, что ожидаемо.

Торсионный угл мы будем менять от -180° до 180° с шагом 12°. Дальше все аналогично.

In [ ]:
text = '''!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 %i
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
for i in range(-180,181, 12):
    #l = 111.2 + i*0.02
    name = 'C'+str(i)+'_CC.inp'
    out_file = open (name, 'w')
    out_file.write(text %(i))
    out_file.close()
In [ ]:
x_o=np.arange(-180,181, 12)
y=[]
for i in range(-180,181, 12):
    name = 'C'+str(i)+'_CC.inp'
    en = run_orca(name)
    y.append(float(en))
y_o = np.array(y)
In [ ]:
#Plot it
plt.plot(x_o, y_o, "ro", c='blue',alpha=0.5)
#plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.xlim(-180,180)
#plt.show()
plt.savefig("pr6_gr3")


Рис 3. График зависимости энергии соединения от величины торсионного угла H-C-C-H.

Мы видим более одного минимума: -180, -60, 60 и 180, что, по-видимому, соответствует торсионным углам H-C-C-H для различных пар водородов в молекуле этана.