Gaussian型cube文件简介及读、写方法和简单应用
Gaussian型cube文件简介及读、写方法和简单应用
Introduction to Gaussian cube files, reading and writing methods and simple applications
文/Sobereva @北京科音
First release: 2012-Feb-8 Last update: 2022-Feb-18
1 前言
格点文件一般是指记录了三维空间中均匀分布的各个点的特定函数的数值的文件。在计算化学界最通用的格点文件格式是Gaussian型cube文件（以下简称cube文件）。它不仅可以由Gaussian程序产生，也可以由很多其它计算化学软件生成，如Multiwfn、CP2K、GROMACS的gmx spatial等等；它也可以被很多可视化软件所读取以显示等值面或进一步分析处理，如Multiwfn、VMD、ChemCraft、Chimera、Molekel、GaussView、VESTA等等。OpenDX格式（.dx）、.vti也是被一些化学软件所支持的格点文件格式，但用得远没cube格式普遍。本文将介绍cube文件的格式，通过一个简单的Fortran程序的代码介绍如何利用、读取和输出cube文件。熟练、灵活运用cube文件，可以解决不少计算化学上不方便处理的问题。
本文的代码不含注释的版本可在此下载：/usr/uploads/file/20150610/20150610023600_28748.f90。
如果读者仅仅是要获得某个体系的特定函数的cube文件，可以用Multiwfn（http://sobereva.com/multiwfn）载入体系的波函数文件，然后用主功能5，按提示选择被计算的函数、设置格点分布，算完后选择导出cube文件即可。Multiwfn的相关知识看《Multiwfn入门tips
http://sobereva.com/167》和《Multiwfn FAQ》（http://sobereva.com/452）。Multiwfn手册4.5.1节有个完整的计算格点数据和产生cube文件的例子。波函数文件的产生方法看《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）。
2 cube文件格式介绍
注意很多程序输出的cube文件并不很规范，本文严格按照Gaussian程序输出的cube文件格式进行说明。
cube文件包含了标题、平移矢量、原子坐标、格点数据这四个段落。记录的是标量数据，每个格点位置上只有一个数值。但还有一种特殊的cube文件，它多了分子轨道编号段落，允许在格点数据段落同时记录多条分子轨道波函数数值。
下面通过一个水分子的cube文件的开头部分来介绍具体格式。注意所有坐标、平移矢量单位都是bohr而不是埃。
Generated by Multiwfn //标题行，内容任意
Total 531440 grids //标题行，内容任意
3 -6.000000 -7.424912 -6.867160 //原子数（如果是记录分子轨道的cube文件则原子数为负值）；格点数据原点的X/Y/Z坐标
73 0.165751 0.000000 0.000000 //第一个平移矢量方向上有73个数据点；平移矢量x,y,z分量
91 0.000000 0.165751 0.000000 //第二个平移矢量方向上有91个数据点；同上
80 0.000000 0.000000 0.165751 //第三个平移矢量方向上有80个数据点；同上
8 8.000000 0.000000 0.000000 0.216790 //第一个原子的原子序号（氧）；原子核的有效电荷数（等于实际计算时原子的电子数）；X/Y/Z坐标。注意在使用赝势的时候由于一部分内核电子和核电荷的效果被赝势所代替，因此原子核的有效电荷数将小于原子序数
1 1.000000 0.000000 1.424912 -0.867160 //第二个原子的原子序号（氢）；同上
1 1.000000 0.000000 -1.424912 -0.867160 //第三个原子的原子序号（氢）；同上
//若是分子轨道cube文件（原子数为负值），这里会有分子轨道信息。第一个数字为此cube文件包含的分子轨道的数目，接下来是每个分子轨道的序号。比如此处有一行3 1 5 7就代表cube文件含3个分子轨道，分别是1、5、7号MO。分子轨道信息内容中每10个数字换一行
8.97452E-19 1.68963E-18 3.12519E-18 5.67888E-18 1.01380E-17 1.77805E-17 //每个格点位置上的数值，每六个换一行，1PE13.5格式，后同。
3.06365E-17 5.18606E-17 8.62458E-17 1.40910E-16 2.26177E-16 3.56662E-16
5.52546E-16 8.40976E-16 1.25748E-15 1.84723E-15 2.66590E-15 3.77982E-15
5.26502E-15 7.20496E-15 9.68649E-15 1.27939E-14 1.66014E-14 2.11635E-14
2.65054E-14 3.26125E-14 3.94217E-14 4.68156E-14 5.46196E-14 6.26051E-14
...（直到格点数据写完为止）
假设我们将三个平移矢量的x,y,z分量分别写为v1x,v1y,v1z、v2x,v2y,v2z、v3x,v3y,v3z，把三个平移矢量方向的数据点编号用i,j,k表示，把坐标原点写为orgx,orgy,orgz，那么
(i,j,k)点的x坐标=orgx+(i-1)*v1x+(j-1)*v2x+(k-1)*v3x
(i,j,k)点的y坐标=orgy+(i-1)*v1y+(j-1)*v2y+(k-1)*v3y
(i,j,k)点的z坐标=orgz+(i-1)*v1z+(j-1)*v2z+(k-1)*v3z
通常使用的cube文件是立方型的，也就是三个平移矢量分别平行于x,y,z轴，正如上面这个cube文件的例子。那么就可以简化写成
(i,j,k)点的x坐标=orgx+(i-1)*v1x
(i,j,k)点的y坐标=orgy+(j-1)*v2y
(i,j,k)点的z坐标=orgz+(k-1)*v3z
其中i的范围是从1~73，j是1~91，k是1~80。因此，(1,1,1)的位置就是(orgx,orgy,orgz)
格点数据部分记录顺序是先循环k，再循环j，之后循环i。每次将k循环完一遍后，即便那行未满6个数据，也会换一行。
对于记录多条分子轨道的cube文件，其循环顺序和上面一致，但是每个数据值被扩展为nmo个，nmo是记录的分子轨道数目。比如记录了三条分子轨道，则第一次循环k的段落可示意为
MO1(1,1,1) MO2(1,1,1) MO3(1,1,1) MO1(1,1,2) MO2(1,1,2) MO3(1,1,2)
MO1(1,1,3) MO2(1,1,3) MO3(1,1,3) MO1(1,1,4) MO2(1,1,4) MO3(1,1,4)
...
3 读取cube文件的代码
本文的例子程序叫做cubelite，由Fortran90语言编写，它可以载入并输出cube文件。主体代码如下：
!cube文件各种信息作为全局变量在defvar这个module里定义，还定义两种数据类型，atomtype用来记录原子信息，content用来记录格点的坐标和相应数据值
module defvar
type atomtype
integer index
real*8 x,y,z,charge
end type
type content
real*8 x,y,z,value
end type
integer n1,n2,n3,ncenter !三个方向的格点数，以及原子总数
real*8 orgx,orgy,orgz,v1x,v1y,v1z,v2x,v2y,v2z,v3x,v3y,v3z
type(content),allocatable :: cubmat(:,:,:) !记录格点数据
type(atomtype),allocatable :: a(:) !记录原子信息
end module
!主程序
program cubelite
use defvar
implicit real*8 (a-h,o-z)
character inpfile*200,outfile*200
write(*,*) "cubelite: Read, manipulate and write cube file"
write(*,*) "Programmed by Sobereva 2012-Feb-8"
write(*,*)
write(*,*) "Input the cube filename to be loaded"
read(*,"(a)") inpfile
call readcube(inpfile) !载入cube文件
!这里可以填入对cube文件的各种操作的代码，见第5节
write(*,*)
write(*,*) "Input the cube filename to be outputted"
read(*,"(a)") outfile
call outcube(outfile) !输出cube文件
pause
end program
下面介绍载入cube文件的子程序。注意在读取信息时都以自由格式而非固定格式来载入，这可以避免某些程序产生的cube文件的数据格式不规范而导致载入失败。
subroutine readcube(cubname)
use defvar
implicit real*8 (a-h,o-z)
character(len=*) cubname
integer,allocatable :: mo_serial(:)
real*8,allocatable :: temp_readdata(:)
open(10,file=cubname,access="sequential",status="old")
read(10,*)
read(10,*)
read(10,*) ncenter,orgx,orgy,orgz !载入总原子数，原点坐标
read(10,*) n1,v1x,v1y,v1z !第一个矢量方向的总格点数和矢量的x,y,z分量
read(10,*) n2,v2x,v2y,v2z
read(10,*) n3,v3x,v3y,v3z
nmo=0
if (ncenter<0) then
nmo=1 !由于总原子数为负值，必然此cube文件包含一条或一条以上分子轨道信息。在这里将分子轨道数nmo暂且设为1
ncenter=abs(ncenter) !将总原子数还原为正值
end if
allocate(a(ncenter)) !分配内存
allocate(cubmat(n1,n2,n3))
!输出格点文件的基本信息
write(*,"('Total number of atoms is',i8)") ncenter
write(*,"('Translation vector: x y z (Bohr)')")
write(*,"(a20,3f12.6)") "Vector 1: ",v1x,v1y,v1z
write(*,"(a20,3f12.6)") "Vector 2: ",v2x,v2y,v2z
write(*,"(a20,3f12.6)") "Vector 3: ",v3x,v3y,v3z
endx=orgx+(n1-1)*v1x+(n2-1)*v2x+(n3-1)*v3x
endy=orgy+(n1-1)*v1y+(n2-1)*v2y+(n3-1)*v3y
endz=orgz+(n1-1)*v1z+(n2-1)*v2z+(n3-1)*v3z
write(*,"('The range of x from ',f12.6,' to ',f12.6,' bohr,' i5,' points')") ,orgx,endx,n1
write(*,"('The range of y from ',f12.6,' to ',f12.6,' bohr,',i5,' points')") ,orgy,endy,n2
write(*,"('The range of z from ',f12.6,' to ',f12.6,' bohr,',i5,' points')") ,orgz,endz,n3
write(*,"('Total number of grid points is ',i10)") n1*n2*n3
do i=1,ncenter !载入原子信息
read(10,*) a(i)%index,a(i)%charge,a(i)%x,a(i)%y,a(i)%z
end do
if (nmo==1) then !此cube文件记录了分子轨道数据，需特殊处理
read(10,"(i5)",advance="no") nmo !获得实际包含的分子轨道数目
if (nmo>1) then !含有多个轨道，需要读入轨道编号
allocate(mo_serial(nmo))
allocate(temp_readdata(n3*nmo))
read(10,*) mo_serial !载入各个轨道编号
write(*,"('There are ',i5,' MOs in this grid file, the serial number are: ')") nmo
do i=1,nmo !输出各个轨道编号和选项，让用户选择载入其中哪个轨道的数据
write(*,"('Number ',i5,' : MO= ',i5)") i,mo_serial(i)
end do
write(*,*) "Which MO you want to load? Input the serial number"
read(*,*) mo_select !mo_select代表要载入第几个轨道。比如记录了12 57 88 102共四条轨道，输入3就表明要载入MO=88的轨道数据
else !只有1个轨道，不用读取分子轨道编号，因此空过此行
read(10,*)
end if
end if
write(*,*)
write(*,*) "Loading cube file, please wait..."
do i=1,n1 !循环第一个向量方向
do j=1,n2 !循环第二个向量方向
!一次性将第三个向量方向上的n3个数据都载入内存
if (nmo==0.or.nmo==1) then !如果记录的不是分子轨道数据或只含一个分子轨道，就直接将数据存到cubmat里即可
read(10,*) cubmat(i,j,:)%value
else !含有多个分子轨道，因此每次循环第三个向量方向时将有nmo*n3个数据，先载入到临时为位置temp_readdata
read(10,*) temp_readdata
cubmat(i,j,:)%value=temp_readdata(mo_select:size(temp_readdata):nmo) !从第mo_select号数据开始每经过nmo个数据就将数据从临时数组挪到cubmat里一次，直到nmo*n3号数据。这样就把指定的分子轨道数据载入了
end if
end do
end do
!将每个格点都赋予坐标信息
do i=1,n1
do j=1,n2
do k=1,n3
cubmat(i,j,k)%x=orgx+(i-1)*v1x+(j-1)*v2x+(k-1)*v3x
cubmat(i,j,k)%y=orgy+(i-1)*v1y+(j-1)*v2y+(k-1)*v3y
cubmat(i,j,k)%z=orgz+(i-1)*v1z+(j-1)*v2z+(k-1)*v3z
end do
end do
end do
write(*,*) "Loading completed!"
close(10)
end subroutine
4 输出cube文件的代码
这里介绍如何将内存中的格点数据输出到cube文件中，主程序中相应的outcube子程序如下。在输出的时候用的都是固定格式，与Gaussian程序输出的格式是完全一致的。我建议各种程序输出的cube文件格式都严格按照Gaussian的格式，以免发生不兼容。
subroutine outcube(cubname)
use defvar
implicit real*8 (a-h,o-z)
character cubname*200
open(10,file=cubname,access="sequential",status="replace")
write(10,"(' Generated by cubelite')")
write(10,"(' Totally ',i12,' grid points')") n1*n2*n3
write(10,"(i5,3f12.6)") ncenter,orgx,orgy,orgz
write(10,"(i5,3f12.6)") n1,v1x,v1y,v1z
write(10,"(i5,3f12.6)") n2,v2x,v2y,v2z
write(10,"(i5,3f12.6)") n3,v3x,v3y,v3z
do i=1,ncenter
write(10,"(i5,4f12.6)") a(i)%index,a(i)%charge,a(i)%x,a(i)%y,a(i)%z
end do
write(*,*) "Outputting cube file..."
do i=1,n1
do j=1,n2
do k=1,n3
write(10,"(1PE13.5)",advance="no") cubmat(i,j,k)%value !Gaussian输出时都是1PE13.5格式，即比如3.06365*10^(-17)这个数会输出为3.06365E-17，而不是E13.5格式时输出的0.30636E-16。因此相同宽度下1PE13.5比E13.5格式的有效位数上增加了一位
if (mod(k,6)==0.or.k==n3) write(10,*) !每满六个数据，或者第三个向量方向数据每次输出完毕时都换一行
end do
end do
end do
close(10)
write(*,*) "Outputting finished"
end subroutine
5 几个简单的使用格点数据的例子
注：以下的格点数据操作例子用Multiwfn的格点数据处理功能（主功能13）都能直接实现。启动Multiwfn，载入cube文件，按提示操作（不懂的话参考Multiwfn手册3.16节的相关部分的说明）即可，都不需要自己专门写程序。Multiwfn在格点数据处理方面极为灵活、强大。
这里给出几个使用格点数据的例子。注意编译器往往预留的堆栈内存空间不够，这会导致此cubelite程序在载入较大格点文件或者在数据处理时出现栈溢出的错误，因此需要适当加大栈内存空间。对于Intel visual fortran，方法是在“项目”-“cubelite属性”-Linker-System里将Stack Reserve Size设一个稍大的数值，比如设为268435456（即256MB）。对于Compaq visual fortran，则进入"Project"-"Settings..."-"Link", 在"Category"中选择"Output"，将"Reserve"填入0x10000000（转换为十进制为268435456，也是256MB）。
第一个例子是将记录了分子轨道波函数的格点文件变换为分子轨道电子密度的格点文件。方法很简单，就是将原先每个格点数值求平方，也就是将readcube和outcube之间插入这行代码即可：cubmat%value=cubmat%value**2
第二个例子是屏蔽掉格点文件的一部分区域。比如我们只想在可视化程序中显示出等值面的某个感兴趣的部分，比如x<=0的部分，那么可以将所有x>0的区域的数值设得非常大（比如1000），远大于原先格点数据数值范围，那么在显示的时候x>0的等值面区域将不会被显示出来，相当于被屏蔽掉了。在readcube和outcube之间插入这行代码即可：where (cubmat%x>0) cubmat%value=1000
第三个例子是将格点数据中指定平面上的数据伴随着坐标以文本文件方式输出出来，这样就可以在比如sigmaplot、surfer等软件中绘制出平面图了。此处的例子是让用户输入Z值，找出离输入Z值最近的XY平面，然后将这个平面上的数据输出到当前目录下的output.txt文件里。注意必须是立方格点的cube文件才能直接用此功能，否则若格点平移向量与笛卡尔坐标不平行就必须做额外的插值处理，会很麻烦而且精度低。具体做法是在readcube和outcube之间插入下面的代码：
b2a=0.529177249D0 !用于bohr和埃相互转换。在本程序内部一律以bohr为单位，在用户输入以及数据输出时变换为多数用户习惯的以埃为单位。
open(10,file="output.txt",status="replace")
write(*,*) "Input Z (in angstrom) to define a XY plane"
read(*,*) posZ
posZ=posZ/b2a !转换为bohr
rmindist=abs(cubmat(1,1,1)%z-posZ)
k=1 !先假设k=1的XY平面的z值与输入的z值最接近
do icycz=2,n3
disttmp=abs(cubmat(1,1,icycz)%z-posZ) !相同k值的格点数据的z值都是相同的，因此i,j值任取。这里就用i=1、j=1的那个数据点的z值作为当前XY平面的z值
if (disttmp<rmindist) then !k=icycz的XY平面的z值与输入的z值离得比当前最小值还小，故更新最小值和k
rmindist=disttmp
k=icycz
end if
end do
write(*,"(a,f14.6,a)") "The X-Y plane closest to your input is Z=",cubmat(1,1,k)%z*b2a," angstrom"
do i=1,n1
do j=1,n2
!输出离指定z值最近的XY平面的数据，前三列是数据坐标，以埃为单位，最后一列是数值
write(10,"(3f11.6,f22.15)") cubmat(i,j,k)%x*b2a,cubmat(i,j,k)%y*b2a,cubmat(i,j,k)%z*b2a,cubmat(i,j,k)%value
end do
end do
write(*,*) "The data have been exported to output.txt in current folder"
close(10)
将output.txt导入sigmaplot程序，做平面图，将第一、二列数据作为X、Y轴，将第四列作为Z轴，就能做出对应的截面图了。