NBO程序的.31至.40文件格式及轉換為.wfn文件的方法

我曾經寫過《高斯fch文件與wfn波函數文件的介紹及轉換方法》一文，見http://www.shanxitv.org/55，本文是對那篇文章內容的拓展，重復的內容就不再介紹了，建議先回顧一下。

1 .31至.40文件簡介

想獲得某種軌道的解析表達式，進而研究軌道特點，需要兩類信息：1 基函數信息 2 軌道的組合系數。在.wfn和.fch文件中都有這兩類信息。NBO程序產生的一系列plot文件也包含這樣的信息，后綴由.31到.40，這兩類信息不在同一個文件中。其中.31文件只包含基函數信息，是繪制各種軌道圖形所必需的，.32至.40只儲存了各種軌道向.31記錄的基函數的展開系數，各個后綴名記錄的軌道為：
32 PNAO
33 NAO
34 PNHO
35 NHO
36 PNBO
37 NBO  //最重要
38 PNLMO
39 NLMO
40 MO  //這就是普通的分子軌道

只要將.31文件的基函數信息和.32至.40文件中的組合系數信息相結合，就能計算相應軌道的各種屬性，最一般的自然就是獲得空間各點波函數數值，由此可以做出軌道圖形（故曰plot文件）。觀看NBO產生的各類軌道可以通過官方指定軟件NBOview（收費）以及筆者開發的Multiwfn 1.5（免費）、Chemcraft（收費）等，它們都需要載入這樣的文件。

plot文件可以由NBO或GENNBO生成。NBO程序分為兩種，一種是掛著ESS(electronic structure system，即各種量化程序)的，作為子程序來用，輸入是通過內部調用傳過去的，比如高斯的L607.exe。而GENNBO(General NBO)是可以獨立使用的NBO，不需要掛著ESS，輸入靠.47文件。GENNBO的功能由于沒有ESS提供的一些信息或功能而不得不受限，比如自然能量分解(NEDA)、刪除軌道、一些與單、雙電子積分或自洽場迭代相關的功能等等。

如果用高斯的NBO模塊生成plot文件，在route section里面寫上pop=nboread，說明要從輸入文件末尾讀入傳給NBO模塊的指令。然后分子坐標部分末尾空一行寫上比如$nbo plot file=c:\sob $end，這樣運行后就得到了plot文件，從c:\sob.31到c:\sob.40（盡管還有.41，但此文不涉及）。

如果用GENNBO生成plot文件，需要先獲得GENNBO輸入文件.47。在route section里面寫上pop=nboread，然后分子坐標部分末尾空一行上比如$nbo archive file=c:\sob $end，運行后就得到了c:\sob.47，在此文件開頭的$NBO和$END之間寫上plot，用GENNBO運行之就得到了plot文件。

2 .31的文件格式

.31文件的內容是.47文件的子集，在NBO程序手冊上雖然沒對.31文件格式有明確解釋，但通過對.47文件的介紹容易猜測到.31文件的格式。

.31文件開頭的三行沒用，
第4行：原子數 殼層數 原始層數
這里說的原始層是指《高斯fch文件與wfn波函數文件的介紹及轉換方法》一文的圖中橘紅色水平細線，每個殼層包括了多個原始層，每個原始層對應一種軌道指數和收縮系數。

第6行開始：元素序號，XYZ坐標（埃）

之后是類似下面的內容，有多少個殼層就有多少個類似的內容
      1     6    11     1
    201   204   206   202   203   205
其中的第一行是層所在的原子序號、層中基函數的數目、層中的第一個原始層的序號 收縮度
其中的第二行是此層里面每個基函數的類型的label，由于此例的層中有6個基函數，所以這里有6項。
各種label的含義見此圖


可見s的label為1，p、d、f、g軌道以100、200、300、400（g的情況在圖中沒列出）為基準然后依次遞增序號，對于球諧型d、f高斯函數則是從250、350開始遞增。此文中只考慮笛卡爾型GTF的情況。

.31文件的科學記數法記錄的內容分為五大部分，內容分別是原始殼層的指數、s型GTF的收縮系數、p型GTF的收縮系數、d型GTF的收縮系數、f型GTF的收縮系數，較新NBO版本生成的.31文件還會有g型GTF的收縮系數。每個部分的項數顯然就是總原始層數。注意收縮系數是已經乘了GTF的歸一化系數的，不像.fch里記錄的是沒乘歸一化系數的收縮系數。如果用的基函數不含SP殼層，同一原始層內將不會有不同類型的基函數，此時同一項不可能在多個部分同時不為0；如果包含SP殼層，同一原始層會同時包括s和p型基函數，所以此時同一項可以在s、p收縮系數部分同時不為0。

有一點要特別注意，對于d函數，相同原始層內（即具有相同的指數）的XX、YY、ZZ的歸一化系數是相同的，XY、XZ、YZ的歸一化系數也是相同的，然而第一種歸一化系數卻不同于第二種。對于.31文件里的d殼層的收縮系數，其中包含的歸一化系數是XX、YY、ZZ的，因此對于XY、XZ、YZ來說，想獲得正確的收縮系數，就必須將收縮系數除以XX的（或YY、ZZ的）然后再乘以XY的（或XZ、YZ）。對于f函數，.31文件中收縮系數中的歸一化系數是XXX、YYY、ZZZ的，獲得其它f函數的正確的收縮系數也需要像d函數那樣處理，即除一次再乘一次。

3 .32至.40的文件格式

.32至.40文件記載了每個軌道向.31文件記錄的基函數的展開系數。每5個數據換行一次，達到基函數數目后，即此軌道信息記錄完畢，即便此行沒到5個數據也換行，然后開始記錄下一個軌道。軌道數目與基函數數目一致。對于NBO(.37)、NLMO(.39)，在末尾還有些關于軌道特點的信息，包括占據數。

開殼層時，.34至.40文件內的不同自旋的軌道分別記錄，先記錄完所有ALPHA SPIN的再記錄BETA SPIN的，對于PNAO/NAO(.32,.33)不同自旋不分開記錄。注意開殼層體系Gaussian03的NBO3.1模塊輸出的.32和.33文件有bug，缺失了前三行，而GENNBO5.0生成的.32和.33文件格式是正確的。

4 讀入plot文件的代碼

這里介紹read31子程序的寫法，實際上這就是Multiwfn 1.5的sub.f90中的一個子程序，它能讀入.31以及.32至.40文件中的一個。和之前文章介紹的readfch子程序目的相似，要把讀入的文件內容轉換為.wfn文件里的信息的形式，這樣也方便計算屬性。read31的工作也就是要完成給以下變量賦值的任務：
nmo 總軌道數
nprims 總GTF數目
ncenter 總原子數
CO(:,:) 系數矩陣，第i,j個元素代表第i個軌道在第j個未歸一化的GTF上的展開系數
a(:) 儲存原子信息的數組，它是自定義類型，其中a%name,a%index,a%charge分別是相應原子的名稱、在周期表中的序號、核電荷數（使用贗勢時數值小于在周期表中的序號），a%x,a%y,a%z是原子XYZ坐標
b(:) 儲存每個GTF的信息的數組，b%exp,b%center,b%functype分別是相應GTF的指數、所在中心序號、函數類型標識。
MOocc(:) 記錄軌道的占據數。只有.37、.39文件中才記錄了這樣的信息而能被讀取。
還有幾個信息在plot文件中并沒有，也沒法算，但為了與.wfn文件中的信息一致，也要賦值，但內容是隨意的，包括totenergy=總能量；virialratio=維里系數；nelec=總電子數；MOene(:)=記錄軌道能量的數組；
另外MOtype(:)數組用來記錄相應軌道是什么類型，這在編寫其它代碼時會用到，0代表雙占據軌道，1是alpha自旋軌道，2是beta自旋軌道。

read31子程序代碼如下：
subroutine read31(name)
use define !此module包括如上所示的要被賦值的內容，各數組都是allocatable的狀態。
character(len=*) name !傳入的.31文件的文件名
character :: name2*200=" ",chartemp*80 !name2用來裝.32至.40文件名，chartemp是作為臨時判斷之用的字符串
logical alive
integer i,j,k,nbasis,nshell,nprimshell !nbasis=總基函數數目，nshell=總殼層數，nprimshell=總原始層數
integer bastype2func(415) !將NBO的基函數的label轉換為.wfn文件中的類型編號，415是最大的軌道序號，即ZZZZ型的g函數
integer,allocatable :: shellcon(:),shell2atom(:),shellnumbas(:),shell2prmshl(:),bastype(:) !殼層收縮度、殼層所在原子、殼層內基函數數目、殼層中第一個原始層的序號、基函數的label
real*8,allocatable :: orbcoeff(:,:),prmshlexp(:),cs(:),cp(:),cd(:),cf(:),cg(:) !以.31中基函數為基的軌道系數矩陣、原始層的指數、原始層中GTF類型為s,p,d,f,g時的收縮系數（已乘了歸一化系數）

!對照上面圖中label與GTF類型的對應關系，以及.wfn文件中如下所示的各類型的標識：
!1/2/3/4/5/6/7/8/9/10=S/X/Y/Z/XX/YY/ZZ/XY/XZ/YZ
!11/12/13/14/15/16/17/18/19/20=XXX/YYY/ZZZ/XXY/XXZ/YYZ/XYY/XZZ/YZZ/XYZ
!可以得到label到wfn類型標識的轉換關系
bastype2func(1)=1 !s
bastype2func(101)=2 !x
bastype2func(102)=3 !y
bastype2func(103)=4 !z
bastype2func(201)=5 !xx
bastype2func(202)=8 !xy
bastype2func(203)=9 !xz
bastype2func(204)=6 !yy
bastype2func(205)=10 !yz
bastype2func(206)=7 !zz
bastype2func(301)=11 !xxx
bastype2func(302)=14 !xxy
bastype2func(303)=15 !xxz
bastype2func(304)=17 !xyy
bastype2func(305)=20 !xyz
bastype2func(306)=18 !xzz
bastype2func(307)=12 !yyy
bastype2func(308)=16 !yyz
bastype2func(309)=19 !yzz
bastype2func(310)=13 !zzz
!.wfn文件沒有對g函數標識有統一的定義，以下20~35的序號順序設為與Gaussian的fch文件一樣。
bastype2func(401)=35 !XXXX
bastype2func(402)=34 !XXXY
bastype2func(403)=33 !XXXZ
bastype2func(404)=32 !XXYY
bastype2func(405)=31 !XXYZ
bastype2func(406)=30 !XXZZ
bastype2func(407)=29 !XYYY
bastype2func(408)=28 !XYYZ
bastype2func(409)=27 !XYZZ
bastype2func(410)=26 !XZZZ
bastype2func(411)=25 !YYYY
bastype2func(412)=24 !YYYZ
bastype2func(413)=23 !YYZZ
bastype2func(414)=22 !YZZZ
bastype2func(415)=21 !ZZZZ

write(*,*) "Input filename with suffix from .32 to .40 (e.g. ltwd.35)" !讀入.32至.40文件
do while(.true.)
    read(*,*) name2
    inquire(file=name2,exist=alive)
    if (alive.eqv..true.) exit
    write(*,*) "File not found, input again"
end do

tmplen=len_trim(name2) !輸出當前讀入的軌道類型
if (name2(tmplen-1:tmplen)=="32") write(*,*) "Loaded .32 file(PNAO)"
if (name2(tmplen-1:tmplen)=="33") write(*,*) "Loaded .33 file(NAO)"
if (name2(tmplen-1:tmplen)=="34") write(*,*) "Loaded .34 file(PNHO)"
if (name2(tmplen-1:tmplen)=="35") write(*,*) "Loaded .35 file(NHO)"
if (name2(tmplen-1:tmplen)=="36") write(*,*) "Loaded .36 file(PNBO)"
if (name2(tmplen-1:tmplen)=="37") write(*,*) "Loaded .37 file(NBO)"
if (name2(tmplen-1:tmplen)=="38") write(*,*) "Loaded .38 file(PNLMO)"
if (name2(tmplen-1:tmplen)=="39") write(*,*) "Loaded .39 file(NLMO)"
if (name2(tmplen-1:tmplen)=="40") write(*,*) "Loaded .40 file(MO)"

open(10,file=name,access="sequential",status="old")
read(10,*) !.31的前三行沒用
read(10,*)
read(10,*)
read(10,*) ncenter,nshell,nprimshell
!讀入原子數、殼層數、原始層數，之后就可以對數組分配內存了。但目前并不知道到底有多少基函數，但bastype由于馬上要用，所以取上限來分配內存，即假設最高角量子數的類型為g，又由于假設為笛卡爾型(15g)，故分配nshell*15。
allocate(a(ncenter),shellcon(nshell),shell2atom(nshell),shellnumbas(nshell),shell2prmshl(nshell))
allocate(prmshlexp(nprimshell),cs(nprimshell),cp(nprimshell),cd(nprimshell),cf(nprimshell),cg(nprimshell))
allocate(bastype(nshell*15))
read(10,*)
do i=1,ncenter !讀入原子信息
    read(10,*) a(i)%index,a(i)%x,a(i)%y,a(i)%z
    a(i)%charge=a(i)%index !.31并沒將核電荷數與原子在周期表中序號分別記錄，故直接將它們設為相同
end do
a%name=name2ind(a%index) !根據原子序號給原子賦上元素符號，name2ind數組在define module內定義。
a%x=a%x/b2a !.31內的坐標以埃為單位，轉換為波爾，b2a在define module內定義，為0.529177249D0
a%y=a%y/b2a
a%z=a%z/b2a
read(10,*)
j=1
do i=1,nshell !讀入每個殼層的信息
    read(10,*) shell2atom(i),shellnumbas(i),shell2prmshl(i),shellcon(i)
    read(10,*) bastype(j:j+shellnumbas(i)-1) !讀入此層中每個基函數的label
    j=j+shellnumbas(i)
end do
read(10,*)
read(10,*) prmshlexp !讀入每個原始層的軌道指數。用自由格式讀入。
read(10,*) !空行
read(10,*) cs !讀入每個原始殼層中GTF為s型時的收縮系數，接下來為p,d,f時的。
read(10,*)
read(10,*) cp
read(10,*)
read(10,*) cd
read(10,*)
read(10,*) cf
!較新的NBO程序輸出的.31還有cg段落記錄GTF為g型時的收縮系數，有則讀入。
read(10,"(a80)",iostat=ierror) chartemp
if (ierror==0) then
    backspace(10)
    read(10,*) cg
end if
close(10)

totenergy=0.0D0 !由于plot文件未記錄總能量、維里系數，故都設為0
virialratio=0.0D0
nbasis=sum(shellnumbas) !每個殼層內基函數數目加起來為總基函數
nprims=0
do i=1,nshell
    nprims=nprims+shellcon(i)*shellnumbas(i) !殼層內基函數乘以殼層收縮度并加和就是總GTF數。
end do

open(10,file=name2,access="sequential",status="old") !開始讀.32至.40文件
read(10,*)
read(10,*)
read(10,*)
read(10,"(a80)") chartemp !進行試探，如果此行內容是 ALPHA SPIN，說明是開殼層體系（但.32和.33并不區分體系類型）
if (chartemp(1:11)==" ALPHA SPIN") then
    wfntype=4 !代表這是開殼層
    nmo=2*nbasis !開殼層時軌道數為基函數數目二倍
    allocate(orbcoeff(nmo,nbasis),MOocc(nmo),MOene(nmo),MOtype(nmo),b(nprims),co(nmo,nprims))
    MOocc=1.0D0 !占據數這么設當然毫無道理，僅是示意
    MOtype(1:nbasis)=1 !前nbasis個軌道是alpha軌道
    read(10,*) ((orbcoeff(iorb,ibasis),ibasis=1,nbasis),iorb=1,nbasis) !把alpha軌道組合系數讀入系數矩陣的前nbasis行。
    if (name2(tmplen-1:tmplen)=="37".or.name2(tmplen-1:tmplen)=="39") read(10,*) MOocc(1:nbasis) !如果是.37/.39文件，還有占據數信息可讀入
    MOtype(nbasis+1:nmo)=2 !后nbasis個軌道是beta軌道
    call loclabel(10," BETA  SPIN") !loclabel見之前文章的介紹，此子程序用于將當前讀寫位置定位到文件中含有" BETA  SPIN"的行的行首
    read(10,*) !跳過" BETA  SPIN"字樣的行。
    read(10,*) ((orbcoeff(iorb,ibasis),ibasis=1,nbasis),iorb=nbasis+1,nmo) !把beta軌道組合系數讀入系數矩陣的后nbasis行。
    if (name2(tmplen-1:tmplen)=="37".or.name2(tmplen-1:tmplen)=="39) read(10,*) MOocc(nbasis+1:nmo)
else
    wfntype=3 !代表這是閉殼層
    nmo=nbasis !閉殼層時軌道數等于基函數數目
    allocate(orbcoeff(nmo,nbasis),MOocc(nmo),MOene(nmo),MOtype(nmo),b(nprims),co(nmo,nprims))
    MOocc=2.0D0
    MOtype=0 !軌道類型接設為無自旋型
    backspace(10) !回退一行，彌補讀入chartemp造成讀寫位置的下移。
    read(10,*) ((orbcoeff(iorb,ibasis),ibasis=1,nbasis),iorb=1,nmo)
    if (name2(tmplen-1:tmplen)=="37".or.name2(tmplen-1:tmplen)=="39) read(10,*) MOocc
end if
close(10)
MOene=0.0D0 !plot文件無軌道能量信息，都設為0
nelec=0
if (name2(tmplen-1:tmplen)=="37".or.name2(tmplen-1:tmplen)=="39) nelec=sum(MOocc) !為.37/.39文件時，MOocc有有意義值，故計算總電子數。否則設為0。

!下面把以.31內的基函數為基的系數矩陣orbcoeff變換到以GTF為基的系數矩陣CO，循環過程與fch轉wfn的過程類似。并且給記錄GTF信息的b向量賦值。
iGTF=1 !用于記錄當前GTF序號
ibasis=1 !用于記錄當前基函數序號
do i=1,nshell !循環每一殼層
    b(iGTF:iGTF+shellcon(i)*shellnumbas(i)-1)%center=shell2atom(i) !同一殼層中GTF所屬中心相同，故把這一殼層中所有GTF所屬中心都設好
    do j=1,shellnumbas(i) !循環殼層內每個基函數
        b(iGTF:iGTF+shellcon(i)-1)%functype=bastype2func(bastype(ibasis)) !同一基函數中的GTF類型都一樣，一次性設好。bastype(ibasis)是當前基函數的label，通過bastype2func數組轉換到wfn文件中的類型標識。
        do k=1,shellcon(i) !循環組成當前基函數的每個GTF
            iprmshlpos=shell2prmshl(i)+k-1 !代表當前GTF所在原始層的序號
            b(iGTF)%exp=prmshlexp(iprmshlpos) !給當前GTF以相應的指數
            if (bastype(ibasis)==1) then !根據當前基函數的label的不同，選擇類型對應的收縮系數向量，將其中對應的原始層的收縮系數取出。
                contract=cs(iprmshlpos)
            else if (bastype(ibasis)<=200) then !p函數
                contract=cp(iprmshlpos)
            else if (bastype(ibasis)<=300) then !d函數
                contract=cd(iprmshlpos)
            !202,203,205對應的是XY,XZ,YZ型GTF。上一行得到的收縮系數contract對XX,YY,ZZ是正確的（見本文第二節最后一段），下面將之轉換為對于XY,XZ,YZ正確的收縮系數
                if (bastype31(ibasis)==202.or.bastype31(ibasis)==203.or.bastype31(ibasis)==205) then
                    valnorm31=normgau(5,prmshlexp(iprmshlpos)) !normgau函數用于計算歸一化系數，這里獲得XX的歸一化系數
                    valnormnew=normgau(8,prmshlexp(iprmshlpos)) !獲得XY的歸一化系數
                    contract=contract/valnorm31*valnormnew !將XX,YY,ZZ的收縮系數變成XY,XZ,YZ的
                end if
            else if (bastype(ibasis)<=400) then !f函數
                contract=cf(iprmshlpos)
            !類似對d函數的處理，上一行得到的contract只是對XXX,YYY,ZZZ正確，下面將之轉換成其它類型f函數的收縮系數
                if (bastype31(ibasis)/=301.and.bastype31(ibasis)/=307.and.bastype31(ibasis)/=310) then !XXX,YYY,ZZZ以外的函數
                    valnorm31=normgau(11,prmshlexp(iprmshlpos)) !XXX,YYY,ZZZ的歸一化系數
                    if (bastype31(ibasis)==302.or.bastype31(ibasis)==303.or.bastype31(ibasis)==304&
                    .or.bastype31(ibasis)==306.or.bastype31(ibasis)==308.or.bastype31(ibasis)==309) then
                        valnormnew=normgau(14,prmshlexp(iprmshlpos)) !XXY,XXZ,XYY,XZZ,YYZ,YZZ的歸一化系數
                    else if (bastype31(ibasis)==305) then
                        valnormnew=normgau(20,prmshlexp(iprmshlpos)) !XYZ的歸一化系數
                    end if
                    contract=contract/valnorm31*valnormnew
                end if
            else if (bastype(ibasis)<=500) then
                contract=cg(iprmshlpos)
            !類似對d,f函數的處理，對收縮系數進行轉換。上一行得到的contract只是對XXXX,YYYY,ZZZZ正確
                if (nt/=401.and.nt/=411.and.nt/=415) then
                    valnorm31=normgau(21,prmshlexp(iprmshlpos)) !XXXX,YYYY,ZZZZ的歸一化系數
                    nt=bastype31(ibasis)
                    if (nt==402.or.nt==403.or.nt==407.or.nt==410.or.nt==412.or.nt==414) then
                        valnormnew=normgau(22,prmshlexp(iprmshlpos)) !XXXY,XXXZ,XYYY,XZZZ,YYYZ,YZZZ的歸一化系數
                    else if (nt==404.or.nt==406.or.nt==413) then
                        valnormnew=normgau(23,prmshlexp(iprmshlpos)) !XXYY,XXZZ,YYZZ的歸一化系數
                    else if (nt==405.or.nt==408.or.nt==409) then
                        valnormnew=normgau(27,prmshlexp(iprmshlpos)) !XYZZ,XYYZ,XXYZ的歸一化系數
                    end if
                    contract=contract/valnorm31*valnormnew
                end if
            end if
            CO(:,iGTF)=orbcoeff(:,ibasis)*contract !軌道向當前基函數的展開系數乘以當前GTF在此基函數中的收縮系數就是軌道向當前GTF的展開系數
            iGTF=iGTF+1 !更新當前的GTF的序號
        end do
        ibasis=ibasis+1 !更新當前的基函數的序號
    end do
end do
!輸出匯總信息
if (nelec==0) write(*,"(' There are',i6,' atoms,',i6,' basis,',i6,' orbitals')") ncenter,nbasis,nmo
if (nelec/=0) write(*,"(' There are',i6,' atoms,',i6,' basis,',i6,' orbitals,',f10.4,' electrons')") ncenter,nbasis,nmo,nelec
if (wfntype==4) then
    write(*,*) "This is open-shell system"
    write(*,"(' Orbitals from 1 to',i6,' are alpha spin, from',i6,' to',i6,' are beta spin')") nbasis,nbasis+1,nmo
end if
end subroutine