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

import subprocess
In [2]:
def run_orca(inp, fileout):
    with open(fileout, 'w') as outfile:
        outfile.write(inp)
    p = subprocess.Popen("/srv/databases/orca/orca " + fileout, 
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]
    outlines = out.splitlines()
    for line in outlines:
        if "FINAL SINGLE" in line:
            k = line.split()
            energy = float(k[-1])
            return energy
In [3]:
def give_plot(x_o, y_o, start, end, p0):
    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(start,end)
    plt.show()
In [4]:
def give_plot_(x_o, y_o, start, end):

    #Plot it
    plt.plot(x_o, y_o, "ro", c='blue',alpha=0.5)
    plt.xlim(start,end)
    plt.show()
In [5]:
#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

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

In [6]:
lengths = []
energies = []
for n in range(-10, 11):
    length = 1.52986 + 0.02*n
    lengths.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))
    energies.append(run_orca(inp, 'CC_bond_outfile.inp'))
    
print(energies)

[-79.045965590542, -79.053603047463, -79.060094932942, -79.06554485476, -79.070047523843, -79.073689505281, -79.076549916355, -79.078701061159, -79.080208995167, -79.081134076598, -79.081531431378, -79.081451404936, -79.080939963093, -79.080039062387, -79.078786990807, -79.077218675992, -79.075365973146, -79.073257919932, -79.070920984592, -79.068379281351, -79.065654773152]
In [7]:
p0 = [1,1, -79]
give_plot(lengths, energies, 1, 2, p0)

Optimized params: [  0.63991792   1.55699294 -79.08219399]

Зависимость энергии молекулы этана от размера угла H-C-C:

In [8]:
angles = []
energies = []
for n in range(21):
    angle = 109.2 + 0.2*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))
    energies.append(run_orca(inp, 'HCC_angle_outfile.inp'))
    
print(energies)    

[-79.081417669214, -79.081443881466, -79.081466814993, -79.081486450148, -79.081502772458, -79.0815157829, -79.081525485378, -79.081531884287, -79.081534984087, -79.081534789216, -79.081531304089, -79.081524533069, -79.081514480488, -79.081501150634, -79.081484547759, -79.081464676074, -79.081441539752, -79.081415142944, -79.081385489745, -79.081352584204, -79.081316430358]
In [9]:
give_plot(angles, energies, 109, 114, p0)

Optimized params: [  4.10675329e-05   1.10890567e+02  -7.90815352e+01]

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

In [10]:
angles = []
energies = []
for n in range(31):
    angle = - 180 + 12*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 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(angle))
    energies.append(run_orca(inp, 'CC_angle_outfile.inp'))
    
print(energies)        

[-79.197572400974, -79.197137699071, -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.195996509903, -79.197137698941, -79.197572392766]
In [11]:
give_plot_(angles, energies, -185, 185)

У этой функции 3 минимума, соответствующих заторможенным конформациям.

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

In [12]:
lengths = []
energies = []
for n in range(-10, 11):
    length = 1.52986 + 0.1*n
    lengths.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))
    energies.append(run_orca(inp, 'CC_bond_outfile2.inp'))
    
print(energies)

[-73.376057637917, -75.600038672809, -76.986925190759, -77.840020402069, -78.364137030277, -78.683902686411, -78.875495245128, -78.986233920144, -79.045965568718, -79.07368949118, -79.081531419723, -79.077218664924, -79.065654757297, -79.049930677904, -79.031982697254, -79.013023521422, -78.993820247811, -78.974866462702, -78.956483580267, -78.938878296987, -78.922175247055]
In [13]:
give_plot_(lengths, energies, 0.25, 2.5)

Эту функцию можно было бы аппроксимировать гиперболической функцией.

CC_angle_outfile.inp
CC_bond_outfile.inp
CC_bond_outfile2.inp
HCC_angle_outfile.inp

In [ ]: