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

Суть этой работы состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов. Нам дана оптимизированная структура этана в виде z-matrix (bond_length,bond angle, torsion 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 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-C в диапазоне 1.33-1.71Å с шагом 0.2Å, рассчитываем энергию с помощью ORCA.

In [3]:
import subprocess, os
def run_orca(inp):
    with open('orca.inp', 'w') as outfile:
        outfile.write(inp)
    my_env = os.environ.copy()
    p = subprocess.Popen("/srv/databases/orca/orca orca.inp", 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=my_env)
    p.wait()
    out=p.communicate()[0]
    lines=out.splitlines()
    for line in lines:
        if "FINAL SINGLE POINT ENERGY" in line:
            line=line.split()[4]
            return float(line)
In [4]:
start=1.52986
bondls=[]
energies=[]
for i in range(-10,10):
    bondl=start+i*0.02
    bondls.append(bondl)
    inp = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 %6.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 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
    cur_energ=run_orca(inp %bondl)
    energies.append(cur_energ)
print bondls
print energies
[1.32986, 1.34986, 1.36986, 1.38986, 1.4098600000000001, 1.42986, 1.44986, 1.46986, 1.48986, 1.50986, 1.52986, 1.54986, 1.56986, 1.58986, 1.60986, 1.62986, 1.6498599999999999, 1.66986, 1.68986, 1.70986]
[-79.16463106888, -79.17188407523, -79.178023208669, -79.183150205081, -79.187357994011, -79.190731446504, -79.19334807093, -79.195278640475, -79.196587755386, -79.197334397576, -79.197572385189, -79.197350823782, -79.196714506477, -79.195704278345, -79.194357395094, -79.192707753728, -79.190786285168, -79.188621137567, -79.186237935046, -79.183659992885]

Получили различныe значения энергии молекулы в зависимости от длины С-С связи. Строим график этой завиcимости в matplotlib.

In [5]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [34]:
x_o=np.array(bondls)
y_o=np.array(energies)
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, "bo", x_o,fitfunc(p1,x_o),"g-",alpha=0.5)
plt.xlim(1.3,1.73)
plt.ylabel('Energy, eV')
plt.xlabel('C-C bond length, angstrem')
plt.show()
Optimized params: [  0.64460658   1.54848048 -79.19820325]

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

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

In [7]:
start=111.200
angles=[]
ang_energies=[]
for n in range(-10,10):
    angle=start+n*0.2
    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 %6.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 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
    ang_energy=run_orca(inp %angle)
    ang_energies.append(ang_energy) 
print angles
print ang_energies
[109.2, 109.4, 109.60000000000001, 109.8, 110.0, 110.2, 110.4, 110.60000000000001, 110.8, 111.0, 111.2, 111.4, 111.60000000000001, 111.8, 112.0, 112.2, 112.4, 112.60000000000001, 112.8, 113.0]
[-79.197410484523, -79.19744131329, -79.197468918281, -79.197493268922, -79.197514349558, -79.197532160197, -79.197546703804, -79.197557983845, -79.197566003856, -79.197570767364, -79.197572277885, -79.197570538888, -79.197565553819, -79.197557326093, -79.197545859096, -79.197531156182, -79.197513220676, -79.197492055893, -79.197467665094, -79.197440051511]
In [33]:
x_o=np.array(angles)
y_o=np.array(ang_energies)
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, "bo", x_o,fitfunc(p1,x_o),"g-",alpha=0.5)
plt.xlim(109,113.2)
plt.ylabel('Energy, eV')
plt.xlabel('HCC angle, angstrem')
plt.show()
Optimized params: [  4.06310384e-05   1.11194916e+02  -7.91975723e+01]

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

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

In [9]:
start=180
tors_angles=[]
tors_energies=[]
for n in range(-30,1):
    tors_angle=start+n*12
    tors_angles.append(tors_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 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
*
'''    
    tors_energy=run_orca(inp %tors_angle)
    tors_energies.append(tors_energy) 
print tors_angles
print tors_energies
[-180, -168, -156, -144, -132, -120, -108, -96, -84, -72, -60, -48, -36, -24, -12, 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180]
[-79.197572396846, -79.197137699075, -79.195996511441, -79.194579710958, -79.193428553474, -79.192987702649, -79.193428552912, -79.194579709509, -79.195996509904, -79.197137698941, -79.197572392766, -79.197137699058, -79.195996511441, -79.194579710958, -79.193428553474, -79.192987702649, -79.193428552912, -79.194579709509, -79.195996509904, -79.197137698941, -79.197572392766, -79.197137699058, -79.195996511441, -79.194579710958, -79.193428553474, -79.192987702649, -79.193428552912, -79.194579709509, -79.195996509904, -79.197137698941, -79.197572392766]
In [32]:
x_o=np.array(tors_angles)
y_o=np.array(tors_energies)

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
 
#Plot it
plt.plot(x_o, y_o, "bo", x_o, fitfunc(p1,x_o),"g-",alpha=0.5)
plt.ylabel('Energy, kJ/mol')
plt.xlabel('Bond angle, deg')
plt.xlim(-190,190)
plt.show()
Optimized params: [  9.94834578e-01   2.29205075e-03  -2.08588956e-06  -7.91952842e+01]

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