Вычисление параметров для молекулярной механики
In [2]:
import subprocess
Зависимость энергии молекулы этана от длины С-С связи
In [5]:
def run_orca(inp, fileout):
with open(fileout, 'w') as outfile:
outfile.write(inp)
p = subprocess.Popen("/srv/databases/orca/orca " + fileout,
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out=p.communicate()[0]
outlines = out.splitlines()
for line in outlines:
if "FINAL SINGLE" in line:
energy = float(line.rstrip().split()[-1])
return energy
In [36]:
lengths = []
energies = []
step = 0.02
for n in range(-10, 11):
length = 1.52986 + step*n
lengths.append(length)
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))
energies.append(run_orca(inp, "orca_bond_%s.inp" % (str(n))))
In [10]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [17]:
#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 [18]:
def plot(x, y, startlim, finlim, p0):
p1, success = optimize.leastsq(errfunc, p0[:], args=(x, y))
print "Optimized params:", p1
#Plot it
plt.plot(x, y, "ro", x,fitfunc(p1,x),"r-",c='purple',alpha=0.5)
plt.xlim(startlim,finlim)
plt.show()
In [20]:
initguess = [1, 1, -79]
plot(lengths, energies, 1.2, 1.8, initguess)
Зависимость энергии молекулы этана от величины валентного угла НСС
In [22]:
valangles = []
valangle_energies = []
step = 0.2
for n in range(21):
valangle = 109.2 + step*n
valangles.append(valangle)
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(valangle))
valangle_energies.append(run_orca(inp, "orca_hcc_%s.inp" % (str(n))))
In [23]:
plot(valangles, valangle_energies, 108, 114, initguess)
Зависимость энергии молекулы этана от величины торсионного угла СС
In [14]:
torangles = []
torangle_energies = []
step = 12
for n in range(31):
torangle = -180 + step*n
torangles.append(torangle)
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(torangle))
with open('orca_cc_%s.inp' % (str(n)), 'w') as outfile:
outfile.write(inp)
p = subprocess.Popen("/srv/databases/orca/orca orca_cc_%s.inp" % (str(n)),
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out=p.communicate()[0]
out = out.splitlines()
for outline in out:
if ('FINAL SINGLE' in outline):
torangle_energies.append(float(outline.strip('\n').split()[-1]))
In [16]:
plt.plot(torangles, torangle_energies, "ro",c='purple',alpha=0.5)
plt.xlim(-190,190)
plt.show()
Функция периодическая с тремя минимумами, соответствующими трем возможным заторможенным конформациям этана.
Зависимость энергии молекулы этана от длины С-С связи (* увеличим шаг до 0.1)
In [35]:
new_lengths = []
new_energies = []
step = 0.1
for n in range(-10, 11):
new_length = 1.52986 + step*n
new_lengths.append(new_length)
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(new_length))
new_energies.append(run_orca(inp, "orca_bond_newstep_%s.inp" % (str(10 + n))))
In [34]:
plt.plot(new_lengths, new_energies, "ro",c='purple',alpha=0.5)
plt.xlim()
plt.show()
