我想生成一个像下面的链接图
http://en.wikipedia.org/wiki/Reaction_coordinate
该图是根据安装的python库的计算生成的。
我希望使用cspline gnuplot类型的线是平滑的
值E_ads = 234.4211,E_dis = 0.730278和E_reac = -0.8714
谁能帮我
from ase import *
from ase.calculators.jacapo import *
import Gnuplot as gp
# -- Read in all energies
datadict = {'H2O' :'water.nc',
'Pt' :'out-Pt.nc',
'H2OPt' :'H2O.Pt.nc',
'OHPt' :'OHPt.nc',
'HPt' :'HPt.nc',
}
E = {}
for label, file in datadict.items():
print 'Reading energy for %5s from file %20s' % (label, file),
atoms = Jacapo.read_atoms(file)
E[label] = atoms.get_potential_energy()
print '\tE = %14.6f eV'% E[label]
print
# -- Calculate adsorption and disassociation energies
E_ads = (E['H2OPt'] - 2*E['H2O'] - E['Pt'])/2
print 'H2O adsorption energy on Pt:'
print 'E_ads =', E_ads, 'eV\n'
E_dis = E['HPt'] - E['Pt'] + E['OHPt'] - E['Pt'] - E['H2O']
print 'H2O -> OH + H disassociation energy on Pt:'
print 'E_dis =', E_dis, 'eV\n'
E_reac = E['H2OPt'] - E['HPt'] - E['OHPt'] + E['Pt']
print 'H2O@Pt -> OH@Pt +H@Pt reaction energy on Pt:'
print 'E_reac =', E_reac, 'eV\n'
# -- Collect reaction path
Epath = np.asarray([1.0, E_ads, E_dis, E_reac])
PathLabels= ['']
# -- Plot the reaction path
import pylab as p
import numpy as np
import matplotlib.path as mpath
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
from scipy.interpolate import spline
import matplotlib.pyplot as plt
from numpy import array, linspace
from scipy.interpolate import spline
fig = p.figure(1)
sp = p.subplot(1,1,1)
p.plot(Epath, color='black', linestyle=':', linewidth=2.0, alpha=0.8)
p.text(0.37, 10.05, 'Free H$_2$O',fontsize=12, color='black',ha='right', va='bottom', alpha=1.0)
p.text(1.1, 238, 'H$_2$O + Pt',fontsize=12, color='black',ha='right', va='bottom', alpha=1.0)
p.title('H$_2$O disassociation')
p.ylabel('Energy [eV]')
p.xlabel('Reaction path')
#p.xlim([-0.5, 2.5])
#p.ylim([-0.5, 1.5])
sp.get_xaxis().set_ticks([]) # Turn off ticks on xaxis
#p.savefig('Teste.png')
p.show()
最佳答案
您可以通过执行以下操作来绘制数据的常规三次样条曲线版本:
plot(np.linspace(0,3),spline([0,1,2,3],Epath,np.linspace(0,3)))
这将产生类似以下内容:
但是我怀疑那不是您想要的。您可能需要求助于Monotone splines或shape preserving splines之类的形状,才能获得与Wikipedia链接中显示的曲线相同的形状。我不相信这些插值方法中的任何一种目前都在scipy中实现。
如果您对这些曲线的数学形式有一个大概的了解,则可以始终将自己的近似函数拟合为连续部分,而只需将函数限制在该范围之外即可。例如:
plot(np.linspace(0,3),np.maximum(E_react,spline([0,1,2,3],Epath,np.linspace(0,3))))
将产生:
即使它不是正确的拟合,它至少也“看起来”像链接到的曲线。