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

In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize 
import subprocess
In [8]:
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
*
'''
In [9]:
def run_orca(inp, suffix):
    with open('orca_%s.inp' % suffix, 'w') as outfile:
        outfile.write(inp)
    p = subprocess.Popen(str('/srv/databases/orca/orca orca_%s.inp' % suffix), 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]
    lines = out.splitlines()
    
    # extract energy: FINAL SINGLE POINT ENERGY'
    # and return it as float
    
    for s in out.splitlines():
        if s.startswith("FINAL SINGLE POINT ENERGY"):
            return(float(s[len("FINAL SINGLE POINT ENERGY"):]))

Варьирование длины одной из связей (C-C)

In [10]:
cc_length = []
cc_energy = []
for n in range(20):
    length = (1.52986 - 0.02 * 10) + 0.02 * n
    cc_length.append(length)
    energy = run_orca(inp % length, '_cc_' + str(length))
    cc_energy.append(energy)
In [74]:
def plot(x_o, 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), "r-", c='blue', alpha=0.5)
#     plt.xlim(1,2)
    plt.show()
In [12]:
plot(cc_length, cc_energy)

Optimized params: [  0.66824617   1.55389274 -79.08233322]

Варьирование валентного угла HCС

In [13]:
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
*
'''
In [14]:
hcc_length = []
hcc_energy = []
for n in range(41):
    length = 109.2 + 0.1 * n
    hcc_length.append(length)
    energy = run_orca(inp % length, '_hcc_' + str(length))
    hcc_energy.append(energy)
In [15]:
plot(hcc_length, hcc_energy)

Optimized params: [  4.10772398e-05   1.10890640e+02  -7.90815354e+01]

Варьирование торсионного угла CC

In [66]:
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 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
In [67]:
tcc_length = []
tcc_energy = []
for n in range(-180, 181, 12):
    length = n
    tcc_length.append(length)
    energy = run_orca(inp % length, '_tcc_' + str(length))
    tcc_energy.append(energy)
In [68]:
plt.plot(tcc_length, tcc_energy, "ro", c='blue', alpha=0.5)
print(tcc_energy)

[-79.081531423212, -79.034437767115, -78.964316116676, -78.771764535876, -77.829036111883, None, -77.834350139803, -78.772948682768, -78.947574373287, -78.985918028233, -78.99271833053, -78.98591802822, -78.947574376835, -78.772948680674, -77.834350139692, None, -77.829036111985, -78.771764535885, -78.964316116673, -79.034437767115, -79.081531423234, -79.120871238173, -79.15240032985, -79.175042279978, -79.188523815737, -79.192987718045, -79.188523815638, -79.175042280007, -79.15240032985, -79.120871238218, -79.081531423212]

Функция имеет три минимума. Для некоторых значений угла энергия не посчиталась.

Дополнительное задание

Увеличим шаг до 0.1 ангстрема при расчёте связи

In [79]:
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
*
'''
In [80]:
ecc_length = []
ecc_energy = []
for n in range(20):
    length = (1.52986 - 0.1 * 10) + 0.1 * n
    ecc_length.append(length)
    energy = run_orca(inp % length, '_ext_cc_' + str(length))
    ecc_energy.append(energy)
In [81]:
plot(ecc_length, ecc_energy)

Optimized params: [  2.95586958   1.759108   -79.56411653]

Наблюдаемую зависимость можно аппроксимировать гиперболой.

Файлы с результатами