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

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

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

def orcan(inp, filename):
    with open(filename, 'w') as outfile:
        outfile.write(inp)
    p = subprocess.Popen("/srv/databases/orca/orca %s"%(filename), 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]
    lines = out.splitlines()
    result = 0
    for i in lines:
        if 'FINAL SINGLE' in i:
            k = i.split()
            result = float(k[-1])
    return result

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

lens = []
energ = []
step = 0.02
for n in range(-11, 11):
    length = 1.52986 + step*n
    lens.append(length)
    inp = '''!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(length))
    energ.append(orcan(inp, "orca_bond_%s.inp" % (str(n))))

#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

def plot(x, y, startlim, finlim, p0):   ##р1 - эррор функция, р0 - initial guess
    p1, success = optimize.leastsq(errfunc, p0[:], args=(x, y))
    print "Optimized params:", p1
    #Plot it
    plt.plot(x, y, "ro", x,fitfunc(p1,x),"r-",c='olive',alpha=1)
    plt.xlim(startlim,finlim)
    plt.show()

Оценим длину связи С-С:

plot(lens, energ, 1.2, 1.9, p0)

Optimized params: [  0.6772355    1.55759154 -79.08260198]

Теперь попробуем изменить шаг изменения длины связи в 5 раз (0.1).

lens_new = []
energ_new = []
step = 0.1
for n in range(-10, 11):
    length = 1.52986 + step*n
    lens_new.append(length)
    inp = '''!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(length))
    energ_new.append(orcan(inp, "orca_bond2_%s.inp" % (str(n))))

plt.plot(lens_new, energ_new, "ro", c='maroon',alpha=1)
plt.xlim()
plt.show()

Зависимость гиперболическая.

Значение энергии этана от величины валентного угла НСС

angles = []
ang_energies = []
step = 0.2
for n in range(21):
    angle = 109.2 + step*n
    angles.append(angle)
    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 %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(angle))
    ang_energies.append(orcan(inp, "orca_hcc_%s.inp" % (str(n))))

plot(angles, ang_energies, 107, 115, p0)

Optimized params: [  4.10773253e-05   1.10890768e+02  -7.90815353e+01]

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

tor_ang = []
tor_energ = []
step = 12
for n in range(31):
    torangle = -180 + step*n
    tor_ang.append(torangle)
    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 %s
    H 2 1 6 1.08439 111.200 120
    H 2 1 6 1.08439 111.200 -120
    *
    '''%(str(torangle))
    with open('orca_cc_%s.inp' % (str(n)), 'w') as outfile:
        outfile.write(inp)
    p = subprocess.Popen("/srv/databases/orca/orca orca_cc_%s.inp" % (str(n)), 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]
    out = out.splitlines()
    for outline in out:
        if ('FINAL SINGLE POINT' in outline):
            tor_energ.append(float(outline.strip('\n').split()[-1]))

plt.plot(tor_ang, tor_energ, "ro",c='maroon',alpha=1)
plt.xlim(-190,190)
plt.show()

Зависимость периодическая, а не параболическая, поэтому фиттить лучше функцией косинуса. Хорошо различимы 3 минимума (4, но 2 являются по сути одним), соответствующие заторможенной конформации этана.