Вычисление параметров для молекулярной механики
Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов. Работать будем со структурой молекулы этана. Будем варьировать параметры структуры и смотреть, какие структуры наиболее выгодны по энергии, посчитанной с помощью ORCA.
Имеем оптимизированную структуру этана в формате z-matrix. Создаем 20 файлов с этой z-matrix, меняя значение длины связи С-С с шагом 0.02 ангстрем (по 10 в каждую сторону от оптимальной длины связи 1.52986):
# make files
lin='''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 %f 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
*
'''
for i in range(-10,10):
j=0.02*i
bond=1.52986+j
name="bond/bond"+str(i)
fi=open(name,'w')
fi.write(lin %(bond))
fi.close()
i=i+1
Будем запускать ORCA для каждого из этих файлов и получать на выходе значение энергии структуры (FINAL SINGLE POINT ENERGY). Затем будем строить график зависимости энергии от варьируемого параметра:
import subprocess
def run_orca(infile):
p = subprocess.Popen("/srv/databases/orca/orca %s" %(infile),shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out=p.communicate()[0]
e = float
for line in out.splitlines():
if "FINAL SINGLE POINT ENERGY" in line:
e=float(line.split(" ")[-1])
# extract energy: FINAL SINGLE POINT ENERGY'
# and return it as float
if e != 0:
return e
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
def plot(infile,label,x_o,p0,outfile):
open(outfile, 'w').close() #erase previous numbers in outfile
es=[]
for i in range(-len(x_o)/2,len(x_o)/2):
e=run_orca("%s%s" %(infile,i))
es.append(e)
#write energies to file
fi=open(outfile,'a')
fi.write(str(e)+"\n")
i=i+1
fi.close()
y_o=np.array(es)
# print x_o
# print y_o
print "min energy: %s" %min(y_o)
#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
# Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print "Optimized parameters:", p1
#Plot it
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
# plt.xlim(1.2,1.8)
plt.xlabel(label)
plt.ylabel('FINAL SINGLE POINT ENERGY')
plt.show()
Посмотрим зависимость энергии структуры от длины С-С связи:
plot("bond/bond","C-C bond length",np.arange(1.329860,1.729860,0.02),[0.5,1, -79],"es_bond.txt")
Зависимость параболическая, но точки не четко совпадают с параболой. Видимо, межатомный потенциал на самом деле описывается более сложными функциями.
Теперь выполним аналогичные действия для валентного угла НСС. Варьируем угол от 109.2 до 113.2 и создаем 20 файлов с такими структурами этана.
lin='''!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 %f 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
*
'''
for i in range(-10,11):
ang=111.2+i*0.2
name="val/val"+str(i)
fi=open(name,'w')
fi.write(lin %(ang))
fi.close()
i=i+1
Смотрим зависимость энергии от значения угла НСС. Она параболическая:
plot("val/val","HCC angle",np.arange(109.2,113.2,0.2),[1,1, -79],"es_val.txt")
Теперь будем работать с торсионным углом НССН. Создаем 30 файлов с его значениями от -180 до 180:
lin='''!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 %f
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*
'''
for i in range(-15,16):
ang=0+i*12
print ang
name="tors/tors"+str(i)
fi=open(name,'w')
fi.write(lin %(ang))
fi.close()
i=i+1
Получем энергию этих структур с разными углами НССН:
open("es_tors.txt", 'w').close() #erase previous numbers in outfile
es=[]
for i in range(-15,16):
e=run_orca("%s%s" %("tors/tors",i))
print e
es.append(e)
#write energies to file
fi=open("es_tors.txt",'a')
fi.write(str(e)+"\n")
i=i+1
fi.close()
y_o=np.array(es)
И строим график зависимости. Как оказалось, зависимость лучше описывается синусоидой:
x_o=np.arange(-180,192,12)
#function is f(x)=k*cos(nx) + a
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
# Initial guess for the parameters
p0 = [1,1, -79]
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print "Optimized parameters:", p1
#Plot it
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.xlabel("HCCH torsion angle")
plt.ylabel('FINAL SINGLE POINT ENERGY')
plt.show()
Наблюдаем три минимума энергии при изменении торсионного угла HCCH от -180 до 180. Видимо, они соответствуют трем заторможенным конформациям молекулы этана. Атомы водорода одной СН3-группы располагаются "в промежутках" между водородами другой СН3-группы - это называется заторможенной конформацией.