DFT交換相關泛函庫的使用方法

DFT交換相關泛函庫的使用方法

文/Sobereva @北京科音
First release: 2013-Nov-10    Last update: 2019-Feb-10

1 前言

DFT大行其道的今天，寫量子化學程序總是免不了要涉及DFT計算，而DFT的核心就是交換相關泛函。現有的交換相關泛函數目甚巨，每年都有新的泛函被提出。交換相關泛函形式往往比較復雜，要想求其導數更是件麻煩事（計算交換相關勢時要用到），自行手寫不僅費時間還容易出錯。想在自己的程序里加入對各種各樣的泛函的支持，最方便的方法就是直接用現成的交換相關泛函庫。本文就介紹五個泛函庫，并對其中兩個的使用進行著重介紹，給出示例代碼，但對于零碎的細節并不會一覽無余的介紹，用戶應在閱讀本文后去閱讀手冊和相關網頁。這些庫都是免費的。

注意以下介紹的五個庫對于雜化泛函都只能產生純泛函的部分。也就是說，要想做雜化泛函計算，HF交換項部分得自己想辦法計算，這要涉及到雙電子積分，這些庫是沒法給出的。雙電子積分是量化中最復雜的積分，圖省事的話可以考慮用調用免費、現成的電子積分庫libint來實現（http://sourceforge.net/p/libint/home），ORCA、PSI、CP2K也都用的是這個庫。至于雙雜化泛函，那還得自己整MP2的代碼。

2 Density functional repository庫

用這個庫之前建議閱讀一下CPL,199,557，這篇文章涉及到求解kohn-Sham方程的一些細節，特別是KS矩陣的交換相關部分。閱讀這篇文章有助于了解這個庫的設計思想。

這個庫的網址是ftp://ftp.dl.ac.uk/qcg/dft_library/index.html。這是個比較老的庫，到2006年就不更新了。支持的泛函不算很多，是LDA/GGA/hybrid_GGA級別的，meta-GGA的都沒有。這個庫完全使用Fortran77寫的，一個優點是每個泛函的代碼都是獨立的。比如說自己的程序只想用PW91，那么就把PW91的代碼從相應頁面里拷貝出來插到自己的程序里并調用即可。

這個庫的網址里還給出了每種泛函的定義以及計算原子時的參考數據，可以用來檢驗自己的代碼算得對不對。實際上沒必要用它的參考數據，直接和Gaussian的結果對比一下就行了，詳見本文附錄A。

這個庫的子程序命名規則比較明確，例如uks_x_pbe，x就代表這個子程序用來計算開殼層形式的PBE交換泛函。計算相關部分就是c，如果是交換相關合在一起就是xc。閉殼層形式的子程序和開殼層是不同的，把uks換為rks即可。例如rks_xc_b97就是用來計算閉殼層形式下B97交換相關泛函。為了免得大家自己一個個去從網頁上拷貝代碼，筆者直接把這個庫里常用的泛函的代碼都合并在了同一文件DFTxclib.F里，可以從這里下載：/usr/uploads/file/20150609/20150609211428_67937.rar。其中包括這些子程序
uks_x_lda
uks_x_b88
uks_x_pbe
uks_x_pw91
uks_c_vwn5
uks_c_p86
uks_c_lyp
uks_c_pw91
uks_c_pbe
uks_xc_b97
uks_xc_hcth407
相應的rks的版本也在DFTxclib.F這個文件里。

值得一提的是這個庫的代碼并不是手寫的，而是通過Knowles等人的dfauto程序通過計算機自動生成的。因此查看這個庫里的代碼，會發現很不符合人類的邏輯，代碼很長，凈是一大堆數值和搞不懂的變量名。之所以作者不手寫而靠計算機自動生成，原因是這個庫不僅產生泛函的數值，還產生泛函對電子密度和對sigma（見后文）的導數，以及混合導數，對于開殼層又分為兩種自旋，情況會更為復雜。如果對各種泛函各個導數都一個一個地手推公式會非常困難，不僅很難保證不出錯，代碼寫起來也極其勞神。關于dfauto和Maple的介紹見本文附錄B。

這個庫的子程序的調用方式都是標準、統一的，可參見ftp://ftp.dl.ac.uk/qcg/dft_library/design.html#calling_convention里的說明。對于閉殼層形式，需要提供的信息有
rhoa1：總密度
sigmaaa1：總密度梯度模的平方，稱為σ，也即|▽ρ|^2=▽ρ·▽ρ
返回的信息包括
zk：泛函的數值。設交換相關能寫為E_xc=∫f(ρ(r),σ(r))dr，zk就是指f的值
vrhoa：泛函對總密度的一階導數，即df/dρ（d代表偏導，后同）
vsigmaaa：泛函對σ的一階導數，即df/dσ
v2rhoa2：泛函對總密度的二階導數
v2rhoasigmaaa：泛函同時對總密度的一階導數和對σ的一階導數
v2sigmaaa2：泛函對σ的二階導數

下面是個這個庫里的典型的子程序的開頭部分
      subroutine rks_x_pbe
     & (ideriv,npt,rhoa1,sigmaaa1,
     &  zk,vrhoa,vsigmaaa,
     &  v2rhoa2,v2rhoasigmaaa,v2sigmaaa2)
...略
      implicit real*8 (a-h,o-z)
      integer ideriv,npt
      real*8 rhoa1(npt)
      real*8 sigmaaa1(npt)
      real*8 zk(npt),vrhoa(npt),vsigmaaa(npt)
      real*8 v2rhoa2(npt),v2rhoasigmaaa(npt),v2sigmaaa2(npt)
...略
ideriv=0時代表只計算泛函數值zk；=1代表把一階導數vrhoa、vsigmaaa也計算出來；=2代表把二階導數v2rhoa2、v2rhoasigmaaa、v2sigmaaa2也計算出來。
npt代表輸入多少個點的信息。從子程序定義可見，輸入、輸出的信息并不是標量，而是長度為npt的數組。也就是說，比如傳進來5個點的密度值和σ值，程序就會計算并返回這5個點的泛函值和導數。之所以子程序設計成一次計算一批點的形式，是因為實際DFT計算時要用到格點積分方法，需要算數量很多的點，調用一次子程序就算完N個點比起調用N次每次只算一個點的效率明顯要高。不過當然也可以讓npt=1，一次就只計算一個點了。

特別要注意的是，這個庫里的閉殼層形式的子程序算出來的導數不對。想得到正確結果，vsigmaaa和v2rhoasigmaaa都必須除以4，v2rhoa2必須除以2，v2sigmaaa2必須除以16。大家可以通過有限差分方法獲得zk對密度和對σ的導數來驗證這一點。

下面是一個例子，用戶每次輸入密度和σ，然后返回泛函值和導數
program xctest
implicit real*8 (a-h,o-z)
real*8 rho(1),sigma(1),d1rho(1),d1sig(1),d2rho(1),d2rhosig(1),d2sig(1)
do while(.true.)
 write(*,*) "Select functional:"
 write(*,*) "0 LDA exchange functional"
 write(*,*) "1 Becke 88 exchange functional"
 write(*,*) "2 LYP correlation functional"
 write(*,*) "3 HCTH407 exchange-correlation functional"
 read(*,*) isel
 write(*,*) "Input rho"
 read(*,*) rho
 write(*,*) "Input sigma"
 read(*,*) sigma
 if (isel==0) then
  call rks_x_lda(2,1,rho,sigma,value,d1rho,d1sig,d2rho,d2rhosig,d2sig)
 else if (isel==1) then
  call rks_x_b88(2,1,rho,sigma,value,d1rho,d1sig,d2rho,d2rhosig,d2sig)
 else if (isel==2) then
  call rks_c_lyp(2,1,rho,sigma,value,d1rho,d1sig,d2rho,d2rhosig,d2sig)
 else if (isel==3) then
  call rks_xc_hcth407(2,1,rho,sigma,value,d1rho,d1sig,d2rho,d2rhosig,d2sig)
 end if
 write(*,"(' Value=',f17.12)") value
 write(*,"(' d1rho=',f17.12,'  d1sig=',f17.12)") d1rho,d1sig/4D0
 write(*,"(' d2rho=',f17.12,'  d2rhosig=',f17.12,'  d2sig=',f17.12,/)") d2rho/2D0,d2rhosig/4D0,d2sig/16D0
end do
end program
對于LDA級別的泛函，由于不依賴于密度的梯度，所以輸入σ的時候隨便提供一個數就行了。輸出的信息中也只有value、d1rho、d2rho有意義。
像是B97和HCTH407都是交換和相關部分結合在一起的。如果想計算BLYP這樣交換和相關部分彼此獨立的泛函，就要分別調用rks_x_b88和rks_c_LYP并把結果加在一起。

有興趣的話可以寫一個基于.wfn或.fch文件計算整個體系的交換相關能的小程序，并不困難。筆者在《利用wfn文件計算電子密度的代碼的編寫方法》（http://www.shanxitv.org/182）里面已經詳細介紹了怎么讀取.wfn文件獲得波函數信息，并由此計算指定點的電子密度和梯度。在《密度泛函計算中的格點積分方法》（http://www.shanxitv.org/69）中筆者介紹了如何通過Becke提出的多中心格點積分方法對特定函數在全空間進行積分。因此，把這兩個帖子和本文內容結合在一起就可以寫出代碼了。也就是先載入波函數信息，在格點積分過程中對每一個點的位置計算電子密度和梯度，并調用泛函庫的子程序求得這個點的泛函值，然后乘上積分權重進行累加，最終積分值便是交換相關能。

交換相關勢就是交換相關能對密度的變分，即Vxc=δExc/δρ，計算它是求解KS方程中重要的一環。KS矩陣元中交換相關部分計算公式為

其中χ代表基函數。這種矩陣元也是通過格點積分方式來計算的，其中要用到的f對ρ和對σ的導數，這在前面的例子中已經通過調用庫里的子程序得到了，因此矩陣元可以比較容易地計算出來。

計算交換相關勢在某個點的數值的計算公式為

▽▽ρ就是電子密度的Hessian矩陣。可見計算某個點的交換相關勢需要用到泛函的二階導數（庫里的子程序也都輸出了）以及密度的二階導數。Multiwfn程序（3.2.1版及以后。http://www.shanxitv.org/multiwfn）的function.f90中DFTxcpot_close函數就是用來根據這個式子計算交換相關勢的，有興趣者可參閱相應代碼。

上面只討論了閉殼層形式。對于開殼層情況，alpha和beta部分不同，故要輸入和返回的信息都明顯增加，例如這是uks_x_pbe開頭的定義
      subroutine uks_x_pbe
     & (ideriv,npt,rhoa1,rhob1,sigmaaa1,sigmabb1,sigmaab1,
     &  zk,vrhoa,vrhob,vsigmaaa,vsigmabb,vsigmaab,
     &  v2rhoa2,v2rhob2,v2rhoab,
     &  v2rhoasigmaaa,v2rhoasigmaab,v2rhoasigmabb,
     &  v2rhobsigmabb,v2rhobsigmaab,v2rhobsigmaaa,
     &  v2sigmaaa2,v2sigmaaaab,v2sigmaaabb,
     &  v2sigmaab2,v2sigmaabbb,v2sigmabb2)
要輸入的包括alpha密度ρ_α（rhoa1）、beta密度ρ_β（rhob1）、▽ρ_α·▽ρ_α（sigmaaa1）、▽ρ_β·▽ρ_β（sigmabb1）、▽ρ_α·▽ρ_β（sigmaab1）。由于輸入的變量就有5個，所以輸出的一階導數、二階導數和混合導數數量巨多，這里就不一一說明了，每一項的含義見ftp://ftp.dl.ac.uk/qcg/dft_library/design.html#calling_convention。開殼層形式的子程序產生的導數沒有bug，可以直接用，不需要像閉殼層的情況那樣還得除以特定的因子來修正。

3 Libxc庫

Libxc庫的網址是http://www.tddft.org/programs/octopus/wiki/index.php/Libxc，介紹這個庫的文章見Comput. Phys. Commun. 183, 2272 (2012)。

這個庫是目前最為全面、最強大的泛函庫，一直在不斷維護中，力圖將所有已被提出的泛函全都支持，包括近幾年大量提出的meta級別的泛函（依賴于動能密度或電子密度的拉普拉斯值），另外還支持大量動能泛函。目前這個庫支持的泛函約200種。這個庫已經被很多程序所使用，很多都是第一性原理的程序，包括Octopus，Abinit、CP2K、BigDFT、GAPW、Atomistix ToolKit、AtomPAW等。實際上這個庫的前身是Octopus程序里面的計算泛函的部分，后來才獨立出來成為單獨的庫。

這個庫是由C語言編寫的，但也帶有Fortran的接口。雖然這個庫支持的泛函數量極多，而且也能產生各種一階、二階導數，但是代碼并不是機器自動生成的，而是完全靠手寫的，主要是開發者嫌自動生成的代碼太冗長、效率低、不易讀。

這個庫有一點我個人不太喜歡的地方是數據結構比較復雜，代碼間相互依賴，沒法像Density functional repository那樣想用哪個泛函就把哪個泛函的代碼拷出來放到自己的程序里直接調用就行，而必須先得把這個庫的代碼編譯成庫文件然后才能在自己的程序中調用，而且使用時還得做初始化之類的麻煩事。

Libxc的編譯方式是下載后并解壓，進入目錄后運行
autoreconf -i
./configure FC=ifort
make
make install
然后.lib庫文件和.h頭文件，以及Fortran用的.mod文件就會被安裝到/opt/etsf下面的lib和include目錄中。./configure這一步也可以自己指定安裝目錄。

下面是個Fortran90調用libxc的例子，和上一節的例子的用途和結果完全一致，即自行輸入密度和σ，選擇泛函，然后程序給出泛函的數值和導數。比如此文件叫做xctset.f90，編譯的時候就用ifort xctest.f90 -I/opt/etsf/include /opt/etsf/lib/libxc.a -o xctest
program xctest
use xc_f90_types_m       //必須寫
use xc_f90_lib_m       //必須寫
implicit real*8 (a-h,o-z)
TYPE(xc_f90_pointer_t) :: xc_func,xc_info    //必須寫，定義libxc用到的指針

do while(.true.)
   write(*,*) "Input rho"
   read(*,*) rho
   write(*,*) "Input sigma"
   read(*,*) sigma
   write(*,*) "Select functional:"
   write(*,*) "0 LDA exchange functional"
   write(*,*) "1 Becke 88 exchange functional"
   write(*,*) "2 LYP correlation functional"
   write(*,*) "3 HCTH407 exchange-correlation functional"
   read(*,*) isel
   if (isel==0) ifunc=XC_LDA_X
   if (isel==1) ifunc=XC_GGA_X_B88
   if (isel==2) ifunc=XC_GGA_C_LYP
   if (isel==3) ifunc=XC_GGA_XC_HCTH_407

   call xc_f90_func_init(xc_func,xc_info,ifunc,XC_UNPOLARIZED)  //初始化

   select case (xc_f90_info_family(xc_info))
   case(XC_FAMILY_LDA)
     call xc_f90_lda(xc_func,1,rho,value,d1rho,d2rho,d3rho)
   case(XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
!      call xc_f90_gga_exc(xc_func,1,rho,sigma,value)
!      call xc_f90_gga_exc_vxc(xc_func,1,rho,sigma,value,d1rho,d1sig)
      call xc_f90_gga(xc_func,1,rho,sigma,value,d1rho,d1sig,d2rho,d2rhosig,d2sig)
   end select

   write(*,"(' value=',f17.12)") value*rho
   write(*,"(' d1rho=',f17.12,'  d1sig=',f17.12)") d1rho,d1sig
   write(*,"(' d2rho=',f17.12,'  d2rhosig=',f17.12,'  d2sig=',f17.12)") d2rho,d2rhosig,d2sig
   if (ifunc==XC_LDA_X) write(*,"(' d3rho=',f17.12)") d3rho
   call xc_f90_func_end(xc_func)  //清理
   write(*,*)
end do
end program xctest

例子中諸如XC_LDA_X、XC_GGA_C_LYP這都是在module libxc_funcs_m里面定義的整型變量，打開libxc的src子目錄下的libxc_funcs.f90（或xcfuncs.h）就能看到定義，例如
  integer, parameter :: XC_GGA_XC_B97_GGA1   =  96  !  Becke 97 GGA-1                          
  integer, parameter :: XC_GGA_XC_HCTH_A     =  97  !  HCTH-A                                  
  integer, parameter :: XC_GGA_X_BPCCAC      =  98  !  BPCCAC (GRAC for the energy)
  integer, parameter :: XC_GGA_C_REVTCA      =  99  !  Tognetti, Cortona, Adamo (revised)
  integer, parameter :: XC_GGA_C_TCA         = 100  !  Tognetti, Cortona, Adamo
  integer, parameter :: XC_GGA_X_PBE         = 101  !  Perdew, Burke & Ernzerhof exchange   
因此例子中ifunc的值就是所選取的泛函在libxc當中的內部編號。

xc_f90_func_init子程序依據ifunc的值和選擇的自旋形式，返回相應的xc_func和xc_info。例子中調用xc_f90_func_init時傳入XC_UNPOLARIZED說明接下來計算泛函時是閉殼層形式的，如果用開殼層形式計算，就傳入XC_POLARIZED。XC_UNPOLARIZED和XC_POLARIZED實際上是在module libxc_funcs_m里面定義的整型變量，對應的數值分別是1和2。

libxc在調用泛函的時候是根據類別來進行的，而不是由用戶直接去調用相應的泛函的子程序。利用xc_info我們可以得到所選的泛函的一些信息，例子中將xc_info傳入xc_f90_info_family函數返回的就是我們所選的泛函所屬的類別編號，然后通過分支語句去調用相應類別的子程序。例子中XC_FAMILY_LDA、XC_FAMILY_GGA、XC_FAMILY_HYB_GGA都是module libxc_funcs_m里面定義的整型變量，數值分別為1、2、32。

如果所選的泛函是LDA類型的，就通過
call xc_f90_lda(xc_func,1,rho,value,d1rho,d2rho,d3rho)
形式調用，若是GGA或者雜化GGA級別，就通過
call xc_f90_gga(xc_func,1,rho,sigma,value,d1rho,d1sig,d2rho,d2rhosig,d2sig)
形式調用。可見不同類別的泛函的輸入輸出參數是不同的，但是相同類別的泛函的輸入輸出參數是一致的。調用對應類別的子程序，再給它們傳入xc_func，就確切指定了到底要計算的是哪一個泛函的值。

libxc調用泛函時所需要輸入的信息以及輸出的信息和Density functional repository里的子程序一樣。諸如GGA級別，就是輸入密度、σ，然后返回泛函對密度、對σ以及對二者的混合導數。但對于LDA，libxc庫還額外可以計算對密度的三階導數。若計算時用的是開殼層形式，則也需要把alpha和beta的密度信息分別傳入，詳見libxc網頁上的manual里的解釋。libxc也可以同時計算很多點的值，但是上例每次只計算一個點，所以調用xc_f90_gga的時候第二個參數為1。

對于每類泛函其實libxc提供了不同形式的調用。例子中我們用xc_f90_gga算出了泛函值和一階、二階導數。如果只想獲得泛函值，改用call xc_f90_gga_exc即可，顯然比xc_f90_gga更省時。如果只需要泛函值和一階導數，調用xc_f90_gga_exc_vxc就夠了。

計算最后，要通過call xc_f90_func_end(xc_func)來進行清理工作。

特別注意的是libxc定義交換相關能的形式為E_xc=∫ρ*ε(ρ(r),σ(r))dr，給出的泛函值是ε項而非積分號里的整體，因此還需要再乘以ρ才和Density functional repository給出的泛函值一致。但libxc對泛函的導數的定義和Density functional repository則完全一樣，數值能直接對應。

和上一節的例子相比，可見libxc在使用的時候還是略微繁瑣的，架構有點復雜，筆者還是更喜歡Density functional repository的子程序的簡單直接的調用方式。

4 MFM庫（Minnesota Functional Module）

MFM庫的網址為http://comp.chem.umn.edu/mfm/index.html，只支持明尼蘇達系列泛函，隨著這個泛函系列的不斷發展這個庫也在時不時更新，目前支持M05, M05-2X, M06-L, M06-HF, M06, M06-2X, M08-HX, M08-SO, M11, M11-L, MN12-L, MN12-SX, SOGGA, SOGGA11, SOGGA11-X, N12, N12-SX。

MFM庫完由Fortran77編寫，是手寫的。輸入的信息是alpha和beta電子的密度，它們的密度梯度和動能密度。這個庫只能計算泛函值和對密度、σ和動能密度的一階導數。

壓縮包里面每種泛函對應一個.F文件，像Density functional repository庫一樣，每個泛函代碼之間彼此獨立，調用方便。由于這個庫用處不大所以不再舉例。

5 funclib庫

funclib庫下載地址為http://www.theochem.kth.se/~pawsa/，在JCC,28,2569(2007)有這個庫的介紹。這個庫后來沒怎么更新，沒有meta級別泛函。由于funclib支持的泛函遠不及libxc全面，且同為C語言編寫，所以funclib沒什么優勢和使用價值，這里不舉例了。

funclib的代碼和Density functional repository一樣是通過計算機自動推導、生成的，但是用的是比maple歷史更為悠久的著名計算機代數系統Maxima（http://maxima.sourceforge.net）。雖然Maxima的界面遠不及maple，不過最大的好處是完全免費。

6 XCFun庫

XCFun庫的地址為http://dftlibs.org/xcfun/，由Ulf Ekstrom開發，Dalton、Dirac用的就是這個庫。支持的泛函不多，主要優點是可以計算無限階導數，另外也直接提供了計算LDA/GGA級別的交換相關勢的函數。這個庫是C++開發的，但是也提供了Fortran接口。這個庫的介紹文檔及其簡陋，調用方式也不是一個泛函對應一個子程序，而和funclib稍微有點類似。

7 自寫泛函代碼

實際上自己手寫泛函代碼也不太困難，下面給出筆者自己寫的計算指定點的B88和LYP泛函值的函數。其中gendensgradab是Multiwfn程序中用來計算給定位置的ρ_α、ρ_β、總密度ρ、|▽ρ_α|、|▽ρ_β|、|▽ρ|和▽ρ_α·▽ρ_β的子程序。
!!!----- Integrand of Becke88 exchange functional
real*8 function xBecke88(x,y,z)
implicit real*8 (a-h,o-z)
real*8 x,y,z
call gendensgradab(x,y,z,adens,bdens,tdens,agrad,bgrad,tgrad,abgrad)
adens4d3=adens**(4D0/3D0)
bdens4d3=bdens**(4D0/3D0)
slatercoeff=-3D0/2D0*(3D0/4D0/pi)**(1D0/3D0)
slaterxa=slatercoeff*adens4d3
slaterxb=slatercoeff*bdens4d3
slaterx=slaterxa+slaterxb
redagrad=agrad/adens4d3
redbgrad=bgrad/bdens4d3
arshredagrad=log(redagrad+dsqrt(redagrad**2+1))
Beckexa=adens4d3*redagrad**2/(1+6*0.0042D0*redagrad*arshredagrad)
arshredbgrad=log(redbgrad+dsqrt(redbgrad**2+1))
Beckexb=bdens4d3*redbgrad**2/(1+6*0.0042D0*redbgrad*arshredbgrad)
Beckex=-0.0042D0*(Beckexa+Beckexb)
xBecke88=slaterx+Beckex
end function

!!!----- Integrand of LYP corelation functional
real*8 function cLYP(x,y,z)
implicit real*8 (a-h,o-z)
real*8 x,y,z
call gendensgradab(x,y,z,adens,bdens,tdens,agrad,bgrad,tgrad,abgrad)
parma=0.04918D0
parmb=0.132D0
parmc=0.2533D0
parmd=0.349D0
densn1d3=tdens**(-1D0/3D0)
parmw=exp(-parmc*densn1d3) / (1+parmd*densn1d3) * tdens**(-11D0/3D0)
parmdel=parmc*densn1d3+parmd*densn1d3/(1+parmd*densn1d3)
parmCf=3D0/10D0*(3*pi*pi)**(2D0/3D0)
tmp1=-parma*4D0/(1+parmd*densn1d3)*adens*bdens/tdens
tmp2a=2**(11D0/3D0)*parmCf*(adens**(8D0/3D0)+bdens**(8D0/3D0))
tmp2b=(47D0/18D0-7D0/18D0*parmdel)*tgrad**2
tmp2c=-(2.5D0-parmdel/18D0)*(agrad**2+bgrad**2)
tmp2d=-(parmdel-11D0)/9D0*(adens/tdens*agrad**2+bdens/tdens*bgrad**2)
tmp2=adens*bdens*(tmp2a+tmp2b+tmp2c+tmp2d)
tmp3=-2D0/3D0*tdens**2*tgrad**2+(2D0/3D0*tdens**2-adens**2)*bgrad**2+(2D0/3D0*tdens**2-bdens**2)*agrad**2
cLYP=tmp1-parma*parmb*parmw*(tmp2+tmp3)
end function

附錄 A：

在Gaussian中，使用IOp(5/33=1)后程序就會在每一次迭代中輸出DFT交換能和相關能，例如B3LYP下計算水的最后一輪迭代的輸出：
One-electron energy=-0.122935559078D+03（電子動能+核吸引勢能）
 1/2 <PG(P)>=          44.994396247255（由雙電子積分模塊得到的電子間作用能。即電子間經典庫侖互斥能加上HF交換能）
 Calling FoFDFT with Acc= 0.10D-09 ICntrl=  400011 IRadAn=       4.
 Ex=          -7.116238484455 Ec=          -0.438937943288（DFT純泛函的交換能和相關能。對于雜化泛函Ex不包含HF交換能）
 1/2 <PG(P)> + E(ex) + E(corr)=          37.439219819512（即電子間總相互作用能）
 E= -76.4089533399206     Delta-E=       -0.000000000005 Rises=F Damp=F
總能量-76.4089533399206即是-0.122935559078D+03與37.439219819512之和再加上當前體系的核互斥能。

附錄 B：

介紹dfauto程序的文章是Computer Physics Communications,136,310(2001)，實際上就是個bash腳本。用戶通過Maple程序的語言格式寫好泛函的定義，dfauto在讀取后就可以自動調用當前系統中安裝的Maple程序產生計算泛函數值以及各種導數的Fortran代碼，導數的推導完全通過Maple的符號運算功能自動地實現。dfauto會控制Maple產生Fortran代碼的過程并添油加醋，以使得代碼的輸入輸出接口變得標準化。不過遺憾的是dfauto御用的是非常老的Maple 5，筆者嘗試在Maple16上使用卻通不過。

Maple的符號運算功能很強大而且能產生Fortran代碼，在量化上很有用，這里就以LDA交換泛函作為例子簡單說明一下。這個泛函形式極簡單，就是常數乘以密度的4/3次方。在Maple中首先定義此泛函，即輸入
LDAx:=-cf*rho^(4/3):
然后輸入
with(CodeGeneration):
Fortran([value=LDAx,vrho=diff(LDAx,rho),v2rho=diff(LDAx,rho$2)]);
程序就會輸出計算此泛函數值以及計算它對密度的一階和二階導數的Fortran代碼：
value = -cf * rho ** (0.4D1 / 0.3D1)
vrho = -0.4D1 / 0.3D1 * cf * rho ** (0.1D1 / 0.3D1)
v2rho = -0.4D1 / 0.9D1 * cf * rho ** (-0.2D1 / 0.3D1)
其中vrho正是LDA交換勢。