
import subprocess
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize


def run_orca(inp):
    with open('orca.inp', 'w') as outfile:
        outfile.write(inp)
    with subprocess.Popen("/srv/databases/orca/orca orca.inp", shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) as p:
        out=p.communicate()[0].decode()
    for line in out.splitlines():
        if 'FINAL SINGLE POINT' in line:
            return float(line.strip().split()[4])
    return None

bond_length_template = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 {:.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 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
angle_template = '''!MP2 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 {:.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
*'''
torsion_template = '''!MP2 6-311G**
* 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 {:d}
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*'''


bond_lengths = np.arange(1.3, 1.7, 0.02)
with open('bond_length_energies.txt', 'w') as outfile:
        outfile.write(str([run_orca(bond_length_template.format(value)) for value in bond_lengths]))
        
long_bond_lengths = np.arange(1, 2.6, 0.1)
with open('long_bond_length_energies.txt', 'w') as outfile:
        outfile.write(str([run_orca(bond_length_template.format(value)) for value in long_bond_lengths]))
        
angle_values = np.arange(109.2, 113.21, 0.2)
with open('angle_energies.txt', 'w') as outfile:
        outfile.write(str([run_orca(angle_template.format(value)) for value in angle_values]))
        
torsion_values = np.arange(-180, 181, 12)
with open('torsion_energies.txt', 'w') as outfile:
        outfile.write(str([run_orca(torsion_template.format(value)) for value in torsion_values]))


bond_lengths = np.arange(1.3, 1.7, 0.02)
with open('bond_length_energies.txt', 'r') as infile:
    bond_length_energies = np.array([float(i) for i in infile.readline()[1:-1].split(', ')])


def length_energy(length, k, l0, a):
    return k / 2 * (length - l0) ** 2 + a

popt_l, pcov = optimize.curve_fit(length_energy, bond_lengths, bond_length_energies)
popt_l

array([  1.54918132,   1.54783145, -79.08299494])


plt.plot(bond_lengths, bond_length_energies, 'o', label='orca')
plt.plot(bond_lengths, length_energy(bond_lengths, *popt_l), label='fitted', color='red')
plt.legend()
plt.show()


def morze_potential(length, d, a, l0, b):
    return d * (1 - np.exp(- a * (length - l0))) ** 2 + b

popt_m, pcov = optimize.curve_fit(morze_potential, bond_lengths, bond_length_energies)
popt_m

array([  0.19079668,   1.74329552,   1.53584299, -79.0815584 ])


plt.plot(bond_lengths, bond_length_energies, 'o', label='orca')
plt.plot(bond_lengths, morze_potential(bond_lengths, *popt_m), label='fitted', color = 'red')
plt.legend()
plt.show()


from functools import partial

morze_for_optimization = partial(morze_potential, **dict(zip(['d', 'a', 'l0', 'b'], popt_m)))
optimization_results = optimize.minimize_scalar(morze_for_optimization, bounds=(1.3, 1.6))
optimization_results.x

1.5358429414347123


angle_values = np.arange(109.2, 113.21, 0.2)
with open('angle_energies.txt', 'r') as infile:
    angle_energies = np.array([float(i) for i in infile.readline()[1:-1].split(', ')])


def angle_potential(x, k, f0, a):
    return k / 2 * (x - f0) ** 2 + a

popt_a, pcov = optimize.curve_fit(angle_potential, angle_values, angle_energies)
popt_a

array([ 7.23517961e-05,  1.11037057e+02, -7.94918725e+01])


plt.plot(angle_values, angle_energies, 'o', label='orca')
plt.plot(angle_values, angle_potential(angle_values, *popt_a), label='fitted', color = 'red')
plt.legend()
plt.show()


angle_for_optimization = partial(angle_potential, **dict(zip(['k', 'f0', 'a'], popt_a)))
optimization_results = optimize.minimize_scalar(angle_for_optimization, bounds=(110, 112))
optimization_results.x

111.03705069866203


torsion_values = np.arange(-180, 181, 12)
with open('torsion_energies.txt', 'r') as infile:
    torsion_energies = np.array([float(i) for i in infile.readline()[1:-1].split(', ')])


def torsion_potential(x, v, n, g, a):
    return v / 2 * (1 + np.cos(np.radians(n * np.pi * x - g))) + a

popt_t, pcov = optimize.curve_fit(torsion_potential, torsion_values, torsion_energies)
popt_t

array([ 5.19874183e-03,  9.55055359e-01,  6.34421050e-03, -7.95705882e+01])


plt.plot(torsion_values, torsion_energies, 'o', label='orca')
plt.plot(torsion_values, torsion_potential(torsion_values, *popt_t), label='fitted', color = 'red')
plt.xticks(np.arange(-180, 181, 60))
plt.legend()
plt.show()


torsion_for_optimization = partial(torsion_potential, **dict(zip(['v', 'n', 'g', 'a'], popt_t)))
optimization_results = optimize.minimize_scalar(torsion_for_optimization, bounds=(0, 120))
optimization_results.x

59.99421061538492


long_bond_lengths = np.arange(1, 2.6, 0.1)
with open('long_bond_length_energies.txt', 'r') as infile:
    long_bond_length_energies = np.array([float(i) for i in infile.readline()[1:-1].split(', ')])


popt_long, pcov = optimize.curve_fit(morze_potential, long_bond_lengths, long_bond_length_energies)
popt_long

array([  0.23994803,   1.65275979,   1.53227029, -79.08500184])


plt.plot(long_bond_lengths, long_bond_length_energies, 'o', label='orca')
plt.plot(long_bond_lengths, morze_potential(long_bond_lengths, *popt_long), label='fitted', color = 'red')
plt.legend()
plt.show()


morze_long_for_optimization = partial(morze_potential, **dict(zip(['d', 'a', 'l0', 'b'], popt_long)))
optimization_results = optimize.minimize_scalar(morze_long_for_optimization, bounds=(1.4, 1.6))
optimization_results.x

1.5322702064517515

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

Расчёт минимума энергии для длины C-C связи¶

Расчёт минимума энергии для валентного угла H-C-C¶

Расчёт минимума энергии для торсионного угла C-C¶

Расчёт минимума энергии для длины C-C связи (шаг = 0.1 А)¶