Вычисление параметров для молекулярной механики¶
In [5]:
import numpy as np
import subprocess
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import optimize
from IPython.display import Image, display
Функция по запуску ORCA с kodomo¶
In [7]:
def run_orca(inp):
with open('orca.inp', 'w') as outfile:
outfile.write(inp)
with subprocess.Popen("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
Варьируем длину связи¶
In [20]:
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
*
'''
In [21]:
bond_lengths = np.arange(1.3, 1.7, 0.02)
energies = []
for length in bond_lengths:
energy = run_orca(bond_length_template.format(length))
energies.append(energy)
energies
Out[21]:
In [23]:
plt.scatter(bond_lengths, energies)
Out[23]:
Варьируем валентный угол HCС¶
In [ ]:
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
*'''
Шаг не был указан и был взят по аналогии с предыдущим, поэтому получилось достаточно много значений.
In [26]:
angle_HCC = np.arange(109.2, 113.2, 0.02)
energies = []
for angle in angle_HCC:
energy = run_orca(angle_template.format(angle))
energies.append(energy)
energies
Out[26]:
In [27]:
plt.scatter(angle_HCC, energies)
Out[27]:
Варьируем торсионный угол CC¶
In [25]:
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
*'''
In [28]:
torsion_CC = np.arange(-180, 180, 12)
energies = []
for angle in torsion_CC:
energy = run_orca(torsion_template.format(angle))
energies.append(energy)
energies
Out[28]:
In [29]:
plt.scatter(torsion_CC, energies)
Out[29]:
Минимумы энергии соответствуют конформации этана с максимальным удалением атомов водорода друг от друга. Таких различных случаев возможно 3.