In [1]:
import numpy as np
import subprocess
import re
import os
import matplotlib.pyplot as plt
from scipy import optimize
In [2]:
inp = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 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 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
In [3]:
regex = re.compile("FINAL\sSINGLE\sPOINT\sENERGY\s+(?P<energy>[\d\-\.]+)")
def process_file(in_name):
p = subprocess.Popen("/srv/databases/orca/orca {}".format(in_name),
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out=p.communicate()[0]
with open(in_name + "out.txt", "w") as out_file:
#print out
out_file.write(out)
for line in out.splitlines():
match = regex.match(line)
if match:
print match.group("energy")
return float(match.group("energy"))
break
In [4]:
init_value = 1.0
result_bond = np.zeros(100)
values_bond = np.arange(100) * 0.02 + init_value
for ind, val in enumerate(values_bond):
in_name = "mol_" + str(val) + "bond.zmat"
with open(in_name , "w") as in_file:
in_file.write(inp.format(val))
result_bond[ind] = process_file(in_name)
In [5]:
plt.plot(values_bond, result_bond)
plt.show()
Проделайте аналогичные операции для валентного угла HCH, его значения должны изменяться от 109.2 до 113.2. Сохраните полученные коэффициенты в отчёт. Примечание: не перезаписывайте файл со значениями энергий они Вам нужны для отчёта. Сравните полученные константы с данными из статьи о GAFF. Укажите возможные причины расхождений ваших результатов и публикацией.
In [6]:
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 {0:.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
*
'''
In [7]:
init_value = 109.2
step = 0.2
num_of_steps = (113.2 - init_value) * 5 + 1
values_angle = np.arange(num_of_steps) * step + init_value
result_angle = np.empty_like(values_angle)
for ind, val in enumerate(values_angle):
in_name = "mol_" + str(val) + "angle.zmat"
with open(in_name , "w") as in_file:
in_file.write(inp.format(val))
result_angle[ind] = process_file(in_name)
In [8]:
plt.plot(values_angle, result_angle)
plt.show()
"Проделайте аналогичные операции для торсионного угла CC, его значения должны изменяться от -180 до 180 c шагом 12. Сохраните изображение и укажите в отчёте количество минимумов функции."
In [9]:
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 {}
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
In [10]:
init_value = -180
step = 12
num_of_steps = (180 - init_value) / step + 1
values_tors = np.arange(num_of_steps) * step + init_value
result_tors = np.empty_like(values_tors, dtype=np.float64)
for ind, val in enumerate(values_tors):
in_name = "mol_" + str(val) + "angle.zmat"
with open(in_name , "w") as in_file:
in_file.write(inp.format(val))
result_tors[ind] = process_file(in_name)
In [11]:
import matplotlib.ticker as ticker
import matplotlib.pyplot as plt
ax = plt.axes()
ax.xaxis.set_major_locator(ticker.FixedLocator(np.arange(len(values_tors) )))
ax.xaxis.set_major_formatter(ticker.FixedFormatter((values_tors)))
#ax.yaxis.set_major_locator(ticker.FixedLocator(np.arange(len(result_tors) )))
#ax.yaxis.set_major_formatter(ticker.FixedFormatter((result_tors)))
plt.xticks(rotation=70)
plt.title("")
plt.yticks(result_tors)
plt.plot(result_tors)
plt.show()
