Gaussian型cube文件簡介及讀、寫方法和簡單應用

Gaussian型cube文件簡介及讀、寫方法和簡單應用

文/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://www.shanxitv.org/multiwfn）載入體系的波函數文件，然后用主功能5，按提示選擇被計算的函數、設置格點分布，算完后選擇導出cube文件即可。Multiwfn的相關知識看《Multiwfn入門tips
http://www.shanxitv.org/167》和《Multiwfn FAQ》（http://www.shanxitv.org/452）。Multiwfn手冊4.5.1節有個完整的計算格點數據和產生cube文件的例子。波函數文件的產生方法看《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/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軸，就能做出對應的截面圖了。