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


In [1]:
import subprocess

Предоставлена оптимизированная структура этана в виде z-matrix :

In [2]:
inp = '''!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 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''

Надо определить константы ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов.

Меняя значения длины связи C-C в диапазоне 1.33-1.71 ангстрем с шагом 0.2 ангстема рассчитываем энергию с помощью ORCA. Получаем различныe значения энергии молекулы в зависимости от длины С-С связи. Строим график этой завиcимости в matplotlib.

In [35]:
inps = []
dists = []
do = 1.52986
for i in range(-10,11):
#    print i, i*0.02, do + i*0.02
    j = i*0.02
    di = do + j
    inpi = inp[0:47] + "%.5f" % di + inp[54:]
#    print inpi
    inps = inps + [inpi]
    dists = dists + [di]
In [36]:
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]
    d = inp[47:54]
    with open('orca'+d+'.out', 'w') as outfile:
        outfile.write(out)
    for i in out.splitlines():
        if i.startswith('FINAL SINGLE POINT ENERGY'):
            energy = float(i.strip().split()[-1])
            return energy
In [37]:
energies = []
for i in inps:
    en = run_orca(i)
    energies = energies + [en]
In [7]:
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [38]:
x_o=np.asarray(dists)
y_o=np.asarray(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

Optimized params: [  0.61699538   1.55137226 -79.1980466 ]
In [39]:
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.show()

Полученный график отражает зависимость значения энергии от длины С-С связи в молекуле этана. Он был аппроксимирован параболой. Минимум графика соответствует значению длины связи в оптимизированной структуре.

Проделываем аналогичные действия для валентного угла HCC. Его значения меняем от 109.2° до 113.2° с шагом 0.2°.

In [11]:
def run_orca2(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]
    d = inp[76:82]
    with open('orca'+d+'.out', 'w') as outfile:
        outfile.write(out)
    for i in out.splitlines():
        if i.startswith('FINAL SINGLE POINT ENERGY'):
            energy = float(i.strip().split()[-1])
            return energy
In [12]:
inps2 = []
vals = np.arange(109.2,113.4,0.2)
for i in vals:
    inpi = inp[0:76] + "%.3f" % i + inp[82:]
    inps2 = inps2 + [inpi]
In [14]:
energies2 = []
for i in inps2:
    en = run_orca2(i)
    energies2 = energies2 + [en]
In [16]:
x_o=vals
y_o=np.asarray(energies2)

#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

Optimized params: [  4.05910929e-05   1.11195324e+02  -7.91975723e+01]
In [17]:
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.show()

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

Проделываем аналогичные действия для торсионного угла, значения которого меняем от -180° до 180° с шагом 12°.

In [19]:
def run_orca3(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]
    d = inp[167:170]
    with open('orca'+d+'.out', 'w') as outfile:
        outfile.write(out)
    for i in out.splitlines():
        if i.startswith('FINAL SINGLE POINT ENERGY'):
            energy = float(i.strip().split()[-1])
            return energy
In [28]:
inps3 = []
tors = np.arange(-180,180,12)
#print tors
for i in tors:
#    print i, i*0.02, do + i*0.02
    inpi = inp[0:167] + "%i" % i + inp[170:]
#    print inpi
    inps3 = inps3 + [inpi]
In [29]:
energies3 = []
for i in inps3:
    en = run_orca3(i)
    energies3 = energies3 + [en]
In [33]:
x_o=tors
y_o=np.asarray(energies3)

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_o)/(2**0.5)
guess_phase = 0
guess_b = np.mean(y_o)

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

p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print "Optimized params:", p1

Optimized params: [  9.94834904e-01   2.29231331e-03  -1.70484762e-04  -7.91952844e+01]
In [34]:
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.show()

График отражает зависимость энергии от значения торсионного угла в молекуле этана. Этот график имеет 4 минимума (при -180, -60, 60 и 180 °). Эти четыре значения угла соответствуют торсионным углам H-C-C-H для различных пар атомов водорода в оптимизированной молекуле этана (заторможенная конформация).