高斯fch文件與wfn波函數文件的介紹及轉換方法
了解高斯的格式化檢查點文件fch和波函數文件wfn的格式和每項的意義對編寫第三方量化輔助程序是很重要的，前者蘊含的信息豐富，后者更為簡明而且通用。本文將解讀wfn文件的格式和讀、寫代碼，對fch中的一些要點進行介紹，介紹如何寫出將fch轉化為wfn文件的程序fch2wfn，以幫助大家認識這兩種文件并利用fch和wfn寫出相關程序。
1. wfn(Wavefunction)文件的格式和意義
wfn文件最早作為Bader的AIMPAC程序的輸入文件出現，AIMPAC是歷史最久的AIM分析軟件，后來其它的AIM分析軟件如Multiwfn (http://www.shanxitv.org/multiwfn)、AIM2000、AIMALL、Morphy、Xaim等等也都需要wfn文件。wfn文件包含的是所有高斯函數的信息和波函數向它們的展開系數，結構簡潔，計算空間中波函數值、電子密度等屬性很方便。Gaussian、GAMESS等許多主流量化程序都可以輸出wfn文件，例如Gaussian中只需要在route section加上output=wfn，并在分子結構部分末尾空一行寫上wfn文件輸出路徑即可，如c:\divokej_bill.wfn。
這里以高斯中HF/3-21G計算的CH4分子輸出的wfn文件為例對格式進行介紹
第1行
Title Card Required
wfn中第一行沒有任何意義，只是文件注釋
第2行
GAUSSIAN 5 MOL ORBITALS 27 PRIMITIVES 5 NUCLEI
開頭的GAUSSIAN說明文件中記錄的函數類型是高斯類型，一般wfn中都是高斯類型函數。
5 MOL ORBITALS代表文件中的分子軌道數，wfn中只記錄占據軌道，虛軌道不記錄，因為AIM分析的對象是電子密度，而虛軌道對計算電子密度沒有貢獻。
27 PRIMITIVES代表原始高斯函數的總數目，可獨立變分的基函數是由它們收縮成的。
5 NUCLEI是原子數
第3-7行
C 1 (CENTRE 1) 0.00000000 0.00000000 0.00000000 CHARGE = 6.0
H 2 (CENTRE 2) 1.16740622 1.16740622 1.16740622 CHARGE = 1.0
H 3 (CENTRE 3) -1.16740622 -1.16740622 1.16740622 CHARGE = 1.0
H 4 (CENTRE 4) -1.16740622 1.16740622 -1.16740622 CHARGE = 1.0
H 5 (CENTRE 5) 1.16740622 -1.16740622 -1.16740622 CHARGE = 1.0
這些是原子的信息。例如其中的第一行，C是原子符號，后面的1是原子序號。CENTRE后面的也是原子序號，與前面的序號是重復的，所以取其一即可。接下來是原子坐標，單位為波爾，最后是原子的核電荷數。使用贗勢時核電荷數不等于原子序數，而等于被明確表達的價層電子數。
第8、9行
CENTRE ASSIGNMENTS 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 3 3
CENTRE ASSIGNMENTS 3 4 4 4 5 5 5
每項代表每個高斯函數所在的原子序號。由于一共有27個高斯函數，這里也就有27項。
第10、11行
TYPE ASSIGNMENTS 1 1 1 1 1 2 2 3 3 4 4 1 2 3 4 1 1 1 1 1
TYPE ASSIGNMENTS 1 1 1 1 1 1 1
這是每個高斯函數的類型。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
生成wfn文件時應注意應使用6d/10f型高斯函數而不要用球諧型。wfn文件對于g及更高角動量函數沒有統一的序號定義，而且由于涉及這樣的函數機會不多，所以多數AIM軟件還不支持這些函數，如果體系中有這樣的函數，高斯也無法輸出wfn文件。
第12-17行
EXPONENTS 0.1722560D+03 0.2591090D+02 0.5533350D+01 0.3664980D+01 0.7705450D+00
EXPONENTS 0.3664980D+01 0.7705450D+00 0.3664980D+01 0.7705450D+00 0.3664980D+01
EXPONENTS 0.7705450D+00 0.1958570D+00 0.1958570D+00 0.1958570D+00 0.1958570D+00
EXPONENTS 0.5447178D+01 0.8245472D+00 0.1831916D+00 0.5447178D+01 0.8245472D+00
EXPONENTS 0.1831916D+00 0.5447178D+01 0.8245472D+00 0.1831916D+00 0.5447178D+01
EXPONENTS 0.8245472D+00 0.1831916D+00
這是每個高斯函數的指數，即N*x^ix*y^iy*z^iz*exp(-α*r^2)中的α。顯然也是27項。
第18-52行
......
MO 3 MO 0.0 OCC NO = 2.0000000 ORB. ENERGY = -0.548298
0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00
0.00000000D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.64363172D+00
0.33350218D+00 0.00000000D+00 0.00000000D+00 0.00000000D+00 0.57302486D-01
0.63659475D-01 0.89429130D-01 0.29385384D-01 0.63659475D-01 0.89429130D-01
0.29385384D-01 -0.63659475D-01 -0.89429130D-01 -0.29385384D-01 -0.63659475D-01
-0.89429130D-01 -0.29385384D-01
......
這部分是分子軌道信息，這里只列出第3條以說明。
MO 3代表是第3條軌道，也就是它在高斯中的軌道序號。對于閉殼層，這個序號和此軌道在本文件中的序號一致。
MO 0.0沒有任何意義，只是個標簽。
NO = 2.0000000是占據數，這里就是代表有兩個電子占據，若是非限制性波函數即為1。若是限制性開殼層計算，高斯會誤把單占據軌道占據數寫為2，需手動修改為1。
ORB. ENERGY = -0.548298是軌道能，單位為Hartree。
接下來就是分子軌道向這27個高斯函數的展開系數，一定要注意，每個展開系數已經包含了高斯函數的歸一化系數，也就是分子軌道向沒有歸一化的高斯函數（即x^ix*y^iy*z^iz*exp(-α*r^2)形式）展開的系數。
這里的軌道并非一定是分子軌道，如果用后HF方法，并且加了density=current表明用后HF密度，則輸出的每個軌道就是自然軌道，其占據數一般不是整數，輸出的軌道數目與可獨立變分的基函數一致。因為后HF密度的自然軌道沒有是否占據之分，只有占據多少之分，每條軌道都對電子密度有貢獻，所以所有的自然軌道都會列出。如果是開殼層計算，還應當加上pop=NOAB，這樣就會記錄自然自旋軌道。
軌道記錄順序是能量由低往高。注意對于自旋軌道，wfn中先按這個次序記錄完所有alpha軌道再記錄beta軌道。由于第一個beta軌道的能量肯定比最后的alpha軌道的能量低，所以通過檢驗軌道能量與前一個軌道能量的關系就能判斷beta軌道從哪里開始記錄。如果是自然自旋軌道，則要判斷它與前一個軌道的占據數關系。
wfn中軌道的軌道序號，記錄的前后順序是隨意的。高斯函數的順序也是隨意的，只要在類型、所在中心、指數、軌道展開系數段落中的位置對應就行，不會影響結果。
第53行
END DATA
其中END DATA代表標準wfn文件要求的全部信息的說明已結束，在此之后允許附加一些信息。
第54行及之后
THE HF ENERGY = -39.976405845082 THE VIRIAL(-V/T)= 2.00033480
這些是附加信息，是隨意的。一般量化程序輸出的wfn都會按照這個格式輸出能量和維里值。THE HF ENERGY是體系HF或DFT總能量，如果是后HF方法，這里仍然是HF能量。有的程序還會附上其它一些信息。
2.讀、寫wfn文件的fortran代碼解析
首先根據wfn所擁有的數據，定義一些變量，以及兩種自定義類型用于儲存原子和高斯函數的數據。
type atomtype
character*2 name !元素名稱
real*8 x,y,z,charge !原子坐標和電荷
end type
type basistype
integer center,functype !高斯函數所在原子序號，函數類型
real*8 exp !指數
end type
character*8 mode !用來記錄函數類型，一般wfn里都是高斯函數，所以可以不定義。
integer nmo,nprims,ncenter !軌道、高斯函數、原子的數目
real*8 totenergy,virialratio !總能量、維里值
type(atomtype),allocatable :: a(:) !儲存原子信息
type(basistype),allocatable :: b(:) !儲存基函數信息
real*8,allocatable :: MOocc(:),MOene(:) !軌道占據數和能量
real*8,allocatable :: CO(:,:) !系數矩陣，CO(i,j)代表第i個分子軌道向第j個高斯函數的展開系數，已包括歸一化系數。
還定義兩個變量，并不是直接從wfn文件中讀取，需要用代碼來判斷，在寫一些相關程序時會用得到它們。
integer,allocatable :: MOtype(:) !軌道類型，0代表雙占據的空間軌道，1是alpha自旋軌道，2是beta自旋軌道
integer wfntype !用來說明波函數類型，0/1/2代表R/U/ROHF波函數，3/4代表是閉殼層/開殼層后HF波函數，分別代表儲存的是空間自然軌道和自旋自然軌道。
下面是讀取wfn文件的子程序，是作為主程序的內部函數使用的，能共享前面在主程序中定義的變量，主程序中已使用了隱式聲明implicit real*8(a-h,o-z)。
subroutine readwfn
open(10,file=filename,access="sequential",status="old")
read(10,*) !跳過wfn開頭
read(10,"(a8,i15,13x,i7,11x,i9)") mode,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,"(2x,a2,20x,3f12.8,10x,f5.1)") a(i)%name,a(i)%x,a(i)%y,a(i)%z,a(i)%charge
end do
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,*)
read(10,"(18x,f19.12,20x,f12.8)") totenergy,virialratio !總能量和維里值
close (10)
write(*,"('There are',i6,' MOs,',i6,' atoms,',i7,' electrons,',i7,' Gaussian functions')") nmo,ncenter,int(sum(moocc)),nprims
! 下面分析波函數的類型
if (sum(MOocc)==2*nmo) then !總電子數是軌道數二倍（再次注意wfn只記錄占據軌道），為HF閉殼層
wfntype=0
MOtype=0 !軌道類型都設雙占據
write(*,"('This is restricted wavefunction')")
else if (sum(MOocc)==nmo) then !總電子數等于軌道數，為HF開殼層
wfntype=1
MOtype=1 !先都設為alpha
do i=2,nmo
if (MOene(i)<=MOene(i-1)) exit !確定哪里是alpha軌道與beta軌道的分界處
end do
MOtype(i:nmo)=2 !將分界處后面的軌道都設成beta
write(*,"('This is unrestricted wavefunction,',i6,' alpha orbitals',i6,' beta orbitals')") i-1,nmo-i+1
else if (MOocc(nmo)/=int(MOocc(nmo))) then !最后一個軌道占據數不為整數，為后HF波函數
wfntype=3 !先假設是閉殼層
MOtype=0
do i=2,nmo
if (MOocc(i)>MOocc(i-1)) then !尋找alpha與beta軌道分界位置，如果能找到，說明是開殼層
wfntype=4
MOtype(1:i-1)=1 !給自旋自然軌道賦予alpha/beta類型
MOtype(i:nmo)=2
exit
end if
end do
if (wfntype==3) write(*,"('This is restricted Post-HF wavefunction')")
if (wfntype==4) write(*,"('This is unrestricted Post-HF wavefunction')")
if (wfntype==4) write(*,"(i6,' alpha electrons',i6,' beta electrons')") int(sum(MOocc(1:i-1))),int(sum(MOocc(i:nmo)))
else
wfntype=2 !以上情況都不是，就是限制性開殼層波函數
MOtype=0
j=0
do i=1,nmo
if (MOocc(i)==1) then
MOtype(i)=1 !單占據軌道都設成alpha類型
j=j+1
end if
end do
write(*,"('This is restricted open-shell wavefunction,',i6,' single occupied orbitals')") j
end if
write(*,"('Total energy:',f19.12,' hartree, Virial ratio:',f12.8)") totenergy,virialratio
write(*,*)
end subroutine
了解了wfn讀入方法，輸出wfn也就會了，以下是輸出子程序。
subroutine outwfn
open(11,file=outfilename,status="replace")
write(11,*) "Generated by fch2wfn"
write(11,"('GAUSSIAN',i15,' MOL ORBITALS',i7,' PRIMITIVES',i9,' NUCLEI')") nmo,nprims,ncenter
do i=1,ncenter
write(11,"(2x,a2,i4,4x,'(CENTRE',i3,')',1x,3f12.8,' CHARGE =',f5.1)") a(i)%name,i,i,a(i)%x,a(i)%y,a(i)%z,a(i)%charge
end do
write(11,"('CENTRE ASSIGNMENTS ',20i3)") b%center
write(11,"('TYPE ASSIGNMENTS ',20i3)") b%functype
write(11,"('EXPONENTS ',5D14.7)") b%exp
do i=1,nmo
write(11,"('MO',I5,' MO 0.0 OCC NO = ',f12.7,' ORB. ENERGY =', f12.6)") i,MOocc(i),MOene(i)
write(11,"(5D16.8)") (co(i,j),j=1,nprims)
end do
write(11,"('END DATA',/,' THE HF ENERGY = ',f19.12,' THE VIRIAL(-V/T)= ',f12.8)") totenergy,virialratio
close(11)
end subroutine
3. fch(Formatted check point)文件簡介
任何平臺下高斯執行任務生成的二進制的chk文件都可以用formchk生成fch/fchk文件，fch文件是ASCII文件，是人可讀的。它獨立于平臺，可用于在不同平臺上交換chk文件，也可方便地從中提取計算過程中有用的信息。fch文件在Gaussian03與09中內容大體不變，但是記錄項目的順序有所改變。
fch中的信息遠比wfn豐富，包含了獲得wfn文件所需要的一切信息，里面也有不少冗余的信息。本文主要目的是介紹fch向wfn的轉換方法，如何寫fch2wfn程序，下面只談一些相關要點，一些條目靠fch自帶的解釋就可以明白就不談了。
第一行是標題，是任意的，可以自行寫一些附加信息在此告訴fch2wfn程序對它做哪些特殊處理。
第二行是任務類型、方法、基組，格式是A10,A30,A30。看第11列是R、U還是RO就能判斷波函數類型。
每一類數組型數據前面都有一行標簽，讀取格式都是A40,3X,A1,5X,I12。比如
Primitive exponents R N= 18
18代表下面有18個數據。R代表下面的是實數型數據。若是I、C、L則分別代表是整型、字符型、邏輯型數據。實數型數據讀取格式是5(1PE16.8)，整型是6I12。
fch中與wfn顯著不同的一點是基函數的儲存方法。wfn中是按高斯函數順序儲存的，而fch中是一層一層地儲存的。Shell types代表了每個殼層的類型，0=s, 1=p, -1=sp, 2=6d, -2=5d, 3=10f, -3=7f, 4=15g, -4=9g。Number of primitives per shell代表的是每一殼層內的基函數的收縮度。
fch中每個殼層里笛卡爾型基函數的類型是按如下順序排列S,X,Y,Z, XX,YY,ZZ,XY,XZ,YZ, XXX,YYY,ZZZ,XYY,XXY,XXZ,XZZ,YZZ,YYZ,XYZ，對比wfn中的函數類型編號，會發現f軌道的順序不一致，這在轉換時需要特別注意。如前所述，wfn文件里不能記錄g及更高角動量函數，但是在fch文件中可以，g函數的順序為ZZZZ,YZZZ,YYZZ,YYYZ,YYYY,XZZZ,XYZZ,XYYZ,XYYY,XXZZ,XXYZ,XXYY,XXXZ,XXXY,XXXX。
Coordinates of each shell其實是多余的，因為已經有Shell to atom map（殼層所屬原子）和Current cartesian coordinates（每個原子xyz坐標），可以直接推出。
Total SCF Density就是SCF密度矩陣，用了后HF方法并用density=current還會有后HF的密度矩陣被記錄，如Total MP2 Density。由于是對稱的，所以只記錄了下三角部分（含對角元），順序是(1,1)(2,1)(2,2)(3,1)(3,2)(3,3)(4,1)...，
限制性方法計算的軌道能量和展開系數記錄在Alpha Orbital Energies/Alpha MO coefficients里（盡管雙占據稱為Alpha軌道并不妥），如果是非限制性的計算還會多出Beta的。用density=current并不會使fch文件像wfn一樣使記錄的軌道成為自然軌道，它仍然是SCF的系數。如果想把自然軌道系數寫入到Alpha MO coefficients當中，則先用后HF方法帶著density=current關鍵字跑一次，然后再運行Guess=(Save,Only,NaturalOrbitals) ChkBasis即可。此時fch里軌道能量部分就不再是能量，而是占據數。
加了pop=saveNBO使高斯調用NBO 3.1進行NBO分析后，會把NBO軌道保存到fch的分子軌道系數部分。如果用的SCF密度，則分子軌道能量部分保存的是NBO軌道能量，讀入下文的fch2wfn程序前應把第一行改為saveNBOene；如果用的是后HF密度，這部分保存的則是NBO軌道占據數，讀入fch2wfn前應把第一行改為saveNBOocc。fch中倒數的幾個NBO軌道的能量或占據數誤為1000，應根據NBO分析輸出的結果手動修改。
Primitive exponents、Contraction coefficients與P(S=P) Contraction coefficients的關系：
這里舉HF/STO-3G算CH4的fch的這部分作為例子，這三類只列出第一行來說明，最開頭的就是C的S和SP殼層。
Primitive exponents R N= 18
7.16168373E+01 1.30450963E+01 3.53051216E+00 2.94124936E+00 6.83483096E-01
...
Contraction coefficients R N= 18
1.54328967E-01 5.35328142E-01 4.44634542E-01 -9.99672292E-02 3.99512826E-01
...
P(S=P) Contraction coefficients R N= 18
0.00000000E+00 0.00000000E+00 0.00000000E+00 1.55916275E-01 6.07683719E-01
...
下面是一些linux版高斯自帶的sto3g.gbs文件，
-C
S 3 1.00
0.7161683735D+02 0.1543289673D+00
0.1304509632D+02 0.5353281423D+00
0.3530512160D+01 0.4446345422D+00
SP 3 1.00
0.2941249355D+01 -0.9996722919D-01 0.1559162750D+00 注：第一列是高斯函數指數，后兩列分別是S和P型的收縮系數
0.6834830964D+00 0.3995128261D+00 0.6076837186D+00
0.2222899159D+00 0.7001154689D+00 0.3919573931D+00
將fch的內容豎著看，而將sto3g.gbs的內容橫著看，可見fch中與sto3g.gbs的內容是對應的。對于S殼層，P(S=P) Contraction coefficients對應的項都是0，到了SP殼層就不為0了，代表的就是P的收縮系數。如果用的是例如Dunning相關一致性基組，由于S和P殼層是獨立的，沒有SP殼層，則這部分就不會出現。
4.讀取fch文件的fortran代碼
為了完成fch->wfn的轉換，不需要讀入所有fch的數據，只需要讀入其中有用的數據并經過轉換，把將前面outwfn子程序要用的變量全都填好就可以了。讀取fch文件的代碼中要用到在fch中定位標簽的子程序，例如call loclabel(10,'Alpha MO coefficients')就會使當前讀取位置處在Alpha MO coefficients這行開頭，如果沒找到要找的內容，則全局變量noentryinfch=1，否則為0。loclabel子程序代碼如下：
subroutine loclabel(fileid,label)
integer fileid
integer error
character*80 :: c80=' '
CHARACTER(LEN=*) label
rewind(fileid) !每次都從頭開始找，因為fch里內容的順序與高斯版本有關，從頭開始找比較保險。
error=0
noentryinfch=0
do while(index(c80,trim(label))==0.and.error==0) !index函數的用處是在c80中找label代表的字符串所在位置，返回0說明沒找到
read(fileid,"(a80)",iostat=error) c80 !每次讀一整行
if (error/=0) noentryinfch=1
end do
backspace(fileid)
end subroutine
下面是讀入fch文件的子程序readfch：
subroutine readfch
character*80 c80
real*8 temp
integer :: i,j,k,l,nbasis !nbasis是基函數數目
integer :: type2norb(-3:3)=(/ 7,5,4,1,3,6,10 /) !將殼層類型轉化為殼層所含軌道數，0=s->1,1=p->3,-1=sp->4,2=6d->6,-2=5d->5,3=10f->10,-3=7f->7
integer,allocatable :: shelltype(:),shellcon(:),shell2atom(:) !殼層類型、收縮度、所屬原子
integer :: s2f(-3:3,30)=0 !s2f(i,j)就是第i類殼層中第j個基函數的類型編號
real*8,allocatable :: primexp(:),concoeff(:),SPconcoeff(:) !每一殼層中高斯函數的指數、收縮系數和SP殼層的P型的收縮系數
real*8,allocatable :: amocoeff(:,:),bmocoeff(:,:) !fch文件中alpha/beta分子軌道系數矩陣
s2f(0,1)=1
s2f(-1,1:4)=(/ 1,2,3,4 /)
s2f(1,1:3)=(/ 2,3,4 /)
s2f(2,1:6)=(/ 5,6,7,8,9,10 /)
s2f(3,1:10)=(/ 11,12,13,17,14,15,18,19,16,20 /) !這是10f殼層，11~20不是按順序排列，用來起到fch中f軌道順序與wfn中順序的轉換。
mode="GAUSSIAN" !Gaussian程序里用的都是高斯函數
open(10,file=infilename,access="sequential",status="old")
read(10,*) fchtitle !fch第一行
!因為wfn只含占據軌道，故一般情況readfch只需讀入占據軌道。如果fch第一行含有saveNO，說明存的是自然軌道。若第一行是saveNBOocc或saveNBOene，說明是NBO軌道。這三種情況所有軌道占據數皆不為0，都應輸出到wfn，故應讀入所有軌道。
isaveNO=0
isaveNBOocc=0
isaveNBOene=0
if (index(fchtitle,'saveNBOocc')/=0) isaveNBOocc=1
if (index(fchtitle,'saveNBOene')/=0) isaveNBOene=1
if (index(fchtitle,'saveNO')/=0) isaveNO=1
if (isaveNBOocc==1.or.isaveNBOene==1) write(*,*) "The file contains NBO information"
if (isaveNO==1) write(*,*) "The file contains natural orbitals information"
read(10,"(a80)") c80
if (c80(11:11)=="R") wfntype=0
if (c80(11:11)=="U") wfntype=1
if (c80(11:12)=="RO") wfntype=2
call loclabel(10,'Number of electrons') !將讀入位置移動到Number of electrons這行的開頭
read(10,"(49x,i12)") ntotelec !總電子數
read(10,"(49x,i12)") nalphaelec !總alpha電子數
read(10,"(49x,i12)") nbetaelec !總beta電子數
read(10,"(49x,i12)") nbasis
call loclabel(10,'Virial Ratio')
read(10,"(49x,1PE22.15)") virialratio
call loclabel(10,'Total Energy')
totenergy=0.0D0
if (noentryinfch==0) read(10,"(49x,1PE22.15)") totenergy !使用guess=(save,only,naturalorbitals)保存自然軌道的fch文件無Total Energy項，此時不讀
call loclabel(10,'Atomic numbers')
read(10,"(49x,i12)") ncenter
allocate(a(ncenter))
read(10,"(6f12.0)") (a(i)%charge,i=1,ncenter)
do i=1,ncenter
a(i)%name=name2ind(int(a(i)%charge)) !根據原子的質子數轉換到元素名稱
end do
call loclabel(10,'Current cartesian coordinates')
read(10,*)
read(10,"(5(1PE16.8))") (a(i)%x,a(i)%y,a(i)%z,i=1,ncenter) !讀入每個原子的坐標，用隱式循環很方便
call loclabel(10,'Shell types')
read(10,"(49x,i12)") nshelltype
allocate(shelltype(nshelltype))
read(10,"(6i12)") (shelltype(i),i=1,nshelltype)
call loclabel(10,'Number of primitives per shell')
read(10,"(49x,i12)") nshellcon
allocate(shellcon(nshellcon))
read(10,"(6i12)") (shellcon(i),i=1,nshellcon)
call loclabel(10,'Shell to atom map')
read(10,"(49x,i12)") nshell2atom
allocate(shell2atom(nshell2atom))
read(10,"(6i12)") (shell2atom(i),i=1,nshell2atom)
call loclabel(10,'Primitive exponents')
read(10,"(49x,i12)") nprimexp
allocate(primexp(nprimexp))
read(10,"(5(1PE16.8))") (primexp(i),i=1,nprimexp)
call loclabel(10,'Contraction coefficients')
read(10,"(49x,i12)") nconcoeff
allocate(concoeff(nconcoeff))
read(10,"(5(1PE16.8))") (concoeff(i),i=1,nconcoeff)
read(10,"(a80)") c80
if (index(c80,"P(S=P) Contraction coefficients")/=0) then
backspace(10)
read(10,"(49x,i12)") nSPconcoeff
allocate(SPconcoeff(nSPconcoeff))
read(10,"(5(1PE16.8))") (SPconcoeff(i),i=1,nSPconcoeff)
end if
nprims=0
!計算總高斯函數數目，用于分配內存給b(:)。每個殼層的基函數數目乘以殼層的收縮度就是每個殼層的高斯函數數目
do i=1,nshelltype
nprims=nprims+type2norb(shelltype(i))*shellcon(i)
end do
allocate(b(nprims))
call loclabel(10,'Alpha Orbital Energies') !定位到Alpha Orbital Energies處
read(10,*)
if (wfntype==0.or.wfntype==2) then !閉殼層或限制性開殼層，很多代碼可以共用所以放到一起。
if (wfntype==0) nmo=ntotelec/2 !根據電子數和波函數類型判定占據軌道數
if (wfntype==2) nmo=nalphaelec !開殼層情況alpha電子>=beta電子數，故限制性開殼層取alpha電子數為占據軌道數
if (isaveNO==1.or.isaveNBOocc==1.or.isaveNBOene==1) nmo=nbasis !fch中為自然軌道或NBO軌道時讀取所有軌道
allocate(MOene(nmo))
allocate(MOocc(nmo))
allocate(MOtype(nmo))
allocate(amocoeff(nmo,nbasis)) !給alpha系數矩陣分配內存
MOtype=0
MOocc=2.0D0
if (wfntype==2) then
MOtype(nbetaelec+1:nmo)=1 !限制性開殼層時設定單占據軌道
MOocc(nbetaelec+1:nmo)=1.0D0
write(*,"('This is restricted open-shell wavefunction,',i6,' single occupied orbitals')") nalphaelec-nbetaelec
else if (wfntype==0) then
write(*,"('This is restricted wavefunction')")
end if
read(10,"(5(1PE16.8))") (MOene(i),i=1,nmo)
call loclabel(10,'Alpha MO coefficients')
read(10,*)
read(10,"(5(1PE16.8))") ((amocoeff(imo,ibasis),ibasis=1,nbasis),imo=1,nmo) !雙隱式循環讀入軌道系數
else if (wfntype==1) then !開殼層情況
nmo=ntotelec !占據軌道數等于總電子數
if (isaveNO==1.or.isaveNBOocc==1.or.isaveNBOene==1) then
nmo=2*nbasis
nalphaelec=nbasis !nalphaelec/nbetaelec控制著讀取范圍、如何設定軌道類型。含自然軌道/NBO信息時軌道全都讀入，故將alpha軌道和beta軌道數目都設成基函數數。
nbetaelec=nbasis
end if
allocate(MOene(nmo))
allocate(MOocc(nmo))
allocate(MOtype(nmo))
allocate(amocoeff(nalphaelec,nbasis))
allocate(bmocoeff(nbetaelec,nbasis))
MOocc=1.0D0
MOtype(1:nalphaelec)=1 !設為alpha軌道
MOtype(nalphaelec+1:nmo)=2 !設為beta軌道
write(*,"('This is unrestricted wavefunction,',i6,' alpha orbitals',i6,' beta orbitals')") nalphaelec,nbetaelec
read(10,"(5(1PE16.8))") (MOene(i),i=1,nalphaelec) !讀alpha軌道能量
call loclabel(10,'Beta Orbital Energies')
read(10,*)
read(10,"(5(1PE16.8))") (MOene(i),i=nalphaelec+1,nmo) !讀beta軌道能量。alpha/beta軌道能量都一起存在MOene里。
call loclabel(10,'Alpha MO coefficients')
read(10,*)
read(10,"(5(1PE16.8))") ((amocoeff(imo,ibasis),ibasis=1,nbasis),imo=1,nalphaelec)
call loclabel(10,'Beta MO coefficients')
read(10,*)
read(10,"(5(1PE16.8))") ((bmocoeff(imo,ibasis),ibasis=1,nbasis),imo=1,nbetaelec)
end if
if (isaveNBOocc==1.or.isaveNO==1) MOocc=MOene !fch開頭為saveNBOocc或saveNO時說明軌道能量部分其實是自然軌道/NBO的占據數，所以調換過來。
if (isaveNBOocc==1.or.isaveNO==1) MOene=0.0D0
if (isaveNBOene==1) MOocc=0.0D0 !fch開頭為saveNBOene時fch里沒有NBO占據數信息，故都設0。可以轉換到wfn后手動補上占據數信息。
where (MOocc==1000) !保存了NBO信息后，軌道能量/占據數最后幾個數誤為1000，都統一設為0。
MOocc=0.0D0
end where
where (MOene==1000)
MOene=0.0D0
end where
allocate(co(nmo,nprims))
!!!!!!
!現在就剩wfn的系數矩陣co(:,:)沒有填好了，這部分用于生成它，代碼及解釋見下節
!!!!!!
write(*,"('There are',i6,' MOs,',i6,' atoms,',i7,' electrons,',i7,' Gaussian functions')") nmo,ncenter,int(dnint(sum(MOocc))),nprims
write(*,"('Total energy:',f19.12,' hartree, Virial ratio:',f12.8)") totenergy,virialratio
write(*,*)
close(10)
end subroutine
5.由fch的系數矩陣生成wfn的系數矩陣
后面的代碼要用到計算高斯函數歸一化系數的函數，如下所示。函數需要輸入高斯函數類型代表的數字和指數。其中ft()代表階乘函數。
real*8 function normgau(type,exp)
real*8 exp
integer type,ix,iy,iz
!類型順序是S/X/Y/Z/ XX/YY/ZZ/XY/XZ/YZ/ XXX/YYY/ZZZ/XXY/XXZ/YYZ/XYY/XZZ/YZZ/XYZ，與wfn的類型順序一致，與fch中的順序不一致
integer :: type2ix(20)=(/ 0,1,0,0, 2,0,0,1,1,0, 3,0,0,2,2,0,1,1,0,1 /)
integer :: type2iy(20)=(/ 0,0,1,0, 0,2,0,1,0,1, 0,3,0,1,0,2,2,0,1,1 /)
integer :: type2iz(20)=(/ 0,0,0,1, 0,0,2,0,1,1, 0,0,3,0,1,1,0,2,2,1 /)
ix=type2ix(type)
iy=type2iy(type)
iz=type2iz(type)
normgau=(2*exp/pi)**0.75*sqrt( (8*exp)**(ix+iy+iz)*ft(ix)*ft(iy)*ft(iz)/(ft(2*ix)*ft(2*iy)*ft(2*iz)) )
end function
從fch的系數矩陣amocoeff和bmocoeff轉化到輸出wfn文件要用的系數矩陣CO過程略微復雜，故單獨說明。下面的圖對3-21G計算的CH4的高斯函數循環過程進行了說明，每一排球是一個殼層，每個球是一個基函數，球中每條線是一個高斯函數，箭頭指明了高斯函數循環次序。要把fch中通過層層循環描述的基函數的系數轉換到wfn中的每個高斯函數的系數，就需要先把每個高斯函數都循環一遍，循環過程中賦予每個高斯函數系數、所在中心、指數和類型。
下面的代碼可以進行這個工作。最外層循環每一殼層（桔色箭頭），在殼層中循環每個基函數（粉色箭頭），在基函數中循環每個高斯函數（綠色箭頭），在每個高斯函數中循環每個分子軌道。
k=1 !初始化累加變量
iexp=1
ibasis=1
do i=1,nshelltype !循環每一殼層
b(k:k+shellcon(i)*type2norb(shelltype(i))-1)%center=shell2atom(i) !同一殼層中高斯函數所屬中心唯一，故把這一殼層中所有高斯函數所屬中心都設好
do j=1,type2norb(shelltype(i)) !循環殼層中每個基函數，通過type2norb獲得每個殼層所屬的殼層類型含有基函數的數量
b(k:k+shellcon(i)-1)%functype=s2f(shelltype(i),j) !同一基函數中高斯函數所屬類型唯一，故把這一基函數中所有高斯函數類型設好
do l=1,shellcon(i) !循環基函數中每個高斯函數，其數目由這一殼層收縮度shellcon(i)確定
b(k)%exp=primexp(iexp+l-1) !給高斯函數賦予指數
tnormgau=normgau(b(k)%functype,b(k)%exp) !獲得高斯函數歸一化系數
temp=concoeff(iexp+l-1) !確定高斯函數的收縮系數
if (shelltype(i)==-1.and.j/=1) temp=SPconcoeff(iexp+l-1) !如果是SP殼層且是殼層中第2、3、4個基函數（即P型），則改用P(S=P) Contraction coefficients里的收縮系數
do imo=1,nmo !循環每條分子軌道
if (wfntype==0.or.wfntype==2) then !閉殼層或限制性開殼層
co(imo,k)=amocoeff(imo,ibasis)*temp*tnormgau !wfn文件中分子軌道向高斯函數的展開系數=分子軌道向基函數的展開系數*高斯函數在此基函數中的收縮系數*高斯函數的歸一化系數
else if (wfntype==1) then !開殼層情況根據當前軌道序號判斷所屬自旋類型，使用fch中不同的系數矩陣生成co(imo,k)
if (imo<=nalphaelec) co(imo,k)=amocoeff(imo,ibasis)*temp*tnormgau
if (imo>nalphaelec) co(imo,k)=bmocoeff(imo-nalphaelec,ibasis)*temp*tnormgau !imo-nalphaelec就是當前軌道在beta軌道中的序號
end if
end do
k=k+1 !當前高斯函數在全部高斯函數中的序號
end do
ibasis=ibasis+1 !當前基函數在全部基函數中的序號，以確定取amocoeff中的哪個值
end do
iexp=iexp+shellcon(i) !當前殼層第一個收縮系數在收縮系數列表concoeff中的位置
end do
將上述內容組合在一起，就成為了fch2wfn程序，完整的代碼及編譯好的程序見/usr/uploads/file/20150609/20150609184044_26809.rar，在CVF6.5、ifort編譯通過。可以運行后輸入fch文件名和輸出的wfn文件名，也可以直接用命令行模式，如fch2wfn c:\sob.fch ..\saint.wfn。再次提醒，含自然軌道的fch開頭一行應改為saveNO。含NBO的fch若用了后HF密度要把開頭改為saveNBOocc，若用了SCF密度要把開頭改為saveNBOene，以便fch2wfn正確處理。
轉換的結果和直接用高斯輸出的wfn文件幾乎完全一致，一般只有末位由于數值精度問題相差1或2。但有時會發現含自然軌道信息的fch轉化出來的wfn中有幾條軌道系數和高斯輸出的wfn不一致，但實際上都是正確的，只是自然軌道取向不同而已。
AIMALL中的模塊aimqb也能將fch向wfn的轉換，但aimqb轉換出的wfn雖然內容正確，與fch2wfn一致，但格式與標準格式不符，可能造成一些AIM程序無法讀取。而且高斯函數的儲存順序與一般規則有異（當然這完全不影響結果），比如3-21G，一般SP殼層的高斯函數的類型儲存順序是X X Y Y Z Z | X Y Z，其中|代表分裂價基的分隔，而aimqb轉換出來的順序是X Y Z X Y Z X Y Z。fch2wfn的轉換不僅結果正確，也完整保留了wfn文件的習俗。fch2wfn的功能也已經整合進了Multiwfn中，在Multiwfn中讀入fch后選功能6，再選保存wfn文件即可。