利用Multiwfn令Dalton計算時使用其它程序產生的軌道作為初猜

利用Multiwfn令Dalton計算時使用其它程序產生的軌道作為初猜
Using Multiwfn to take orbitals generated by other programs as initial guess in Dalton calculations

文/Sobereva@北京科音 2025-Apr-5

1 前言

Dalton這個量子化學程序雖然用起來復雜，但有很多獨特的亮點，例如可以在SOPPA(CCSD)和CASSCF級別下計算各種磁屬性（筆者最近在Chem. Eur. J., e202404138 (2025)文中就用Dalton在CASSCF下算了NICS），可以用響應方式較準確地計算磷光壽命，可以算雙光子吸收（北京科音高級量子化學培訓班http://www.keinsci.com/KAQC里專門有一節講了)，可以用高級別耦合簇方法算激發態（見《使用Dalton通過CC3方法極高精度計算激發態》http://www.shanxitv.org/463），等等。可惜Dalton的SCF收斂方面做得一般，比不上Gaussian，而且單靠Dalton無法產生一些重要的軌道，比如UHF計算產生的UNO軌道經常被用于CASSCF的初猜，但Dalton甚至連UHF都做不了，因為Dalton里開殼層只支持RO形式。

為解決以上問題，從2025-Apr-5更新的Multiwfn開始，Multiwfn產生Dalton輸入文件的功能支持了將Multiwfn中當前的軌道波函數寫入到Dalton的輸入文件.dal里作為初猜波函數。因此借助Multiwfn，Dalton用戶可以用Gaussian、ORCA、GAMESS-US等諸多主流的量子化學程序產生的波函數當初猜，從而解決Dalton遇到SCF不收斂的問題。而且，Dalton用戶還可以直接用Multiwfn產生的UNO、定域化分子軌道（LMO，見http://www.shanxitv.org/380）等當CASSCF初猜。

下面就通過例子簡單演示一下具體實現方法，本文中涉及的文件都可以在http://www.shanxitv.org/attach/740/file.zip下載。Multiwfn可以在http://www.shanxitv.org/multiwfn免費下載，相關常識看《Multiwfn FAQ》（http://www.shanxitv.org/452）。如果讀者還不知道Dalton怎么安裝、運行、產生輸入文件的話，看《量子化學程序Dalton的編譯方法和運行方式簡介》（http://www.shanxitv.org/463）。本文用的Multiwfn是2025-Apr-5更新的版本，Gaussian是G16 B.01，Dalton是2022版。

如果本文介紹的Multiwfn的功能給你的研究帶來了幫助，請在論文中引用Multiwfn剛啟動時提示的Multiwfn程序的原文，這是對Multiwfn開發和維護最好的支持！

2 把Gaussian收斂的HF波函數當Dalton做CCSD(T)的初猜

這一節以CH3NH2為例演示把Gaussian收斂的HF/cc-pVTZ波函數當Dalton做CCSD(T)/cc-pVTZ的初猜。首先用Gaussian做一個普通的單點計算，輸入文件是本文文件包里的CH3NH2.gjf，內容如下所示。算完后把得到的CH3NH2.chk用formchk轉成CH3NH2.fch。

%chk=D:\CH3NH2.chk
# HF/cc-pVTZ int=NoBasisTransform
[空行]
test
[空行]
0 1
C                  0.05159500    0.70381800    0.00000000
H                  0.59439800    1.06209400    0.88166000
H                  0.59439800    1.06209400   -0.88166000
H                 -0.94293100    1.18451300    0.00000000
N                  0.05159500   -0.76078400    0.00000000
H                 -0.45830300   -1.10306000    0.81258800
H                 -0.45830300   -1.10306000   -0.81258800
[空行]
[空行]

之后啟動Multiwfn，載入CH3NH2.fch，依次輸入
100 //其它功能（Part 1）
2 //產生輸入文件
19 //產生Dalton輸入文件
CCSDpT.dal //要產生的輸入文件名
[回車] //要產生的.mol文件名用默認的CH3NH2.mol
y //由于當前內存里有基函數信息，Multiwfn問你是否把軌道展開系數寫入到.dal文件里
1 //使用Dalton標準格式（4F18.14）方式寫入。即四個值換一行，每個值占18列、小數位數14個

此時當前目錄下就有了CCSDpT.dal和CH3NH2.mol。把CH3NH2.mol里的基組名都替換為cc-pVTZ_emsl。目前CCSDpT.dal里的計算設置對應B3LYP算單點，自行把關鍵詞改為CCSD(T)計算用的，此時CCSDpT.dal的內容如下
**DALTON INPUT
.RUN WAVE FUNCTIONS
**WAVE FUNCTIONS
.CC
*CC INPUT
.CC(T)
*ORBITAL
.MOSTART //代表從**MOLORB部分讀取軌道展開系數作為初猜
FORM18
**MOLORB (Punched by Multiwfn)
-0.00001750284370 0.00060444251600 0.00013849165900 -0.00001024308040
0.00002442752580 -0.00001726827680 -0.00014237439200 0.00004058047020
...略
**END OF INPUT

用Dalton進行計算，可以看到SCF一輪就收斂了，達到了我們的目的：
@ *** DIIS converged in   1 iterations !
@     Converged SCF energy, gradient:    -95.252488149534    1.44D-07
    - total time used in SIRFCK :              0.00 seconds
之后CCSD(T)計算也順利跑完了。

一些注意事項和相關信息：
? 對于廣義或部分廣義收縮基組，都應當用int=NoBasisTransform關鍵詞，拿不準就始終帶上。詳見《將Gaussian等程序收斂的波函數作為ORCA的初猜波函數的方法》（http://www.shanxitv.org/517）里對此關鍵詞的解釋。
? 把Gaussian的cc-pVDZ、cc-pVTZ的波函數當Dalton的初猜時，Dalton要用cc-pVDZ_emsl、cc-pVTZ_emsl基組，而不要用不帶_emsl后綴的，否則還是需要迭代很多次才能收斂（甚至可能不收斂），但帶不帶_emsl的版本收斂后的結果是相同的。而如果是def、def2系列基組，沒什么特別要注意的，比如Gaussian里def-TZVP和def2-TZVP做的計算，在Dalton里基組名分別寫成Turbomole-TZVP、def2_tzvp即可，注意對大小寫有要求，必須對應Dalton目錄下的basis子目錄里的基組文件名。
? Dalton默認是基于球諧型基函數做的計算，因此如果在Gaussian里使用的是比如6-31G*等D殼層默認是笛卡爾型基函數的基組，應當寫上5d關鍵詞強行要求用球諧型，詳見《談談5d、6d型d殼層基函數與它們在Gaussian中的標識》（http://www.shanxitv.org/51）。雖然Dalton也可以用笛卡爾型基函數，但Multiwfn往Dalton輸入文件里寫入初猜軌道系數的功能不支持此情況。
? 只要一個量子化學程序能輸出Multiwfn可以識別的帶有基函數信息的波函數文件，并且用的是球諧型基函數，則波函數都可以通過以上方式作為Dalton的初猜。支持的文件詳見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。例如ORCA做計算產生的molden文件也可以載入Multiwfn后像上面那樣產生帶初猜波函數的Dalton輸入文件。

3 把Gaussian收斂的DFT波函數當Dalton做DFT的初猜

此例以順鉑體系為例，把Gaussian收斂的PBE0波函數當Dalton做PBE0計算的初猜，對Pt用Stuttgart小核贗勢，對配體用6-311G*。Gaussian輸入文件如下（本文文件包里的Ptcoord.gjf）：

%chk=D:\Ptcoord.chk
#P PBE1PBE/genecp int=NoBasisTransform
[空行]
b3lyp/def2TZVP opted //眾所周知這是純給自己看的標題行，別問為什么寫成這個
[空行]
0 1
[坐標略...]
[空行]
N H Cl
6-311G*
****
Pt
SDD
****
[空行]
Pt
SDD
[空行]
[空行]

按照和上一節例子相同的過程，把chk轉成fch，再用Multiwfn產生.dal和.mol文件。.dal里的泛函名改成PBE0，然后設置.mol里的基組和贗勢成為下面這樣，當前對Pt用的贗勢和贗勢基組和Gaussian里的情況一樣（若不確定的話，可以Gaussian里用gfinput關鍵詞輸出定義，并和Dalton的基組、贗勢目錄下的相應文件對照）。

ATOMBASIS
test molecule
Generated by Multiwfn
Atomtypes=4 Angstrom Nosymmetry charge=0
Charge=1.0 Atoms=6 Basis=6-311G*
H1   -0.82596600    1.64390200    2.14978600
H2    0.00000000    2.40772900    0.93575300
H3    0.82596600    1.64390200    2.14978600
H4   -0.82596600   -1.64390200    2.14978600
H5   -0.00000000   -2.40772900    0.93575300
H6    0.82596600   -1.64390200    2.14978600
Charge=7.0 Atoms=2 Basis=6-311G*
N1    0.00000000    1.59755500    1.56108400
N2   -0.00000000   -1.59755500    1.56108400
Charge=17.0 Atoms=2 Basis=6-311G*
Cl1    0.00000000    1.70827400   -1.36819100
Cl2   -0.00000000   -1.70827400   -1.36819100
Charge=78.0 Atoms=1 Basis=stuttgart_rsc_1997_ecp ECP=stuttgart_rsc_1997_ecp
Pt1   -0.00000000   -0.00000000    0.18195700

用Dalton計算后，經過4輪SCF就輕易收斂了，結果為-1152.597056，和Gaussian給出的-1152.596983非常接近。而如果此例Gaussian和Dalton都用HF的話，會和上例一樣，結果精確相同且SCF一輪可收斂。

4 把UNO軌道當Dalton做CASSCF的初猜

之前在《CASSCF計算雙自由基以及雙自由基特征的計算》（http://www.shanxitv.org/264）里示例過用Gaussian基于UNO初猜軌道做CASSCF(2,2)計算的方法，不了解這種計算和UNO概念的話建議先看看。這一節演示使用Gaussian產生UHF波函數，讀入Multiwfn并產生UNO軌道，然后寫入.dal文件用于Dalton做CASSCF(2,2)計算的初猜軌道的流程。示例體系還是264號博文里的C4H8雙自由基，基組還是此文里的6-31G*。

首先用以下Gaussian輸入文件（本文文件包里的C4H8.gjf）對C4H8做UHF計算得到對稱破缺波函數。不了解為什么用guess=mix的話看《談談片段組合波函數與自旋極化單重態》（http://www.shanxitv.org/82）。

%chk=D:\C4H8.chk
# UHF/6-31G* guess=mix 5d
[空行]
ub3lyp/6-31g(d) opted
[空行]
0 1
C     -0.74742092    1.77656753    0.00000000
H     -0.62438907    2.32965189    0.92333358
H     -0.62438907    2.32965189   -0.92333358
C     -0.74742092    0.30933780    0.00000000
H     -1.24978876   -0.09122321    0.88294562
H     -1.24978876   -0.09122321   -0.88294562
C      0.74742092   -0.30933780    0.00000000
H      1.24978876    0.09122321   -0.88294562
H      1.24978876    0.09122321    0.88294562
C      0.74742092   -1.77656753    0.00000000
H      0.62438907   -2.32965189   -0.92333358
H      0.62438907   -2.32965189    0.92333358
[空行]
[空行]

計算完畢后把C4H8.chk轉成C4H8.fch，載入Multiwfn，然后依次輸入
200 //其它功能（Part 2）
16 //基于fch文件里的密度矩陣產生自然軌道。這個功能在《在Multiwfn中基于fch產生自然軌道的方法與激發態波函數、自旋自然軌道分析實例》（http://www.shanxitv.org/403）里有專門的介紹
SCF //讀取fch里的SCF類型的密度矩陣
1 //產生空間自然軌道
此時從屏幕上看到如下信息，這便是UNO軌道的占據數。其中最高占據自然軌道（HONO）和最低非占據自然軌道（LUNO）的占據數分別是1.123225和0.876775，都偏離整數顯著，很滿足期望。

Occupation numbers:
    2.000000    2.000000    1.999999    1.999999    1.999904    1.999898
    1.999881    1.999835    1.999309    1.999120    1.999004    1.998977
    1.996829    1.996536    1.996053    1.123225    0.876775    0.003947
    0.003464    0.003171    0.001023    0.000996    0.000880    0.000691
...略

接著輸入
n //不把自然軌道導出成new.mwfn文件并重新載入，因為當前不需要如此（如果你想看一下當前的UNO軌道圖像的話，應選y，之后用Multiwfn主功能0觀看軌道）
0 //返回主菜單
100 //其它功能（Part 1）
2 //產生輸入文件
19 //產生Dalton輸入文件
CASSCF.dal //要產生的輸入文件名
[回車] //要產生的.mol文件名用默認的C4H8.mol
y //把當前內存里的軌道（即UNO）的展開系數寫入到.dal文件里
1 //使用Dalton標準格式（4F18.14）方式寫入

當前C4H8.mol里的基組默認就是6-31G*，不用改。修改CASSCF.dal，把以下兩行
.DFT
B3LYPg
替換為CASSCF(2,2)計算的設置：
.MCSCF
*CONFIGURATION INPUT
.CAS SPACE
2
.ELECTRONS
2
.INACTIVE
15

Dalton計算完畢后，活性空間內自然軌道占據數為1.2679和0.7320，能量為-156.024938 Hartree，和Gaussian在同級別做CASSCF精確一致。你若想看Dalton做CASSCF產生的自然軌道，把計算產生的.tar.gz文件里的molden.inp解壓出來載入Multiwfn，并效仿《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）的方式觀看即可。

5 關于線性依賴基函數的問題

較大基組（尤其是帶彌散函數的），計算原子分布較致密的體系很容易出現線性依賴問題，Gaussian和Dalton此時都會自動去除線性依賴基函數，導致最終求解出來的軌道數會少于基函數數目。這種情況下，Gaussian輸出信息中的NBsUse（實際用的基函數數目，等同于實際求解出來的軌道數）小于NBasis，Dalton輸出文件中的Total number of orbitals小于Number of basis functions。這給上文的傳軌道帶來了麻煩。有兩個解決方法：
(1)不去除線性依賴基函數：Gaussian計算時帶著IOp(3/32=2)關鍵詞、Dalton計算時在**WAVE FUNCTIONS里面加上以下內容，從而讓兩個程序都不去除線性依賴基函數
*ORBITAL
.AO DELETE
1E-100
這種做法的缺點是可能導致數值不穩定性。
(2)自動去除線性依賴基函數，但讓Dalton少讀取軌道：比如Gaussian計算顯示NBsUse=50、NBasis=54，因此自動去除了4個線性依賴基函數。因此當Multiwfn載入fch文件后，最后4個軌道的展開系數都會為0，導出到.dal文件里也是如此。此時Dalton若照常讀取全部軌道，計算一開始會報錯*** ORTHO-FATAL ERROR ***。解決方法是在*ORBITAL里寫上
.DELETE
4
這樣Dalton就不會讀取最后4個軌道了，計算可以正常進行。

6 .dal里的軌道展開系數存在星號而無法讀取的問題

有時候Multiwfn寫入到.dal里的軌道展開系數會存在一堆星號，導致Dalton讀取軌道時提示input conversion error錯誤。這不是Multiwfn的問題，而是Dalton自身設計不周。Dalton定義的**MOLORB字段中記錄展開系數的格式在前面說了，每個數值都是F18.14格式，因此當有的軌道展開系數特別大，超過了其整數位記錄上限，就會被記錄成一堆星號。這種問題在基組較彌散而又不讓自動去除線性依賴基函數的情況下出現概率大。這是為什么Multiwfn導出Dalton輸入文件的時候特意讓你選擇記錄的格式，前面的例子都是選1用4F18.14格式，還可以選2用4E20.12格式，這是使用科學計數法輸出，因此就不會出現數值過大記錄成星號的問題了。相應地，必須自己修改Dalton的讀取軌道展開系數的源代碼。對于Dalton 2022版，需要修改Dalton源代碼目錄下的DALTON/sirius/sirgp.F文件，把PFMT = '(4F18.14)'替換為PFMT = '(4E20.12)'，之后重新編譯即可。

7 其它

Multiwfn就相當一個大熔爐，能讀取諸多量子化學程序產生的波函數，又能導出波函數作為諸多量子化學程序的初猜，這在《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）里專門說了。對于傳軌道到ORCA的情況我之前還專門寫過《將Gaussian等程序收斂的波函數作為ORCA的初猜波函數的方法》（http://www.shanxitv.org/517）。

MOKIT也是很好的可以實現不同程序間傳軌道的工具，fch到.dal和.mol的轉換可以用里面的fch2dal子程序。相對于MOKIT，本文介紹的用Multiwfn傳軌道有幾個額外的好處：
(1)Multiwfn不需要安裝MOKIT所需的Python環境。尤其是對于Windows用戶，Multiwfn解壓即用明顯更為方便
(2)Multiwfn特意考慮了.dal里記錄軌道展開系數的格式的問題并給了解決辦法，上一節說了
(3)Multiwfn自身可以產生諸多類型的軌道（定域化分子軌道、空間/自旋自然軌道、NAdO、AdNDP、NTO、雙正交化軌道等等）并輸出作為Dalton的初猜
而fch2dal也有額外的優點，它會把fch里的基組、贗勢的具體定義直接寫到輸出的.mol里，因此如果Gaussian計算時用的基組、贗勢在Dalton里沒有自帶也能直接計算而不需要自定義。并且fch2dal支持用笛卡爾型基函數的情況。