In [68]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

функция по запуску ORCA

In [69]:
import subprocess
def run_orca(inp):
    with open('orca.inp', 'w') as outfile:
        outfile.write(str(inp))
    p = subprocess.Popen(str('/srv/databases/orca/orca orca.inp'  ), 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]
    out = str(out)
    lines = out.split("\\n")

    
    for s in lines:
        
        if 'FINAL SINGLE POINT ENERGY' in s:
            return float(s.split()[4])

варьирование длины одной из связей (C-C). Меняем значения длины связи C-C в диапазоне 1.4-1.78Å с шагом 0.02Å, рассчитываем энергию с помощью ORCA.

In [70]:
initial_bond=1.52986 #написано в input-файле
bonds=[]
energies=[]
bonds =  np.linspace(1.4, 1.78, 20)

#for i in range(-10,11):
 #   bond=initial_bond+i*0.02
 #   bonds.append(round(bond,5))
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
*
'''
for bond in bonds:
    cur_energ=run_orca(inp % bond)
#    cur_energ=run_orca(inp)
    energies.append(cur_energ)
print (bonds)

[1.4  1.42 1.44 1.46 1.48 1.5  1.52 1.54 1.56 1.58 1.6  1.62 1.64 1.66
 1.68 1.7  1.72 1.74 1.76 1.78]
In [71]:
print (energies)

[-79.185392866975, -79.189167575782, -79.19214806606, -79.194408414768, -79.196016199469, -79.197033104503, -79.197515438223, -79.197514577159, -79.197077392427, -79.196246631013, -79.195061288756, -79.193556853986, -79.191765711092, -79.189717344045, -79.187438590084, -79.184953877768, -79.182285428981, -79.179453454225, -79.176476328448, -79.173370740933]
In [72]:
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.ylabel('Energy, kJ/mol')
    plt.xlabel('Bond length, angstrem')
#     plt.xlim(1,2)
    plt.show()
In [73]:
plot(bonds, energies)

Optimized params: [  0.44993239   1.54424109 -79.1967562 ]

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

Расчет зависимости энергии от значения валентного угла HCC:

In [74]:
inp_2='''!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.5f 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
    *
    '''
bonds1=[]
energies1=[]
bonds1 =   np.linspace(109.2,113.2,21)
for bond1 in bonds1:
    cur_energ1=run_orca(inp_2 % bond1)
#    cur_energ=run_orca(inp)
    energies1.append(cur_energ1)
print (bonds1)

[109.2 109.4 109.6 109.8 110.  110.2 110.4 110.6 110.8 111.  111.2 111.4
 111.6 111.8 112.  112.2 112.4 112.6 112.8 113.  113.2]
In [75]:
plot(bonds1, energies1)

Optimized params: [ 4.10664822e-05  1.10890548e+02 -7.90815352e+01]

Расчет зависимости энергии от значения торсионного угла CC:

In [94]:
inp_3 = '''!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 %6.5f
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
In [95]:
bonds2=[]
energies2=[]
bonds2 =   np.linspace(-180,180,31)
for bond2 in bonds2:
    cur_energ2=run_orca(inp_3 % bond2)
#    cur_energ=run_orca(inp)
    energies2.append(cur_energ2)
print (bonds2)

[-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.]
In [102]:
plot(bonds2, energies2) #при использовании функции plot выдается странная аппроксимация. нужно зафиттить иначе

Optimized params: [-3.03437616e-08  6.38359545e-01 -7.91950087e+01]

не будем использовать заданную функцию,а зафиттим с помощью косинуса

In [105]:
fitfunc = lambda p, x: p[0]*(np.cos(p[1]*x)) + p[2]
errfunc = lambda p, x, y: fitfunc(p, x) - y
p0 = [1, 1, 1]

p1, success = optimize.leastsq(errfunc, p0[:], args=(bonds2, energies2))
print ("Optimized params:", p1)
plt.ylabel('Energy, kJ/mol')
plt.xlabel('Bond length, angstrem')
#Plot it
plt.plot(bonds2, energies2, "ro", bonds2, fitfunc(p1,bonds2),"r-",c='blue',alpha=0.5)
plt.xlim(-190,190)
plt.show()

Optimized params: [ 2.29205146e-03  9.94834578e-01 -7.91952842e+01]

функция имеет четыре минимума

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

In [106]:
bonds4=[]
energies4=[]
bonds4 =  np.linspace(0.4, 4.3, 40)

#for i in range(-10,11):
 #   bond=initial_bond+i*0.02
 #   bonds.append(round(bond,5))
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
*
'''
for bond4 in bonds4:
    cur_energ4=run_orca(inp % bond4)
#    cur_energ=run_orca(inp)
    energies4.append(cur_energ4)
print (bonds4)

[0.4 0.5 0.6 0.7 0.8 0.9 1.  1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.  2.1
 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.  3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9
 4.  4.1 4.2 4.3]
In [107]:
plot(bonds4, energies4)

Optimized params: [  0.9281547    2.7764879  -79.80581609]

наблюдаемую зависимость нельзя аппроксимировать гиперболой (хотя при маленьком диапазоне казалось, что можно)