Практикум 6: вычисление параметров для молекулярной механики¶
26.04.2018
In [3]:
import subprocess
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
original bond length: 1.52986
In [4]:
def get_fin_energy_orca(inp, filename):
with open(filename, 'w') as outfile:
outfile.write(inp)
p = subprocess.Popen("/srv/databases/orca/orca %s"%(filename),
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out=p.communicate()[0]
lines = out.splitlines()
final_energy = 0
for i in lines:
if 'FINAL SINGLE' in i:
k = i.split()
final_energy = float(k[-1])
return final_energy
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]:
def give_me_a_plot(x_range, y_range, start, end, initial_guess):
p0 = initial_guess # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_range, y_range))
print "Optimized params:", p1
#Plot it
plt.plot(x_range, y_range, "ro", x_range,fitfunc(p1,x_range),"r-",c='blue',alpha=0.5)
plt.xlim(start,end)
plt.show()
Часть 1. С-С связь¶
In [2]:
CC_bond_lens = []
x_range = []
for i in range(-10, 11):
length = 1.52986 + 0.02*i
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))
CC_bond_lens.append(inp)
x_range.append(length)
In [9]:
energies = []
for k in CC_bond_lens:
en = get_fin_energy_orca(k, 'orca_CC_bond.inp')
energies.append(en)
print(energies)
In [29]:
initial_guess_1 = [1,1, -79]
give_me_a_plot(x_range, energies, 1, 2, initial_guess_1)
Часть 2. Валентный угол HCC¶
In [72]:
HCC_val_angles = []
x_range_va = []
for i in range(21):
val_angle = 109.2 + 0.2*i
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(val_angle))
HCC_val_angles.append(inp)
x_range_va.append(val_angle)
In [73]:
energies_val_angl = []
for k in HCC_val_angles:
en = get_fin_energy_orca(k, 'orca_HCC_val_angle.inp')
energies_val_angl.append(en)
print(energies_val_angl)
In [75]:
initial_guess_1 = [1,1, -79]
give_me_a_plot(x_range_va, energies_val_angl, 108, 114, initial_guess_1)
Осталось найти в статье параметры
Файлы с результатами¶
Часть 3. Торсионный угол СС¶
In [13]:
CC_torsion = []
x_range_to = []
for i in range(31):
torsion = - 180 + 12*i
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(torsion))
CC_torsion.append(inp)
x_range_to.append(torsion)
In [14]:
energies_torsion = []
for k in CC_torsion:
en = get_fin_energy_orca(k, 'orca_CC_torsion_angle.inp')
energies_torsion.append(en)
print(energies_torsion)
In [16]:
plt.plot(x_range_to, energies_torsion, "ro", c='blue', alpha = 0.5)
plt.xlim(-180,180)
plt.show()
Получилось 3 минимума
Часть 4. Дополнительное задание¶
In [9]:
CC_bond_lens2 = []
x_range2 = []
for i in range(-10, 11):
length = 1.52986 + 0.1*i
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))
CC_bond_lens2.append(inp)
x_range2.append(length)
In [16]:
energies2 = []
for k in CC_bond_lens2:
en = get_fin_energy_orca(k, 'orca_CC_bond2.inp')
energies2.append(en)
print(energies2)
In [17]:
plt.plot(x_range2, energies2, "ro", c='blue', alpha = 0.5)
plt.xlim()
plt.show()
Это похоже на гиперболу.