layout: post
title: 绘制反应机理势能面的小工具
categories:
2015-07-15 10:03:33
缘起
在使用量子化学计算反应机理的论文中, 经常需要绘制反应的势能面(或称势能剖面图).
我以前也做过从头算方面的工作, 曾经写了一个Fortran小程序用来帮助绘制势能面.
主要思路是给出各个反应物种的能量及相对位置, 然后使用程序生成Origin的数据文件, 将数据文件导入Origin就可以直接绘出势能面了.
现在我基本已经不再做从头算方面的工作了, 所以这个小程序已经很长时间未用了, 一直躺在我电脑的某个文件夹中.
可时不时地, 我还会看到有人请教有了计算结果之后如何画反应的势能面. 网上给出的答案有ChemDraw, Origin甚至PowerPoint.
当然, 如果你有闲情逸致, 用Windows自带的画图软件也是可以做出很漂亮的图的. 但那不是我喜欢的方式.
既然我已经不再画势能面了, 这个小程序放在我这里也就没多大用处了. 所以, 我就花了点时间, 将其改造成了一个网上的在线工具,
你可以在直接在浏览器中编辑控制文件, 然后看到绘制效果. 程序还会给出用于Origin的数据文件, 将其导入Origin就可以使用Origin作图了.
我现在将这个小工具放在这里, 希望能对那些需要绘制势能面的人有所帮助.
如果你觉得在线使用不方便, 你可以直接将本页面保存在你的电脑上, 这样不需要上网, 在本机上也可以使用.
使用
在线工具的使用很简单, 在控制文件中编辑好每个翻译物种的编号, 名称, 相对位置, 能量, 与其连接的物种的编号, 然后点击绘图就可以看到绘制效果了.
下面的Origin单列数据和Origin多列数据两个文本框内给出了用于导入Origin的数据.
单列数据绘图时只绘制能量与连接线, 多列数据绘图时可将每一物种在图例中显示出来.
根据需要选择使用哪种数据. 具体区别见下面的说明.
对单列数据, 将其保存为文本文件后, 导入Origin后, 将Pos和adjPos列指定为X, 绘制Eng和adjEng列就可以了.
操作步骤如下:
对多列数据, 操作方法类似, 导入Origin后, 将Pos和adjPos列指定为X, 绘制每个物种对应的列和adjEng列就可以了.
在上面的图中, 使用了绝对能量(Hartree为单位)作图, 由于数值差异过小, Origin确定Y轴坐标范围出现问题, 所以建议使用以kcal/mol为单位的相对能量作图,
这样得到的相对能量数值上差异较大, 导入Origin后容易设置显示范围.
只要将基准能量设为反应物是能量-567.9303804, 能量因子设为627.5095重新生成数据就可以了. 这样得到的图就好多了.
你可以对得到的图进行修饰, 要改变反应物种的相对位置, 也可以直接在Origin中修改.
在线工具
<table>
<tr>
<td><div id="echarts" style="height:600px; width:600px"></div></td>
<td>控制文件(编号, 名称, 位置, 能量, 连接编号)<BR>
<textarea id="input" style="width:600px; height:550px;"></textarea><BR>
基准能量: <input type="text" id="subEng" value="0">
能量因子: <input type="text" id="facEng" value="1">
<input type="button" id="btn" value="绘图" onClick="plot()">
</td>
</tr>
<tr>
<td>Origin单列数据<BR><textarea id="singCol" style="width:600px; height:400px;"></textarea></td>
<td>Origin多列数据<BR><textarea id="multCol" style="width:600px; height:400px;"></textarea></td>
</tr>
</table><script src="http://echarts.baidu.com/build/dist/echarts.js"></script><script>
var myChart, option
require.config({ paths: {echarts: 'http://echarts.baidu.com/build/dist'} });
require( ['echarts', 'echarts/chart/line'],
function (ec) {
myChart = ec.init(document.getElementById('echarts'));
option = {
title: { text: '反应势能面' },
legend: { data:['相对能量'] },
tooltip: {trigger:'axis'},
toolbox: {
show: true,
feature: {
mark: {show: true},
dataZoom: {show: true},
dataView: {show: true, readOnly: false},
magicType: {show: true, type: ['bar','line']},
restore: {show: true},
saveAsImage: {show: true}
}
},
dataZoom: { show: true, realtime: true, start: 0, end: 100 },
xAxis: [ {type:'value', axisLine:{show: false}, axisLabel: {formatter: '{value}'}} ],
yAxis: [ {type:'value', axisLabel: {formatter: '{value}'} } ],
series: [
{ name:'相对能量', type:'line', data: [ [0,0] ] }
]
};
myChart.setOption(option);
}
);
var $=function(id){return document.getElementById(id)};
$('input').value=
"1 Rea 0 -567.9303804 2 5 8 11"
+"\n2 IM-1_1 1 -567.9351594 3"
+"\n3 TS-1 2 -567.9289591 4"
+"\n4 IM-1_61 3 -567.9852427 "
+"\n5 IM-2_1 1 -567.9365541 6"
+"\n6 TS-2 2 -567.9268951 7"
+"\n7 IM-2_61 3 -567.9480452 "
+"\n8 IM-3_61 1 -567.9364975 9"
+"\n9 TS-3 2 -567.9364548 10"
+"\n10 IM-3_1 3 -567.9568412 "
+"\n11 IM-4_1 1 -567.9309362 12"
+"\n12 TS-4 2 -567.920802 13"
+"\n13 IM-4_61 3 -567.94818 "
function plot() {
$('btn').value='正在绘图...'
var txt=$('input').value.replace(/^\s*\n*/,"").replace(/\s*\n*$/,"").replace(/\s+[\n|$]/g,"\n"),
txt=txt.split("\n"),
molNum=txt.length,
subEng=$('subEng').value,
facEng=$('facEng').value
var i, j, idx, jdx, data, adjNum, eng, minEng=1E99, maxEng=-1E99,
molIdx=[], molEng=[], molLab=[], minPos=[], maxPos=[], molAdj=[], txtEng=[], txtAdj=[]
option.series=[]
for(i=0; i<molNum; i++) {
data=txt[i].split(/\s+/)
idx=data[0]
molIdx[i]=idx
molLab[idx]=data[1]
minPos[idx]=parseFloat(data[2])
eng=(parseFloat(data[3])-subEng)*facEng
molEng[idx]=eng
molAdj[idx]=data.slice(4)
maxPos[idx]=minPos[idx]+.5
if(eng>maxEng) maxEng=eng
if(eng<minEng) minEng=eng
option.series.push({name: molLab[idx], type:'line', data: [ [minPos[idx],molEng[idx]],[maxPos[idx],molEng[idx]] ]})
}
option.yAxis[0].min=minEng
option.yAxis[0].max=maxEng
for(i=0; i<molNum; i++) {
idx=molIdx[i]
for(j=0; j<molAdj[idx].length; j++) {
jdx=molAdj[idx][j]
if(jdx) option.series.push({name: molLab[idx]+"-"+molLab[jdx], type:'line', itemStyle:{normal:{ lineStyle:{ color: '#000000' } }}, data: [ [maxPos[idx],molEng[idx]],[minPos[jdx],molEng[jdx]] ]})
}
}
$('btn').value='绘图'
for(i=0; i<molNum; i++) {
idx=molIdx[i]
txtEng[2*i-1] = molLab[idx]+','+minPos[idx]+','+molEng[idx]
txtEng[2*i ] = molLab[idx]+','+maxPos[idx]+','+molEng[idx]
}
adjNum=0
for(i=0; i<molNum; i++) {
idx=molIdx[i]
for(j=0; j<molAdj[idx].length; j++) {
jdx=molAdj[idx][j]
txt = molLab[idx]+'->'+molLab[jdx]+','
txtAdj[2*adjNum-1] = txt+maxPos[idx]+','+molEng[idx]
txtAdj[2*adjNum ] = txt+minPos[jdx]+','+molEng[jdx]
adjNum++
}
}
$('singCol').value=" "
+"\n Lab,Pos,Eng,adjLab,adjPos,AdjEng"
+"\n --,--,--,--,--,--"
+"\n --,--,--,--,--,--"
+"\n --,--,--,--,--,--"
+"\n --,--,--,--,--,--"
j=Math.min(molNum, adjNum)
for(i=0; i<j; i++) {
$('singCol').value +=
"\n"+txtEng[2*i-1]+","+txtAdj[2*i-1]
+"\n"+txtEng[2*i ]+","+txtAdj[2*i ]
+"\n--,--,--,--,--,--"
}
if(molNum>adjNum) {
for(i=j; i<molNum; i++) {
$('singCol').value +=
"\n"+txtEng[2*i-1]+",,,"
+"\n"+txtEng[2*i ]+",,,"
+"\n--,--,--,--,--,--"
}
} else {
for(i=j; i<adjNum; i++) {
$('singCol').value +=
"\n,,,"+txtAdj[2*i-1]
+"\n,,,"+txtAdj[2*i ]
+"\n--,--,--,--,--,--"
}
}
var coma=[]
i=molNum+10; while(i--) coma[i]=''
txt = ''
for(i=0; i<molNum; i++) { txt += molLab[molIdx[i]]+',' }
$('multCol').value = "\nLab,Pos,"+txt+'Lab,adjPos,adjEng'
txt="\n"+coma.slice(0, molNum+6)
$('multCol').value += txt+txt+txt+txt
for(i=0; i<molNum; i++) {
idx=molIdx[i]
txt=coma.slice(0,i+1)+molEng[idx]+coma.slice(0, molNum-i)
txtEng[2*i-1] = molLab[idx]+','+minPos[idx]+","+txt
txtEng[2*i ] = molLab[idx]+','+maxPos[idx]+","+txt
}
txt="\n"+coma.slice(0, molNum+5)
j=Math.min(molNum, adjNum)
for(i=0; i<j; i++) {
$('multCol').value +=
"\n"+txtEng[2*i-1]+","+txtAdj[2*i-1]
+"\n"+txtEng[2*i ]+","+txtAdj[2*i ]
+txt
}
if(molNum>adjNum) {
for(i=j; i<molNum; i++) {
$('multCol').value +=
"\n"+txtEng[2*i-1]+",,,"
+"\n"+txtEng[2*i ]+",,,"
+"\n"+coma.slice(0, molNum+2)
}
} else {
for(i=j; i<adjNum; i++) {
$('multCol').value +=
txt+txtAdj[2*i-1]
+txt+txtAdj[2*i ]
+txt
}
}
require('echarts').init(document.getElementById('echarts')).setOption(option);
myChart.setOption(option);
}
</script>相应的Fortran程序
如果你需要一个编译好的程序, 可以点击这里下载一个压缩包, 里面有Fortran的源代码, 编译好的程序, 两个示例控制文件及其输出文件.
Fortran源码是很久前写的了, 有点粗糙.
如果有什么问题, 你可以在下面留言.
2016-02-15 附记
最近看到一篇类似的文章在Origin中绘制能量折线图的方法, 里面提供了一些方法, 可供参考.