• 從Gaussian和GAMESS-US中提取電子積分

    從Gaussian和GAMESS-US中提取電子積分

    文/Sobereva @北京科音  2016-Jul-31


    獲得電子積分是自己寫量化程序的關鍵。在網上能找到一些零散的從常用的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)的上面一行。
    久久精品国产99久久香蕉