分子轨道处理懒人脚本
One-Key Research Scripts for Molecular Orbitals
Jiawei Xu
Released: 2022-11-18 / Updated: 2022-11-18
当有大量分子轨道数据需要处理和作图时，懒人脚本是十分实用的。类似的做法也可以用于处理其他格点数据的作图。
1. 提取分子轨道能量
使用方法：python OrbEne.py Gaussian16.log [orblist.txt] [orbene.txt]
程序默认从当前目录下的orblist.txt获取需要输出能量的分子轨道序号，从Gaussian16.log获取结果输出到orbene.txt。orblist.txt的语法格式示例：
15 18 to 24 30 to 28 # Get MOs numbered with 15, 18-24 and 28-30.
h to h-5 l+1 # HOMO-5 to HOMO and LUMO+1
a h-3 to l+3 # Alpha-MOs from HOMO-3 to LUMO+3程序会自动识别是否为开壳层体系。可在orblist.txt中写入多行，程序会全部读取、整理排序后输出，并同时输出HOMO-LUMO gap信息。
import sys
def is_number(s):
try:
float(s)
return True
except ValueError:
pass
try:
import unicodedata
unicodedata.numeric(s)
return True
except (TypeError, ValueError):
pass
return False
def orb_ene(output,spin=0):
flag=False
occ,vir=[],[]
ab=['Alpha','Beta']
with open(output,'r') as read:
for line in read.readlines():
l=line.split()
if len(l)>4:
if l[0]=='The' and l[1]=='electronic' and l[2]=='state':
flag=True
occ,vir=[],[]
continue
if flag and l[0]==ab[spin] and l[1]=='occ.' and l[2]=='eigenvalues':
for i in range(len(l)-4):
occ.append(l[4+i])
continue
elif flag and l[0]==ab[spin] and l[1]=='virt.' and l[2]=='eigenvalues':
for i in range(len(l)-4):
vir.append(l[4+i])
continue
if l[0]=='Condensed' and l[1]=='to' and l[2]=='atoms':
flag=False
return occ,vir
def parse(occ,src):
res=0
mo1=['h','l']
mo2=['-','+']
for i in range(2):
if is_number(src):
res=int(src)
elif src[0]==mo1[i]:
if src==mo1[i]:
res=len(occ)+i
elif len(src)>2 and src[1]==mo2[i]:
dis=int(src[2:])
res=len(occ)+dis*(2*i-1)+i
return res
def relist(occ,listfile):
orblist=[[],[]]
tmp1,tmp2=0,0
spin=0
with open(listfile,'r') as read:
for line in read.readlines():
l=line.split()
if l[0]!='b':
pass
else:
spin=1
if len(l)!=0:
for i in range(len(l)):
if is_number(l[i]) or l[i][0]=='h' or l[i][0]=='l':
tmp1=parse(occ[spin],l[i])
orblist[spin].append(tmp1)
elif l[i]=='to':
tmp1,tmp2=parse(occ[spin],l[i-1]),parse(occ[spin],l[i+1])
for j in range(abs(tmp1-tmp2)-1):
orblist[spin].append(1+j+min(tmp1,tmp2))
elif i==0 and (l[i]=='a' or l[i]=='b'):
pass
else:
print('Error in input! Exit program')
exit()
read.close()
aorblist,borblist=orblist[0],orblist[1]
aorblist.sort(),borblist.sort()
return aorblist,borblist
def list_ene(output,listfile,writefile):
aocc,avir=orb_ene(output,0)
aorb=aocc+avir
bocc,bvir=orb_ene(output,1)
borb=bocc+bvir
occ=[aocc,bocc]
vir=[avir,bvir]
orb=[aorb,borb]
smallab,bigab=['a','b'],['Alpha','Beta']
orblist=relist(occ,listfile)
with open(writefile,'w') as ene:
for spin in range(2):
if len(bocc)==0 and spin==1:
break
orbs=orb[spin]
ene.write(bigab[spin]+'-HOMO is orbital '+str(len(occ[spin]))+smallab[spin]+'\n')
gap=float(orbs[len(occ[spin])])-float(orbs[len(occ[spin])-1])
ene.write(bigab[spin]+'-HOMO-LUMO gap is '+format(gap,'.5f')+' a.u. or '+format(gap*27.2114,'.5f')+' eV\n')
for i in range(len(orblist[spin])):
orbene_au=orbs[orblist[spin][i]-1]
orbene_ev=format(float(orbene_au)*27.2114,'.5f')
ene.write(str(orblist[spin][i])+smallab[spin]+'\t'+orbene_au+'\t'+orbene_ev+'\n')
ene.write('\n')
ene.close()
return
output=sys.argv[1]
listfile='orblist.txt'
try:
listfile=sys.argv[2]
except:
pass
writefile='orbene.txt'
try:
writefile=sys.argv[3]
except:
pass
list_ene(output,listfile,writefile)
输出示例：
Alpha-HOMO is orbital 217a
Alpha-HOMO-LUMO gap is 0.06212 a.u. or 1.69037 eV
213a -0.25040 -6.81373
214a -0.24144 -6.56992
215a -0.23823 -6.48257
216a -0.23769 -6.46788
217a -0.14421 -3.92416
218a -0.08209 -2.23378
219a -0.07925 -2.15650
220a -0.07452 -2.02779
221a -0.03748 -1.01988
222a -0.02886 -0.78532
Beta-HOMO is orbital 215b
Beta-HOMO-LUMO gap is 0.05099 a.u. or 1.38751 eV
211b -0.24697 -6.72040
212b -0.24529 -6.67468
213b -0.23739 -6.45971
214b -0.23734 -6.45835
215b -0.23716 -6.45346
216b -0.18617 -5.06595
217b -0.07925 -2.15650
218b -0.07677 -2.08902
219b -0.07473 -2.03351
220b -0.07442 -2.025072. 批量提取分子轨道能量
搭配1中脚本使用，一次性获取当前目录下所有.log文件记录的分子轨道信息。对于优化任务，默认获取最后一帧结构的信息。
for i in ./*.log
do
fname=`basename $i .log`
orbene $i ./orblist.txt ./orbene_$fname.txt
done
mkdir orb_ene
mv ./orbene_*.txt orb_ene/3. 批量导出分子轨道格点数据
前线轨道的分布是时常会关心的问题，以下脚本会调用Multiwfn（假定通过mwfn命令运行）产生.cub文件，接受一个参数N用于定义轨道范围。事实上也可以与1中脚本结合，使轨道范围定义更佳灵活，不过一般我们只关心前线轨道，使用一个参数定义已经足够。
使用方法：./FMO.sh N（产生HOMO-(N-1)至LUMO-(N-1)范围的分子轨道格点数据）
function judge_u(){
mwfn $1 > ./null << EOF
q
EOF
set `grep 'determinant' ./null`
rm -f ./null
if [ $4 = 'unrestricted' ]
then
return 1
else
return 0
fi
}
function get_orb(){
mwfn $1 > ./null << EOF
200
3
$2
3
0
q
EOF
rm -f ./null
}
mkdir fmos
for i in ./*fchk
do
fname=`basename $i .fchk`
mkdir fmos/$fname
judge_u $i
if [ $? = 0 ]
then
get_orb $i h
get_orb $i l
j=1
while ((j<$1))
do
get_orb $i h-$j
get_orb $i l+$j
((j++))
done
else
for s in {a,b}
do
get_orb $i h$s
get_orb $i l$s
j=1
while ((j<$1))
do
get_orb $i h$s-$j
get_orb $i l$s+$j
((j++))
done
done
fi
mv ./*.cub fmos/$fname
done4. 批量渲染分子轨道等值面图
该脚本需要安装VMD1.9.3和vcube2.0，可从http://bbs.keinsci.com/thread-18150-1-1.html下载。运行后自动调用VMD和vcube，使用Tachyon渲染器进行批量渲染。需要当前目录下有一个show.tcl文件用作可视化模板，使用时先对一个分子进行调整，效果满意后将操作过程通过VMD导出即可。
使用方法：./MO_Isosurface.sh
for i in ./*.cub
do
fname=`basename $i .cub`
vmd << EOF
vmol $i
`cat ./show.tcl`
vrender
n
quit
EOF
rm -f ./VCUBE/renderall.bat ./VCUBE/renderall.sh
mv ./VCUBE/renderall.dat $fname.dat
done
cd ./VCUBE
for i in ./*.dat
do
/usr/local/lib/vmd/tachyon_LINUXAMD64 $i -format BMP -o $i.bmp -trans_raster3d -res 2000 1500 -numthreads 16 -aasamples 24 -mediumshade
done
cd ..