Вычисление параметров для молекулярной механики¶
Рассчитаем 20 разных длинн связи с шагом 0.02 ангстрема для расчёта энергии в ORCA с разными значениями по длине одной из связей C-C, валентного угла H-C-С (его значения должны изменяться от 109.2 до 113.2), орсионного угла H-C-C-H (его значения должны изменяться от -180 до 180 c шагом 12). На месте # в z-matrix расположены те значения длин и углов связей, которые будут меняться.
In [1]:
import numpy
import scipy.special
import scipy.misc
from IPython.display import display,Image
import os, sys
In [2]:
import subprocess
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
Оптимизированная структура этана в виде z-matrix:
!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 180
H 2 1 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*Связь C-C¶
Расчет 20 разных длинн связи с шагом 0.02 ангстрема
'''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 # 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
*
'''Преобразуем z-matrix для удобства использования при модификации связей и углов:
In [3]:
temp = '''!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 %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 %s
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
In [4]:
values = ['1.52986', '109.2', '-180'] # bond, angle, dihedral angle
In [5]:
def run_orca(inp, filename):
with open(filename, 'w+') as wf:
wf.write(inp)
p = subprocess.Popen("/home/domain/data/prog/orca_4_0_1_2_linux_x86-64_openmpi202/orca " + filename,
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 [21]:
lengths = []
energies = []
step = 0.02
name = "C-C"
for n in range(-10, 11):
length = float(values[0]) + step * n
lengths.append(length)
inp = temp % (str(length), values[1], values[2])
energies.append(run_orca(inp, name + "_%s.inp" % (str(n))))
Посмотрим содержание одного из полученных файлов:
In [7]:
! cat C-C_10.inp
Построим график:
In [24]:
def plot(x0, y0, morse = "False"):
if morse == "True":
# f(x) = D*(1-e^(-a(x - x0)))^2
fitfunc = lambda p, x: p[0]*((1 - np.exp(-p[1]*(x-p[2])))**2) # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [-79, 1, -1.5] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x0,y0))
print "Optimized params:", p1
#Plot it
plt.plot(x0,y0, "ro", x0 ,fitfunc(p1,x0),"r-",c='blue',alpha=0.5)
plt.show()
elif morse == "False":
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
p0 = [1,1, -79] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x0,y0))
print "Optimized params:", p1
#Plot it
plt.plot(x0,y0, "ro", x0 ,fitfunc(p1,x0),"r-",c='blue',alpha=0.5)
plt.show()
In [25]:
plot(lengths[:15], energies[:15], morse = "True")
In [26]:
plot(lengths, energies)
Сделаем шаг изменения длины связи в 5 раз больше:
In [27]:
lengths = []
energies = []
step = 0.1
name = "0.1_C-C"
for n in range(-10, 11):
length = float(values[0]) + step * n
lengths.append(length)
inp = temp % (str(length), values[1], values[2])
energies.append(run_orca(inp, name + "_%s.inp" % (str(n))))
In [28]:
plot(lengths[:15], energies[:15], morse = "True")
In [44]:
plot(lengths, energies)
При шаге в 0.1 потенциал распределение точек лучше воспроизводится потенциалом Морзе, чем обычным параболическим потенциалом.
Валентного угла H-C-С¶
Изменение угла от 109.2 до 113.2
'''!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
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 [14]:
angles = []
energies = []
step = 0.2
name = "H-C-C"
for n in range(21):
angle = float(values[1]) + step * n
angles.append(angle)
inp = temp % (values[0], str(angle), values[2])
energies.append(run_orca(inp, name + "_%s.inp" % (str(n))))
In [15]:
! cat C-H-C_20.inp
In [16]:
plot(angles, energies)
Торсионный угол H-C-C-H¶
Изменение угла от -180 до 180 c шагом 12
'''!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 #
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 [17]:
dihedrals = []
energies = []
step = 12
name = "H-C-C-H"
for n in range(31):
dihedral = float(values[2]) + step * n
dihedrals.append(dihedral)
inp = temp % (values[0], values[1], str(dihedral))
energies.append(run_orca(inp, name + "_%s.inp" % (str(n))))
In [18]:
def plot_cosine(x0,y0):
x0 = np.array(x0)
y0 = np.array(y0)
# f(x) = D*(1-e^(-a(x - x0)))^2
fitfunc = lambda p, x: p[0]*(np.cos(p[1]*x)) + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1, 1, 1] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x0,y0))
print "Optimized params:", p1
#Plot it
plt.plot(x0,y0, "ro", x0 ,fitfunc(p1,x0),"r-",c='blue',alpha=0.5)
plt.xlim(-180,180)
plt.show()
In [19]:
plot_cosine(dihedrals, energies)
Ссылки¶
Визуализация: Matplotlib
Моделирование: z-matrix, ORCA
http://iopenshell.usc.edu/howto/zmatrix/ http://kodomo.fbb.msu.ru/~golovin/pdf/OrcaManual_2_9.pdf
https://sites.google.com/site/orcainputlibrary/home