6. Changing molecular mechanics parameters.

In [31]:
import os
import subprocess 
import numpy as np
import pandas as pd
from scipy import optimize
import matplotlib.pyplot as plt

pd.set_option('display.max_rows', 10)
In [93]:
def run_orca(name, path, refval, startval, step, iter_number):
    result={}
    init_val = refval - startval
    for i in range(iter_number):
        with open(os.path.join(os.getcwd(), f"{path}", '{}_{}.inp'.format(i, path)), 'w') as outfile:
            value = init_val + i*step
            newinp = name.format(value)
            outfile.write(newinp)
        p = subprocess.Popen(". ~/.bash_profile; orca {0}/{1}_{0}.inp".format(path, i), shell=True,
                            stdout=subprocess.PIPE, stderr=subprocess.PIPE)
        out = p.communicate()[0]
        out = out.splitlines()
        for outline in out:
            if "FINAL SINGLE" in str(outline):
                dict_value = float(outline.split()[-1])
                result[value] = dict_value
    return result
In [95]:
def fit_func_plot(x_0, y_0, criterium):
    # Target function for dihedral angle
    if criterium == "dih":
        fitfunc = lambda p, x: p[0]*(np.cos(p[1]*x)) + p[2]
        p0 = [1,1, 1]
        errfunc = lambda p, x, y: fitfunc(p, x) - y 
        p1, success = optimize.leastsq(errfunc, p0[:], args=(x_0, y_0))
        print ("Optimized params:", p1)
    # Target function for cc_bond with big step
    elif criterium == "bond_ds":
        fitfunc = lambda p, x: np.multiply(p[0],pow(x,p[1])) +  np.multiply(p[2], pow(x,p[3])) - p[4]
        p0 = [1, -12, -10, -2, 79]
        errfunc = lambda p, x, y: fitfunc(p, x) - y  
        p1, success = optimize.leastsq(errfunc, p0[:], args=(x_0, y_0))
        print("Optimized params:", p1)
    # Target function for hcc angle and cc_bond
    else:
        fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] 
        p0 = [1,1, -79]
        errfunc = lambda p, x, y: fitfunc(p, x) - y 
        p1, success = optimize.leastsq(errfunc, p0[:], args=(x_0, y_0))
        print ("Optimized params:", p1)

    plt.plot(x_0, y_0, "ro", x_0,fitfunc(p1,x_0),"r-",c='blue',alpha=0.5)
    return plt.show()

6.1 Different bond length.

In [91]:
# int 0 1 -- neutral charge (0), 1 multiplicity (2S + 1, where S -- total spin)
# S = (# of alpha electrons)(+1/2) + (# of beta electrons)(-1/2)

# Three first digits -- cartesian coordinates:
# first -- the number of the atom with which the atom forms a bond;
# second -- the number of the atom with which two atoms (they form a bond) form the angle;
# third -- the number of the atom with which three atoms (two -- common) form the dihedral angle;
# Last three values -- values of a bond, angle and dihedral angle.

orca_inp_bond = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 {} 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
*
'''
In [94]:
energies_bond = run_orca(orca_inp_bond,"bond",1.52986,0.4,0.02,41)
In [16]:
df_bond = pd.DataFrame.from_dict(list(energies_bond.items()))
df_bond.rename(columns={0:"bond_length", 1:"energy"}, inplace=True)
df_bond
Out[16]:
bond_length energy
0 1.12986 -78.875495
1 1.14986 -78.902957
2 1.16986 -78.927518
3 1.18986 -78.949432
4 1.20986 -78.968934
... ... ...
36 1.84986 -79.046472
37 1.86986 -79.042938
38 1.88986 -79.039339
39 1.90986 -79.035684
40 1.92986 -79.031983

41 rows × 2 columns

In [17]:
# Looks like a parabola.

x_bond = df_bond.bond_length
y_bond = df_bond.energy

fit_func_plot(x_bond, y_bond, "bond")

Optimized params: [  0.77567117   1.60806717 -79.08937121]

6.2 Different valent angles.

In [14]:
orca_inp_angle = '''!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 {} 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
*
'''
In [15]:
energies_angle = run_orca(orca_inp_angle,"angle",111.200,2,0.1,41)
In [18]:
df_angle = pd.DataFrame.from_dict(list(energies_angle.items()))
df_angle.rename(columns={0:"hcc_angle", 1:"energy"}, inplace=True)
df_angle
Out[18]:
hcc_angle energy
0 109.2 -79.081418
1 109.3 -79.081431
2 109.4 -79.081444
3 109.5 -79.081456
4 109.6 -79.081467
... ... ...
36 112.8 -79.081386
37 112.9 -79.081370
38 113.0 -79.081353
39 113.1 -79.081335
40 113.2 -79.081317

41 rows × 2 columns

In [19]:
# Looks like a parabola.

x_angle = df_angle.hcc_angle
y_angle = df_angle.energy

fit_func_plot(x_angle, y_angle, "angle")

Optimized params: [ 4.10769603e-05  1.10890633e+02 -7.90815353e+01]

6.3 Different dihedral angles.

In [57]:
orca_inp_dih = '''!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 {}
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
In [58]:
energies_dih = run_orca(orca_inp_dih,"dih",-180,0,12,31)
In [59]:
df_dih = pd.DataFrame.from_dict(list(energies_dih.items()))
df_dih.rename(columns={0:"dih_angle", 1:"energy"}, inplace=True)
df_dih
Out[59]:
dih_angle energy
0 -180 -79.197572
1 -168 -79.197138
2 -156 -79.195997
3 -144 -79.194580
4 -132 -79.193429
... ... ...
26 132 -79.193429
27 144 -79.194580
28 156 -79.195997
29 168 -79.197138
30 180 -79.197572

31 rows × 2 columns

In [60]:
# Looks like a cosinus.

x_dih = df_dih.dih_angle
y_dih = df_dih.energy

fit_func_plot(x_dih, y_dih, "dih")

Optimized params: [ 2.29205061e-03  9.94834579e-01 -7.91952842e+01]

There are 3 minima corresponding to the staggered conformation of ethane: one in 180 (a.k.a. -180) and two in -60, 60.

The 3 maxima correspond to the eclipsed conformation (120, -120 and 0).

6.4 Different bond length (step 0.1 A.)

In [66]:
energies_bond_ds = run_orca(orca_inp_bond,"bond_ds",1.52986,0.8,0.1,44)
In [86]:
df_bond_ds = pd.DataFrame.from_dict(list(energies_bond_ds.items()))
df_bond_ds.rename(columns={0:"bond_length", 1:"energy"}, inplace=True)
df_bond_ds
Out[86]:
bond_length energy
0 0.72986 -76.986925
1 0.82986 -77.840020
2 0.92986 -78.364137
3 1.02986 -78.683903
4 1.12986 -78.875495
... ... ...
39 4.62986 -78.739000
40 4.72986 -78.735691
41 4.82986 -78.732704
42 4.92986 -78.730007
43 5.02986 -78.727573

44 rows × 2 columns

In [96]:
# Looks like a sum of power functions.

x_bond_ds = df_bond_ds.bond_length
y_bond_ds = df_bond_ds.energy

fit_func_plot(x_bond_ds, y_bond_ds, "bond_ds")

Optimized params: [ 31.70403902  -2.16402165 -31.75938601  -2.07625653  78.56270824]