自寫Link生成Gaussian的IRC任務中每個點的波函數文件

1 前言

在Gaussian的IRC任務，以及scan、優化等任務中，只有最后一步波函數才允許在任務結束時被l9999輸出為PROAIM波函數文件。如果想要用外部工具，比如Multiwfn繪制這個過程中實空間函數（電子密度、ELF等）的動態變化過程，就必須獲得每一步的波函數文件。為實現這個目的，最笨的方法是寫個腳本，提取每一幀結構再做單點計算，顯然這要多費很多時間。另一個辦法是自行編寫、修改link，然后篡改link的執行順序，使波函數文件在每一步中都輸出，這就避免了重算波函數。本文將首先介紹Gaussian的link的基本代碼結構、編譯方法，然后再談怎么自寫link來實現使IRC任務的每一步都能直接輸出波函數文件。

建議閱讀本文前先閱讀《Gaussian的Link、IOp與非標準計算路徑》一文，見http://www.shanxitv.org/57。并且下載《Gaussian 09 Programmer's Reference》以備用，最好把前58頁讀一遍。本文會直接提供編譯好的link，如果讀者想自行編譯，必須擁有Gaussian09源代碼，最好是A02版本，并且自己已經成功編譯過一遍，編譯方法可參考《Gaussian09 A02、C01 64bit編譯方法》，見http://www.shanxitv.org/81，本文假設$g09root環境變量被設為/sob，用戶是root。如果對波函數文件不熟悉，可以閱讀《高斯fch文件與wfn波函數文件的介紹及轉換方法》http://www.shanxitv.org/55。

本文涉及的源代碼文件、例子輸入文件都可以在這里下載：/usr/uploads/file/20150609/20150609201400_25520.rar

2 Gaussian的link的基本結構

Link，也就是Gaussian目錄下一個個lxxxx.exe，基本的模板是這樣：
      Program mytest
      Implicit Real*8(A-H,O-Z)
      Dimension Work(1)
      Call InitSCM(1,O,O,Work,IOff,MDV)
      Call ltwd(Work(IOff+1),MDV)
      End

      Subroutine ltwd(V,MDV)
      Implicit Real*8(A-H,O-Z)
      Integer Gen
      Dimension V(MDV),Gen(1000)
      Common /IO/ In, IOut, IPunch
      Common /IOp/ IOp(200)
      COMMON /PHYCON/ PHYCON(30)
#include "commonmol.inc"
      Call Drum(V,MDV)
      Call FileIO(2,-501,1000,Gen,0)
      Write(*,*) "Yoooo",nbasis,phycon(5)  !示例
      Write(*,*) Gen(1),Gen(32)  !示例
      Call ChainX(O)
      End
這里mytest程序是一個殼，也稱dummy程序，而它所調用的ltwd子程序才是本體。Gaussian的link基本都是這種組織模式，各個link的本體都在g09/lxxx.F里，而相應的dummy程序是g09/bsd/main.F里的mlxxx。比如NBO模塊l607的本體是GauNBO子程序，當Gaussian運行到l607的時候會先經過ml607程序才調用它。

外殼中InitSCM函數用來分配一塊空內存給這個link里本體使用，這塊內存起始地址是Work數組的IOff+1元素所在地址，MDV是這塊內存的大小，這樣本體就會接收到長度為MDV的名為V的數組，本體中各種工作都利用這個數組做為臨時空間。實際上帶著#P執行任務時，輸出中的MaxMem=后面的數字就是這個MDV的值（單位為word），%mem的值與MDV相對應。

/IO/是輸入輸出位置與實際文件編號的映射（Drum子程序會將之填充為In=5, IOut=6, IPunch=7）。在Gaussian中對屏幕輸出用IOut，從輸入文件中讀取信息用In，而punch關鍵詞所對應的輸出文件就對應于IPunch，比如punch=MO關鍵詞會將分子軌道信息輸出到fort.7里。
對每一組link所賦予的IOp信息都是通過讀取IOp數組被識別的。比如計算路徑中有一行為6/4=2,7=80/65;，那么在l665的代碼里，IOp(4)就等于2，IOp(7)就等于8。
/PHYCON/里是一些常用的物理常數。
commonmol.inc是一個頭文件，它就在g09目錄下，里面定義了common /mol/，在其中記錄了很多重要信息，比如基函數數(nbasis)、alpha/beta電子數(NAE/NBE)、原子坐標(C數組)、原子序號(IAn數組)等等。

除了Link 0和1以外的所有link在正式工作開始前都要調用Drum子程序，Drum起到初始化的作用，包括打開文件，將上述common中目前空著的數組里填充應有的數值等等。

Call FileIO(2,-501,1000,Gen,0)代表從rwf文件（讀寫文件）中編號501的子文件讀取長度為1000的信息到Gen數組里。rwf文件包含眾多子文件，在手冊里有介紹，其中501子文件對應于/Gen/的信息。Gen數組記錄了計算過程中重要的標量信息，比如維里值Gen(1)、SCF能量Gen(32)等

ChainX必須用在link的末尾，它用來做一些掃尾工作，比如更新rwf文件，關閉已打開的文件，確定接下來運行哪個link等等。參數0代表對執行順序不進行修改，這樣運行完當前link就會按照執行路徑轉到下一個link。

這個模板程序中實際有用的內容就是輸出一些信息。將它編譯成比如叫l666.exe，然后將它插入到單點任務中的某一處，比如在l601之后執行，那么屆時就會看到程序輸出比如：
 Yoooo                      23  6.0221417899999999E+023
    1.998327281295026       -40.19463972178421
說明基函數數(nbasis)為23，阿伏伽德羅常數，即phycon(5)為6.0221417899999999E+023，后兩個分別是維里值和能量，和Gaussian在SCF完畢后輸出的值是一致的。

3 編譯自己寫的link和修改的link

再強調一下，編譯link前必須先完整編譯過一遍Gaussian程序，一些中間剩下的文件要用到。

這里以編譯上面的模板程序為例，假設此文件名叫666lite.F。首先確認$g09root已經設好，確認已經進入了csh，然后輸入source /sob/g09/bsd/g09.login。將/sob/g09/bsd/link.make、666lite.F都拷到當前目錄下，將link.make改名為Makefile，并將其中的MAIN502=后面的內容刪掉，將OBJ502=后面寫上666lite.o。然后輸入mk，這是Gaussian自帶的編譯腳本，就得到了l502.exe。如果想讓這個link叫l666.exe，就直接將它改名為l666.exe即可。

還有一種情況，我們只想對現成的某個Link中某個子程序做一點點修改，而不想重寫，比如這里想讓l601中輸出的偶極矩小數點位數更多些。首先確認$g09root已經設好，確認已經進入了csh，然后輸入source /sob/g09/bsd/g09.login。運行gau-get DQ l601，這里DQ是與輸出多極矩相關的子程序，此時DQ子程序就會被gau-get從l601.F中提取出來成為dq.F，將其中2002號format的格式從F20.4改為F20.8。然后將/sob/g09/bsd/link.make、dq.F都拷到一起，將link.make改名為Makefile，將其中502都替換為601，將MAIN601=后面填上ml601.o，這表明使用Main.F中ml601作為殼。將OBJ601=后面填上dq.o。然后運行mk，編譯腳本會把修改過的DQ、ml601以及l601中其它子程序（在g09目錄下l601.a里，是編譯Gaussian過程剩下的）連接到一起成為l601.exe。用這個l601.exe代替原先的執行一遍，會發現已經生效了，偶極矩小數位數從原先4為變成了8位。

若自己寫一個新link，但又要調用某個link里的子程序，編譯方法見下文的例子。如果Gaussian任務調用自己編譯的模塊出現...World accessible. This must be fixed.字樣，用chmod修改文件權限即可，root的話就改為700。

4 輸出每一步的波函數文件

為了簡單起見，本文只考慮生成SCF波函數的波函數文件。在IRC、優化、scan等任務的每一次循環中，l502做完之后SCF波函數信息就有了，所以只需要寫一個link，假設叫l555（源代碼文件名隨意，這里假設叫555.F），將它插在l502之后，就能讓每一步都輸出波函數文件。原先，波函數文件由l9999產生，AllDun子程序是l9999的本體。打開l9999.F查看其代碼可知，AllDun是通過調用DoWfn，然后DoWfn再調用WrtWfn來輸出波函數文件的，所以要寫的link中只要調用DoWfn即可。然而，DoWfn需要接收大量參數，在調用DoWfn之前必須做很多鋪墊，才能給予DoWfn足夠的信息來輸出波函數文件。這些鋪墊代碼自己寫比較麻煩也容易出錯，所以這里用的策略是將l9999的AllDun全部內容拷過來，然后將與DoWfn無關的代碼刪掉，比如讀取IOp、初始化無關變量、更新Z矩陣、對導數進行旋轉、保存check文件等等內容。有用的信息則必須保留，尤其是計算IEnd變量（記錄的是臨時數組V中已被占用的長度）的過程不要動，避免引起麻煩。由于改動的地方比較繁多細碎，無法一一說明，完整的555.F已經提供在壓縮包里了，和AllDun、WrtWfn原始代碼對照一下就能知道哪些被改動過，其中插入了一些向屏幕輸出內部變量的語句便于調試。判斷哪些變量對DoWfn有用，主要方法是將其接收的各個參數往上搜，直接或間接相關的語句最好進行保留。

在l9999.F中只有IOp(6)的百位數不為0的時候才會調用DoWfn，這里則把條件判斷語句去掉，因為l555存在的目的就是輸出波函數文件。原先l9999.F中在初始化一些變量時用到MOp函數，它是一個在AllDun內部直接定義的函數，但是在我這里編譯不過，可以將之單獨寫成一個function，但由于與MOp相關的變量對輸出波函數文件都沒用，就把MOp聲明和那些變量都刪了。讀取IOp數組、初始化變量部分雖然多數也沒用，但是留著也不礙事，為避免潛在問題，就沒去掉。

還有個問題必須要考慮，也就是每一步輸出的波函數文件名字必須不同才行，但是從輸入文件中讀取的文件名只有一個。分析WrtWfn會發現輸出的文件名最終是在這里控制，所以也要把WrtWfn改了才行，于是首先將WrtWfn的代碼也拷進555.F準備進行修改。

WrtWfn開頭的Read(In,'(A)') FilNam代表從輸入文件讀入字符串到FilNam，FilNam儲存的就是要輸出的波函數文件名。Gaussian假設這個模塊只被調用一遍，但是我們需要將它調用許多遍。每次都讀入一行的話，第二次調用時讀入的文件名就可能為空或者因為讀到文件末尾而報錯，為了避免這種情況，在這行下面插一行backspace(In)。

假設從輸入文件中讀入的文件名叫/sob/H2.wfn，我們想讓每一步的編號順延，即調用l555的時候依次生成/sob/H2_0000.wfn、/sob/H2_0001.wfn...這就需要將這樣內容加到進去
      inamelen=LenFil-4
      inamefullen=LenFil+5
      do itestname=0,9999
           write(namedigit,"i4.4") itestname  ! namedigit是character*4型新增的變量
           FilNam(inamelen+1:inamelen+1)='_'
           FilNam(inamelen+2:inamelen+5)=namedigit
           FilNam(inamefullen-3:inamefullen)=".wfn"
           Inquire(file=FilNam(1:inamefullen),exist=wfnalive)
           if (.not.wfnalive) then    ! wfnalive是新增的邏輯變量
             goto 111  ! 111加在隨后的Write(IOut,1000) FilNam(1:inamefullen)語句開頭
           endif
      end do
LenFil是讀入的文件名的長度，inamelen是不含擴展名的文件名長度，inamefullen是加上編號的文件名長度。這段代碼也就是說，程序會先檢查/sob/H2_0000.wfn是否存在，若不存在這次就用這個名字作為波函數文件名，如果已經存在，就再檢查/sob/H2_0001.wfn是否存在，反復如此，這樣文件名編號就會依次累加了。然后將后面幾行中ISt:LenFil都替換成1:inamefullen。ISt是指讀入的文件名的起始位置（之前是空格），這里不想搞得這么麻煩，也假設用戶不會在開頭加空格，所以令ISt成了1。輸出波函數文件title的部分也進行了修改，title內容會和步數編號一致，便于以后核查。

現在開始編譯。編譯方法和第三節介紹的大體相似，但有些不同，這里link完全是新寫的，但是又要里利用l9999的其它子程序（其實也可以只修改AllDun和WrtWfn而仍用ml9999的殼）。首先進入csh，source一下g09.login，將link.make拷貝為Makefile并與555.F放到一起。將Makefile里502都改為9999，然后把MAIN9999=后面的內容去掉，OBJ9999=后面寫為555.o。運行mk，編譯器編譯l9999.exe時相當于把原本l9999里的外殼、AllDun、WrtWfn都用555.F里所定義的替換，l9999內其余子程序依然沿用。編譯好后將l9999.exe改名為l555.exe。

現在開始測試，這里用三重態H2拉伸scan為例，輸入文件在壓縮包里，共20步。輸入文件中原始的route section仍被保留，但已經用嘆號注釋掉了。注意由于要從很近的距離0.2埃開始拉，必須用geom=nocrowd避免距離檢測。此任務的執行路徑如下
#P nonstd <-----使用非標準路徑
 1/38=1,60=1/1,8;
 2/12=1,17=6,18=5,29=3,40=1/2;
 3/5=16,8=22,11=2,16=1,25=1,30=1,74=-5/1,2,3;
 4//1;
 5/5=2,38=5/2,55;
 5//55;  <-----插入了l555
 6/7=2,8=2,9=2,10=2,18=1,28=1/1;
 1/60=1/8(1);
 99/9=1/99;
 2/12=1,29=3/2;
 3/5=16,8=22,11=2,16=1,25=1,30=1,74=-5/1,2,3;
 4/5=5,16=3/1;
 5/5=2,38=5/2,55;  <-----插入了l555
 1/60=1/8(-4);
 99/9=1/99;  <-----由于每步都已經輸出了波函數文件，就不需要末尾再次輸出了，所以沒有傳入IOp(6)=100/200/300的信息。
在這個scan任務中，從執行路徑可以看到是先對初始結構做一次SCF，然后再把每個距離也都做SCF，每次遇到(-4)標簽后回退開始新一輪循環，直到到達步數上限。這里用了兩種插入link的寫法，一種是直接寫一行5//55;代表只調用l555且不傳入任何IOp，第二種是在5/5=2,38=5/2;的分號前插入",55"，兩種做法效果是一致的。盡管后者會把IOp(5)=2、IOp(38=2)這樣的信息也傳給l555，但是AllDun里并沒利用到它們，等于沒寫。如果沒有將l555.exe放進g09目錄下而是別的地方，比如當前目錄，為了讓Gaussian能找到它，要在輸入文件開頭寫上%subst l555 .（小點代表當前目錄）。執行后就生成了H2_0000.wfn~H2_0021.wfn，其中H2_0000.wfn對應初始結構。

還有一種插入l555的方法，雖然沒有非標準路徑方法那么靈活，但更為方便，免得列出冗長的執行路徑，也就是在route section里寫上ExtraLinks=l555，這樣在每次調用完所有overlay 5的link后（此例僅為l502），就會調用l555。

5 輸出IRC任務每個點的波函數文件

上一節看似已經達到了我們的目的，但是對于IRC任務還需要考慮更多。附件里包含一個HCN的氫轉移IRC任務，執行路徑應當這樣寫：
[略] （第一次的l502后面也插入l555）
 2/9=110,29=1/2;
 3/5=1,6=6,7=101,11=2,16=1,25=1,30=1,71=2,74=-5/1,2,3;
 4/5=5,16=3/1;
 5/5=2,38=5/2;
 5//55; <----插入了l555
 8/6=4,10=90,11=11/1;
 11/6=1,8=1,9=11,15=111,16=1/1,2,10;
 10/6=1,7=6,13=1/2;
 7/10=1,18=20,25=1/1,2,3,16;
 1/14=-1,18=20,39=30,44=3/23(-8);  <---- 更新坐標，開始一輪新循環
[略]
執行后會發現直接這樣做行不通。因為IRC的每一步實際上執行的是限制性優化，這個限制性優化需要幾步才能收斂。假設IRC包含20步，平均每步需要三次優化，那么l555將被調用約60次，即最后我們將得到約60個波函數文件！但是其中只有21個（其中1個是初始的TS結構）波函數文件是有用的。因此必須想辦法只讓每次限制性優化收斂的步才輸出波函數文件，有兩種實現辦法：

第一種方法是每次循環都輸出波函數文件，而在最后根據輸出文件信息將波函數文件中有用的分揀出來。使用shell腳本比較方便，內容如下：
#!/bin/bash
dir=HCN             # The directory where all wavefunction files were stored
outname=HCN_IRC.out # Gaussian output filename
i=0
j=0
res=(pass `grep -E "Delta-x Convergence NOT Met|Delta-x Convergence Met" $outname|awk '{print $3}'`)
for filnam in `ls $dir/*.wfn`
do
echo ${res[$i]},$filnam
if [ "${res[$i]}" == "NOT" ] && [ $i!=0 ]; then
 rm -f $filnam
else
 mv -f $filnam $dir/`printf "%4.4i" $j`.wfn
 j=$(($j+1))
fi
i=$(($i+1))
done
首先執行IRC任務，輸出文件為當前目錄下HCN_IRC.out，波函數文件都輸出到HCN/目錄下。執行這個腳本后（別忘了先切換回bash），HCN_IRC.out會被讀取，每一步是否是收斂步會在屏幕上顯示，"Not"代表這個波函數文件對應的是未收斂步，會被刪掉；"Met"代表這個波函數文件是收斂步，將被保留，并且文件名將會被整理得連貫。過渡態，也就是第一步也這么處理，但顯示為"Pass"。這樣，最終剩下的文件名為0000.wfn~0020.wfn，0是過渡態，1~10和11~20對應兩個方向。

第二種辦法這里只給出梗概，怎么具體實現就略過了。通過分析，會發現l123.F的PreDWI子程序與這個問題相關，當cordon變量為true時說明限制性優化已收斂。所以可以改寫PreDWI，讓收斂時傳遞一個信息給隨后的l555。傳遞信息的方式多種多樣，比如可以通過/Gen/里面的空余位置，也可以用call system("touch conv")命令，即僅當收斂的時候在當前目錄建立conv空文件，WrtWfn中通過Inquire(file="conv",exist=ircconv)就能靠ircconv邏輯變量判斷出是否已經收斂了，這就能使只有收斂的時候才輸出波函數文件。注意還有一些細節問題要考慮，因為當WrtWfn得知限制性優化已收斂的時候，此時的波函數已經是下一步的了，需要靠一點trick來解決這個矛盾，請自行嘗試。

有了IRC過程每個點的波函數文件，我們就可以用它們做“這樣或那樣”的事情了。筆者在隨后的文章中將介紹怎么繪制IRC過程中各種實空間函數變化的動畫，這比起觀看靜態的圖像更有趣，也能獲得更豐富的信息。由于本文的l555里面刪掉了一些信息，有可能導致特殊情況下所得波函數文件不正確，屆時要以l9999輸出的為準。