利用wfn文件計算電子密度的代碼的編寫方法
利用wfn文件計算電子密度的代碼的編寫方法
文/Sobereva @北京科音 2013-Apr-11
在《雜談Multiwfn從1.0到3.0版的開發經歷》（http://www.shanxitv.org/180）一文中曾經提到過，寡人三年多以前開發Multiwfn最初的動機是為了寫個帖子介紹怎么通過wfn文件里的信息計算電子密度，然而那帖子一直沒寫。正好最近有人在其代碼里需要計算電子密度，于是借這個機會寫個帖子介紹下怎么實現。
本文的代碼來自Multiwfn 3.0，為了簡明易懂進行了大量精簡，用于一般的wfn文件完全沒問題，但是無法像Multiwfn一樣兼容一些格式特殊、非標準的wfn文件，并且計算速度會略低于Multiwfn。如果已經讀懂了本文的代碼，并且想要更強的功能、更快的速度和更好的兼容性，可以再把完整的子程序從Multiwfn的源碼里提取出來。本文的代碼的Fortran90源文件可從此處下載：/usr/uploads/file/20150609/20150609205057_32281.f90。編譯在CVF6.5、Intel visual fortran上通過。
程序主要包含以下幾個部分：
module defvar：這個Module定義了一些重要的全局變量，諸如總電子數、原子信息、原始高斯函數(GTF)信息。另外也定義了兩種類型，一個叫atomtype，包含原子的名字、序號、坐標和核電荷數；另一個叫primtype，包含每個GTF的所屬中心、類型和指數。
module function：包含一個名為orbderv的子程序，用于計算指定位置的軌道波函數的數值及其導數；另一個是名為fdens的函數，用于計算指定位置的電子密度值。
program wfnprop：主程序，是一個很簡單的界面，用戶輸入wfn文件名字，以及要計算的點的位置，然后就返回電子密度。
readwfn子程序：用于從wfn文件中載入波函數信息，包括原子信息、GTF的定義、展開系數、軌道占據數等，并且由此確定體系中電子數、波函數類型等。
uc2lc和lc2uc子程序：用于進行字符的大小寫轉換。因為不同程序輸出的wfn中原子符號有的是大寫有的是小寫，這兩個字程序用來將大小寫統一，使得readwfn子程序總能夠根據原子符號正確辨認出原子序號。
以下是程序代碼，進行了適量注釋
module defvar
type atomtype !自定義的記錄原子信息的類型
character*2 name !原子名
integer index !原子序號。當使用了贗勢時原子核電荷將小于原子序號
real*8 x,y,z,charge !原子坐標(Bohr)及原子核電荷
end type
type primtype !自定義的記錄GTF信息的專用類型
integer center,functype !GTF所屬原子以及GTF的類型
real*8 exp !指數
end type
character*2 :: name2ind(0:109)=(/ "Bq","H ","He", & !原子名和原子序號的轉換表。0號是虛原子
"Li","Be","B ","C ","N ","O ","F ","Ne", & !3~10
"Na","Mg","Al","Si","P ","S ","Cl","Ar", & !11~18
"K ","Ca","Sc","Ti","V ","Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr", & !19~36
"Rb","Sr","Y ","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I ","Xe", & !37~54
"Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu", & !55~71
"Hf","Ta","W ","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn", & !72~86
"Fr","Ra","Ac","Th","Pa","U ","Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr", & !87~103
"Rf","Db","Sg","Bh","Hs","Mt" /) !104~109
!下面三行是GTF類型與GTF的x,y,z上的指數的轉換表。例如XY型GTF的x,y,z的指數分別是1,1,0，和下面三行的第8列對應，所以這樣的GTF的functype是8。
integer :: type2ix(35)=(/ 0,1,0,0, 2,0,0,1,1,0, 3,0,0,2,2,0,1,1,0,1, 0,0,0,0,0,1,1,1,1,2,2,2,3,3,4 /)
integer :: type2iy(35)=(/ 0,0,1,0, 0,2,0,1,0,1, 0,3,0,1,0,2,2,0,1,1, 0,1,2,3,4,0,1,2,3,0,1,2,0,1,0 /)
integer :: type2iz(35)=(/ 0,0,0,1, 0,0,2,0,1,1, 0,0,3,0,1,1,0,2,2,1, 4,3,2,1,0,3,2,1,0,2,1,0,1,0,0 /)
integer :: nmo=0,nprims=0,ncenter=0 !軌道數、GTF數、原子數
integer wfntype !0/1/2代表波函數是R/U/ROHF波函數，3/4分別代表閉殼層和開殼層后HF波函數
real*8 :: nelec=0,naelec=0,nbelec=0 !總電子數、alpha和beta電子數
type(atomtype),allocatable :: a(:) !記錄原子信息的數組
type(primtype),allocatable :: b(:) !記錄GTF信息的數組
real*8,allocatable :: MOocc(:),MOene(:) !軌道占據數和軌道能量
integer,allocatable :: MOtype(:) !記錄軌道類型。0/1/2分別代表無自旋/alpha/beta型軌道
real*8,allocatable :: CO(:,:) !軌道展開系數矩陣。CO(i,j)代表第j個GTF在第i號軌道中的展開系數，系數中已經把收縮系數、歸一化系數全都包含進去了。
module function
use defvar
implicit real*8 (a-h,o-z)
contains
! 當orbderv子程序的runtype==1，表明只計算x,y,z處的軌道波函數的數值，結果存在wfnval數組里，其第i個元素代表第i個軌道波函數的數值，并且此時此子程序的最后一個參數grad可以空缺。如果runtype==2，則代表不僅計算數值，還計算其導數，grad(1/2/3,i)代表第i號軌道波函數在x/y/z方向的導數。istart和iend代表計算的軌道序號范圍，一般都是1,nmo，也就是計算所有的軌道（顯然計算軌道的數目越少計算時間越短）。Multiwfn中完整的orbderv子程序還可以計算出軌道波函數的二、三階導數。
subroutine orbderv(runtype,istart,iend,x,y,z,wfnval,grad)
real*8 x,y,z,wfnval(nmo)
real*8,optional :: grad(3,nmo)
integer runtype,istart,iend
wfnval=0D0
if (present(grad)) grad=0D0
!雖然可以依次循環每個軌道計算它的波函數值及其導數，但是這樣計算量比較大。效率最好的方法是下面這樣循環每個GTF，得到它對各個軌道波函數及其導數的貢獻。
!GTF的表達式如下，式中的x,y,z是當前點相對于它所屬原子的原子核位置的坐標。式中r^2=x^2+y^2+z^2。ε代表GTF的指數。
!下面的代碼算出來的GTFval就是第j個GTF在輸入的x,y,z處的函數值
do j=1,nprims
ix=type2ix(b(j)%functype)
iy=type2iy(b(j)%functype)
iz=type2iz(b(j)%functype)
ep=b(j)%exp
sftx=x-a(b(j)%center)%x
sfty=y-a(b(j)%center)%y
sftz=z-a(b(j)%center)%z
rr=sftx**2+sfty**2+sftz**2
expterm=exp(-ep*rr)
GTFval=sftx**ix *sfty**iy *sftz**iz *expterm
!將GTF的數值乘上它在軌道中的展開系數就是它對軌道波函數的貢獻。這次把j號GTF對各個軌道的貢獻都算出來并累加進去。等全部GTF循環完之后，wfnval的各個元素就是各個軌道在當前點的波函數數值了
do imo=istart,iend
wfnval(imo)=wfnval(imo)+co(imo,j)*GTFval
end do
if (runtype==2) then !計算軌道波函數的導數。GTFdx、GTFdy、GTFdz就是當前GTF在x,y,z方向的導數
tx=0.0D0
ty=0.0D0
tz=0.0D0
if (ix/=0) tx=ix*sftx**(ix-1)
if (iy/=0) ty=iy*sfty**(iy-1)
if (iz/=0) tz=iz*sftz**(iz-1)
GTFdx=sfty**iy *sftz**iz *expterm*(tx-2*ep*sftx**(ix+1))
GTFdy=sftx**ix *sftz**iz *expterm*(ty-2*ep*sfty**(iy+1))
GTFdz=sftx**ix *sfty**iy *expterm*(tz-2*ep*sftz**(iz+1))
!將GTF導數的數值乘上它在軌道中的展開系數就是它對軌道波函數導數的貢獻
do imo=istart,iend
grad(1,imo)=grad(1,imo)+co(imo,j)*GTFdx
grad(2,imo)=grad(2,imo)+co(imo,j)*GTFdy
grad(3,imo)=grad(3,imo)+co(imo,j)*GTFdz
end do
end if
end do
end subroutine
!計算指定位置的電子密度，有了軌道波函數之后計算密度十分容易
real*8 function fdens(x,y,z)
real*8 x,y,z,wfnval(nmo)
call orbderv(1,1,nmo,x,y,z,wfnval) !先調用orbderv得到軌道波函數
fdens=0.0D0
do i=1,nmo
fdens=fdens+MOocc(i)*wfnval(i)**2 !軌道波函數的模方（由于這里用的都是實函數，所以等于平方）乘上軌道的占據數就是此軌道對總密度的貢獻。累加起來就是總電子密度。
end do
end function
end module
!====== 主程序
program wfnprop
use function
implicit real*8 (a-h,o-z)
character c80*80
write(*,*) "Input .wfn file name"
read(*,"(a)") c80
call readwfn(c80)
do while(.true.)
write(*,*) "Input the coordinate in Bohr e.g. 1.0,1.2,-3.5"
read(*,*) x,y,z
write(*,"(' The density is',f16.10,' e/Bohr^3',/)") fdens(x,y,z)
end do
end program
!從wfn中載入波函數信息的子程序。wfn文件的格式在《高斯fch文件與wfn波函數文件的介紹及轉換方法》（http://www.shanxitv.org/55）里面有詳細介紹，這里不再累述。
subroutine readwfn(name)
use defvar
implicit real*8 (a-h,o-z)
CHARACTER(LEN=*) name
character*80 wfntitle,c80tmp*80
open(10,file=name,access="sequential",status="old")
read(10,"(a80)") wfntitle
read(10,"(8x,i15,13x,i7,11x,i9)") nmo,nprims,ncenter
allocate(a(ncenter))
allocate(b(nprims))
allocate(co(nmo,nprims))
allocate(MOocc(nmo))
allocate(MOene(nmo))
allocate(MOtype(nmo))
!讀取原子信息
do i=1,ncenter
read(10,"(a24,3f12.8,10x,f5.1)") c80tmp,a(i)%x,a(i)%y,a(i)%z,a(i)%charge
read(c80tmp,*) a(i)%name
call lc2uc(a(i)%name(1:1)) !對原子名進行轉換，以確保第一個字符是大寫，第二個字符是小寫
call uc2lc(a(i)%name(2:2))
do j=0,109 !根據原子名確定原子序號
if (a(i)%name==name2ind(j)) then
a(i)%index=j
exit
end if
end do
end do
!讀取GTF定義
read(10,"(20x,20i3)") (b(i)%center,i=1,nprims)
read(10,"(20x,20i3)") (b(i)%functype,i=1,nprims)
read(10,"(10x,5D14.7)") (b(i)%exp,i=1,nprims)
!讀取軌道信息
do i=1,nmo
read(10,"(35X,f12.7,15x,f12.6)") MOocc(i),MOene(i)
read(10,"(5D16.8)") (CO(i,j),j=1,nprims)
end do
read(10,*)
close (10)
!根據軌道信息確定波函數類型和各個軌道的類型
if (sum(MOocc)==2*nmo) then !全部軌道占據數之和恰好是軌道數目的二倍，說明是閉殼層HF/DFT波函數
wfntype=0
MOtype=0 !所有軌道設為閉殼層類型
else if (sum(MOocc)==nmo) then !全部軌道占據數之和恰好是軌道數目，表明每個軌道占據數都是1，說明是開殼層HF/DFT波函數
wfntype=1
MOtype=1 !先假定所有軌道都是alpha軌道
do i=2,nmo !wfn中軌道能量是由低到高排序的。先記錄alpha再記錄beta。所以掃描各個軌道，如果發現當前軌道能量比上個軌道能量還小，表明從當前軌道起都是beta軌道
if (MOene(i)<MOene(i-1)) exit
end do
MOtype(i:nmo)=2
else if (any(MOocc/=int(MOocc))) then !只要有任何一個軌道占據數不為整數，就表明波函數是后HF波函數，軌道是自然軌道。接下來判斷是閉殼層還是開殼層后HF波函數
if (nint(maxval(MOocc))==2) then !占據數最大的自然軌道占據數接近于2，說明是閉殼層后HF波函數
wfntype=3
MOtype=0
else
wfntype=4 !開殼層后HF波函數
MOtype=0
do i=2,nmo
if (MOocc(i)>MOocc(i-1)) then !wfn文件中自然軌道占據數是從高往低排，先記錄alpha再記錄beta，如果發現當前軌道占據數比上個軌道的還大，表明從當前軌道起都是beta軌道
MOtype(1:i-1)=1
MOtype(i:nmo)=2
exit
end if
end do
end if
else
wfntype=2 !上述條件都不滿足，說明是ROHF/RODFT波函數
MOtype=0 !假設所有軌道都是雙占據
do i=1,nmo
if (MOocc(i)==1) MOtype(i)=1 !單占據軌道設為alpha軌道
end do
end if
!統計電子數
nelec=0
naelec=0
nbelec=0
do imo=1,nmo
if (MOtype(imo)==0) then
naelec=naelec+MOocc(imo)/2D0
nbelec=nbelec+MOocc(imo)/2D0
else if (MOtype(imo)==1) then
naelec=naelec+MOocc(imo)
else if (MOtype(imo)==2) then
nbelec=nbelec+MOocc(imo)
end if
end do
nelec=naelec+nbelec
!輸出匯總信息
write(*,*)
write(*,"('Total/Alpha/Beta electrons:',3f12.4)") nelec,naelec,nbelec
write(*,"('Net charge:',f12.5,' Expected multiplicity:',i5)") sum(a(:)%charge)-nelec,nint(naelec-nbelec)+1
write(*,"('The number of orbitals:',i6,', Atoms:',i7,', GTFs:',i7)") nmo,ncenter,nprims
if (wfntype==0) write(*,"('This is restricted close-shell single-determinant wavefunction')")
if (wfntype==1) write(*,"('This is unrestricted single-determinant wavefunction')")
if (wfntype==2) write(*,"('This is restricted open-shell wavefunction')")
if (wfntype==3) write(*,"('This is close-shell post-HF wavefunction')")
if (wfntype==4) write(*,"('This is open-shell post-HF wavefunction')")
if (wfntype==1.or.wfntype==4) then
do i=1,nmo
if (MOtype(i)==2) exit
end do
write(*,"('Orbitals from 1 to',i6,' are alpha type, from',i6,' to',i6,' are beta type')") i-1,i,nmo
end if
write(*,"(/,'Title line of this file:',/,a79,/)") wfntitle
end subroutine
!--------- 將字符轉換為小寫
subroutine uc2lc(inc)
character*1 inc
itmp=ichar(inc)
if (itmp>=65.and.itmp<=90) itmp=itmp+32
inc=char(itmp)
end subroutine
!--------- 將字符轉換為大寫
subroutine lc2uc(inc)
character*1 inc
itmp=ichar(inc)
if (itmp>=97.and.itmp<=122) itmp=itmp-32
inc=char(itmp)
end subroutine
將上述代碼進行編譯，啟動程序后，輸入wfn文件名，就會看到波函數的匯總信息（應確認信息無誤）。然后輸入坐標，就立刻輸出電子密度值了。只要把上述module和子程序移植到自己的程序里，就可以在自己的程序里調用fdens函數計算指定位置的密度值了。
利用orbderv子程序的信息，計算許多屬性都很容易。這里再給出兩個例子，將它們塞進module function里就可以在主程序中調用它們了。
!計算自旋密度的函數，也就是alpha密度減去beta的密度
real*8 function fspindens(x,y,z)
real*8 :: x,y,z,wfnval(nmo)
call orbderv(1,1,nmo,x,y,z,wfnval)
adens=0.0D0
bdens=0.0D0
do i=1,nmo
if (MOtype(i)==1) then !Alpha軌道
adens=adens+MOocc(i)*wfnval(i)**2
else if (MOtype(i)==2) then !Beta軌道
bdens=bdens+MOocc(i)*wfnval(i)**2
end if
end do
fspindens=adens-bdens
end function
!計算電子密度導數的函數，公式非常簡單，不妨自己推一下。label=="x","y","z","t"分別代表返回密度在x,y,z方向的導數和梯度的模
real*8 function fgrad(x,y,z,label)
real*8 x,y,z,wfnval(nmo),wfnderv(3,nmo),gradrho(3)
character label
call orbderv(2,1,nmo,x,y,z,wfnval,wfnderv)
gradrho=0D0
do i=1,nmo
gradrho(:)=gradrho(:)+MOocc(i)*wfnval(i)*wfnderv(:,i)
end do
gradrho=2*gradrho
if (label=='x') fgrad=gradrho(1)
if (label=='y') fgrad=gradrho(2)
if (label=='z') fgrad=gradrho(3)
if (label=='t') fgrad=dsqrt( sum(gradrho(:)**2) )
end function
Multiwfn的function.f90里還有很多計算實空間函數的子程序可以挪來使用。在Multiwfn程序手冊附錄2里有個列表。注意如果讀者自行開發的程序中使用了上述代碼或Multiwfn里的代碼，并且打算發表，請引用Multiwfn的原文Tian Lu, Feiwu Chen, J. Comp. Chem. 33, 580-592 (2012)。