使用Dalton通過CC3方法極高精度計算激發態

使用Dalton通過CC3方法極高精度計算激發態

文/Sobereva@北京科音   2023-Dec-30

0 前言

耦合簇方法可以通過線性響應（linear response, LR）計算激發態。在常見范疇內，精度和耗時關系是LR-CC3>LR-CCSDR(3)≈EOM-CCSD(T)>LR-CCSD=EOM-CCSD>=CC2。LR-CC3如今被普遍視為是激發態計算的金標準（多參考特征很強以及雙電子激發特征占主導的情況除外），激發能能精確到零點零幾eV。但LR-CC3結合像樣基組只能算得動個位數原子，對更大體系則可以用明顯更便宜的LR-CCSDR(3)，精度也已經非常好了。

Dalton程序在耦合簇方面功能豐富，可以做LR-CC3、LR-CCSDR(3)、LR-CCSD、LR-CC2激發能計算，對于LR-CC2/CCSD/CC3還可以得到包括振子強度在內的諸多躍遷屬性。Dalton的輸入文件對于初學者來說比較復雜抽象，本文用盡可能簡單易懂的語言完整介紹一下用免費的Dalton做這些激發態耦合簇計算的方法，使用有C2v對稱性的非常典型的分子甲醛為例。讀者如果對Dalton不熟悉，務必先閱讀我之前寫的《量子化學程序Dalton的編譯方法和運行方式簡介》（http://www.shanxitv.org/463），我假定讀者已經認真看過此文。

值得一提的是，雖然Dalton做LR-CC3計算如今在文獻中很常見，但Dalton在這方面不是最快的，而且代碼還沒有并行化。如果追求更好的效率，可以用專門做耦合簇計算的e^T程序（https://www.etprogram.org），它的LR-CC3是所有程序中明顯最快的。另外，LR-CCSD激發能和振子強度沒任何必要非得用Dalton計算，常用的Gaussian和ORCA的EOM-CCSD也都做得很好，結果和Dalton的LR-CCSD是等同的。Dalton做耦合簇的激發能計算不支持考慮溶劑模型，而G16的EOM-CCSD則可以。還值得一提的是，Dalton的衍生程序LSDalton還額外支持RI-CC2，可以比Dalton的LR-CC2更高效率地計算較大體系，這不屬于本文的范疇。

本文涉及的輸入輸出文件都可以在http://www.shanxitv.org/attach/693/file.zip里獲得。計算用的是2023-Dec-22通過git下載的Dalton最新的開發版。

本文下面的例子首先要重復知名的TBE1激發能測試集（J. Chem. Phys., 128, 134110 (2008)）里的LR-CC3/def-TZVP算的甲醛的1-1A2、1-1B1、2-1A1三個單重態的垂直激發能。如文章表1所示，結果分別為3.95、9.18、10.45 eV。此文用的幾何結構是MP2/6-31G*級別優化的，對應的結構文件是file文件包里的H2CO.xyz，本文的計算也將基于這個結構。

1 確定不可約表示順序

為了指定各個不可約表示的激發態各算幾個態，需要先知道C2v點群的各個不可約表示在Dalton程序中的順序。最簡單的方法是做個DFT單點計算，看一開始輸出的不可約表示的順序。為此，我們用Multiwfn（http://www.shanxitv.org/multiwfn）創建這個任務的輸入文件。啟動Multiwfn，然后依次輸入
H2CO.xyz
100   //其它功能（Part 1）
2   //導出各種文件
19   //Dalton
DFT.dal  //在當前目錄下產生DFT.dal，對應B3LYP/6-31G*計算設置
H2CO.mol   //在當前目錄下產生Dalton格式的體系定義文件H2CO.mol

由于當前計算要考慮對稱性，因此把H2CO.mol里的Nosymmetry刪掉。然后基于DFT.dal和H2CO.mol做計算，從輸出文件DFT_H2CO.out里可以看到在SCF開始之前就輸出了當前自動判斷的點群，以及相應的各個不可約表示的順序

@    Full point group is: C(2v)         
@    Represented as:      C2v

@  * The irrep name for each symmetry:    1: A1     2: B1     3: B2     4: A2

2 做LR-CC3激發能計算

寫一個文本文件LRCC3.dal，內容如下（這是最簡單寫法，默認設置就已經適合的選項就沒再做額外設置）

**DALTON
.RUN WAVE FUNCTIONS
**WAVE FUNCTION
.CC
*CC INP
.CC3
*CCEXCI
.NCCEXCI
2 1 0 1
**END OF DALTON

對當前體系默認是做閉殼層HF計算得到參考態波函數。**DALTON控制任務類型，.RUN WAVE FUNCTIONS要求算單點并得到波函數。**WAVE FUNCTION設置波函數計算類型，.CC要求做耦合簇計算。*CC INP設置做什么形式耦合簇計算，.CC3要求做CC3計算。*CCEXCI要求做耦合簇的激發能計算，.NCCEXCI下面按照不可約表示的順序指定各個不可約表示的能量最低激發態各計算幾個。當前為了重復TBE1測試集的數據，要求A1算兩個（由此得到1-1A1和2-1A1），B1算1個（得到1-1B1），B2不算，A2算1個（得到1-1A2）。

然后把前述的H2CO.mol復制為H2CO_TZVP.mol，手動把里面的6-31G*替換為Turbomole-TZVP，這是Dalton內置的def-TZVP基組的寫法，對應Dalton目錄下的basis子目錄中的Turbomol-TZVP這個文件的名字。

然后基于LRCC3.dal和H2CO_TZVP.mol做計算，得到LRCC3_H2CO_TZVP.out。本節的文件都已經提供在了前述的file文件包里。打開輸出文件，從里面可以看到如下激發能信息，2-1A1為10.44842 eV，1-1B1為9.18343 eV，1-1A2為3.94711 eV，可見和TBE1原文里的10.45、9.18、3.95 eV完全一致。表中||T1||是單激發貢獻，和TBE1原文里給出的也是一致的。

 +=============================================================================+
 |  sym. | Exci.  |        CC3        Excitation energies            | ||T1||  |
 |(spin, |        +------------------------------------------------------------+
 | spat) |        |     Hartree    |       eV.      |     cm-1       |    %    |
 +=============================================================================+
 | ^1A1  |    1   |     0.3502215  |       9.53001  |     76864.734  |  91.00  |
 | ^1A1  |    2   |     0.3839723  |      10.44842  |     84272.170  |  91.32  |
 +-----------------------------------------------------------------------------+
 | ^1B1  |    1   |     0.3374848  |       9.18343  |     74069.360  |  90.93  |
 +-----------------------------------------------------------------------------+
 | ^1A2  |    1   |     0.1450537  |       3.94711  |     31835.610  |  91.16  |
 +=============================================================================+

如果想要計算三重態激發態，就在.NCCEXCI下面的第二行定義。比如下面的寫法，代表除了計算如上那些單重態激發態以外，還計算1-3A2和1-3A1三重態激發態
.NCCEXCI
2 1 0 1
1 0 0 1
如果不需要計算單重態激發態，就寫
.NCCEXCI
0 0 0 0
1 0 0 1

筆者試了幾個Dalton版本，包括2016、2018、2022開發版，算出來的三重態激發能都是錯的，和TBE1表2里的1-3A2和1-3A1三重態激發能明顯不符，而且||T1||嚴重偏低，我認為是bug。如果不利用對稱性，即.mol文件里帶著Nosymmetry，并且.NCCEXCI下面直接指定計算的激發態的總態數，則問題輕得多，可以看到||T1||都顯著高于90%。

值得一提的是LR-CC3是對激發態一個一個獨立計算的，故算的態數越多耗時明顯越高。

如果只需要計算CC3基態能量，去掉*CCEXCI及下面的部分即可，相對于LR-CC3算激發態部分的耗時來說可以忽略不計。

3 做LR-CC3激發能+振子強度的計算

如果對上面算的激發態還要計算振子強度，就創建LRCC3mom.dal文件，寫入以下內容，然后進行計算。其中**INTEGRAL指定要算基函數之間哪些類型的積分，.DIPLEN要求算長度形式的偶極積分，這是因為算振子強度要算躍遷偶極矩，會用到這類積分。*CCOPA模塊用來計算耦合簇的基態到激發態的躍遷屬性，對當前用的Dalton版本支持CCS、CC2、CCSD、CC3，其中CC3僅限單重態激發態。.DIPOLE要求計算長度形式的躍遷偶極矩及振子強度。

**DALTON
.RUN WAVE FUNCTIONS
**INTEGRAL
.DIPLEN
**WAVE FUNCTION
.CC
*CC INP
.CC3
*CCEXCI
.NCCEXCI
2 1 0 1
*CCOPA
.DIPOLE
**END OF DALTON

從輸出文件（file文件包里的LRCC3mom_H2CO_TZVP.out）可以看到各個激發態的躍遷偶極矩和振子強度，比如2-1A1的如下所示，振子強度為0.34867845。這個值很接近TBE1原文表VII里這個態的LR-CC2和LR-CCSD振子強度（這篇文章里沒給CC3的振子強度），分別為0.368和0.374。

     Transition from ground state to:
     number, multiplicity, symmetry :    2  ^1A1
     frequency :   0.3839722612 a.u.    10.44842 e.V.     84272.2 cm^-1

     +-----------+-----------------+-----------------+---------------------+
     | operator  |   left moment   |  right moment   | transition strength |
     +-----------+-----------------+-----------------+---------------------+
     | XDIPLEN   |      0.00000000 |      0.00000000 |       0.00000000    |
     | YDIPLEN   |      0.00000000 |      0.00000000 |       0.00000000    |
     | ZDIPLEN   |     -0.83921304 |     -1.62309631 |       1.36212358    |
     +-----------+-----------------+-----------------+---------------------+
       oscillator strength (length gauge)   :      0.34867845

上面表格里的ZDIPLEN（偶極矩的Z分量算符）的transition strength值1.36212358對應的是躍遷偶極矩Z分量的平方。

如果你只需要躍遷偶極矩的特定分量，比如Z分量，可以把.DIPOLE改為
.OPERATOR
ZDIPLEN
這樣耗時會比計算所有分量時稍微低一點。

4 LR-CC2/CCSD/CCSDR(3)激發態計算

LR-CC2、LR-CCSD、LR-CCSDR(3)計算只需要把前面例子里的.CC3分別改成.CC2、.CCSD、.CCR(3)即可。注意LR-CCSDR(3)只能算激發能而無法算包括振子強度在內的躍遷屬性。

這些級別的激發能計算的輸入輸出文件在file文件包里的other目錄下都給了，耗時跟LR-CC3比都不值得一提。整體來說按LR-CC2、LR-CCSD、LR-CCSDR(3)的順序，結果和金標準LR-CC3的越來越接近。比如它們的1-1B1激發能分別為9.348、9.257、9.186 eV，LR-CC3的為9.183 eV。所以LR-CC3激發能很難算得動的時候用LR-CCSDR(3)是很好的平替。

5 凍核設置

Dalton默認是不凍核的。顯然恰當凍核可以節約耗時而幾乎不影響精度。這里演示對H2CO做LR-CC3計算時凍結能量最低兩個分子軌道（對應C和O的1s軌道）的方式。首先得得知它們的不可約表示，從之前B3LYP/6-31G*的輸出文件中可以看到各個占據軌道能量的負值（Koopmans定理對應的電離能，相關知識參考《正確地認識分子的能隙(gap)、HOMO和LUMO》http://www.shanxitv.org/543）。

@ Summary of DFT KT ionization potentials [eV]:

@ Symmetry 1 (A1 ) :   521.511     279.798      28.621      17.429      12.139

@ Symmetry 2 (B1 ) :    10.711

@ Symmetry 3 (B2 ) :    13.454       7.303

@ Symmetry 4 (A2 ) :  No IPs in this symmetry

顯然能量最低兩個軌道是A1。

寫一個LR-CC3激發能計算輸入文件LRCC3_FC.dal，內容如下。其中.FROIMP下面第一行指定對各個不可約表示凍結多少個能量最低占據軌道，下面第二行指定對各個不可約表示凍結多少個能量最高空軌道。當前凍結的是能量最低的兩個A1軌道。

**DALTON
.RUN WAVE FUNCTIONS
**WAVE FUNCTION
.CC
*CC INP
.CC3
.FROIMP
2 0 0 0
0 0 0 0
*CCEXCI
.NCCEXCI
2 1 0 1
**END OF DALTON

從file文件包里的此任務的輸出文件LRCC3_FC_H2CO_TZVP.out可看到，現在凍核后耗時是1分49秒，1-1B1激發能是9.18872 eV，之前沒凍核時是3分32秒，1-1B1激發能是9.18343 eV。明顯凍核對激發能的影響可完全忽略不計，而耗時大為降低，因而有重要實際意義。