從Gaussian和GAMESS-US中提取電子積分
獲得電子積分是自己寫量化程序的關鍵。在網上能找到一些零散的從常用的Gaussian或GAMESS-US程序中提取單、雙電子積分的帖子,但是往往語焉不詳,說得不確切、不全面,有的做法不雅,有的還會令人白繞彎路。這里把從Gaussian和GAMESS-US中提取電子積分的方式清楚、詳細、準確地說一下。
無論輸出的是單電子積分還是雙電子積分,都是基函數已歸一化后的值。
雙電子積分形式是(IJ|KL),I,J基函數是r1坐標,K,L基函數是r2坐標。由于積分有對稱性,即IJKL=IJLK=JIKL=JILK=KLIJ=LKIJ=KLJI=LKJI,所以程序輸出雙電子積分的時候都只輸出對稱唯一的部分。
1 從Gaussian中提取積分
以下內容對Gaussian09 D.01親測有效。由于Gaussian會自動調整朝向,會影響積分值,建議加上nosymm,這樣也避免因對稱性等價的積分不被輸出。
1.1 單電子積分
計算時加上IOp(3/33=1)即可,輸出的包括*** Overlap ***:重疊積分
*** Kinetic Energy ***:動能積分
***** Potential Energy *****:核吸引勢積分(沒考慮電子帶負電荷,所以數值皆為正)
****** Core Hamiltonian ******:核哈密頓矩陣。其數值就是動能積分矩陣元減去核吸引勢積分矩陣元
Multipole matrices:輸出三次,IX= 1、IX= 2、IX= 3分別對應X,Y,Z方向的偶極矩積分矩陣。
1.2 雙電子積分
Gaussian可以通過L316模塊雙電子積分。用這些關鍵詞即可:IOp(3/33=3) extralinks(L316) scf=conventional noraff。同時上一節的單電子積分也會輸出。輸出信息諸如:
IntCnt= 4238596 ITotal= 4238596 NWIIB= 262144 ISym2E=0
I= 44 J= 29 K= 43 L= 10 Int= 0.147872132379D-07
I= 44 J= 43 K= 29 L= 10 Int= -0.396977105787D-08
I= 45 J= 29 K= 43 L= 10 Int= -0.403826333553D-07
...
這里IntCnt是實際輸出的積分總數。注意并行計算時每次輸出的積分順序可能不同。
Gaussian計算并輸出雙電子積分的值閾值可以用IOp(3/27=N)調節,小于10^-N的積分不會被計算和輸出。默認N=10,精度已經很高了,若要計算并輸出數值更小的可以比如設N=12。
用32bit版本時這種方式輸出積分時基函數最多127個。雖然64bit版沒這個限制,但基函數多于這個時輸出文件會巨大。
2 從GAMESS-US中提取積分
以下方法對2014-Dec版親測有效。相同基組下(Dunning基組不算,因Gaussian會自動會去掉其中冗余的GTF)輸出的結果和Gaussian中的一致。
輸出單電子積分和雙電子積分都有改代碼和不改代碼兩種做法,結果一樣,前者省事,但后者可以自由控制輸出格式,而且還能同時用NPRINT=-5來避免計算時輸出一堆翔一樣的大量無用的信息。
輸出積分時絕對不要用并行運算,否則輸出內容會癲狂。
2.1 單電子積分
2.1.1 方法1(不用改源碼)
$CONTRL里寫上NPRINT=3即可在計算時輸出基函數間的各種積分矩陣,會看到這樣的段落********************
1 ELECTRON INTEGRALS
********************
下面輸出的包括:
OVERLAP MATRIX:重疊矩陣
BARE NUCLEUS HAMILTONIAN INTEGRALS (H=T+V):核哈密頓矩陣
KINETIC ENERGY INTEGRALS:動能積分矩陣
程序沒直接給出核吸引勢積分,這只需自行把核哈密頓矩陣減去動能積分矩陣即可。
2.1.2 方法2(需改源碼)
在gamess/source/int1.src中----- OPTIONAL DEBUG PRINTOUT -----有兩處,在第一處的下面加上write(*,"(/,' Overlap integrals, num:',i8)") LL2
write(*,"(4D18.10)") S
write(*,"(/,' Kinetic energy integrals, num:',i8)") LL2
write(*,"(4D18.10)") T
write(*,"(/,' Core Hamiltonian matrix, num:',i8)") LL2
write(*,"(4D18.10)") H
write(*,"(/,' Potential energy integrals, num:',i8)") LL2
write(*,"(4D18.10)") H-T
write(*,*)
可見,改過之后可以直接把核吸引勢積分輸出出來,已經考慮電子是帶負電荷。
重新編譯此文件并鏈接成新的可執行文件使之生效,即在GAMESS-US目錄下執行./comp int1;./lked gamess 00。瞬間就完事,用不著把整個GAMESS-US都重頭編譯一遍。
之后一般計算時屏幕上就會輸出各種單電子積分,上述幾種積分矩陣都是對稱矩陣,這里依次輸出的是下三角部分的元素。如重疊積分:
Overlap integrals, num: 28
0.1000000000D+01 0.2367039205D+00 0.1000000000D+01 0.0000000000D+00
0.0000000000D+00 0.1000000000D+01 0.0000000000D+00 0.0000000000D+00
0.0000000000D+00 0.1000000000D+01 0.2005800301D-16 0.2537204019D-17
0.0000000000D+00 0.0000000000D+00 0.1000000000D+01 0.5362008406D-01
0.4729709355D+00 0.0000000000D+00 -0.3205127531D+00 0.2265160217D+00
0.1000000000D+01 0.5362008406D-01 0.4729709355D+00 0.0000000000D+00
0.3205127531D+00 0.2265160217D+00 0.2327636225D+00 0.1000000000D+01
這里num后面是輸出的元素數目。
2.2 雙電子積分
一定要用$SCF DIRSCF=.F. $END使用conventional的SCF方式,否則不會輸出這些積分。2.2.1 方法1(不用改源碼)
$CONTRL里寫上NPRINT=4即可在計算時輸出雙電子積分,會看到如II,JST,KST,LST = 1 1 1 1 NREC = 1 INTLOC = 1
II,JST,KST,LST = 2 1 1 1 NREC = 1 INTLOC = 2
1 1 1 1 1.0 4.785065752 2 1 1 1 1.0 0.741380321
2 2 1 1 1.0 1.118946840 3 3 1 1 1.0 1.115813819
4 4 1 1 1.0 1.115813819 5 5 1 1 1.0 1.115813819
2 1 2 1 1.0 0.136873367 3 1 3 1 1.0 0.024477411
...
TOTAL NUMBER OF NONZERO TWO-ELECTRON INTEGRALS = 228
最后一行就是實際輸出的雙電子積分數,當前是228個。默認是小于10^-9的積分不被輸出,可以在$CONTRL里用ICUT=N來把輸出閾值調為10^-N,顯然N越大輸出得越多。
對(IJ|KL)積分,有8種等價形式,程序自動把指標調換成I>=J、K>=L、I>=K的形式。
2.2.2 方法2(需改源碼)
在gamess/source/int2a.src的I4 = LOCL+L下面加上write(*,"(4i6,D25.16)") I1,I2,I3,I4,VAL
重新編譯此文件并鏈接成新的可執行文件使之生效,即在GAMESS-US目錄下執行./comp int2a;./lked gamess 00。
之后一般計算時屏幕上就會輸出雙電子積分,諸如
1 1 1 1 0.4785065752035100D+01
2 1 1 1 0.7413803207944081D+00
2 2 1 1 0.1118946840438162D+01
3 3 1 1 0.1115813818563188D+01
4 4 1 1 0.1115813818563188D+01
5 5 1 1 0.1115813818563188D+01
...
如上一節所述,非常小的積分不會輸出,閾值可以用ICUT來設。如果想輸出所有積分,直接把int2a.src中IF(ABS(VAL).LT.CUTOFF) GO TO 200這一行注釋掉即可。
如果想輸出NPRINT=4那樣的將序號調換為I>=J、K>=L、I>=K后的積分,把前述輸出語句加到IF (OUT) CALL INTOUT(I1,I2,I3,I4,QQ4,IJKL_INDEX,VAL)的上面一行。