http://www.shanxitv.org/category/QC/ 量子化學 http://www.shanxitv.org/702 2024-02-29T10:22:00+08:00 使用SYSMOIC程序繪制磁感生電流圖和計算鍵電流強度 文/Sobereva@北京科音  2024-Feb-29 0 前言 筆者在《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）中介紹過，磁感生電流是一種重要、直觀的考察分子芳香性的方式。SYSMOIC是一個基于continuous transformation of the origin of the current density (CTOCD)方法考察化學體系的磁感生電流相關問題的程序。本文將結合苯分子為例，介紹這個程序的最關鍵用法。SYSMOIC程序和相關原理的詳細介紹可以看原文J. Chem. Inf. Model., 61, 270 (2021)。 AICD和GIMIC也是兩個重要的考察磁感生電流的程序，筆者在下文有專門的介紹。這些文章里介紹的背景知識對理解本文很有幫助 ? 使用AICD 2.0繪制磁感應電流圖（http://www.shanxitv.org/294） ? 使用AICD程序研究電子離域性和磁感應電流密度（http://www.shanxitv.org/147） ? 考察分子磁感生電流的程序GIMIC 2.0的使用（http://www.shanxitv.org/491） ? Utilizing AICD and GIMIC programs to study magnetically induced current density（http://www.shanxitv.org/148） 之所以筆者寫此文介紹SYSMOIC，在于它對以上兩個程序有明顯互補性。AICD程序是將感生電流以箭頭方式繪制在AICD的等值面上，而等值面以外區域的感生電流就無法在圖上體現出來了，而且AICD沒法給出感生電流的鍵截面積分值（鍵電流強度）來定量考察電流大小。SYSMOIC則能夠給出三維空間或二維平面上的感生電流向量場圖，從而能展現出AICD圖表現不出來的信息，而且SYSMOIC還能非常方便地給出鍵電流強度。GIMIC程序雖然能做感生電流積分、結合Paraview程序能給出漂亮的向量場圖，但是GIMIC卻無法將感生電流分解成不同分子軌道的貢獻，而SYSMOIC則能夠方便地指定計算感生電流時只考慮或不考慮哪些分子軌道。可見SYSMOIC有獨特的價值。 本文介紹的情況對應SYSMOIC 2024.1版，對其它版本不一定適用，請結合相應版本程序提示和手冊隨機應變。 1 SYSMOIC的基本特征 SYSMOIC是免費而不開源的程序，Windows、Linux、Mac版都有，可以在http://SYSMOIC.chem.unisa.it/DISTRIB下載。諸如STABIN_LNX-2024.1.tar.gz就是Linux的2024.1版可執行文件包。Linux版需要用較新的系統，筆者用Rocky Linux 9可以正常使用，CentOS 7.4用不了。 SYSMOIC是由一堆子程序組成的，對應壓縮包里不同可執行文件，Windows版有.exe后綴，其它平臺的沒有。例如JBMAP子程序用于計算感生電流格點數據，BOCUST子程序用于計算鍵電流強度，等等。v3d子程序用來對其它子程序計算出來的.3d文件進行可視化，如繪制感生電流矢量圖。各個程序的用法在SYSMOIC的在線手冊http://sysmoic.chem.unisa.it/MANUAL/里都有說明。 SYSMOIC是解壓即用的程序。解壓后只需要把SYSMOIC目錄加入到操作系統的PATH環境變量里（不會的話自行google），然后就可以在任何目錄下輸入子程序名直接使用了。 SYSMOIC有不少官方例子文件，是http://sysmoic.chem.unisa.it/WORKED_EXAMPLES/里的WORKEXA.tar.gz。例子按照體系/理論方法/基組對目錄層級進行命名，里面RUN_TESTS文件是對例子體系做各種分析的shell腳本，可以從中看到都是以什么命令調用程序的，以及傳入的命令文件是什么。REFERENCES目錄下的文件是運行RUN_TESTS腳本后產生的.3d文件和程序輸出信息，對前者可以直接用v3d作圖。 SYSMOIC必須結合最主流的量子化學程序Gaussian使用，目前只支持閉殼層形式的HF和DFT。只需要用Gaussian做NMR計算的同時產生wfx文件，再用SYSMOIC自帶的unpackwfx將之轉化成fort.3、fort.11、fort.23、fort.28這四個二進制文件，之后就可以用JBMAP、BOCUST等子程序做各種計算了，算完了可以用v3d子程序直接作圖。 2 產生苯的SYSMOIC輸入文件 下面以最典型的芳香性分子苯為例，講解SYSMOIC繪制感生電流矢量圖和計算鍵電流強度的方法。其它的用法我覺得在一般研究中用到的機會很少，所以就不在本文介紹了，感興趣者請看程序手冊、原文并結合程序自帶的例子文件摸索。下文涉及的各種輸入輸出文件都可以在http://www.shanxitv.org/attach/702/file.rar中獲得。 首先寫一個在B3LYP/6-31G*級別下做NMR計算的Gaussian輸入文件benzene.gjf，內容如下，在本文的文件包里已經提供了。其中的坐標是在當前級別下已經優化過的。NMR=CSGT要求以CSGT方法做NMR計算，out(wfx,CSGT)要求產生.wfx波函數文件，并且把CSGT計算過程中產生的擾動的軌道展開系數寫入到wfx文件中，輸出位置為末尾指定的D:\benzene.wfx。 # b3lyp/6-31g(d) NMR=CSGT out(wfx,CSGT) [空行] benzene [空行] 0 1  C                  0.00000000    1.39621600    0.00000000  C                  1.20915900    0.69810800    0.00000000  C                  1.20915900   -0.69810800    0.00000000  C                  0.00000000   -1.39621600    0.00000000  C                 -1.20915900   -0.69810800    0.00000000  C                 -1.20915900    0.69810800    0.00000000  H                  0.00000000    2.48297800    0.00000000  H                  2.15032200    1.24148900    0.00000000  H                  2.15032200   -1.24148900    0.00000000  H                  0.00000000   -2.48297800    0.00000000  H                 -2.15032200   -1.24148900    0.00000000  H                 -2.15032200    1.24148900    0.00000000 [空行] D:\benzene.wfx [空行] [空行] 算完后就有了benzene.wfx，確保它在當前目錄下，并且已經把SYSMOIC目錄加入到了PATH環境變量中，然后運行unpackwfx benzene命令，當前目錄下立刻就出現了fort.3、fort.11、fort.23、fort.28，它們都已提供在了本文的文件包里。后文說的所有操作都是對于這四個fort文件都處于當前目錄下的情況而言的。 注意當前的苯分子是處在Z=0的XY平面上，用Multiwfn載入Gaussian的out文件然后進主功能0，或者用GaussView載入out文件并顯示笛卡爾軸，都可以確認這一點。后面我們要考察的都是外磁場方向垂直于苯環所導致的感生電流的情況。 3 繪制磁感生電流向量場圖的例子 3.1 默認方式繪制 先來看怎么產生一個最簡單的磁感生電流量場圖，之后再解釋細節。直接在命令行窗口輸入JBMAP運行之，屏幕上會出現好多默認用的參數，比如下面這行，就是告訴你默認的外磁場是(0 0 1)矢量，顯然正垂直于當前苯環，和我們期望的一致 magnetic field...........B    0.0000    0.0000    1.0000 程序問你當前顯示出來的參數是否ok，這里就直接按回車代表用默認的y（yes）。然后再依次問你三個問題： ? external grid points y/n? [n]：即是否讀取外部格點，這里直接按回車代表用默認的n（如果要計算的點不是均勻分布在一個矩形區域，而是比如分布在一個球層表面，就需要輸入y并且提供格點定義文件） ? reference arrow y/n? [n]：含義不明。直接按回車用默認的n ? B direction y/n? [n]：即是否把外磁場箭頭畫在圖上，直接按回車用默認的n。如果你想把外磁場箭頭繪制在圖上就輸入y，并且輸入箭頭的位置和長度 馬上當前目錄下就有了JBM.3d文件。輸入v3d JBM.3d通過v3d程序對這個文件作圖，就蹦出了圖形窗口，如下所示，可以看到在一個矩形區域內均勻分布的每個格點位置都有一個箭頭指示此位置的磁感生電流的方向和大小。 想退出的話在圖形窗口按ESC鍵，或者在文本窗口按ctrl+C。 上面的圖非常丑陋，箭頭特別擁擠，分子結構也看不清楚，完全沒法用。下一節講怎么對JBMAP的計算進行自定義。 3.2 自定義JBMAP的運行參數 JBMAP有兩個方式指定運行參數，比如要把外磁場矢量改為0 1 1（程序自動會做歸一化成為0.0000 0.7071 0.7071），既可以運行JBMAP -B 0 1 1命令，也可以先啟動JBMAP然后直接輸入B 0 1 1后按回車，在屏幕上顯示的參數中都會看到外磁場矢量已經變更成自設的值了。 JBMAP命令后面可以直接接的選項可以輸入JBMAP -h查看，其中幾個關鍵的在這里提一下： -o：輸出的.3d文件名，默認是-o JBM -f：原子球大小的倍數，默認為-f 0.2。想讓作出的圖里原子球大一些就把這個改大 -B：外磁場矢量，默認為-B 0 0 1 -j：感生電流類型，默認是-j TOT，即diatropic和paratropic的總和。寫-j PARA代表只考慮paratropic部分，寫-j DIA代表只考慮diatropic部分 -q：考慮的軌道。例如只想考慮17 20 21三個分子軌道對磁感生電流的貢獻，就寫-q +3 17 20 21，只去掉15 16號分子軌道的貢獻就寫-q -2 15 16。默認是考慮所有分子軌道的貢獻（只有占據軌道才有貢獻） -m：設置SYSMOIC計算感生電流用的CTOCD方法的具體形式，參看http://sysmoic.chem.unisa.it/MANUAL/S4.SS12.html的說明。當基組是完備的時候這些方法結果等價。有些選項是僅對于某些方法才需要設，比如用common gauge-origin的時候需要用-g自定義原點坐標（默認為0 0 0）。默認用的CTOCD-DZ2基本上對于所有情況都是合適的選擇，也不用設額外參數。關于CTOCD不同形式的介紹和差異，可以看前述的SYSMOIC原文以及感生電流綜述文章WIREs Comput Mol Sci, 6, 639 (2016)的Continuous Transformations of the Origin of the Current Density一節。 JBMAP運行后可以修改的參數當中較重要的如下： ? B：外磁場矢量，默認B 0 0 1 ? JTERM：感生電流類型，默認是JTERM 0，即diatropic和paratropic的總和。1和2分別代表paratropic和diatropic ? RI和RF：分別指定格點所均勻分布的矩形盒子的X、Y、Z最小笛卡爾坐標和最大笛卡爾坐標。默認是根據體系坐標自動確定以覆蓋整個體系。對上面的苯分子的例子，從屏幕上默認的參數可見，當前對應的是RI -6.0635 -6.6921 -2.0000和RF 6.0635 6.6921 2.0000 ? FATT：給感生電流乘的系數。比如大部分箭頭在圖上都顯得很不顯著，就可以設大這個值以更便于觀看 ? FMM：磁感生電流過濾范圍。例如FMM 0.01 1就代表忽略感生電流的模小于0.01和大于1的格點（是對乘上FATT之前而言的），免得干擾觀看。默認0 0.1 ? STEP：格點間距，數值越小箭頭越密、計算越耗時。默認0.4 Bohr ? LUNA：0和1分別代表箭頭的長度和面積正比于感生電流的模，默認是后者 ? METH：計算方法，默認的2對應CTOCD-DZ2 為了方便，我推薦以命令行方式運行JBMAP。交互式程序怎么以命令行方式運行，我在《詳談Multiwfn的命令行方式運行和批量運行的方法》（http://www.shanxitv.org/612）里有很詳細說明。下面以這種方式結合自定義參數用JBMAP再次對苯進行計算。創建一個文本文件叫CDMAP.txt，內容如下，每一行對應一條傳給JBMAP程序的命令。此例RI和RF的Z值都是1，X和Y范圍都是-6到6 Bohr，因此格點會分布在苯分子上方1 Bohr的XY平面上的正方形區域內。此文件后面的y、n、n、n對應交互式界面里程序向你提問時依次輸入的四個指令。 B   0 0 1 RI -6 -6 1 RF  6  6 1 FMM 0.01 1 FATT 3 STEP 0.4 LUNA 1 JTERM 0 y n n n 運行JBMAP -f 0.5 < CDMAP.txt，然后再運行v3d JBM.3d，看到下圖，可見清楚地把苯分子上方的磁感生電流描繪了出來。 再來看只考慮pi軌道時的磁感生電流圖的繪制。把benzene.wfx文件載入Multiwfn，進入主功能0，按照《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）說的，從HOMO開始往下挨個看軌道，就能找到三個pi軌道，序號為17、20、21。原子數多的時候也可以用《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）介紹的做法讓Multiwfn自動找出來所有pi軌道序號。運行JBMAP -q +3 17 20 21 -f 0.5 < CDMAP.txt，然后再用v3d繪圖，得到下圖，這便是苯的pi電子的感生電流圖了。相對于之前的所有電子貢獻的圖，此圖少了sigma電子造成的出現在內環的逆時針（paratropic電流）特征，而且也沒有了氫原子區域的感生電流，因為pi電子不在氫上。 3.3 v3d的使用細節 下面專門細談一下v3d界面的操作。v3d圖形窗口出現后，可以用鼠標拖動體系旋轉視角。還可以在圖形窗口處于激活狀態下按h鍵，此時文本窗口會出現可以用的所有按鍵的功能說明。其中較重要的按鍵在這里說一下。 先說視角調節按鍵： a、b、g：視角由alpha、beta、gamma值定義，按a、b、g鍵會分別增加它們1.0。如果按住shift鍵再按，也即按A、B、G，會分別增加它們10 ? f和F：調節縮放因子。按一下f圖像就會放大一點，按一下F就會放大很多 ? d和D：調節視角距離。按一下d就會增加0.1，按一下D就會增加1。視角距離越遠，近大遠小的透視效果越弱、越接近正交視角的效果（據我所知v3d沒法完全用正交視角） ? n和N：設置霧化距離閾值。超過視角一定距離閾值開始就會用霧化效果使得遠處物體被屏蔽。按n增加閾值0.1，按N增加1.0 ? C：令圖像居中 每次按以上按鍵調節設置的時候在文本窗口都會看到當前值。如果你想把這些鍵的操作效果反過來，比如按一次f就令圖像縮小一點，那就先按一下c鍵進入反轉操作狀態（文本窗口里會看到c變成了c-1），然后再按f就行了。再次按c則還原到普通操作狀態。 按p鍵可以把當前圖像保存成當前目錄下的.ps文件。這是矢量圖像格式，在Windows下可以用裝了ghostscript插件的irfanview或photoshop等程序打開，也可以用Acrobat distiller轉成pdf再打開。在Rocky Linux、CentOS等系統里也自帶了ps文件觀看程序。 默認情況下產生的ps文件畫出來的感生電流的箭頭會有較厚的白邊，當箭頭多的時候白邊會導致后面的箭頭被遮擋，比較難看。解決辦法是按一下c鍵先進入反轉操作狀態，然后按住k鍵一會兒使得白邊量減小（文本窗口不會提示減小到了多少），然后再導出ps文件，然后就會發現白邊非常細了。 如果想保存當前的作圖設置，按S，當前目錄下就出現了SpecialValues.txt文件，里面記錄了作圖設置，可以看到f、d等數值都記錄在內。當再次啟動v3d時，程序若發現當前目錄下有這個文件，就會自動載入里面的設置。如果你希望每次產生的ps文件里箭頭的白邊都非常細，可以啟動v3d后先直接按S保存一個記錄默認狀態的SpecialValues.txt文件，然后把里面對應白邊粗細的ridfr前面的值改為0.1。這樣每次啟動v3d就會默認用這個設置了。 4 計算鍵電流強度 對鍵截面做感生電流積分的概念我在《考察分子磁感生電流的程序GIMIC 2.0的使用》（http://www.shanxitv.org/491）中已經專門講過了，這可以定量考察穿越特定化學鍵的電流強度從而便于對比，這在SYSMOIC里叫鍵電流強度。SYSMOIC算鍵電流強度遠比用GIMIC更方便，都省得自己定義積分范圍，結果還能直接可視化。下面我們對苯的一個C-C鍵計算鍵電流強度。 直接在命令行窗口中輸入BOCUST命令運行此程序，程序首先會提示 | plotted    12 atoms | plotted    12 bonds  ====> Pair C1    -  C2     distance=    1.396 Ang 這說明程序發現有12個原子對兒的距離都處在默認的距離下限（0.9埃）和距離上限（1.6埃）之內。其中C1-C2是第一個原子對，現在問你是否對這個C1-C2鍵計算鍵電流強度。這里輸入y然后按回車，很快就會出現結果 CONTRIBUTION  1 -0.6000526E+00 AU -0.1690912E-07 SI CONTRIBUTION  3  0.1784109E+00 AU  0.5027513E-08 SI TOTAL CS ==>    -0.4216417E+00 AU -0.1188161E-07 SI (2) CS/CS0 RELBEN FICS     0.99    3.00  0 上面顯示的兩個CONTRIBUTION分別是diatropic和paratropic對鍵電流強度的貢獻，由于方向相反，顯然符號相反。二者之和取絕對值0.4216417 a.u.就是鍵電流強度大小，換算成常用的nA/T單位就乘以28.179409，即11.88 nA/T。這個值與《考察分子磁感生電流的程序GIMIC 2.0的使用》（http://www.shanxitv.org/491）在B3LYP/6-311G*級別通過GIMIC程序得到的12.086 nA/T非常接近。上面信息中的CS/CS0 RELBEN FICS     0.99代表當前值與參考值的比值是0.99。這里參考值是程序內置的對苯在某個級別計算的C-C鍵電流強度12 nA/T。 接下來程序又問你是否對C1-C6計算鍵電流強度。由于當前所有C-C鍵都是等價的，只需要計算一個就行了，因此輸入N回車，程序就結束運行了。如果你想把所有鍵都進行計算，就接著輸入11次y回車把剩下11個鍵都算完。當前體系一共6個C-C鍵、6個C-H鍵。 當前目錄下出現了BCS.3d文件。運行v3d BCS.3d，就可以看到下圖 此圖里黑色箭頭體現了C1-C2鍵電流的大小和方向。積分截面處在鍵的正中央并與鍵垂直，圖中的曲線是這個截面上的感生電流0.001 a.u.等值線，曲線內部區域是對感生電流積分的區域，紅色和藍色分別對應paratropic和diatropic部分。可見diatropic分布在環的外側，這與之前從向量場圖中看到的順時針環電流的主要分布是一致的。圖中的紅色小點對應截面上感生電流的極值點。鍵上面標的文字99對應于前述的參考比值0.99乘以100。文字出現的位置在鍵的正中央，不容易分辨，可以在圖形窗口激活狀態下按W鍵關閉文字顯示，然后再自行ps上。 如果你想一次性把所有鍵的電流強度都計算，又不想一次次麻煩地輸入y，那么可以寫一個文本文件比如叫ally.txt，里面有12行，每一行都是y。運行BOCUST < ally.txt，這樣所有鍵就都會被計算了。此時對得到的BCS.3d繪制會看到下圖 下面介紹一下BOCUST運行時可以接的選項中比較重要的： -o、-f、-B、-m、-q等：和前面介紹的JBMAP程序的情況完全相同。即BOCUST也可以用-B指定外磁場方向、用-q指定考慮哪些軌道，等等。 -d、-D：確定哪些原子對兒被考慮時用的距離下限和上限（埃），默認-d 0.9 -D 1.6 -e：確定鍵截面上做積分的區域用的感生電流等值線的數值（a.u.），默認-e 0.001 -CS0：參考電流強度（nA/T），默認-CS0 12 -nocd：不顯示積分區域和極值點。如果想讓圖上只有箭頭而沒有曲線和小點，應該加上這個 -l [比值]：在截面上繪制感生電流的不同數值的等值線，后面跟的是相鄰等值線數值的比值 -C [粗細]：繪制截面的邊框 例如這里計算C1-C2的sigma電子的鍵電流強度（即扣除三個pi軌道貢獻），并在截面上顯示等值線、顯示邊框，就運行BOCUST -q -3 17 20 21 -l 3 -C 5，然后輸入y回車，再N回車。之后用v3d對BCS.3d作圖看到下圖。可見整個截面，以及里面的感生電流等值線，都顯示得很清楚。由于這個截面上diatropic和paratropic的大小都差不多（畢竟苯分子沒有sigma芳香性特征），所以鍵電流強度僅為0.007 a.u.。 5 總結 此文對分析磁感生電流的程序SYSMOIC做了介紹并給出了具體例子。由例子可見，通過SYSMOIC結合Gaussian繪制感生電流向量場圖很方便，不需要借助第三方程序。SYSMOIC做鍵電流強度計算也很自動化，不需要像GIMIC那樣謹慎地定義積分區域，而且鍵電流大小和方向還能直接通過箭頭直觀顯示出來。SYSMOIC的計算速度也比較快。 本文只涉及了SYSMOIC最常用的功能，還另外有一些其它功能，比如對感生電流密度做拓撲分析，計算各種感生電流密度相關的場（ACID及其各向異性、磁屏蔽密度、vorticity張量的跡及其各向異性等）并繪制成等值面和填色平面圖、計算任意格點上的感生電流張量和導數等等。 http://www.shanxitv.org/699 2024-02-09T17:09:00+08:00 談談有小虛頻時熱力學量的計算 文/Sobereva@北京科音  2024-Feb-9 0 前言 量子化學研究中在優化完極小點結構然后做振動分析時遇到虛頻是家常便飯的煩人的問題。之前我專門寫過博文《Gaussian中幾何優化收斂后Freq時出現NO或虛頻的原因和解決方法》（http://www.shanxitv.org/278）講了怎么消虛頻。除非程序存在bug，否則原理上不存在按上文處理完全解決不掉的虛頻。一個現實問題是：雖然靠各種嘗試、折騰，最終確實總能消掉虛頻，但這個過程對于大體系來說過于浪費計算時間，也耗費研究者的精力。那么是否帶有虛頻的結構就一定不能拿來算能量相關的問題呢？是否虛頻總是要不惜一切代價消除呢？此文就談談我個人的看法，給出理論分析，并且通過一個例子對我的觀點予以驗證。 本文牽扯到熱力學數據計算的各種常識知識，沒有這些常識的話先看《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》（http://www.shanxitv.org/552），以及Shermo程序手冊的附錄，里面有計算原理的簡潔易懂的概述。在北京科音基礎（中級）量子化學培訓班（http://www.keinsci.com/workshop/KBQC_content.html）里對熱力學量的計算有非常全面深入的介紹，歡迎關注和參加。 本文說的小虛頻和過渡態必然存在的虛頻沒有絲毫聯系，初學者絕對不要把本文內容胡亂套用到過渡態的情況。 1 導致小虛頻的原因 為了避免誤解，首先回顧一下優化極小點后振動分析時出現虛頻的各種常見原因： (1)振動分析時用的影響勢能面的設置與優化時的不同。這通常是沒有基本常識的外行人才會犯的低級錯誤。注：JCTC, 17, 1701 (2021)中提出的single-point Hessian (SPH)方法通過對勢能面添加偏移勢，從而允許振動分析用的計算級別低于優化用的級別來節約時間，但這個做法非主流、經驗性強，本文不涉及 (2)初始結構的點群高于實際極小點的點群，且如果優化后的結構維持住了初始的對稱性，振動分析后就會發現破壞對稱性的虛頻，即按照虛頻方式移動結構會導致結構的對稱性下降。有時候這種有虛頻的結構恰對應于過渡態 (3)程序bug。比如Gaussian 09 D.01的DFT-D3色散校正對解析Hessian的貢獻的計算有bug，容易導致算勢能面平緩體系出現虛頻。這點在《DFT-D色散校正的使用》（http://www.shanxitv.org/210）我專門做了說明 (4)某些方法對勢能面引入明顯數值噪音。如Gaussian中的SMD溶劑模型的這個問題就比較明顯 (5)幾何優化收斂得不精確。幾何優化收斂限過于寬松或Hessian的質量不太好會導致這種問題。 (6)積分格點質量不夠好。這個問題對于M06-2X等對積分格點敏感的泛函尤為顯著。關于量子化學程序用的原子中心的積分格點參考《密度泛函計算中的格點積分方法》（http://www.shanxitv.org/69）。如果是CP2K等利用平面波的第一性原理程序，用的均勻分布的格點依賴于平面波截斷能的設置。 以上前四種情況導致的虛頻是一定要解決掉的，而且本來解決就非常簡單。只要不犯低級錯誤自然就沒有(1)的情況出現，解決(2)按照虛頻模式調節結構后進一步優化即可，解決(3)就換成沒bug的版本或其它計算程序，解決(4)可以換成其它的能在優化+振動分析中實現基本相同目的而沒嚴重數值噪音的方法，如SMD換成IEFPCM，或者用SMD時結合SAS方式定義溶質孔洞表面。對于(5)和(6)情況帶來的虛頻，只要計算能力和精力允許，能消盡量要消，但如果由于計算資源或精力所限實在太難消掉、被迫要容忍虛頻的存在，那么后續的能量相關問題的計算就是后文討論的內容了。 對于柔性大分子、弱相互作用復合物這樣勢能面的某些維度非常平緩的體系，在做完幾何優化和振動分析后很容易發現存在大小非常小的虛頻，如波數為9.8i、27.3i cm^-1這種，普遍小于50，這一般是前述(5)或(6)情況所導致的。此類虛頻的振動模式通常對應于甲基旋轉運動、柔性骨架的變形運動、分子間相對運動。這類(5)或(6)情況所導致的虛頻以下簡稱為“小虛頻”。但要注意也不是所有波數小于50的虛頻都屬于這種情況，前述的其它四種情況也完全可能導致數值很小的虛頻，所以絕對別光拿虛頻大小來判斷情況。 2 有小虛頻時能量相關問題的計算 當遇到小虛頻時，顯然最完美的做法（具體選擇視情況而定）是用諸如更嚴格的幾何優化收斂限、更好的Hessian、更好的積分格點精度進一步做優化和振動分析，使得虛頻完全消除后再計算能量相關問題。但是對于普通泛函算幾百個原子有機體系，很可能算一次Hessian就要花一整天時間，幾何優化一步也可能花費一個小時左右，最終把虛頻完全消掉可能得花幾天甚至一個禮拜，十分痛苦，甚至不可能實現。那么，在后續算能量相關問題時怎么考慮？下面根據被計算的量專門討論。 2.1 計算電子能量 諸如Gaussian這樣的主流量子化學程序默認的幾何優化收斂限是以受力和位移作為判據的。在默認收斂限下優化正常完成后，即便結構還沒有與實際極小點特別精確一致，但通常相差很小。而且當前的最大和平均受力肯定已經很小了（因為收斂限考慮了它們），由于受力對應于能量梯度的負矢量，所以當前結構稍微再挪一點達到精確的極小點的過程中能量變化肯定非常小。即曰，在有小虛頻的情況下，計算的電子能量一般是可以用的，也即結果和消掉虛頻后再算的能量相差極小。 2.2 計算0 K下的內能U(0) U(0)是電子能量與零點振動能（ZPE）之和。在諧振近似模型下，每個振動模式對ZPE的貢獻是1/2*h*ν，h是普朗克常數，ν是振動頻率。頻率很低的振動模式對ZPE的貢獻是非常小的，每波數的振動頻率對ZPE的貢獻僅為5.98 J/mol，因此諸如30波數的振動模式只對ZPE有0.18 kJ/mol（0.04 kcal/mol）的貢獻，數量級遠遠低于理論方法、基組、溶劑模型等因素帶來的誤差。順帶一提，Shermo計算熱力學量的時候如果把settings.ini里的prtvib設為1，每個振動模式對熱力學量的貢獻都會明確輸出。 主流量子化學程序，以及Shermo程序在默認情況下，在計算振動對熱力學量的貢獻時都是不考慮虛頻產生的貢獻的，因此小虛頻對熱力學量的貢獻會被直接忽略掉，這不會對ZPE的計算造成實際問題。因為如果把這個小虛頻消了，從而產生一個實頻，其頻率肯定也是非常低的，對ZPE的貢獻可忽略不計。也就是不管消不消這個虛頻，算出的ZPE都沒多大差別。不過，實際情況遠沒這么簡單，因為在虛頻消掉前后，由于二者所在勢能面的位置的些許不同，或改用不同檔次積分格點導致頻率計算精度的差異，其它的實頻模式也會多多少少受到影響，不過頻率值受影響相對明顯的主要是低頻模式（可能影響幾個甚至10個左右波數，而高頻模式受影響一般小于2波數），它們對ZPE的貢獻恰恰又很小，所以虛頻消掉前后算的ZPE雖然未必差異小到諸如前面說的0.2 kJ/mol左右這種程度，但至少也不會大，比如對柔性大體系頂多也就影響個1-2 kJ/mol。 如上所述，鑒于是否小虛頻對電子能量和ZPE的影響都不太明顯，因此算電子能量和U(0)時一般可以容忍一個或多個小虛頻的存在。 2.3 計算從0 K升溫到當前溫度對內能的影響 研究U(0)沒太大實際意義，畢竟實際化學過程也不可能在0 K發生。計算有限溫度，也就是高于0 K時候的熱力學量，就需要考慮從0 K升溫到當前溫度對內能的改變量，其中振動產生的貢獻下面用Uvib(0→T)表示。首先要知道一點，低頻對Uvib(0→T)的貢獻是明顯高于高頻模式的。這點從下面的J. Phys. Chem., 100, 16502 (1996)的圖1可以明顯看出來，圖中ΔHvib(T)和這里說的Uvib(0→T)是一碼事，圖是對T=298.15K的情況繪制的。   雖然低頻模式的貢獻大，但是由圖中放大的0-100波數的區間可見，在小幾十波數范圍以內，頻率值對Uvib(0→T)的影響雖然不可忽視，但不算多大，都在2.0-2.4 kJ/mol范圍內。這也就是說，由于存在一個小虛頻導致少考慮一個實低頻模式，會給Uvib(0→298.15K)帶來-2.2 kJ/mol左右的誤差，這雖然不算小，但還算可以湊合接受。如前所述消虛頻后還可能會對其它低頻模式產生不可忽視的影響，這對Uvib(0→T)的改變一般也不至于太大。 綜上所述，若帶著小虛頻計算U(T)，會造成一定低估。對于有一個小虛頻的情況，低估程度在0.5 kcal/mol左右。如果有多個小虛頻，那還是別計算U(T)了，除非T遠低于常溫，或者結合后文說的把小虛頻當成實頻的做法。 2.4 計算熵、自由能 低頻對熵的貢獻遠大于高頻，而且在很低頻區間內，頻率越低貢獻明顯越大。在包括Gaussian在內，大多數量子化學程序（ORCA除外）計算熱力學量用的是傳統的rigid-rotor harmonic oscillator (RRHO)模型，振動熵的計算完全基于諧振近似模型。此時在低頻區間內，隨著頻率趨近于0，振動模式對熵的貢獻猛增，如下面的Chem. Eur. J., 18, 9955 (2012)文中圖2的虛線所示。Grimme在這篇文章中還提出了一種quasi-RRHO模型（QRRHO，后來也叫modified RRHO、mRRHO），以下稱為QRRHO(Grimme)，它將自由轉子和RRHO用的諧振近似模型算的振動熵做插值，頻率越接近0自由轉子的權重越大、頻率越高諧振近似的權重越大，如下圖實線所示。QRRHO(Grimme)對高頻模式算的熵貢獻和RRHO完全一樣，而對于低頻模式（特別是波數在100以下的）考慮得比RRHO明顯更合理。由下圖還可見QRRHO(Grimme)算的低頻模式的熵比RRHO整體小得多，沒有RRHO的嚴重高估的問題。此外，QRRHO(Grimme)下低頻模式的振動熵對頻率值的敏感性低于RRHO很多。   直接計算熵的用處不大，算熵的主要目的是用來算自由能，因為自由能里有-T*S項。顯然低頻對熵的貢獻計算的準確度顯著影響到自由能的計算準確度。假設體系有一個小虛頻，并且在消掉之后會多一個20 cm^-1的實頻，并且用的是QRRHO(Grimme)算298.15 K的熵，那么在有虛頻的結構下由于缺失了這個實頻，根據上圖可知會造成熵的誤差約為-18 J/mol/K，對自由能造成的誤差約為-298.15*(-18)/1000=5.4 kJ/mol，略高于1 kcal/mol，這就不小了，都能趕上高質量ωB97M-V泛函結合大基組算普通有機反應的誤差了。不過實際中小虛頻帶來的誤差未必那么大，因為有和沒有小虛頻的結構下實低頻模式的頻率也往往有不可忽視的不同，例如下面第3節的例子在有虛頻的結構下實低頻模式的頻率比精確極小點結構下的整體偏小（勢能面形狀原因所致），導致高估了不少實低頻模式的熵，藉由-T*S項進而對自由能造成低估，于是和少考慮一個低實頻對自由能的高估產生了很大程度的相互抵消，下面管這稱為“熵的誤差抵消”。 根據以上討論，有小虛頻的情況下計算熵和自由能都有一定風險。如果就有一個小虛頻，算常溫的情況往往還說得過去，而如果有多個小虛頻，或計算高溫的情況（注意T越大-T*S越大，熵的誤差造成自由能的誤差越大），我就不建議冒險帶著虛頻了。而且就算非要帶著小虛頻算，也至少應當用QRRHO模型，而切勿用RRHO模型，因為RRHO算的振動熵對低頻頻率值的敏感程度顯著高于QRRHO，因而用有虛頻的結構算的熵的誤差上限比QRRHO的明顯要大。 要注意QRRHO模型不止Grimme提出的這一種。如Shermo手冊附錄部分所介紹的，Truhlar的QRRHO是把所有小于一定閾值（通常為100 cm^-1，可以由Shermo的ravib參數控制）的低頻統一提升為這個閾值再算它們對Uvib(0→T)和熵的貢獻，這個做法的物理意義不甚明確，沒有QRRHO(Grimme)的插值做法優雅。目前還有一種QRRHO是Minenkov等人在J. Comput. Chem., 44, 1807 (2023)中提出的，它對振動熵部分的考慮和QRRHO(Grimme)一樣，但還同時將這種自由轉子與諧振近似做插值的思想同時應用到了計算振動對內能的熱校正量Ucorr上（此時沒法單獨計算ZPE），原理更理想，實際結果還略好一點點。這些QRRHO模型在最新版Shermo中都是支持的，通過ilowfreq參數控制。 2.5 將小虛頻視為實頻的做法 在Angew. Chem. Int. Ed. 2022, e202205735中的3.2節Grimme等人認為存在小虛頻的情況下，只要使用QRRHO模型并且將小虛頻視為實頻處理，比如18i cm^-1直接當做18 cm^-1，則得到的熵就是可以用的、小虛頻的存在就不再是個問題。這種做法表面看起來似乎有點莫名其妙、太任性了，但確實有一定道理，可以這么理解：在消除小虛頻后，必定會得到相同數目的波數也很小的實頻模式，它們都對熵有明顯貢獻。相對于完全忽略掉這些小虛頻的貢獻而導致明顯低估熵，干脆把他們當做實頻處理，從而“湊”出來相應數目的實低頻是更好的做法，而且本身在QRRHO模型下，數值在十幾、幾十波數區間的頻率對熵的貢獻也都相差不懸殊，因此也用不著對這么人為湊的實頻頻率太較真。由于如前所述，Uvib(0→T)也對低頻數值不敏感而低頻的貢獻又沒小到可以隨便忽略的程度，所以計算Uvib(0→298.15K)時也使用這種處理是有道理的。從Shermo 2.5版開始，可以通過settings.ini里的imagfreq選項設置視為實頻的虛頻閾值，比如設為50，就代表大小在50 cm^-1以內的虛頻都當做實頻處理來計算振動對各種熱力學量的貢獻，此時從Shermo程序輸出的頻率信息中也會看到已經沒有虛頻了。 然而，這種將小虛頻視為實頻的做法目前很非主流，其思想在我來看過于主觀和樂觀、把問題想得太簡單了，目前也缺乏系統性的檢驗，我個人不建議輕易使用。這種做法有一個明顯問題：如前所述，有小虛頻和消掉虛頻的兩種情況相差的不僅僅是一個小虛頻轉變成一個小實頻，許多其它低頻模式的頻率也會有不小變化，而且有時是整體明顯增加的。不做這個虛頻到實頻的轉換的情況下能夠有上一節說的熵的誤差抵消效應，而這么做了之后，忽略虛頻模式導致對熵的低估問題確實很大程度解決了，但在有小虛頻結構下其它頻率往往偏低而導致熵的高估帶來的誤差就沒處抵消了。此外，如果當前體系的低頻還有很多都是10 cm^-1 左右這種特別低的，這個范圍附近即便是QRRHO(Grimme)或QRRHO(Minenkov)模型算的熵也對頻率很敏感，問題會更大。 所以，將小虛頻視為實頻的做法雖然確實對一些情況更好地計算熵和自由能是有幫助的，但也很可能起到負面效果。所以這絕不是普適、穩健、理想的做法，更千萬不能一看到存在小虛頻就總是想靠這個做法來完全無視虛頻。 3 有虛頻時算熱力學量的實例 筆者最近通過Gaussian 16在ωB97XD/6-311G*級別研究的一個以pi-pi堆積方式形成的共128個原子的三分子復合物，在優化和振動分析后發現存在6.44i cm^-1的小虛頻，之后進一步優化時用opt=calcfc關鍵詞提供精確的初始Hessian，優化完再做振動分析就發現沒虛頻了，此時最低的頻率成了10.35 cm^-1。正好這個例子可以拿來考察消虛頻前后計算的熱力學量的差異，由此檢驗帶著那個小虛頻時以不同方式算的熱力學量有多大誤差。要強調的是，僅僅這么一個測試，顯然不足以充分證明某種做法理想，因為這必須通過大量體系做統計分析，但至少誤差較大的情況可以用來體現相應做法不夠穩健和普適。 先看一下有虛頻相對于無虛頻時各個實頻的頻率變化情況。先把實頻頻率按照由低到高排列，用有虛頻時候它們的頻率減去無虛頻時候的頻率作為縱坐標值，而兩種情況頻率的平均值作為橫坐標，作的圖如下所示。可見，至少對于此例，有虛頻時數值較低的實頻是普遍明顯偏低不少的，而>200 cm^-1的相對而言的中、高頻則受到的影響很小。這種情形雖然不能說是普遍情況，但至少也是挺常見的一種。 此例有虛頻相對于無虛頻的情況計算的電子能量僅僅相差0.04 kcal/mol，小到完全可以忽略不計。用Shermo程序基于Gaussian振動分析輸出文件計算的有虛頻相對于無虛頻時的ZPE和氣態標況下的熵、內能、自由能熱校正量的差異如以下表格所示。RRHO以及Truhlar、Grimme、Minenkov三種QRRHO分別對應Shermo的ilowfreq=0/1/2/3的情況，imagreal=0和50相當于忽略那個虛頻以及把那個虛頻當大小相同的實頻考慮的情況。 首先看常規的忽略虛頻貢獻的情況。由表格可見： ? ZPE：有虛頻時和無虛頻時ZPE相差僅-0.3 kcal/mol，屬于通常可以接受的差異 ? 熵：由于誤差抵消的巧合，恰好RRHO下有小虛頻時的熵和無虛頻時的差異很小。QRRHO(Truhlar)的情況差異挺大，在于它把實低頻都統一當成了100 cm^-1，這部分在有和沒有小虛頻情況之間幾乎都抵消了，而有小虛頻時相當于少考慮了一個100 cm^-1振動模式對熵的貢獻，因此造成熵的嚴重低估。QRRHO(Grimme)和QRRHO(Minenkov)對熵的部分處理相同，由于前述的熵的誤差抵消效應導致有和沒有小虛頻時熵的差異不算太大 ? 內能的熱校正量Ucorr：等于ZPE+Uvib(0→T)。由于Uvib(0→T)對低頻的具體數值不太敏感，有無小虛頻情況下實低頻的貢獻變化很小，而有小虛頻時由于缺失了一個低實頻對Uvib(0→T)的貢獻，導致有小虛頻時Uvib(0→T)偏低不少，Ucorr因此比ZPE偏低得明顯更多。QRRHO(Minenkov)用的計算Ucorr的方式和其它模型不同，對于此例使得有小虛頻時Ucorr偏低程度比其它模型更小，僅為-0.5 kcal/mol。因此在計算有小虛頻情況下的U(T)時，用QRRHO(Minenkov)可能會比其它模型更有優勢。 ? 自由能的熱校正量Gcorr：相對于其它模型，有小虛頻相對于無虛頻時Gcorr的差異在RRHO下是最大的，這在于有無小虛頻時RRHO算的S的差異對此例恰好頗小，導致-T*S項的差異才0.03 kcal/mol，因此沒能充分對Ucorr的差異產生抵消效應。而各種QRRHO模型下Gcorr的差異明顯較小，這在于這種情況下S的差異不小，而且-T*S項的差異的符號和Ucorr的差異相反，因此相互抵消了很多。其中Gcorr差異最小的是QRRHO(Minenkov)，才-0.11 kcal/mol，體現出它可能適合作為有小虛頻情況下算自由能的優先選擇。 再來看將小虛頻當實頻處理的情況。根據以上表格可見，這種處理使得Ucorr在有和沒有小虛頻情況下計算的結果差異減小很多，因此對于算內能的目的，這種做法有一定意義。但是，除了QRRHO(Truhlar)外，這個做法卻造成了有小虛頻時算的S和Gcorr與無虛頻時的差異變得顯著更大，對RRHO最為嚴重，對QRRHO(Grimme)和QRRHO(Minenkov)問題也挺嚴重，這來自于有小虛頻時實頻模式偏低造成的對振動熵的高估無法與因為少考慮一個實低頻導致的低估一定程度相抵消。把小虛頻當實頻處理時結合QRRHO(Truhlar)對此例倒是很不錯，有小虛頻相對于無虛頻時的熵只相差0.27 cal/mol/K，這在于此種QRRHO將低頻全都上拉到100 cm^-1，此時又把小虛頻轉化成了小低頻，因此無虛頻和有虛頻時相當于有同等數目的100 cm^-1的實頻，自然兩種情況下算的熵很接近。 4 總結 根據以上討論和實際體系的測試可知，對于因(5)和(6)情況造成有小虛頻的情況，算電子能量和U(0)沒大問題，不是非得消虛頻。 對于計算Ucorr或與之相差RT的焓的熱校正量Hcorr，用QRRHO(Minenkov)或QRRHO(Grimme)并同時把小虛頻當實頻處理是較好選擇，不會因為存在小虛頻的問題而帶來太大誤差。 有小虛頻時最難搞的是熵，以及溫度不很低時候Gcorr的計算，難點在于小虛頻消掉前后實低頻模式頻率也可能受到不小影響且影響情況不易估計，然而即便用QRRHO(Grimme或Minenkov)時低頻的具體數值對熵的影響也不算太小，所以明顯不是把小虛頻當成實頻處理以避免少考慮一個實低頻模式就能較好解決的，而且如第3節的測試可見，這種做法還可能造成熵和Gcorr在有小虛頻時計算結果更差。所以當你對熵、Gcorr要求精度高時，小虛頻要盡量消。實在搞不定的話，那就姑且用原理上最好的QRRHO(Minenkov)，并且建議不用小虛頻當實頻的處理，畢竟此做法非主流而且起到正面效果的可能性很有限。雖然把小虛頻當實頻并結合QRRHO(Truhlar)的情況下小虛頻和無虛頻時熵的差異很小，但在J, Comput. Chem., 44, 1807 (2023)的測試中QRRHO(Truhlar)算團簇反應自由能表現得遠不如QRRHO(Minenkov)，甚至比RRHO沒明顯改進，所以我也不推薦。 http://www.shanxitv.org/698 2024-01-29T10:07:00+08:00 談談如何計算環張力能：以CPP和碳單環體系為例 文/Sobereva@北京科音   2024-Jan-29 0 前言 環狀化學體系具有不同程度的環張力，環張力和環狀體系的諸多問題緊密相關，諸如生成焓、穩定性、合成難易程度等。環張力能（strain energy, SE）是衡量環張力大小的最直接的指標，它對應當前環狀結構轉變為假想的沒有環張力的結構過程中能量變化的負值。 北京科音自然科學研究中心（www.keinsci.com）的盧天等人之前受邀發表了一篇專門精確考察碳環（cyclocarbon）體系環張力能的文章，文中對環張力的計算思想有詳細介紹并給出了諸多討論： Accurate theoretical evaluation of strain energy of all-carboatomic ring (cyclo[2n]carbon), boron nitride ring, and cyclic polyacetylene, Chin. Phys. B, 31, 126101 (2022) DOI: 10.1088/1674-1056/ac873a（不方便獲取者也可以看預印版https://doi.org/10.26434/chemrxiv-2022-v8w9h，內容基本一樣） 下面，將基于此論文中的內容簡要說明環張力能的計算原理并給出具體例子。感興趣的讀者請務必仔細閱讀論文原文以了解更多信息，也推薦以類似方式計算其它體系的環張力能時引用此論文。另外，讀者如果對碳環體系感興趣，很建議查看同作者之前發表的此類體系的各種研究文章和相關博文：http://www.shanxitv.org/carbon_ring.html。 1 同聯反應法 設計同聯（homodesmotic）反應是計算環張力能的很常用方法。這個方法設計一個化學反應，既使得環被打開使得環張力被完全釋放，同時又盡可能不影響原本的成鍵情況，顯然這個反應能的負值就對應了環張力能。 例如，要計算氧雜環丁烷的環張力能，就可以設計以下同聯反應。紅色部分對應具有顯著環張力的氧雜環丁烷，將它插入到甲乙醚所示位置就變成了右邊的產物，顯然此時環張力就被完全釋放掉了。與此同時，反應物和產物的成鍵特征并未發生變化，例如下圖紅字的第一個碳在反應前后始終是sp3雜化狀態、始終連著一個氧原子和一個CH2，因此化學環境變化導致的成鍵特征的極輕微改變對這個反應能的影響可以忽略不計，也即這個反應能幾乎只體現出環張力的釋放。 同聯反應法計算環張力能的缺點是反應的設計并不唯一，比如以上反應的甲乙醚也可以是甲醚、乙醚等等，結果或多或少會有點差異。顯然反應應當設計得盡可能令各個鍵的特征在反應過程中不發生改變，但這并非總能很好實現，后文提到的碳環體系就是一例。 2 同聯反應法計算實例：[6]CPP 這一節示例使用同聯反應法計算[6]CPP的環張力能。此體系的幾何結構，以及設計的同聯反應如下所示。可見是把[6]CPP插入到了聯苯中間變成了直線的八聯苯，從而完全釋放了環張力。 上圖涉及的三個分子都在ωB97XD/6-31G*級別下用Gaussian 16程序進行了幾何優化，并做振動分析確認了無虛頻，輸入輸出文件都在http://www.shanxitv.org/attach/698/file.rar里。基于輸出文件里最后一次輸出的電子能量，用反應物能量之和減去產物能量，就得到了環張力能：627.51*(-463.1444911-1385.7232492+1849.0260518) = 99.3 kcal/mol。這個值比C-C單鍵的鍵能還要大一些，可見其環張力是相當強的。 如果想得到更精確結果，可以基于優化完的結構在更好的級別下算高精度單點再求差，比如基組可以改用明顯更好的def2-TZVP。這對結果可能造成幾kcal/mol程度的影響。 3 外推法 如果被考察的環狀體系是由重復單元構成的，比如前面的[6]CPP是苯環作為重復單元，就也可以用這一節介紹的外推法計算其環張力能。這個方法計算過程如下 (1)計算從小到大不同重復單元數（n）的環的每單元的能量E/n，這里E是整體能量。顯然E/n是依賴于n的，n較小時環尺寸較小，自然環張力大，平均到每個單元的環張力能也大 (2)用Origin之類程序通過恰當的函數（一般是二次函數）擬合E/n與n之間的解析關系，并外推得到n無窮大時的E/n，這對應于沒有環張力時候的每單元能量 (3)將當前考察的環狀體系的E/n減去n無窮大時的E/n，就得到了當前體系的每單元環張力能（SE/n） (4)將當前體系的SE/n乘以當前的體系的n即是當前體系的環張力能 相較于同聯反應法，外推法的優點： ?原理非常嚴格 ?沒有同聯反應法設計反應的任意性，并可以用于同聯反應法不適合的體系 ?可以順帶得到環張力能與n的解析關系，由此可以直接預測出任意重復單元數的環張力能 缺點： ?需要計算從小到大不同尺寸的環狀體系，總耗時明顯高于同聯反應法 ?計算步驟略繁瑣，得自己做擬合 如果你只是要考察當前一個環狀體系的環張力能，而且可以設計出很合適的同聯反應，那么就沒必要用昂貴且略麻煩的外推法。 4 外推法計算實例：18碳環 2019年首次在凝聚相實驗中觀測到的18碳環分子由18個sp雜化的碳依次相連構成環狀。由于sp雜化的碳明顯傾向于形成180度鍵角，因此18碳環必然有顯著的環張力，具體環張力能是多少非常值得研究，這正是前述Chin. Phys. B, 31, 126101 (2022)文章研究的主要內容。 此文中指出碳環類體系不適合同聯反應法計算環張力能。碳環類體系是長、短C-C鍵交替構成，兩個碳是一個重復單元，若用同聯反應法計算化學式為C_2n的碳環的環張力能，可以設計如下反應式，相當于把碳環的每個重復單元插入到了一定長度的碳炔里。碳炔是sp雜化的碳形成的鏈狀分子，兩端由氫封端。PS：對碳炔感興趣者建議閱讀《氫封端碳鏈H-(C≡C)n-H (n = 3-9, 15)的電子光譜的尺寸依賴性：性質分析及對碳炔的預測》（http://www.shanxitv.org/679） 然而，如下圖可見，將18碳環的重復單元插入到不同長度（通過m體現）的碳炔里，按照同聯反應算出來的環張力能有明顯區別，因此靠同聯反應法難以得到嚴格的18碳環的環張力能。 碳炔和18碳環都具有C-C鍵長、短交替的特征。文中還發現，不管碳炔的m取多少，處于它中央的較長的C-C鍵和18碳環中的較長的C-C鍵之間的鍵長差異都不可忽視，對較短的C-C鍵也是如此。因此原理上沒法設計一個同聯反應使得18碳環的C-C鍵特征在轉移到產物時幾乎完全不發生改變，這也必然導致靠同聯反應法不可能對此體系得到完全精確的環張力能，也即成鍵特征差異造成的能量改變會不可避免地摻入到計算出的環張力能中。 由于發現了以上問題，文中使用外推法計算碳環的環張力能。為了追求盡可能好的精度，文中用ORCA程序在精度很好的DLPNO-CCSD(T)/cc-pVTZ級別下對n=7到n=24的碳環都計算了單點能（基于ωB97XD/def2-TZVP優化的結構），ORCA輸出文件在http://www.shanxitv.org/attach/699/file.rar里都提供了。每重復單元能量（E/n）與n的關系如下所示。由于n較小的碳環的電子結構特征有特殊性、和較大的碳環存在一定差異，因此文中只對n>=13的數據點做了二次函數的擬合，由下圖的紅色曲線可見擬合質量超極好，R^2近乎精確為1。 令擬合公式E/n = 0.8844/n^2 - 76.00997中的n為無窮大，就得到了無環張力對應的無窮大的碳環的每單元能量-76.00997 Hartree。每單元環張力能SE/n = E/n - 76.00997，于是有SE/n = 0.8844/n^2，也即碳環類體系的環張力能SE = 0.8844/n Hartree。18碳環的n=9，代進去就得到了其環張力能0.8844/9 Hartree = 61.7 kcal/mol。 值得一提的是，如果令前述同聯反應法中碳炔的m設得很大，從而使得環張力能隨m收斂，最終算出來的結果是59 kcal/mol，可見和嚴格的外推法的結果相比誤差并不可忽略。 前述論文中也測試了用ωB97XD/def2-TZVP的能量結合外推法算碳環的環張力能，結果和上面用的昂貴的DLPNO-CCSD(T)/cc-pVTZ能量的結果沒多大差別。這在一定程度上體現出ωB97XD/def2-TZVP級別算碳環類體系很不錯。 文中還基于不同尺寸碳環的環張力能和鍵角之間的關系，推導出了適合碳環體系用的C-C-C鍵角力常數56.23 kcal/mol/rad^2。 此外，文中還對以下所示的18碳環的等電子體硼-氮環，以及具有單套pi共軛特征的環聚乙炔的環張力做了計算，發現硼-氮環類體系的環張力明顯小于碳環類，而環聚乙炔類體系的環張力則略大于碳環類。文中利用Multiwfn程序（http://www.shanxitv.org/multiwfn）做波函數分析從電子結構角度清晰解釋了它們的環張力與碳環之間存在差異的原因。文中還發現用常用的B3LYP泛函算全局pi共軛的環狀聚合物的合理性較差，E/n隨n的變化明顯不合理。這些方面這里就不多提了，請讀者務必閱讀論文原文。 5 總結 環張力是環狀分子的重要特性。本文簡明扼要地介紹了兩種最常用的計算環狀分子的環張力能的方法，并以[6]CPP和18碳環作為具體例子進行了演示。同聯反應法比較方便省事，但對于一些全局pi共軛的環狀體系，諸如18碳環，則更適合外推法，明顯嚴格得多，同時外推法還很有助于探究環尺寸與環張力之間的關系。還應當注意并不是所有環狀體系都能用這兩種方法之一計算環張力能，比如苯分子，它存在顯著的和其環狀結構緊密相關的六中心鍵，因而環張力本來就是難以被定義的。 http://www.shanxitv.org/696 2024-01-26T09:04:00+08:00 18碳環等電子體B6N6C6獨特的芳香性：揭示碳原子橋聯硼-氮對電子離域的關鍵影響 文/Sobereva@北京科音   2024-Jan-26 0 前言 由18個碳原子連接形成的環狀體系18碳環（cyclo[18]carbon）具有十分特殊的幾何和電子結構，自從它于2019年首次在凝聚相中被觀測到后，18碳環及其衍生物就得到了化學界的廣泛關注，并陸續有大量理論研究工作發表。特別是在Carbon, 165, 468-475 (2020)中，北京科音自然科學研究中心（http://www.keinsci.com）的盧天和江蘇科技大學的劉澤玉等人證實18碳環具有明顯的雙pi芳香性特征，這是18碳環非常與眾不同的一點。 18碳環的一種等電子體是B9N9，由B和N原子交替相連構成環狀。此體系早已被研究過，被證實幾乎不具備像18碳環那樣的明顯的芳香性。最近，盧天、劉澤玉和西班牙多諾斯蒂亞國際物理中心的Mesías Orozco-Ic發現，由6個(硼-碳-氮)重復單元依次相連構成的18碳環的另一種等電子體B6C6N6具有和18碳環很接近的芳香性。為什么在無芳香性的B9N9中摻雜進去碳原子后，或者說讓碳原子橋聯硼和氮后，B9N9的芳香性就能得到巨大提升？這是非常有意思的問題。于是以上研究者對B6C6N6的特征做了全面深入的理論研究，詳細分析討論了其芳香性特征和本質，成果近期在美國化學會的知名期刊Inorganic Chemistry上發表，非常推薦閱讀： Exploring the Aromaticity Differences of Isoelectronic Species of Cyclo[18]carbon (C18), B6C6N6, and B9N9: The Role of Carbon Atoms as Connecting Bridges, Inorg. Chem., 62, 19986 (2023) https://doi.org/10.1021/acs.inorgchem.3c02675 此研究的研究意義不僅在于討論了B6C6N6本身，還在于揭示了碳原子橋聯硼和氮原子對電子離域特征的關鍵性影響，這對于通過碳摻雜的方式對純硼、氮構成的體系（如硼-氮納米管、二維層狀h-BN體系等）進行改性提供了重要啟示。此研究通過強大的Multiwfn程序（http://www.shanxitv.org/multiwfn）利用了各種波函數分析手段對B6C6N6進行了充分的考察并和18碳環、B9N9進行了對比，這也是充分運用波函數分析探究新穎體系的很好的范例。 下面本文就對Inorg. Chem., 62, 19986 (2023)這篇文章的主要研究思想、結論進行深入淺出的介紹，并對一些研究細節進行附加說明，以便于讀者更好地理解此文，以及將此文的研究手段用于自己的研究。同作者之前還對18碳環及其衍生物的各方面特征做過充分的理論研究并得到了同行的廣泛關注，成果匯總見http://www.shanxitv.org/carbon_ring.html。 1 B6N6C6的幾何結構 做B6N6C6的各方面研究之前先得得到其可靠的幾何結構。18碳環及衍生物用ωB97XD/def2-TZVP級別做幾何優化是很穩妥的，在Carbon, 165, 468-475 (2020)以及http://www.shanxitv.org/carbon_ring.html里列舉的此類體系的大量研究工作中都已經充分證明了這點。因此此文用Gaussian 16在此級別下優化了B6N6C6的結構。特別需要注意的一點是，B6N6C6基態是單重態，但是其閉殼層單重態波函數存在不穩定性，用DFT計算時需要作為對稱破缺單重態來處理方可得到真正的基態，相關討論見《談談片段組合波函數與自旋極化單重態》（http://www.shanxitv.org/82）。因此，此文對B6N6C6基態的各種研究都是基于對稱破缺單重態的結構和波函數做的。如果誤當成了閉殼層單重態來算，作者發現得到的幾何結構將明顯不合理，而且電子結構也明顯錯誤，比如芳香性嚴重偏高。實際上之前有一篇論文也算過B6N6C6，但由于那些作者沒注意這一點而誤當成了閉殼層計算，導致結果是沒有任何意義的。此例體現出DFT計算B6N6C6這種電子結構不同尋常的體系一定要注意做波函數穩定性測試。而ωB97XD/def2-TZVP級別下的B9N9、18碳環的穩定波函數都是單重態閉殼層的。 ωB97XD/def2-TZVP級別下優化得到的無虛頻的B6N6C6、B9N9和18碳環的結構如下所示，紅/灰/藍分別對應硼/碳/氮，笛卡爾坐標在文章的SI里提供了。可見B6N6C6屬于C6h點群，而B9N9屬于D9h點群。 畢竟B6N6C6是一個新穎的結構，不排除靜態相關很顯著導致ωB97XD描述不理想的可能，故為了絕對穩妥，此文還用ORCA在CASSCF(12,12)/def2-TZVP級別下也做了幾何優化。此活性空間設定下，最終會有12個占據數明顯偏離整數的pi型自然軌道，而若活性空間進一步擴大到CAS(14,14)，則會再多出兩個占據數接近整數的pi型自然軌道，這體現CAS(12,12)恰好足夠充分描述B6N6C6的pi電子的靜態相關效應。CASSCF(12,12)優化出的結構和ωB97XD的非常相近，體現出ωB97XD得到的此體系的幾何結構靠譜。為了進一步體現ωB97XD計算出的波函數也合理，文中也順帶對比了一下這兩個級別得到的單電子分布情況。文章按照《談談自旋密度、自旋布居以及在Multiwfn中的繪制和計算》（http://www.shanxitv.org/353）計算了ωB97XD波函數的自旋密度，而由于CASSCF沒法得到自旋密度，故文章使用《使用Multiwfn計算odd electron density考察激發態單電子分布》（http://www.shanxitv.org/583）介紹的方法對CASSCF(12,12)波函數產生了單電子密度（ODD），等值面圖對比如下所示。在不考慮自旋密度符號（體現在等值面顏色上，不同顏色對應單電子不同自旋方向）的情況下，可見二者展現出的單電子分布是基本吻合的，即單電子主要都分布在碳上。這在很大程度上體現出ωB97XD對稱破缺計算產生的波函數也是合理的，適合用于之后的波函數分析。 2 B6N6C6的電荷分布 原子電荷是定量考察化學體系中電荷分布最常用的指標之一。18碳環由于所有的碳都是等價的，顯然原子電荷都為0，而對B6N6C6和B9N9計算原子電荷，可以直接考察這倆體系中不同原子帶的凈電量、了解二者電荷分布的差異。原子電荷計算方法并不唯一，不同方法基于不同的思想、有不同的特點，在《原子電荷計算方法的對比》（DOI: 10.3866/PKU.WHXB2012281）中有大量介紹和對比。考慮到B6N6C6電子結構的特殊性，為了得到穩妥的結論，文中用Multiwfn計算了很常用的ADCH、Hirshfeld、Hirshfeld-I、CHELPG原子電荷（它們的計算例子看Multiwfn手冊4.7節）、用NBO程序計算了常見的NPA原子電荷，結果如下所示。可見雖然不同原子電荷給出的具體數值不同，但可以得到共同的結論，就是氮明顯帶負電，硼明顯帶正電，這體現出了它們顯著的電負性差異。而在B6N6C6中，硼和氮的凈電荷差異明顯小于它們在B9N9中，這體現出碳元素在硼、氮之間的插入明顯促進了體系整體電荷分布的均衡性。而B6N6C6中碳自己的原子電荷則接近于0，即沒怎么受周圍化學環境的影響。 18碳環具有四類分子軌道：平面內和平面外兩套π軌道，分別稱為π-in和π-out，此外還有σ軌道，以及碳的內核原子軌道組成的內核分子軌道。B6C6N6和B9N9也有這四類軌道。為了深入考察各個原子上的電子是怎么分布在這些軌道上的，文中用Multiwfn基于Hirshfeld-I原子空間劃分對三個體系計算了電子在這些軌道上的布居數，如下所示。具體計算方法我專門寫了篇文章說明：《使用Multiwfn基于Hirshfeld-I劃分計算特定類型電子在各個原子上的分布量》（http://www.shanxitv.org/697）。 從以上數據可以看到，在每個原子上，π-in和π-out電子分布量都是基本一致的，沒有什么偏向性，體現出這些體系里π-in和π-out電子特征的近似等價性。B6C6N6和18碳環中的碳的電子結構特征非常相近，π-in和π-out電子數在兩種體系里都基本為1.0，明顯都處于sp雜化的狀態。 3 B6N6C6的成鍵情況 鍵級是考察化學鍵特征非常常見的方式，參見《Multiwfn支持的分析化學鍵的方法一覽》（http://www.shanxitv.org/471）。文中用Multiwfn對B6C6N6、B9N9和18碳環計算了流行的Mayer鍵級，并且分解為了不同類型軌道的貢獻，如下所示。Mayer鍵級體現了原子間等效的共享電子對數，由數據可見，三個體系里所有鍵都包含典型的單個σ鍵特征（即σ+core電子貢獻的Mayer鍵級接近非常接近于1），同時π-in和π-out電子都構成了一定程度的π鍵特征，且二者貢獻相近。B6C6N6中的π作用程度是N-B ≈ C-N > B-C。B9N9中所有B-N等價。18碳環由長、短兩種C-C鍵構成，后者的π作用顯著強于前者。文中還用Multiwfn計算了模糊鍵級以進一步確認了這些結論，數據見補充材料。 《使用IRI方法圖形化考察化學體系中的化學鍵和弱相互作用》（http://www.shanxitv.org/598）中介紹了盧天提出的IRI-π函數，可以直觀地展現化學體系中π電子相互作用情況。此文對B6C6N6和B9N9繪制了IRI-π=1.2的等值面圖并用sign(λ2)ρ函數著色，如下所示。可見每個鍵上都出現了環狀的等值面，非常直觀地體現出這些鍵都存在明顯的π作用，這和Mayer鍵級展現出的信息相吻合。 4 B6N6C6的芳香性 B6N6C6的芳香性是當前研究的重點。和18碳環一樣，B6N6C6和B9N9的π-in和π-out電子皆為18個，都滿足Hückel的4n+2芳香性規則，因此都有可能存在雙π芳香性。衡量芳香性的方法非常多，在《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）里有充分的介紹。其中基于電子結構和基于磁性質的芳香性衡量方法的可靠性和普適性普遍較好，可以用于考察B6N6C6的芳香性并與18碳環和B9N9的進行對比。 AV1245*1000在《使用Multiwfn計算AV1245指數研究大環的芳香性》（http://www.shanxitv.org/519）里專門介紹過，是基于波函數定義的適合用于定量考察共軛大環體系芳香性的指數。文中首先使用Multiwfn計算了AV1245*1000，結果如下。除了總值外，還分別給出了π-out電子和其它電子的貢獻（注：π-in、σ、core電子對AV1245的貢獻無法精確分離故沒有單獨給出。但顯然σ和core的離域性極低，不會對AV1245有什么貢獻，故“其它電子貢獻”實質上等于π-in電子的貢獻）。由數據可見，B6C6N6的芳香性顯著強于B9N9，并且π-out的芳香性比π-in還略強一點，這可能在于體系的π-in電子的共軛程度相對略低，源自于環的曲率使得平行于環的p原子軌道彼此交疊程度低于垂直于環平面的p原子軌道。單從AV1245*1000的數值來看，B6C6N6的芳香性甚至比18碳環還要更強一點。 NICS是非常流行和穩健的基于磁屬性衡量芳香性的指標，它的不同形式在http://www.shanxitv.org/176里有專門介紹。其中NICS(0)zz對本文研究的這些環狀共軛體系最為適合，結果如下所示，數值越負說明芳香性越強。為了更充分了解芳香性的來源，文中還按照《基于Gaussian的NMR=CSGT任務得到各個軌道對NICS貢獻的方法》（http://www.shanxitv.org/670）介紹的方法分別計算了π-in、π-out和σ+core電子的貢獻。NICS(0)zz體現出B6C6N6的芳香性接近18碳環，并且π-in和π-out的貢獻相仿佛。而B9N9雖然形式上也滿足4n+2規則，但由于其NICS(0)zz基本為0，因此可斷定是非芳香性的，這必然是由于其π電子缺乏整體離域能力所致。雖然B9N9的各個B-N鍵如前述Mayer鍵級和IRI-π所示有明顯π作用，但顯然不意味著其π電子能夠容易地離域。對這些體系，σ和core電子都由于其高度定域性而對芳香性基本沒有影響。 《使用AICD程序研究電子離域性和磁感應電流密度》（http://www.shanxitv.org/147）和《使用AICD 2.0繪制磁感應電流圖》（http://www.shanxitv.org/294）介紹的AICD圖是使用得較普遍的圖形化直觀展現化學體系磁感生電流的方法。為了進一步從磁感生電流的角度展現B6C6N6的芳香性以及和B9N9的差異，文中對這些體系的π-in、π-out以及全部π電子分別繪制了AICD圖，如下所示。外磁場方向垂直于環由下朝上施加。由于等值面上的箭頭太小，為了看得清楚，圖中還用大的橙色箭頭做了清晰的標注。此圖直觀體現出，B6C6N6確實存在顯著的π芳香性，因為感生電流方向符合左手定則并且連貫地環繞整個體系，這正是典型的π芳香性體系具備的特征。而缺乏芳香性的B9N9則明顯不具備這樣的特征，感生電流幾乎只是繞著原子核或原子間區域轉。 此文還按照《考察分子磁感生電流的程序GIMIC 2.0的使用（含24分鐘演示視頻）》（http://www.shanxitv.org/491）介紹的方法繪制了GIMIC磁感生電流的動畫，如下所示，更進一步地體現出芳香性明顯的B6C6N6能夠形成整體感生環電流而非芳香性的B9N9不具備這樣的特征。 B6C6N6：http://www.shanxitv.org/attach/696/B6C6N6_GIMIC.mp4 B9N9：http://www.shanxitv.org/attach/696/B9N9_GIMIC.mp4 文中還對穿越B6C6N6、B9N9和18碳環的化學鍵的感生電流進行了積分，數值分別為21.64、0.70、22.50 nA/T，這定量體現出B6C6N6和18碳環的芳香性相仿佛，而B9N9完全不具備芳香性。 在《通過Multiwfn繪制等化學屏蔽表面(ICSS)研究芳香性》（http://www.shanxitv.org/216）和《使用Multiwfn巨方便地繪制二維NICS平面圖考察芳香性》（http://www.shanxitv.org/682）中介紹了ICSS函數，它直觀地體現了化學體系的磁屏蔽和去屏蔽區域的位置和強度，是NICS的重要擴展，已被非常廣泛應用于討論芳香性問題。此研究使用Multiwfn對B6C6N6和B9N9都繪制了ICSS_ZZ = 10 ppm的等值面圖以及垂直于環的截面的填色等值線圖，如下所示。等值面圖中橙色和藍色分別對應屏蔽和去屏蔽區域。由圖可見，B6C6N6的環內部完全是屏蔽區域，環外面完全是去屏蔽區域，這種特征和苯這樣的典型芳香性體系完全相同，更進一步嚴格地證明了B6C6N6有顯著的芳香性。而且B6C6N6的ICSS的這種特征和Carbon, 165, 468-475 (2020)報道的18碳環的高度一致，體現出二者芳香性的高度共性。而B9N9的ICSS等值面圖的屏蔽和去屏蔽區域彼此交錯，沒有形成遍及整體的連貫的等值面，體現出此體系既沒有芳香性也沒有反芳香性，而是非芳香性。 5 B6N6C6的電子離域性 芳香性的本質來自于電子的離域，因此文中還特意對B6N6C6的電子離域性從多方面進行了考察。在《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）中介紹了Multiwfn可以計算和繪制的LOL-π函數，這對于考察π電子的離域情況極其有用，如今已經非常流行。下圖對比了B6N6C6和B9N9的LOL-π=0.5等值面，明顯可見B6N6C6的等值面比B9N9顯得連貫得多，明顯體現出B6N6C6的π電子有強得多的在18個原子上全局離域的能力，也一定程度解釋了前者的芳香性比后者強得多的原因。 LOL-π函數比較適合直觀定性考察，而《在Multiwfn中單獨考察pi電子結構特征》里介紹的ELF-π函數的二分點數值則適合定量分析對比。下圖是B6N6C6和B9N9的ELF-π函數等值面圖，等值數值設為了令等值面剛好分裂時的數值，這被稱為二分值，數值越大說明電子越容易跨越二分點區域發生離域。在圖上也標注了二分點位置和二分值。由圖可見，B6N6C6的π-in和π-out的二分值都明顯大于B9N9的，這給B6N6C6的整體電子離域性更強又進一步提供了定量證據。 電子的費米穴是考察電子離域特征的一個重要的函數，它包含兩個三維位置坐標，作圖考察時通常將一個位置作為參考點，然后對另一個坐標繪圖。Atoms-in-molecules理論的提出者Bader指出an electron can go where its hole goes and, if the Fermi hole is localized, then so is the electron，也就是假設一個電子出現在某個參考點位置，那么其費米穴函數明顯分布在哪，這個電子就容易離域到哪去。因此，通過繪制費米穴圖，可以更進一步了解體系處于不同位置的電子的離域能力。文中通過Multiwfn對B6N6C6、B9N9和18碳環分別繪制了費米穴的二維圖，如下所示。參考點位置以綠色箭頭標注，設在了不同原子的價層區域距離原子核特定距離的位置，此時圖中的費米穴的分布就體現了這個位置的價電子容易離域到哪去。用Multiwfn繪制這種圖的操作在《制作動畫分析電子結構特征》（http://www.shanxitv.org/86）里有相關說明。 由上圖等值線分布可見，B6N6C6的硼和碳的價電子的費米穴分布范圍較廣，體現出電子有往較遠區域離域的能力，特別是其碳的費米穴分布特征和已知具有較顯著全局離域性的18碳環的很接近。相比之下，B9N9的硼和氮的費米穴分布廣度明顯不及B6N6C6的相應元素的，等值線只分布在距離參考點較近的原子上，體現出B9N9的價電子的離域能力相對較差。 上面的等值線圖適合定性討論，為了能把以上信息轉化成定量形式來更好地對比分析，文中定義了一個新的量叫做原子遠程離域指數（atomic remote delocalization index, ARDI）。離域性指數（DI）在《Multiwfn支持的分析化學鍵的方法一覽》（http://www.shanxitv.org/471）里專門介紹過，是衡量特定兩個原子間等效共享電子對數用的，而ARDI則是對單個原子定義的，定義為距離它兩個化學鍵以上的原子與當前原子的所有DI之和，顯然基于Multiwfn計算的所有原子間的DI可以簡單地手算出來各個原子的ARDI。ARDI體現了特定原子上的電子的遠程離域能力，數值越大說明遠程離域能力越強。環上各個原子上的電子有顯著的遠程離域能力顯然是這個環具備明顯芳香性的前提。各原子ARDI如下所示，可見B6N6C6上的原子的ARDI都較大，平均值與18碳環的碳的相仿佛，而B9N9的硼和氮的ARDI都遠小于B6N6C6相應元素的原子，這更充分體現了B9N9的電子的全局電子離域能力遠不如B6N6C6，進一步展現了B9N9芳香性遠低于B6N6C6的原因。 6 碳摻雜進B9N9使得芳香性巨幅提升的本質 根據前面的討論，文章已經充分證明了芳香性的關系是B6C6N6 ≈ C18 ? B9N9。為什么將碳元素摻雜進入無芳香性的B9N9，就能使其芳香性巨幅提升、電子全局離域能力大幅增加？實際上原因并不難理解。盡管前面的Mayer鍵級體現出B9N9中的硼-氮之間有一定π作用，但由于硼和氮的電負性差異非常大，這無疑會導致電子傾向于在氮上定域而難以往相鄰的硼上離域，至于往更遠的原子上離域就更是難上加難了。而碳的電負性則正好介于硼和氮之間，顯然硼-碳和氮-碳鍵的共價性顯著強于硼-氮鍵，因此電子更容易離域過去，相當于碳給硼和氮之間架設了能夠令電子容易離域過去的橋梁，充分調和了硼和氮之間電負性差異過大導致的電子跨鍵離域太難的矛盾。 為了更好地展現碳原子起到的價值，文中將垂直于環的各原子的2p原子軌道繪制了出來，π-out電子正是分布在這些p軌道上，π-out分子軌道也正是它們線性組合產生的。具體來說，由于NBO程序產生的預自然原子軌道（PNAO）可以較合理地描述分子環境中的原子軌道，因此文中按照《使用Multiwfn繪制NBO及相關軌道》（http://www.shanxitv.org/134）中的做法，用Multiwfn結合VMD將相應的PNAO軌道用等值面圖形式做了展示，如下所示。等值面數值選為了適合對比的0.1。圖中也把PNAO的能量（eV）標注在了上面。由圖可清楚地看出，B9N9中的氮和硼的p軌道的空間延展程度差異很大，而且軌道能量差異也很大，勢必電子難以在它們之間離域。而從B6C6N6的圖中可見，碳的p軌道延展范圍和軌道能量都恰好介于硼和氮之間，它的引入無疑能夠顯著幫助電子離域過去。 7 B6N6C6異構體的芳香性 以上研究的那種以B-C-N作為重復單元排列的B6N6C6結構是否真的是B6N6C6的能量最低結構？是否B6N6C6化學組成的體系還有芳香性更強的異構體？這是個值得考察的有趣的問題。為此，文中還另外考慮了兩種B6N6C6構型，優化后得到的無虛頻結構如下所示，可見B6N6C6'構型讓所有碳盡可能連在一起，其它部分部分保持B9N9那樣硼-氮交錯的結構，而B6N6C6''構型則是由三個C-C-B-N-B-N重復單元構成的結構。這兩個結構的基態是閉殼層的，能量分別比前面研究的B-C-N為重復單元構成的B6N6C6結構能量低280.0和250.2 kcal/mol。這體現出B6N6C6中B-C和C-N的平均鍵能比B-N和C-C是要更小的。雖然前面研究的B6N6C6結構能量不是這種化學組成的能量最低結構，但并不意味著它是不穩定的。根據《使用ORCA做從頭算動力學(AIMD)的簡單例子》（http://www.shanxitv.org/576）中介紹的方法，文中對此體系做了較高溫度（500 K）下20 ps長度的基于wB97X-D3/6-311G*級別的從頭算分子動力學模擬，用VMD每1 ps繪制一次的結構疊加圖如下所示（顏色按照藍-白-紅變化），可見環狀結構始終很好維持著，并沒有發生異構化、解離等現象，體現出B6N6C6具有一定熱力學穩定性。 文中還對B6N6C6'和B6N6C6''構型計算了NICS(0)zz，發現分別僅有1.8和1.4 ppm，可見這倆構型沒有芳香性。這說明只有讓碳插入到硼和氮中間，從而以B-C-N作為重復單元時，碳才能真正起到顯著提升純硼-氮環體系的電子離域性的作用。 8 B6N6C6最低三重態激發態的芳香性 前面考察的都是B6N6C6的單重態基態（S0）。最后，文章還研究了B6N6C6最低三重態激發態（T1態）的芳香性。計算發現S0-T1垂直和絕熱激發能分別為1.62和1.17 eV。下圖(a)是優化后的T1態B6N6C6的結構，可見對稱性相對于S0態發生了下降。下圖(b)是按照《使用Multiwfn作電子密度差圖》（http://www.shanxitv.org/113）繪制的在B6N6C6的S0結構下T1態與S0態的密度差，藍色和綠色分別為電子密度降低和增加部分。由圖可見S0到T1態的電子激發是π-in → π-out的激發。由于垂直激發導致的電子密度變化破壞了B6N6C6基態的C6h對稱性，這也是為什么在T1態勢能面上從Franck-Condon點往極小點弛豫過程中幾何結構對稱性發生了下降。 文中對B6N6C6的T1態計算了NICS(0)zz，并分解為了不同類型電子的貢獻。而且為了考察幾何結構弛豫的影響，在S0極小點和T1極小點分別做了這種計算，結果如下所示。可見無論在哪種結構下，T1態由于NICS(0)zz遠大于0，因此是顯著反芳香性的，這和Baird規則所述情況一致。Baird規則指出單重態基態滿足Hückel 4n+2規則的體系的最低三重態是反芳香性的。從NICS(0)zz還看到，T1態的反芳香性特征在T1極小點結構下比S0極小點結構下更弱，這體現出B6N6C6自發減小T1態的反芳香性可能是其S0→T1垂直激發后結構弛豫的重要驅動力。 9 總結 此文基于量子化學計算和Multiwfn等程序做波函數分析，非常全面系統研究了18碳環的等電子體B6N6C6的幾何結構和電子結構，從不同角度全面考察了其芳香性和電子離域性特征，并且和18碳環及其另外一個重要的等電子體B9N9進行了細致的對比。本文體現出以碳原子橋聯硼和氮原子，可以顯著提升純硼-氮體系的電子離域性。這不僅對電子結構產生重要影響，必然也同時會影響其它性質，諸如光學性質、電子輸運性質等。本文雖然研究的只是一維體系，可以預想電負性介于硼和氮之間的碳元素的摻雜必定也會對二維、三維純硼-氮體系產生顯著影響，而且摻雜方式的選取可能會對這些體系的性質起到一定調控作用。 如果讀者對于18碳環及其相關體系感興趣，非常推薦查看http://www.shanxitv.org/carbon_ring.html里匯總的本文作者的巨量其它研究工作，包含大量深入淺出的論文內容和研究思想介紹文章，以及對計算技術細節的附加說明。 http://www.shanxitv.org/693 2023-12-30T12:35:00+08:00 使用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。明顯凍核對激發能的影響可完全忽略不計，而耗時大為降低，因而有重要實際意義。 http://www.shanxitv.org/685 2023-08-19T20:15:00+08:00 使用sobEDA和sobEDAw方法做非常準確、快速、方便、普適的能量分解分析 文/Sobereva@北京科音 First release: 2023-Aug-31    Last update: 2024-Feb-7 1 前言 計算化學領域的能量分解分析（energy decomposition analysis, EDA）方法很多，思想各有不同，幾種常見的在《Multiwfn支持的弱相互作用的分析方法概覽》（http://www.shanxitv.org/252）中的能量分解部分有簡要介紹，在量子化學波函數分析與Multiwfn程序培訓班（http://www.keinsci.com/workshop/WFN_content.html）里有系統的介紹。雖然也有基于分子力場的能量分解，比如《使用Multiwfn做基于分子力場的能量分解分析》（http://www.shanxitv.org/442）里介紹的EDA-FF，但絕大多數能量分解方法是基于量子化學的，并且其中大多數的目的是將體系中用戶自定義的兩個或多個片段間的相互作用能進行分解，得到不同物理成份的貢獻，從而能夠深入了解片段間相互作用的主要本質（靜電、色散、交換-互斥、軌道相互作用等），以及橫向對比不同情況下相互作用的內在差異，例如《透徹認識氫鍵本質、簡單可靠地估計氫鍵強度：一篇2019年JCC上的重要研究文章介紹》（http://www.shanxitv.org/513）介紹的研究氫鍵的文章中使用了流行的SAPT能量分解方法充分考察不同體系中氫鍵的本質和差異并總結出了規律、《靜電效應主導了氫氣、氮氣二聚體的構型》（http://www.shanxitv.org/209）介紹的文章里用DFT-SAPT考察了不同構型的氫氣/氮氣二聚體相互作用的主導因素。 下文說的能量分解一律指基于片段和量子化學的能量分解方法。這類方法之前已經有很多，但在筆者看來各有各的缺點或局限性，尚沒有一個完全同時滿足耗時低、使用簡單方便、普適性強、精度好、物理意義清晰的能量分解方法。特別是對于最流行的Gaussian量子化學程序的用戶，目前尚沒有一個不需要額外代價就能做能量分解的方法。為了解決這個難題，筆者基于如今非常流行的色散校正的DFT理論提出了名為sobEDA的能量分解方法，同時還提出了一個它專用于考察弱相互作用的變體sobEDAw。在筆者來看，這兩個方法完全掃清了廣大量子化學研究者做能量分解的一切障礙。此方法目前已經正式發表在JPCA上，非常建議讀者仔細閱讀： Tian Lu, Qinxue Chen, Simple, Efficient, and Universal Energy Decomposition Analysis Method Based on Dispersion-Corrected Density Functional Theory, J. Phys. Chem. A, 127, 7023 (2023) https://doi.org/10.1021/acs.jpca.3c04374 同時，此文還給出了筆者開發的名為sobEDA.sh的Bash shell腳本，可以很容易地將非常流行的量子化學程序Gaussian和波函數分析程序Multiwfn（http://www.shanxitv.org/multiwfn）相結合做sobEDA、sobEDAw能量分解。 在下文第2節，筆者將對sobEDA和sobEDAw的特點以盡量簡短的語言進行介紹，使得讀者能夠有最低限度的知識正確理解和使用這兩種方法，完整、全面的介紹和討論請仔細閱讀上述原文。第3節將會對sobEDA.sh腳本的基本用法做簡要介紹。第4節將會給出一系列具體例子充分展示如何用sobEDA.sh對各種類型體系和情況做sobEDA和sobEDAw能量分解。最后在第5節將對本文進行總結。 2 sobEDA和sobEDAw簡介 2.1 sobEDA/sobEDAw的主要特點 首先讓大家快速了解一下sobEDA/sobEDAw方法及實現它們的sobEDA.sh腳本的關鍵優點，從而知道為什么sobEDA和sobEDAw極其值得推薦！ ? 門檻低&方便：只要是裝有Gaussian（G16 A.03及之后）和Multiwfn程序的Linux機子就直接能跑sobEDA.sh腳本。此門檻很低，因為Gaussian是用戶規模最大、使用得最普遍的量子化學程序，雖然收費但不算貴，安裝也很簡單。Multiwfn是用戶最多的波函數分析程序，完全免費，在Linux機子上的安裝也不費事。 sobEDA.sh的輸入文件的創建又容易又簡單。用戶僅僅需要懂得Gaussian算DFT單點能的輸入文件怎么寫就能順利使用，而這屬于幾乎所有量子化學研究者的最基本常識。 ? 快：在像樣的雙路服務器下，算含有幾百個原子的體系都完全沒問題，耗時僅比DFT算整個體系貴一倍多，也即DFT計算單點能能算得動什么體系，就能對什么體系做sobEDA/sobEDAw能量分解。而且計算過程中對內存和硬盤的消耗微乎其微。 ? 普適：兩個和多個片段、閉殼層和開殼層體系、主族和含過渡金屬體系、平衡結構和非平衡結構、化學鍵作用和弱相互作用，全都適用！ ? 合理：各種量子力學效應都充分考慮，結果準確、物理意義清楚，并且與之前廣為接受的很多主流能量分解方法的結果高度吻合。 由于sobEDA、sobEDAw沒有明顯缺點，又有以上眾多優點，因此筆者十分建議把sobEDA、sobEDAw作為做能量分解的默認的方法。此外，這兩種方法的原理和實現很簡單，任何其它量子化學程序也都可以非常容易地將之實現在其中。 已有人在體驗過sobEDA后評論道（http://bbs.keinsci.com/thread-39400-1-1.html）：在修回的文章里想添加EDA部分，本來是想用SAPT做的，結果計算速度慢只算的動SAPT0不說，有些結構還死活不收斂，弄得我心力憔悴。在Sobevera老師的推薦下在截至日期的前幾天轉用了他新開發的sobEDA，試用下來的效果只能用一句話形容，“令人震驚的好”！ 2.2 sobEDA方法簡介 首先要明確一點，sobEDA及下一節介紹的其變體sobEDAw是對于特定幾何結構下片段間相互作用能進行分解，因此不考慮單體的變形能（即片段從其孤立狀態極小點結構變成它在復合物中的結構過程中的能量升高量）、不考慮熱力學效應（比如結合過程中焓或者自由能熱校正量的變化）、不考慮溶劑效應（即能量分解是對真空環境而言的）。如果要明確考慮這些效應，可以自己手動計算并作為單獨的能量項予以討論，也令能量分解給出的信息更為豐富，其過程也一點都不復雜，具體公式在前述的JPCA文章的2.1節末尾給明了。 sobEDA和sobEDAw都是基于考慮了色散校正的DFT方法定義的。如果不了解色散校正的話，看《DFT-D色散校正的使用》（http://www.shanxitv.org/210）和《談談“計算時是否需要加DFT-D3色散校正？” 》（http://www.shanxitv.org/413）。sobEDA和sobEDAw目前只考慮了普通泛函，精度已經足夠滿足實際研究需要，更昂貴的雙雜化泛函目前不納入考慮。 下面來看一下sobEDA的思想、各個能量項是怎么來的。 化學體系的實際狀態的波函數可以視為是由各片段的波函數經歷三個步驟構成的，如下所示，這里為了表示簡單只考慮了A、B兩個片段。ΨA+ΨB對應的是片段之間沒有任何相互作用的狀態，ΨAΨB對應的是兩個片段間已經出現了相互作用，但電子結構沒有發生變化，因此整體的占據軌道是這兩個片段占據軌道的并集，這稱為準分子狀態（promolecule state）。Ψ0AB是進一步考慮這倆片段的電子之間的Pauli互斥令占據軌道變形（滿足正交關系）之后的波函數，也叫做凍結態（frozen state）。最后的ΨAB是真實狀態波函數，也相當于常規DFT計算在SCF收斂后產生的波函數。 如上圖可見，體系狀態的每個變化過程都可能伴隨著動能（T）、經典靜電作用能（J）、交換能（X）、DFT相關能（C）、色散校正能（D）的變化。上圖中彩色字對應的是sobEDA方法分解出來的相互作用中的各個物理成份，顯然它們的加和精確等于色散校正的DFT方法以常規的超分子方式（即整體能量減各個片段能量）算的相互作用能ΔE_int，它可以寫為 ΔE_int = ΔE_els + ΔE_x + ΔE_DFTc + ΔE_dc + ΔE_rep + ΔE_orb sobEDA給出的以上相互作用成份的物理含義如下 ΔE_els：片段間經典靜電作用能，els代表electrostatic。可正可負 ΔE_x：交換作用能，x代表exchange。數值為負，由DFT交換泛函貢獻，對雜化泛函顯然其中還牽扯到一部分HF交換能。這體現交換相關效應對片段間相互作用能的貢獻 ΔE_DFTc：DFT相關能。DFTc代表DFT correlation。一般數值為負，由DFT相關泛函貢獻。體現庫侖相關效應對片段間相互作用能的貢獻 ΔE_rep：Pauli互斥能，rep代表repulsion。數值為正。體現不同片段的電子間為了滿足Pauli互斥原理造成的能量升高 ΔE_orb：軌道相互作用能，orb代表orbital。數值為負。來自于片段內和片段間占據軌道與非占據軌道間混合導致的能量變化，物理上體現片段內電子分布被極化、片段間電子轉移對能量的總影響，共價作用也由此項體現 ΔE_dc：色散校正能。dc代表dispersion correction。數值為負。物理含義取決于具體情況，見下文 為了簡化討論，可以將ΔE_x與ΔE_rep相加作為交換-互斥項ΔE_xrep，這樣的項在很多能量分解方法里都有。如果不需要單獨討論色散作用的話，ΔE_DFTc和ΔE_dc還可以相加作為庫侖相關項ΔE_c，因為色散校正在本質上來說主要是對DFT相關泛函描述不好的長程庫侖相關部分進行彌補。 在目前的sobEDA、sobEDAw實現中，色散校正項ΔE_dc是由極為流行的DFT-D3(BJ)色散校正方法來計算的，原理上也可以用其它色散校正方法，如DFT-D4（http://www.shanxitv.org/464）、VV10等。特別要注意的是，僅當你用的泛函是完全無法描述色散作用的泛函，ΔE_dc才能解釋為色散作用能。什么泛函是完全無法描述色散作用的？在sobEDA原文里用一些方法對Ar二聚體做了能量隨距離的掃描，如下所示。眾所周知Ar和Ar之間是最典型的范德華作用，色散作用就是其中起到吸引作用的部分。可見BLYP、B3LYP、BHandHLYP在勢能曲線上都不存在極小點，即它們完全無法描述色散作用，因此sobEDA結合BLYP-D3(BJ)、B3LYP-D3(BJ)、BHandHLYP-D3(BJ)使用時，ΔE_dc才有明確的物理意義、可以被解釋為色散作用能。 值得一提的是，sobEDA/sobEDAw給出的軌道相互作用能ΔE_orb還可以用Multiwfn進一步通過ETS-NOCV方法進行分解，能得到不同類型軌道相互作用的能量貢獻以及對應的電子密度變化，對于探究片段間軌道相互作用本質極其有用！詳見《使用Multiwfn通過ETS-NOCV方法深入分析片段間的軌道相互作用》（http://www.shanxitv.org/609）。 sobEDA的計算耗時相當低，具體實現的細節和流程在前述的JPCA原文2.1節里說了，總共耗時僅僅是DFT計算整個體系的單點能的兩倍左右。所以但凡DFT算單點沒壓力，sobEDA方法做能量分解也必定輕輕松松。對如今較好的雙路服務器來說，基于Gaussian的DFT功能來實現sobEDA，對于好幾百原子都不費勁。 2.3 sobEDAw方法簡介 在說sobEDAw之前，先專門說一下色散作用項。一般意義上的色散作用指的是隨作用距離r具有-1/r^6衰減行為的長程庫侖相關作用。但是當分子間不遠不近，比如在分子復合物的極小點結構下，分子間庫侖相關作用既有中程的也有長程的。《使用PSI4做對稱匹配微擾理論(SAPT)能量分解計算》（http://www.shanxitv.org/526）中介紹的SAPT方法在對分子間弱相互作用的能量分解方面用得很普遍。SAPT（及其變體DFT-SAPT、SAPT(DFT)等）的色散作用項ΔE_disp包含了分子間全部庫侖相關作用，也就是說，當分子間距離不很遠時，SAPT的ΔE_disp項不僅包括了一般意義的色散作用，實際上還包含了中程庫侖相關作用。因此還有人特意管ΔE_disp叫做dispersion-like（而非dispersion）項。ΔE_disp由此也在一定程度上高估了一般意義的色散作用強度。盡管SAPT的色散項不那么“純粹”，SAPT在分子間弱相互作用的能量分解方面特別流行，而且很多文章都用SAPT給出的ΔE_disp和靜電作用項ΔE_els的比值來對相互作用體系進行分類，說誰是靜電主導、誰是色散主導、誰是混合型作用，比如非常知名的S66弱相互作用測試集的原文里就是這樣做的。前面說了，sobEDA結合比如B3LYP-D3(BJ)時的ΔE_dc對應的是“純”色散部分，因此在分子團簇極小點結構附近的ΔE_dc的大小勢必明顯小于ΔE_disp，也因此沒法把sobEDA的結果和SAPT的直接對比。顯然這不是因為sobEDA不合理，而是在于SAPT的色散項定義的特殊性所致。 考慮到以上問題，為了在sobEDA框架內能得到與SAPT相對比的結果，特別是令色散/靜電作用的比例能與SAPT的基本吻合、消除由內在定義差異帶來的明顯分歧，我提出了sobEDA的變體，名為sobEDAw，其后綴w特意強調此方法和SAPT一樣是專門用于對弱（weak）相互作用能進行分解的。sobEDAw的能量分解方式為 ΔE_int = ΔE_els + ΔE_xrep + ΔE_orb + ΔE_disp 其中ΔE_int、ΔE_els、ΔE_orb都和sobEDA完全一樣，但sobEDAw將交換-互斥作用項ΔE_xrep和色散作用項ΔE_disp定義為 可見，sobEDAw引入了一個介于[0,1]的系數w，用來將sobEDA的ΔE_dc和一部分ΔE_DFTc組合成sobEDAw的ΔE_disp，而剩下的ΔE_DFTc和sobEDA的ΔE_xrep組成了sobEDAw的ΔE_xrep。w的計算除了依賴于sobEDA的ΔE_dc和ΔE_els外，還依賴于針對每種泛函/基組的組合而擬合的一套參數c、a、r。擬合的目標是令sobEDAw的ΔE_disp/ΔE_els比例與SAPT的盡可能一致，實際當中是對S66測試集與其文中的DFT-SAPT的數據盡可能吻合來擬合的。從sobEDA原文的圖2(d)可見這個目的充分達到了，因此sobEDAw可以很好地代替昂貴的SAPT。 由上可見，sobEDAw是對sobEDA的進一步延伸，當sobEDA數據算出來后，通過非常簡單的公式一瞬間就能進一步得到sobEDAw的結果。做sobEDAw分析也因此必然順帶產生sobEDA的數據。 大家都知道計算弱相互作用需要對基組重疊誤差（BSSE）問題做恰當的考慮，相關介紹見《談談BSSE校正與Gaussian對它的處理》（http://www.shanxitv.org/46）。以counterpoise方式校正BSSE問題，絕大多數情況可以在基組層面令弱相互作用能的計算精度進行提升。在sobEDA或sobEDAw能量分解計算中，如果計算每個單體都使用整個體系里所有原子的基函數，則給出的ΔE_int及其各個能量成份就都相當于是counterpoise方法校正后的。 下面的表格是前述JPCA文中的表1。表中體現了非常常用的泛函B3LYP-D3(BJ)結合不同基組對S66測試集計算的相互作用能與極高精度CCSD(T)/CBS參考值的誤差（kcal/mol），MAE是平均絕對誤差。w.CB代表with complex basis functions，也即計算單體時使用整個復合物的所有基函數，理解成在sobEDAw分析時考慮了counterpoise校正也可以。 從以上表格可見6-31+G** w.CB是圖便宜首選，精度已經不錯，MAE才區區0.26 kcal/mol。要求更好精度且可以接受更高耗時的話推薦用6-311+G(2d,p) w.CB。def2-TZVP w.CB也值得推薦，性價比雖然不及6-311+G(2d,p) w.CB，但好處是支持的元素非常多，前六個周期都有定義。嫌def2-TZVP w.CB太貴的話也可以用更便宜的def2-TZVP(-f) w.CB，這是去掉了主族元素的f極化函數的版本，精度基本沒打折扣。平時就用這幾個基組做sobEDAw就足夠了，上表里的其它基組都沒什么值得額外特別推薦的。 上表里也給出了對S66測試集用sobEDAw計算的ΔE_disp/ΔE_els比例與DFT-SAPT的平均絕對誤差，以及對這些基組擬合的c、a、r參數。可見B3LYP-D3(BJ)結合上面推薦的幾個基組時的平均絕對百分比誤差（MAPE）才區區6%左右，表現都很理想，和昂貴的DFT-SAPT數據的吻合度相當高。 以上表格也體現了耗時情況。最后一列的18+224原子的體系是《8字形雙環分子對18碳環的獨特吸附行為的量子化學、波函數分析與分子動力學研究》（http://www.shanxitv.org/674）里介紹的主-客體復合物體系。在《淘寶上購買的雙路EPYC 7R32 96核服務器的使用感受和雜談》（http://www.shanxitv.org/653）介紹的96核機子上，筆者用B3LYP-D3(BJ)/6-31+G** w.CB對這么大的體系做sobEDAw能量分解僅僅花了71分鐘，可見耗時很低。只要有耐心，這樣的條件對500原子體系做能量分解都是可行的。如果是對14+16原子的體系，總共才39秒就完了事。 考慮到B3LYP-D3(BJ)并不總適合所有情況，因此筆者在JPCA文中表2里還給出了其它幾種很常見、適合不同情況的泛函結合前述幾種基組的c、a、r參數，以及誤差統計，可見都表現良好。 從上一節中的Ar-Ar勢能面掃描曲線看到，PBE、PBE0對色散作用的描述較為明顯，而M06-2X對色散作用的描述更是很充分，它們從原理上就不能結合sobEDAw，筆者也發現無法對它們擬合出合理的參數。因此以上表格沒有對這些也很常見的泛函進行考慮。表格里的那些的泛函已經足夠用了。 2.4 sobEDA和sobEDAw搭配泛函和基組的建議 鑒于之前有一些用戶對于sobEDA和sobEDAw分析搭配的泛函和基組有點糊涂，也有些人缺乏量子化學研究經驗、不太懂泛函和基組的選擇，這一節專門強調一下這個問題。 首先要明確一點，sobEDA和sobEDAw都可以用于研究弱相互作用，后者的額外好處是其色散項與流行的SAPT給出的可以相對比。對于化學鍵程度的強相互作用，或者片段間既有強相互作用也有弱相互作用，則只適合sobEDA，而sobEDAw和SAPT都不適合這類情況。 (1) 做sobEDAw用的泛函和基組 sobEDAw計算用的級別（泛函+基組）必須是已經有了現成的c、a、r參數的才行。在上面提到的JPCA原文的表1、2中，已經給了B3LYP-D3(BJ)結合大量基組的參數，以及其它幾種知名泛函TPSS-D3(BJ)、TPSSh-D3(BJ)、BLYP-D3(BJ)、BHandHLYP-D3(BJ)結合3種值得使用的基組的參數，因此它們都可以結合sobEDAw使用。如果你非要用其它級別的，則需要自行擬合參數。擬合用的程序、數據集和做法見http://www.shanxitv.org/soft/sobEDAw_fit.zip。一般完全沒必要自己擬合。 先說泛函方面。非常流行的B3LYP-D3(BJ)適合作為sobEDAw一般情況下默認用的泛函。但如《談談量子化學研究中什么時候用B3LYP泛函優化幾何結構是適當的》（http://www.shanxitv.org/557）所說的，B3LYP-D3(BJ)算pi大共軛體系往往不好，算過渡金屬配合物也很一般。對于pi大共軛體系建議改用離域性誤差較小的BHandHLYP-D3(BJ)，對過渡金屬配合物建議用通常表現更好的TPSSh-D3(BJ)。其它情況可以隨機應變，例如如《簡談量子化學計算中DFT泛函的選擇》（http://www.shanxitv.org/272）所示，TPSS和TPSSh算Cu、Ag、Au團簇較好，因此算這種團簇對小分子的物理吸附做sobEDAw能量分解可以用TPSS-D3(BJ)和TPSSh-D3(BJ)。 再說基組方面。要求準定量準確的話，搭配6-31+G(d,p) w.CB即可，要求更高精度用6-311+G(2d,p) w.CB。若有的元素是這倆沒定義的，改用def2-TZVP w.CB或更便宜的def2-TZVP(-f) w.CB。def2-TZVP(-f) w.CB的耗時和6-311+G(2d,p) w.CB差不太多。 (2) 做sobEDA用的泛函和基組 原理上sobEDA可以結合任意泛函（雙雜化除外）和任意基組使用，因為sobEDA自身不依賴任何經驗參數。實際用的級別只要是計算當前體系片段間相互作用能較好、耗時你能接受的即可。從實際角度來說，在泛函方面，一般用B3LYP-D3(BJ)做sobEDA就可以，如果是計算過渡金屬與配體間相互作用，用TPSSh-D3(BJ)或PBE0-D3(BJ)通常更好。基組方面，對主族元素用6-311G(d,p)，過渡金屬用SDD贗勢結合標配基組就不錯。若想要更好結果，或者對6-311G(d,p)沒定義的元素，可以用def2-TZVP。更多情況的泛函和基組的推薦在本文限于篇幅沒法展開介紹，請參閱這些文章里關于泛函和基組的討論：《簡談量子化學計算中DFT泛函的選擇》（http://www.shanxitv.org/272）、《談談量子化學中基組的選擇》（http://www.shanxitv.org/336）、《談談贗勢基組的選用》（http://www.shanxitv.org/373）。 sobEDAw是沒法用混合基組的，因為不同基組牽扯的擬合參數不同。而sobEDA可以混合基組。 sobEDA不是必須用DFT-D3(BJ)，用其它任何色散校正都可以，也包括DFT-D3(0)。JPCA原文中sobEDAw擬合參數時只考慮了DFT-D3(BJ)的情況因此只能用這個。 無論是sobEDA還是sobEDAw，都不能結合ωB97XD、CAM-B3LYP等范圍分離泛函。不清楚哪些是范圍分離泛函的話看《不同DFT泛函的HF成份一覽》（http://www.shanxitv.org/282）。這不是sobEDA/sobEDAw理論的原因，而在于Gaussian程序（至少截止到目前撰文時最新的G16 C.02為止）的Link 608模塊對范圍分離泛函的情況沒法給出總能量當中的各種成分，也因此下面的sobEDA.sh腳本無法提取需要的信息。如果以后Gaussian改進了這點，sobEDA/sobEDAw也就沒這個限制了。 2.5 sobEDA和sobEDAw的準確性和可靠度 在sobEDA原文的第3節里通過8個涉及截然不同情況的例子（中性和離子、閉殼層和開殼層、平衡構型和非平衡構型、二聚體和三聚體、單重鍵和多重鍵、有機/無機分子和過渡金屬配合物，等等）。充分證明了sobEDA和sobEDAw的準確、可靠、普適和極高的實際應用價值。本文限于篇幅就不再復述一遍了，請讀者務必閱讀原文以對sobEDA和sobEDAw的實際應用有充分的認識，相信這些例子能令很多讀者受到不少啟發。 這里就僅僅展示原文里的一個例子，腺嘌呤-胸腺嘧啶氫鍵和堆積型二聚體。它有氫鍵和pi-pi堆積兩種構型，作用強度和作用本質都不同。下圖是sobEDAw和高階的SAPT的結果的對比，其中SAPT2+3(CCD)/aug-cc-pVTZ雖然精確但極為昂貴。可見盡管二者在原理上差異極大，但不僅兩個方法算的總相互作用能差異很小，而且各項都能相符得很好，ΔE_disp/ΔE_els吻合度相當高，因此便宜的sobEDAw完全可以代替很昂貴的高階SAPT。雖然SAPT里最低階的SAPT0相對便宜，但耗時還是顯著高于sobEDAw結合合適的計算級別，同時精度還差不少。此外，做SAPT的最主流程序PSI4的SCF收斂性比Gaussian有較大差距，改用基于Gaussian+Multiwfn做的sobEDAw就沒這個煩惱了。 3 sobEDA.sh腳本的使用 sobEDA/sobEDAw的原理非常簡單，只要基于DFT的程序的開發者有意向，非常容易地就能在程序里實現。sobEDA.sh是專門結合Gaussian和Multiwfn實現sobEDA/sobEDAw分析的Bash shell腳本，在這一節以盡量簡短的語言介紹一下基本使用方法。sobEDA.sh腳本文件可以在http://www.shanxitv.org/soft/sobEDA_tutorial.zip下載，并且里面的pdf文檔是此腳本最完整、詳細的使用說明，建議讀者有時間時完整看一遍。 3.1 使用要求 sobEDA.sh腳本的運行環境必須滿足以下幾點： (1)必須是Linux環境。不一定非得是實體Linux系統，也可以是比如Windows下用VMware虛擬機創建的Linux客戶機、Windows下的WSL環境等 (2)機子里已正常安裝Gaussian，比如可以通過g16命令正常調用Gaussian 16。必須是G16 A.03及以后版本。不會裝Gaussian的話見《Gaussian的安裝方法及運行時的相關問題》（http://www.shanxitv.org/439） (3)機子里已正常安裝Multiwfn，比如直接輸入Multiwfn命令就可以正常啟動。Multiwfn可以在http://www.shanxitv.org/multiwfn免費下載，必須是2023-Jun-25以后更新的版本（如果不知道當前機子里裝的Multiwfn是幾幾年幾月幾號更新的版本就看啟動程序時明確顯示的release date）。圖安裝省事也可以裝noGUI版，此時不需要額外裝任何圖形庫，但此版本的Multiwfn不具備任何和作圖有關的功能和圖形界面。Multiwfn的安裝方法見《Multiwfn在Linux下安裝的中文說明》（http://www.shanxitv.org/688）。 (4)dos2unix命令必須可以正常使用。機子里沒裝的話就裝上，比如Redhat/Rocky Linux/CentOS系統里用yum install dos2unix即可迅速安裝上。 3.2 sobEDA.sh的運行方法 (1)把sobEDA.sh放在任意目錄，并用chmod +x [文件路徑]命令給它加上可執行權限 (2)根據實際情況恰當修改sobEDA.sh里面的設置 (3)將整個體系的結構保存為system.xyz并放在當前目錄下 (4)在當前目錄下創建fragment.txt，在里面定義各片段的原子序號、凈電荷、自旋多重度 (5)在當前目錄下創建template.gjf，作為sobEDA.sh產生Gaussian輸入文件的模板文件 (6)輸入sobEDA.sh腳本的路徑啟動之 sobEDA.sh啟動后會自動調用Gaussian和Multiwfn進行計算、提取和處理數據，最后把結果輸出在屏幕上。運行中途產生的Gaussian的gjf文件、out文件、chk/fch文件都會自動保留，用戶可以自行查看。用Multiwfn載入fch文件后還可以考察各個任務產生的波函數的狀態，比如根據《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）的做法觀看里面的軌道，根據《談談自旋密度、自旋布居以及在Multiwfn中的繪制和計算》（http://www.shanxitv.org/353）的做法計算自旋密度和自旋布居。 下面說一下具體細節。 3.3 sobEDA.sh的設置 用文本編輯器打開sobEDA.sh腳本，就可以在開頭部分看到各種設置，里面有非常清楚的注釋，這里只說幾個重點： ? gau和mwfn變量分別對應調用Gaussian和Multiwfn用的命令，應根據實際情況設置。比如你的機子上的noGUI版Multiwfn是用Multiwfn_noGUI命令啟動的，顯然設為mwfn="Multiwfn_noGUI" ? sobEDAw=0代表做sobEDA分析，=1代表在sobEDA分析后還輸出sobEDAw的結果 ? iCP=0代表不用counterpoise校正，=1代表用（對應基組帶w.CB后綴的情況）。分解弱相互作用的時候強烈建議用iCP=1以改進結果，分解化學鍵作用時則不建議iCP=1，不僅浪費很多時間，還可能令結果更差，見《計算化學鍵鍵能時考慮BSSE不僅是多余的甚至是有害的》（http://www.shanxitv.org/381）。 用sobEDAw的時候必須在sobEDA.sh里設置c、a、r參數，分別對應parm_c、parm_a、parm_r變量。為了方便用戶，sobEDA.sh腳本里已經自帶了一批常用級別的參數（和JPCA原文里表1、2一致），默認都是被注釋掉的狀態（行首帶#）因此不生效。例如腳本里可以看到 # B3LYP-D3(BJ)/6-31+G(d,p) with iCP=1 #parm_c=0.575;parm_a=0.071;parm_r=2.571 # B3LYP-D3(BJ)/6-311+G(2d,p) with iCP=1: #parm_c=0.550;parm_a=0.037;parm_r=2.750 # BHandHLYP-D3(BJ)/6-31+G(d,p) with iCP=1: #parm_c=0.634;parm_a=-0.065;parm_r=0.702 ...略 你用哪種級別就把哪種情況的參數取消注釋，即去掉開頭的#。比如要用B3LYP-D3(BJ)/6-311+G(2d,p) w.CB，就把以上部分改為 # B3LYP-D3(BJ)/6-31+G(d,p) with iCP=1 #parm_c=0.575;parm_a=0.071;parm_r=2.571 # B3LYP-D3(BJ)/6-311+G(2d,p) with iCP=1: parm_c=0.550;parm_a=0.037;parm_r=2.750 # BHandHLYP-D3(BJ)/6-31+G(d,p) with iCP=1: #parm_c=0.634;parm_a=-0.065;parm_r=0.702 ...略 顯然，其它級別的參數也都可以自行定義在里面。對于sobEDAw=0的情況，就完全不用管這些sobEDAw參數設沒設、怎么設的了。 3.4 system.xyz文件 system.xyz文件包含當前整個體系的坐標。不了解xyz格式的話看《談談記錄化學體系結構的xyz文件》（http://www.shanxitv.org/477）。產生這個文件極為簡單，比如Multiwfn載入pdb、mol、mol2等主流坐標格式的文件（支持什么其它的看http://www.shanxitv.org/379），或載入Gaussian/ORCA優化任務的輸出文件（會讀取最后一幀結構），然后在主界面輸入xyz，再輸入system.xyz，在當前目錄下就有了這個文件了。 3.5 fragment.txt文件 此文件里定義能量分解有幾個片段、各對應哪些原子。比如下面的例子代表sobEDA/sobEDAw分析是對3個片段間相互作用能做分解，第1個片段的凈電荷為0、自旋多重度為1，原子序號為1-3,7,10-12。第2個片段的凈電荷為-1、自旋多重度為1，原子序號為4-6。第3個片段的凈電荷為+1、自旋多重度為1，原子序號為8,9。注意所有片段中的原子的并集必須對應整個體系。 3 0 1 1-3,7,10-12 -1 1 4-6 1 1 8,9 自旋多重度寫為負值時表示片段帶的是beta單電子。比如下面的例子，是兩個片段，片段1包含1-4號原子，是三重態，單電子為alpha自旋。片段2包含5-8號原子，也是三重態，但單電子為beta自旋。 2 0 3 1-4 0 -3 5-8 3.6 template.gjf文件 這是Gaussian輸入文件的模板文件，里面記錄了計算過程要用的關鍵詞。[geometry]部分會被腳本自動替換成實際要算的坐標。此文件里凈電荷和自旋多重度設的是對整個體系計算時用的（算各個片段時用的是fragment.txt里定義的）。#P和nosymm關鍵詞是必須有的。ExtraLinks=L608關鍵詞結合此文件末尾的數字用來讓Gaussian輸出當前總能量中的各個成份，它們是要被sobEDA.sh腳本自動提取和處理的。在此基礎上，若遇到SCF不收斂，可以自行加入幫助SCF收斂的關鍵詞，比如scf(novaracc,noincfock)之類的，參考《解決SCF不收斂問題的方法》（http://www.shanxitv.org/61），但有的關鍵詞不兼容，比如guess=AM1等。 #P b3lyp/6-311+G(2d,p) em=gd3bj ExtraLinks=L608 nosymm   <---空行 title line   <---空行 0 1 [geometry]   <---空行 -5 5   <---空行   <---空行 上面例子末尾的-5對應當前用的泛函在Gaussian中的內部IOp參數，必須根據實際情況設置，需要查閱Gaussian IOps手冊的IOp(3/74)部分，見http://gaussian.com/overlay3/#iop_(3/74)。一些比較常見的泛函的這個值列在這里了：-73 (MN15), -55 (M06-2X), -54 (M06), -53 (M06L), -35 (TPSSh), -13 (PBE0), -5 (B3LYP), -6 (B3PW91), -3 (BHandHLYP), 402 (BLYP), 1009 (PBE), 2523 (TPSS)。上面例子末尾第二個數字（5）不用動它，這是要求Link 608模塊也用Gaussian 16做SCF過程默認用的ultrafine積分格點。 要注意一些泛函的DFT-D3(BJ)參數在Gaussian里沒內置。比如TPSSh的D3(BJ)參數起碼到Gaussian 16 C.02仍然不是內置的。這種情況需要按照《DFT-D色散校正的使用》（http://www.shanxitv.org/210）里的做法自行定義。 為了用戶方便，在http://www.shanxitv.org/soft/sobEDA_tutorial.zip其中的templates目錄下已經提供了幾種常見的泛函的模板文件，比如template_BHandHLYP-D3(BJ).gjf就是對應BHandHLYP-D3(BJ)的情況，把它改名為template.gjf，并且修改里面的基組、電荷、自旋多重度為實際情況，就可以用了。 雖然sobEDA/sobEDAw是基于色散校正的DFT定義的，但也可以不加色散校正，即去掉em=gd3bj，此時輸出的色散校正項將為0。 3.7 常見問題 這一節回答幾個使用sobEDA.sh的常見問題。 (1)怎么讓sobEDA.sh計算時并行？ 《Gaussian的安裝方法及運行時的相關問題》（http://www.shanxitv.org/439）里說了，Default.Rou里設置比如-P- 96就可以默認用96核并行，-M-可以設置用的內存上限。也可以直接在template.gjf里用%nproc和%mem里設置并行核數和內存上限。 (2)為什么sobEDA.sh最后給出的總相互作用能和我自己用整體能量減去各個片段能量得到的不符？ 原理上是必然相符的，不相符有以下常見可能：用的幾何結構不一樣、用的基組不一樣、用的泛函不一樣、影響能量的某些設置不一致（比如溶劑模型、積分格點精度）、counterpoise校正考慮得不一樣、凈電荷和自旋多重度設置有差異、你自己做的或者sobEDA.sh自動做的某些單點計算恰好收斂到了不穩定波函數。仔細對比自己的Gaussian輸入輸出文件和sobEDA.sh產生的Gaussian輸入輸出文件必定能找到導致差異的原因。 (3)sobEDA和sobEDAw能帶溶劑模型么？ 本文2.2節一開始說了，sobEDA/sobEDAw是對真空下相互作用能做的分解，也顯然不能直接在template.gjf寫上scrf關鍵詞帶溶劑模型。想考慮溶劑模型要按JPCA原文2.1節末尾說的，自己計算各個體系在真空和溶劑模型下的單點能，通過簡單運算得到溶劑效應對相互作用能的貢獻。 (4)sobEDA.sh運行中途出錯怎么辦？ 首先要知道，sobEDA.sh按順序計算這些狀態：各個孤立片段（fragment[序號].gjf）、準分子態（promol.gjf）、凍結態（frozen.gjf）、以凍結態波函數為初猜的最終態（final.gjf）。只有fragment[序號].gjf和final.gjf在計算過程中需要做SCF直到收斂。腳本還利用Multiwfn產生準分子態的波函數、創建輸入文件等工作。 sobEDA.sh計算過程中在干什么、做到哪一步了，在輸出信息中顯示得非常清楚，從中可以很容易了解到是在哪個環節上出了錯。如果是調用Gaussian運行時出錯，應仔細看當前目錄下sobEDA.sh產生的相應任務的Gaussian的out文件末尾，結合Gaussian基本使用常識搞清楚原因。如果從提示中看到是Multiwfn執行出錯，>99%的可能性是Multiwfn沒有按手冊2.1.2節寫的正確安裝，導致無法正常啟動，或者由于可調用內存上限太少而無法執行。 絕對不能用太老版本的Multiwfn和Gaussian，sobEDA.sh需要用什么版本在本文3.1節明確說了。 必須確保機子里dos2unix命令能用，在本文3.1節明確說了。 sobEDA.sh中途涉及用Gaussian自帶的formchk和unfchk工具進行chk和fch之間的轉換。如果要算幾百原子的大體系，有可能由于自帶的這倆程序可以調用的內存上限不夠而無法轉換而造成任務失敗。解決方法是在用戶主目錄下的.bashrc里加上比如export GAUSS_MEMDEF=20GB然后重新進入終端，以使得這些Gaussian輔助程序能用的內存上限有20GB，這就足夠大了。 如果死活弄不清楚為什么出錯，又確信自己的輸入文件肯定沒毛病、Multiwfn自帶的例子都能正常跑完，可以在計算化學公社論壇http://bbs.keinsci.com的“量子化學”版發帖求助（選擇[綜合交流]分類），上傳所有給sobEDA.sh用的輸入文件。 順帶一提，sobEDA.sh腳本注釋得特別清楚，有點shell腳本常識就能很容易理解里面的內容，沒常識的話看《詳談Multiwfn的命令行方式運行和批量運行的方法》（http://www.shanxitv.org/612）里面關于shell腳本的相關介紹。sobEDA.sh做的幾個主要計算在sobEDA的JPCA原文7025頁左上角部分都有描述。用戶也可以結合實際需要對sobEDA.sh腳本進行適當的修改和調試。 (5)用Gaussian自帶的基組時sobEDA.sh能正常運行，但是使用自定義基組、混合基組，或者引用外部基組文件時sobEDA.sh中途出錯怎么回事？ 顯然沒定義對。一定記得先仔細看《詳解Gaussian中混合基組、自定義基組和贗勢基組的輸入》（http://www.shanxitv.org/60）了解怎么正確使用gen和genecp、格式是什么樣。sobEDA官方教程里“Cu+與H2O間的相互作用”都是用genecp的具體例子。如果是靠@來引用外部基組/贗勢文件，記得外部基組/贗勢文件末尾不要有空行，否則當此文件被展開后，相當于基組/贗勢部分末尾出現了多余的空行，導致Gaussian無法讀取到最末尾的比如-5 5這樣的信息。 (6)我用的是雙雜化泛函、M06-2X等方法算的弱相互作用能，但是sobEDAw沒有這些泛函的參數怎么辦？ 你可以用你想用的方法算相互作用能，而使用sobEDAw支持的級別做能量分解來得到不同物理成份對相互作用能的貢獻來了解不同物理成份起到的作用，這沒有任何不妥。sobEDAw支持的B3LYP-D3(BJ)算絕大多數弱相互作用足夠好，比如你用M06-2X算的相互作用能是-8.7 kcal/mol，而B3LYP-D3(BJ)算出來可能就是比如-8.1 kcal/mol或9.4 kcal/mol，顯然在B3LYP-D3(BJ)級別下做sobEDAw分析對相互作用能提供deeper insight是完全得當的。當然你也可以直接在文章里就拿B3LYP-D3(BJ)的相互作用能直接說事，就完全沒必要再單獨算M06-2X級別的相互作用能了。實際上如Mol. Phys., 115, 2315 (2017)的圖35所示，B3LYP-D3(BJ)算弱相互作用二聚體的相互作用能比M06-2X還更好。即便碰上比如18碳環這樣B3LYP-D3(BJ)描述不好的情況，還可以用sobEDAw直接支持的BHandHLYP-D3(BJ)，在Phys. Chem. Chem. Phys., 19, 32184 (2017)中證明了它算弱相互作用能也是很好的，如這篇文章補充材料的表20的intermol. NCIs（分子間相互作用能）所示，BHLYP-D3(BJ)（對應BHandHLYP-D3(BJ)）的誤差比B3LYP-D3(BJ)還稍微小點，和M06-2X在伯仲之間。 4 實例 sobEDA/sobEDAw官方教程http://www.shanxitv.org/soft/sobEDA_tutorial.zip里面有非常豐富而且講解得特別清楚易懂的例子，充分涵蓋了各種情況下的能量分解，全都弄懂了就可以把sobEDA.sh用得游刃有余。例子包括： ? 對H2O...NH3相互作用在B3LYP-D3(BJ)/6-311+G(2d,p) w.CB下做sobEDA能量分解 ? 對OC-BH3相互作用在B3LYP-D3(BJ)/6-311G*下做sobEDA能量分解 ? 對Cu+與H2O間的相互作用在TPSSh-D3(BJ)/SDD&6-311G*下做sobEDA能量分解 ? 對CH3-NH2相互作用在B3LYP-D3(BJ)/def2-TZVP下做sobEDA能量分解 ? 對水四聚體的四個分子間相互作用在BHandHLYP-D3(BJ)/6-31+G(d,p) w.CP下做sobEDAw能量分解 ? 利用腳本對CH3-CN相互作用在不同間距下做sobEDA能量分解，并且繪制成“相互作用成分 vs 片段間距”曲線圖 ? 利用腳本對苯...水相互作用在不同間距下做sobEDAw能量分解并且繪制成曲線圖 ? 獲得孤立狀態→準分子態、準分子態→凍結態、凍結態→最終態這三個過程中完整的能量成分變化 此外，教程文件包里還提供了sobEDA/sobEDAw原文里的諸多例子用的原文件，包括： ? 對18碳環...三重態氧氣間的相互作用在BHandHLYP-D3(BJ)/6-311+G(2d,p) w.CB級別下做sobEDAw能量分解 ? 對(CO)5Cr=CH2配位鍵相互作用在TPSSh-D3(BJ)/6-311G*&SDD級別下做sobEDA能量分解 ? 對C6F6...Cl-離子相互作用在BHandHLYP-D3(BJ)/6-311+G(2d,p)級別下做sobEDA能量分解 ? 對HC-CH乙炔中的三重鍵在B3LYP-D3(BJ)/def2-TZVP級別下做sobEDA能量分解 在下面，就只把H2O...NH3相互作用的sobEDAw能量分解這個最簡單的例子說一下，其它的請務必看sobEDA/sobEDAw官方教程里的講解（若在這里把所有英文例子翻譯成中文也沒什么意思）。此體系的結構圖如下所示，是經過幾何優化過的 H2O...NH3在B3LYP-D3(BJ)/6-311+G(2d,p) w.CB下做sobEDAw分析能量分解用到的所有文件在官方教程文件包里的HOH-NH3目錄下都提供了。system.xyz就是上圖對應的結構文件，是經過優化過的水與氨氣的二聚體極小點結構。template.gjf和本文3.6節示例的那個是一致的，不再累述。fragment.txt內容如下，就是定義兩個片段，也沒什么可說的 2 0 1 1-3 0 1 4-7 在此目錄下的sobEDA.sh里可見iCP=1、sobEDAw=1，說明要做sobEDAw分析并且考慮counterpoise。并且可以看到parm_c=0.550;parm_a=0.037;parm_r=2.750這一行開頭的#被去掉了，說明這次sobEDAw計算用這套c、a、r參數，對應于當前要用的B3LYP-D3(BJ)/6-311+G(2d,p) w.CB級別的參數。 把HOH-NH3目錄下的sobEDA.sh、system.xyz、fragment.txt、template.gjf都拷到當前目錄下，先運行chmod +x ./sobEDA.sh加上可執行權限，然后運行./sobEDA.sh |tee result.txt即可開始計算，并很快得到結果。結果既顯示在了屏幕上，也同時輸出到了result.txt里，如下所示 （計算過程輸出的細節信息略...） ************************* ***** Final results ***** ************************* Total interaction energy:     -7.18 kcal/mol Physical components of interaction energy derived by sobEDA: Electrostatic (E_els):    -12.15 kcal/mol Exchange (E_x):     -6.40 kcal/mol Pauli repulsion (E_rep):     20.28 kcal/mol Exchange-repulsion (E_xrep = E_x + E_rep):     13.88 kcal/mol Orbital (E_orb):     -5.38 kcal/mol DFT correlation (E_DFTc):     -2.70 kcal/mol Dispersion correction (E_dc):     -0.83 kcal/mol Coulomb correlation (E_c = E_DFTc + E_dc):     -3.53 kcal/mol 可見總相互作用能為-7.18 kcal/mol，同時sobEDA方法定義的各個物理成份都給出得相當清楚。為了便于討論，也把ΔE_x與ΔE_rep的加和值ΔE_xrep給出了、把ΔE_DFTc和ΔE_dc的加和值ΔE_c給出了。 然后看到的是sobEDAw的結果，如下所示。由于sobEDA和sobEDAw的相互作用能是相同的，差異僅在于物理成份的劃分上，因此總相互作用能就沒再輸出一遍。由提示可見96.26%的ΔE_DFTc與ΔE_dc相組合得到了sobEDAw的ΔE_disp項，即0.9626*(-2.70)-0.83 = -3.43 kcal/mol。其余的ΔE_DFTc與sobEDA的ΔE_xrep組合得到了sobEDAw的ΔE_xrep項，即(1-0.9626)*(-2.70)+13.88 = 13.78 kcal/mol。sobEDAw給出的ΔE_orb和ΔE_els與前面給出的sobEDA相應項完全一致。 Variant of sobEDA for weak interaction (sobEDAw): Note:   96.26% DFT correlation is combined with dispersion correction to yield a SAPT-like dispersion term Electrostatic (E_els):    -12.15 kcal/mol Exchange-repulsion (including scaled DFT correlation):     13.78 kcal/mol Orbital (E_orb):     -5.38 kcal/mol Dispersion (E_disp):     -3.43 kcal/mol 由以上數據可見，HOH...NH3氫鍵作用是靜電相互作用所主導的，ΔE_els占總吸引作用的-12.15/(-12.15-5.38-3.43)*100% = 58%。其次軌道相互作用也有不可忽視的貢獻，占總吸引作用的26%。色散作用對此氫鍵的形成貢獻最小。交換-互斥作用對結合起到了顯著不利因素，抵消掉了總吸引作用的66%，這是普遍現象。 腳本運行完畢后在當前目錄下還可以看到fragment1.gjf/out/fch、fragment2.gjf/out/fch、promol.gjf/out/chk/fch、final.gjf/out/chk這些中間過程產生的文件，現在已經完全沒用了，可以直接刪掉，感興趣的話也可以自己查看以了解腳本都調用Gaussian做了哪些計算。 值得一提的是，《透徹認識氫鍵本質、簡單可靠地估計氫鍵強度：一篇2019年JCC上的重要研究文章介紹》（http://www.shanxitv.org/513）里介紹的筆者的J. Comput. Chem., 40, 2868 (2019)一文中算的42種氫鍵復合物中也包含了HOH...NH3氫鍵復合物，文中在精度很高且非常昂貴的SAPT2+(3)δMP2/aug-cc-pVTZ下得到的相互作用能為-6.47，靜電、交換互斥、誘導（本質同軌道相互作用）、色散的貢獻分別為-11.70、12.63、-4.16、-3.23。當前sobEDAw計算得到的這5個值分別為-7.18、-12.15、13.78、-5.38、-3.43，可見相符程度極好！SAPT2+(3)δMP2給出的色散/靜電比例為0.276，此例sobEDAw給出的為0.282，相符巨好！ 5 總結 本文全面、完整介紹了筆者提出的基于很流行的色散校正DFT理論的sobEDA和sobEDAw能量分解方法，以及基于Gaussian和Multiwfn實現這兩種分析的sobEDA.sh腳本。sobEDA/sobEDAw方法分析準確度高、普適性強、耗時低，而sobEDA.sh用起來十分容易，對于量子化學初學者來說也沒有什么門檻，輸出信息頗為易于理解。sobEDA/sobEDAw結合sobEDA.sh可謂掃清了廣大量子化學研究者做能量分解的各種困難，相信在不久的將來就會成為很流行的方法，非常推薦在大家實際研究中廣泛使用sobEDA/sobEDAw并向同行積極推廣！ 發表的研究工作中使用sobEDA或sobEDAw方法時請記得引用方法原文J. Phys. Chem. A, 127, 7023 (2023)。由于sobEDA.sh腳本利用了Gaussian和Multiwfn程序，因此也請同時按照Multiwfn官方網站上或程序啟動時的提示正確引用Multiwfn的原文，以及引用Gaussian程序。 http://www.shanxitv.org/684 2023-08-14T07:59:00+08:00 理論設計新穎的基于18碳環構成的雙馬達超分子體系 文/Sobereva@北京科音   2023-Aug-14 1 前言 北京科音自然科學研究中心（www.keinsci.com）的盧天和江蘇科技大學的劉澤玉等人近期設計了由18碳環（cyclo[18]carbon, C18）和具有兩個大環的分子（OPP）構成的非常新穎的雙馬達超分子體系，并做了詳細的分子動力學特征的研究，相關工作已發表在知名的Chemical Communications通訊刊物上，文章信息如下，歡迎閱讀和引用： Zeyu Liu,* Xia Wang, Tian Lu,* et al., Theoretical design of a dual-motor nanorotator composed of all-carboatomic cyclo[18]carbon and a figure-of-eight carbon hoop, Chem. Commun., 59, 9770 (2023) DOI: 10.1039/D3CC02262E https://pubs.rsc.org/en/content/articlelanding/2023/CC/D3CC02262E 鏈接：https://pan.baidu.com/s/1H7AcSnDqbjBOE2ajFk-LXA?pwd=tnzv 下面將對此文的主要研究內容和思想進行深入淺出的介紹，以有助于讀者更好地理解文章內容。下圖是被研究的雙馬達超分子體系的結構圖和運動方式示意 值得一提的是，同作者之前還深入細致地考察了上文研究的OPP對18碳環的吸附和釋放的動力學過程及相互作用本質，已發表于Phys. Chem. Chem. Phys., 25, 16707 (2023) DOI: 10.1039/d3cp01896b，并在《8字形雙環分子對18碳環的獨特吸附行為的量子化學、波函數分析與分子動力學研究》（http://www.shanxitv.org/674）中做了詳細介紹，推薦閱讀。同作者之前還對18碳環及衍生物的各方面性質還做了非常廣泛的研究，相關工作匯總見http://www.shanxitv.org/carbon_ring.html。 2 設計思想 18碳環這個奇妙的分子雖然在半個多世紀以前就被理論預測，但直到2019年才首次在凝聚相實驗觀測到，并迅速引發了巨大關注。它的直徑和C60富勒烯很相似，對比見下圖 具有較大雙環結構的OPP分子在Angew. Chem., Int. Ed., 61, e202113334 (2022)中被首次合成出來，實驗上得到了它帶有C60和C70富勒烯的晶體，C60@OPP的晶體結構如下所示，可見C60被嵌入在了大環里面 由以上結構特征可以想到兩個C18也可以分別內嵌到OPP的兩個大環中，成為2C18@OPP。倘若每個C18在OPP的大環里能夠由于熱的驅動而較容易地發生轉動，自然就成了雙馬達超分子體系。這樣的體系在之前的文章中是沒有報道的，此體系或類似物在未來有望成為構造復雜分子機器的關鍵組成部分。本文介紹的Chem. Commun.這篇文章的目的就在于通過量子化學和分子動力學模擬來證實OPP與C18復合作為雙馬達體系的可能性并考察其特點。 3 基于量子化學的能量角度的研究 文中首先使用量子化學方法從能量角度對C18@OPP和2C18@OPP進行了研究。首先使用Gaussian程序通過ωB97XD泛函優化了復合物結構并做了振動分析。2C18@OPP達到了260個原子，為了節約時間，對占體系大部分原子的OPP部分用了6-31G*基組，而關鍵的C18部分用了更好的6-311G*基組（在http://www.shanxitv.org/584中也指出6-31G*無法合理描述C18）。對相應單體也在對應的級別做了優化和振動分析。之后通過ORCA程序使用ωB97X-V/def2-TZVP級別結合counterpoise校正計算了它們的結合能（做法可參考《在ORCA中做counterpoise校正并計算分子間結合能的例子》http://www.shanxitv.org/542），發現C18 + OPP以及C18 + C18@OPP的結合能分別為-18.4和-18.5 kcal/mol，一方面說明C18與OPP有較強的內在結合能力，另一方面說明C18@OPP的已進入的一個C18幾乎不影響OPP對第二個C18的結合能力。之后，文中通過《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》（http://www.shanxitv.org/552）介紹的Shermo程序計算了C18、OPP、C18@OPP和2C18@OPP常溫下自由能熱校正量，并進而得到了常溫下的C18 + OPP和C18 + C18@OPP的結合自由能，分別為-6.0和-4.0 kcal/mol，明顯為負值說明C18的進入是熱力學上可自發的，其大小明顯小于結合能是來自于分子間復合帶來的熵罰效應。 審稿人問及C18是否有可能與OPP存在其它結合方式。為了說明這一點，此文的補充材料里給出了優化出的其它兩個C18與OPP復合物的極小點結構，如下所示，相對于C18平行地嵌入在OPP大環內的結構的能量差也在下圖給出了。可見其它兩種結構都是能量更高、明顯更不穩定的。在《8字形雙環分子對18碳環的獨特吸附行為的量子化學、波函數分析與分子動力學研究》（http://www.shanxitv.org/674）介紹的PCCP文章中實際上也已經通過分子動力學模擬證明了C18沒有其它的能夠與OPP穩定結合的方式。 文中進一步研究了C18的旋轉勢壘。對C18@OPP和2C18@OPP優化的C18旋轉的過渡態結構如下(a)所示，可見虛頻數值非常小，暗示旋轉勢壘肯定很低。圖中的箭頭體現了虛頻方向，可見明顯對應的是C18的整體旋轉。此圖是按照《在VMD中繪制Gaussian計算的分子振動矢量的方法》（http://www.shanxitv.org/567）介紹的方式繪制的。 進一步，文中對C18在OPP大環中的旋轉進行了掃描（計算級別同幾何優化），如上圖(b)所示，其中橫坐標是旋轉角度，0處為過渡態結構。可見旋轉勢壘非常低，只有0.13 kcal/mol，這體現C18在OPP大環中的旋轉必然很容易靠熱運動來驅動。上圖橫坐標范圍對應一個旋轉周期，是40度，這也正對應于18碳環的極小點結構是D9h的特征。圖中看上去每20度有一個極大點，這是因為18碳環是18個原子構成的環狀體系，如果忽視其長-短鍵交替特征的話，就相當于20度一個旋轉周期。 4 分子動力學模擬C18在OPP中的運動行為 分子動力學模擬研究是本文最關鍵的部分，模擬出的動力學行為是C18與OPP復合物能否作為分子馬達的決定性證據。此文使用GROMACS程序通過經典力場進行動力學模擬，使用靈活的sobtop程序（http://www.shanxitv.org/soft/Sobtop）構建拓撲文件，主要基于GAFF力場，部分參數利用OPP的Hessian矩陣由sobtop產生，相關細節在《8字形雙環分子對18碳環的獨特吸附行為的量子化學、波函數分析與分子動力學研究》（http://www.shanxitv.org/674）里有描述，這里就不多提了。 2C18@OPP處于300 K的平衡狀態時的200 ps動畫如下所示（是在文章的補充材料里提供的）。為了能讓C18的旋轉特征看得清楚，其中一個原子以紅色來標記。由動畫可見，的確C18在OPP中能夠非常順利、自如地轉動起來，充分證實了此文提出的雙馬達超分子體系的設想！ http://www.shanxitv.org/attach/684/2C18_OPP_rotator.mp4 （點擊查看動畫） 需要注意的是，由于原子間頻繁碰撞導致的動能交換，單靠熱運動驅動的C18在OPP中的旋轉并非始終是單向的。300 K下模擬的50000幀（0.2 ps保存一幀，故相當于10 ns）中的2C18@OPP中兩個C18的旋轉角速度如下所示（數據已通過Savitzky-Golay算法利用相鄰1000個數據點做了平滑化處理以消除噪音）。由圖可見C18的含符號角速度不斷發生變化，時正時負，時大時小，因此C18的旋轉不能視為是連貫、無摩擦的。從圖上兩條曲線也可以看到2C18@OPP的兩個C18彼此間的運動并沒有什么相關性。由于它們離得遠，自然也不會有什么明顯的直接的耦合或者藉由OPP的結構產生的間接耦合。 5 分子動力學模擬C18在OPP中轉動的統計行為 文中在50、100、200、300、400 K的熱浴下都做了100 ns的分子動力學模擬，并通過自寫的軌跡分析程序考察了C18在OPP中轉動的統計行為，結果如下表所示。在溫度不是很低情況的模擬過程中，偶爾C18會傾斜地倚靠在OPP的內側，此時的體系明顯不算是轉子。這種情況的幀數占總模擬幀數的百分比是下表的p_out，這些幀不納入C18轉動行為的統計。可以看到從200 K開始，隨著溫度上升，處于這種亞穩狀態的比例隨之上升。 上表中ω是平均無符號轉動角速度，可見隨著溫度越高、體系中原子的平均動能越大，C18平均旋轉速度也隨之上升。但是300 K變化到400 K過程中平均旋轉速度增大甚微，這主要在于400 K的時候C18在OPP中的狀態已經很不穩定了，過強的熱運動很大程度影響了C18轉動的有序性。上表中f=ω/(2π)是換算出來的轉動頻率，2C18@OPP屬于亞THz范圍的超快轉子。 上表中αavg、αSD分別是C18環平面與OPP大環平面之間在模擬過程中的夾角的平均值和標準偏差，d_avg和d_SD分別是C18環中心與OPP環中心之間在模擬過程中的距離的平均值和標準偏差。從這些指標來看，隨著模擬溫度提升得越高，2C18@OPP偏離其極小點結構（α和d都近乎為0）的程度越高，越缺乏理想轉子特征。 下圖展現了VMD繪制的2C18@OPP在100 ns動力學過程中的軌跡疊加圖，每1 ns繪制一次，根據幀號從前到后按照紅-白-藍著色。可見在100 K的時候由于熱運動弱，體系結構基本上就是在極小點結構附近很小范圍波動。到了更高的200 K后OPP的運動幅度以及C18的運動范圍都顯著增加了。在300 K的時候結構波動更大，還明顯看到C18都偶爾側貼在OPP大環內側了。到了400 K的時候C18在OPP大環里的運動其實已經很大程度算是無序了，甚至仔細看軌跡都會發現C18在大環內出現了整體翻轉。可以預期如果溫度明顯更高的話，C18甚至就能從OPP的大環中跑出去了，比如臨時性跑到OPP外側乃至飛走。 上圖的(b)和(c)對不同溫度下模擬的C18在OPP中的結構狀態做了統計分析，包括C18與OPP大環間的夾角以及中心距離，可見這兩個衡量C18在OPP中結構狀態的參數的分布都是隨著溫度增加而顯著變寬，也即在模擬過程中C18的朝向和位置的波動程度愈發增大，這和軌跡疊加圖展現的信息是一致的。 6 總結 本文介紹了近期發表于Chem. Commun., 59, 9770 (2023)的理論設計的雙轉子超分子復合物結構，此工作從能量和動力學角度對體系進行了全面的研究，充分證明了此體系作為超快納米轉子的能力，并對其轉子性能以及溫度依賴性做了系統的考察。由于OPP和C18都已經被合成出來，因此此體系或者其變體有望成為有實際價值的分子機器的組成部分，本文的研究思想和發現對于納米轉子類體系在未來的探索也有顯著的啟示作用。 http://www.shanxitv.org/680 2023-08-05T01:25:00+08:00 辨析計算化學中的任務類型和理論方法 文/Sobereva@北京科音 First release: 2023-Aug-5    Last update: 2024-Feb-23 0 前言 在網上答疑時，偶爾看到有初學者搞不清楚任務類型和理論方法，比如今天有人在思想家公社3號群問“md，aimd和qmmm的區別是啥啊？什么需求下會用到這三種呢？”，這三個詞明顯都不是一個層面的東西。此文就把計算化學中的任務類型和理論方法的關系明確一下，簡明扼要地介紹一下相關基本概念，做一個科普，以令計算化學零基礎者一次性搞清楚它們的關系。更多的計算化學的總覽性的知識在《北京科音初級量子化學培訓班》（http://www.keinsci.com/workshop/KEQC_content.html）里有介紹，此培訓也是從零開始入門量子化學研究的最好方式。 1 任務類型 有N個原子的非線型的體系有3N-6個內坐標描述原子間相對位置關系，在Born-Oppenheimer近似下體系的能量是依賴于內坐標的，也即能量是個3N-6維的函數，這被稱為勢能面（potential energy surface, PES）。計算化學領域有很多常見的任務（task）都是基于勢能面做的，即有了勢能面信息（能夠求得勢能面上任意位置的能量及其導數）就能做這些任務。下面羅列一下常見的這種類型的任務： ? 優化極小點：平時說的幾何優化（geometry optimization）一般也是指這個。此任務從一個給定的初猜結構開始，根據特定算法去尋找與之最近的勢能面上極小點的精確位置。在分子動力學程序如GROMACS里這種任務也被叫做能量極小化（energy minimization, em），只不過實際目的不一樣，em的目的主要是在動力學模擬之前釋放體系中可能存在的較大斥力（自行構建的初始模型里往往有原子間距離太近、斥力太大），免得動力學模擬一開始由于過大的原子加速度造成過大的速度而導致動力學模擬異常或崩潰。 ? 優化過渡態：從一個給定的初猜結構開始，根據特定算法去尋找與之最近的勢能面上的一階鞍點的精確位置。過渡態可以視為反應過程中最有代表性的一個結構，可以由此判斷或驗證反應機理，利用它和相鄰極小點的能量求差可以得到反應勢壘，對于討論反應難易有關鍵意義，見《談談如何通過勢壘判斷反應是否容易發生》（http://www.shanxitv.org/506）。優化過渡態的算法介紹見《過渡態、反應路徑的計算方法及相關問題》（http://www.shanxitv.org/44）。 ? 產生反應路徑：用于把基元反應對應的勢能面上兩個相鄰極小點之間最容易相互轉變的路徑產生，也相當于得到了一個軌跡并可以觀看對應的結構變化過程，過渡態是此路徑上能量最高點。在質量權重坐標下產生的反應路徑稱為IRC（intrinsic reaction coordinate），不考慮質量權重時產生的一般稱為MEP（minimum energy path）。具體算法很多，見《過渡態、反應路徑的計算方法及相關問題》（http://www.shanxitv.org/44）。反應路徑可認為是實際化學反應過程最有代表性的路徑，自然對于理解反應機理有重要意義，還可以用Multiwfn對反應路徑上各個點做波函數分析考察反應過程中電子結構的變化以獲得更深入的信息，例如《通過鍵級曲線和ELF/LOL/RDG等值面動畫研究化學反應過程》（http://www.shanxitv.org/200）。 ? 振動分析：通常是在勢能面上的駐點（所有原子受力都為0的點）做的，用于得到相應結構下的振動頻率，可以用來計算熱力學量，公式看《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》（http://www.shanxitv.org/552）介紹的Shermo程序手冊的附錄部分，還可以用來得到振動能級、振動波函數、檢驗幾何優化的準確性（根據虛頻判斷）、繪制振動光譜（參看《使用Multiwfn繪制紅外、拉曼、UV-Vis、ECD、VCD和ROA光譜圖》http://www.shanxitv.org/224）。 ? 分子動力學（molecular dynamics, MD）：上面的任務都是“靜態”的任務，即不含時間這個維度。而分子動力學則引入了時間這個維度，可以模擬體系狀態隨時間的變化，能研究的問題與那些靜態的任務有極大的互補性，在北京科音分子動力學與GROMACS培訓班（http://www.keinsci.com/workshop/KGMX_content.html）里面有全面的介紹。分子動力學模擬過程相當于體系在勢能面上不斷運動，也相當于是對勢能面的一種采樣方式。雖然其名字叫分子動力學，但研究對象并不限于分子，比如無機固體、金屬團簇的動力學模擬也可以叫分子動力學。最常見的分子動力學的實現形式是BOMD（結合量子化學方法時另有CPMD、ADMP等），也相當于根據原子的位置、速度和受力按照經典牛頓運動方程演化原子的坐標和速度。還有其它一些特殊的動力學形式，如朗之萬動力學、耗散粒子動力學等，如今用得很有限。 ? 蒙特卡羅：和分子動力學并列的另一種常見的對勢能面采樣的方法，適用場合相對更少，主要是在多孔材料吸附小分子、計算氣液共存曲線等問題上用得較多。不需要像MD那樣計算受力和速度信息，也沒有明確的時間的概念、無法像MD那樣嚴格地考察含時演化。 ? 構型/構象搜索：勢能面上往往有很多極小點，能量最低的那個是全局最小點，可以視為體系最穩定的構型或構象對應的結構，其它的極小點對應亞穩的構型/構象。前面說的幾何優化只能收斂到與初猜結構最近的極小點，顯然這未必是最小點。如果想確保得到最小點（全局優化，global optimization），或者能量最低的一批構型/構象，就需要做構型/構象搜索，方法很多，比如免費的molclus（http://www.keinsci.com/research/molclus.html）就是很常用的實現工具。 2 理論方法 上一節介紹的那些任務在進行過程中，至少需要計算體系在不同結構下的能量。有的還需要計算受力（勢能面對核坐標一階導數矢量，或者說梯度，的負矢量），如幾何優化、分子動力學。有的還需要計算能量對核坐標的更高階導數，比如振動分析至少需要算二階導數（Hessian矩陣由之構成）。怎么獲得各種任務中涉及的能量及其對核坐標的導數，就是理論方法層面的問題了。下面簡要列舉一些常見的理論方法： ? 分子力場（molecular forcefield）：也稱為經典力場（classical forcefield）或分子力學（molecular mechanics, MM），或者經驗勢函數等等。雖然通常帶著“分子”倆字，但實際中不限于用于分子體系。這類方法一般使用簡單的模型描述體系（通常一個原子作為一個粒子，電子與原子核不分開考慮），并用形式簡單的數學函數（勢函數）描述原子間相互作用。例如計算原子間靜電相互作用時，分子力場普遍是給每個原子核位置分配一個點電荷（數值等于原子電荷），然后基于庫侖公式計算。由于分子力場的形式簡單，因此計算耗時極低。分子力場的一個關鍵的局限性是絕大多數分子力場都無法描述化學反應，因為成鍵關系在計算從始至終是固定不變的，是一開始就定義好的。描述化學反應問題主流的是后面提及的量子化學類型的方法，或者反應力場（遠比普通分子力場要貴、支持的程序要少）。分子力場另一個關鍵局限性是依賴于大量參數，一方面影響模擬精度，一方面決定能描述的體系，尋找和構建合適的參數往往不是易事。 ? 機器學習勢（machine learning potential）：基本思想是通過機器學習的思想構造依賴于原子坐標的分子描述符（如距離矩陣、能量矩陣等）與能量之間的抽象關系，這種關系比上述傳統的分子力場用的勢函數形式明顯更為復雜。 以下提及方法都是量子化學（quantum chemistry, QC）范疇的，也可以叫基于量子力學（quantum mechanism, QM）的方法。 ? 密度泛函理論（DFT）：是目前在量子化學、第一性原理計算領域用得最普遍的理論方法，因為其性價比非常高，交換-相關泛函（影響DFT計算精度的決定性因素）選擇得當的話可以得到很不錯的結果，足夠滿足絕大多數研究需要，而且精度顯著強于分子力場（除非力場是對某類體系專門訓練得特別精妙的力場），普適性強得多得多，但耗時高于分子力場N個數量級，體系越大N越大。DFT應用到實際問題中實際上分為Kohn-Sham DFT（KS-DFT）和orbital free DFT兩種形式，前者是絕對的主流，計算過程牽扯到軌道波函數，本文以及大家平時說的DFT計算都是用的KS-DFT形式；而后者非主流，不牽扯到軌道，應用的局限性非常大，只有很少數程序支持。 ? Hartree-Fock（HF）：在DFT興起之前用得很普遍，如今已經完全過時。耗時不比DFT顯著更低（DFT計算時利用RI技術時反倒還可以明顯更快），而精度差得多。 ? 后HF：是一大類方法的統稱，如CCSD(T)、MP2、CISD等。HF方法精度爛在于同時沒考慮到動態相關和靜態相關。后HF在HF基礎上額外把動態相關在一定程度上考慮了進來，但耗時也高了很多。動態相關、靜態相關的相關概念很有學問，在北京科音基礎量子化學培訓班（http://www.keinsci.com/workshop/KBQC_content.html）里有詳細講授，這里就不多提了。 ? MCSCF：主要彌補HF缺乏對靜態相關的考慮。但由于幾乎沒有或很少考慮動態相關，所以結果也說不上好。MCSCF最常見的實現是CASSCF。 ? 多參考方法：在MCSCF基礎上進一步把動態相關考慮進來，精度整體很好，普適性很強，但使用相對復雜且昂貴。典型的實現如CASPT2、NEVPT2、MRCI。 ? 半經驗：一大類方法的統稱，都是對HF的簡化以巨幅降低耗時，耗時只是HF的微小零頭，整體精度也有不小降低，而且引入了依賴于元素的參數而只能用于有限的元素。典型的如AM1、PM3、PM6。一些較復雜的機器學習勢的耗時已經追平了半經驗。 ? GFN-xTB：相當于納入了一定DFT思想的半經驗級別的方法，耗時和半經驗相仿佛，而整體精度和可靠性>=主流的半經驗方法，但跟DFT比還有很大差距。GFN-xTB從2017年開始發展，和歷史更早的DFTB有極大的共性。 前面介紹的HF、后HF、MCSCF、多參考方法里面都完全沒有引入經驗參數（極個別的除外，諸如CASPT2可以涉及shift參數），計算中的參數只依賴于物理常數，故它們被稱為“從頭算（Ab initio）”類型的方法。分子力場、半經驗、GFN-xTB由于引入了大量經驗參數，因此明顯不屬于從頭算。DFT算不算從頭算模棱兩可。如今尚未找到的理論上嚴格的交換-相關泛函是不含任何參數的，但如今只有近似的交換-相關泛函被提出，其中雖包含經驗參數，但整體相對較少，因此爭論DFT是否屬于從頭算沒意思，只是個說法問題。有的泛函的經驗性較小，有的經驗性則較高。如今也有的文章將HF、后HF、MCSCF、多參考方法等物理思想完全基于波函數的方法統稱為wavefunction theory (WFT)，在討論時使用WFT一詞時DFT就不會和那些方法混淆了。 還有一個詞叫“第一性原理（first-principle）”，其字義和從頭算一樣，但“第一性原理計算”如今的主流指代是“量子化學方法研究周期性體系”，不僅包括嚴格的從頭算方法，也包括DFT，因為這類計算最常用的就是DFT。DFTB雖然依賴很多參數，但DFTB計算周期性體系也往往被不計較地算作第一性原理計算。 由于有了“第一性原理計算”這個專門描述把量子化學方法用于周期性體系的計算的詞，因此平時說的“量子化學計算”、“從頭算方法計算”一般是指計算分子、團簇這樣孤立（非周期性）體系的情況。那些通過量子化學方法主要研究孤立體系的程序被稱為量子化學程序，如Gaussian、ORCA、GAMESS-US等，而通過量子化學方法主要研究周期性體系的程序被稱為第一性原理程序，如CP2K、Quantum ESPRESSO等。注意“量子計算”和量子化學計算又不是一碼事，前者強調的是利用在量子計算機上實現的量子算法，可以但絕不限于做量子化學問題的研究。 接下來再說很知名的QM/MM。這個方法是指使用量子化學方法以量子力學（QM）的思想描述體系的一部分，用分子力場以分子力學（MM）的思想描述體系的另一部分，并且恰當考慮兩部分之間的耦合。對于很大的體系，這顯然比起全都用量子化學方法描述要便宜得多得多（MM部分計算耗時相當于QM部分的九豬一毛），但代價是MM區域描述的精度比QM區要差，而且絕大多數力場不能描述這部分發生的成鍵/斷鍵過程。因此，必須恰當劃分QM和MM區，讓最需要精確描述的部分或者涉及化學反應的原子納入到QM區中，而起到較次要作用的“環境”原子納入到MM區中。 還有一個詞是ONIOM，也是劃分不同區域，不同方法著重描述不同部分，它可以把任意兩種理論方法相組合。如果把量子化學方法和分子力場相組合，特稱為ONIOM(QM:MM)，其用處和QM/MM類似但實現不同，描述時絕對不能搞混。一些相關討論見《要善用簇模型，不要盲目用ONIOM算蛋白質-小分子相互作用問題》（http://www.shanxitv.org/597）。 再說一下對電子激發態的描述。前面介紹的半經驗、HF、后HF、DFT等方法通過delta-SCF方式可以算激發態，但不夠普適。量子化學領域有很多專門的方法可以計算激發態，包括計算激發態的能量及其導數，如最常用的TDDFT、過時的CIS、較高精度的EOM-CCSD等，見《亂談激發態的計算方法》（http://www.shanxitv.org/265）。幾乎所有分子力場都是描述基態電子態的，但如果專門對激發態勢能面擬合力場參數，也可以用力場描述激發態，但已有的這樣的力場都是針對極其具體的體系的激發態擬合的，不具備普適性。機器學習勢同理。QM/MM也可以描述激發態，但僅限于電子激發完全發生在QM區的情況。 再簡單說一下和實際計算關系密切的基函數（basis function）的概念。本節介紹的理論方法，凡是基于量子理論思想提出的（也就是前面說的那些里面除了分子力場和機器學習勢以外的方法），在實際數值求解的過程中一般都要涉及分子軌道，關于分子軌道的概念詳見《正確地認識分子的能隙(gap)、HOMO和LUMO》（http://www.shanxitv.org/543）。絕大多數計算程序中都是把分子軌道展開成基函數的線性組合來描述的。最常見的基函數一類是原子中心基函數（如Gaussian等大多數量子化學程序以及CP2K等部分第一性原理程序用的高斯型基函數），其中心一般位于原子核，還有一類常見的基函數是平面波基函數，是大多數計算周期性體系為主的第一性原理程序用的，它的分布覆蓋整個被計算的晶胞。基組（basis set）是對于原子中心基函數而言的，例如6-31G*、def2-TZVP、cc-pVDZ等都是很常見的基組，它定義了實際計算時對各種元素原子具體用多少、什么參數的基函數。做HF、DFT、后HF、MCSCF、多參考等方法計算時都需要告訴計算程序用什么基組，基組質量越好，也就是越接近于所謂的完備基組極限，這些方法自身的精度發揮得就越充分，但代價就是耗時越高。一個好方法配一個爛基組，以及一個爛方法配一個好基組，結果都不理想，必須好方法配好基組才能得到較準確結果。半經驗、GFN-xTB方法計算時也利用基函數，但是用什么基函數是這類方法自身定義的，就不需要而且也不能再指定基組了。北京科音初級量子化學培訓班（http://www.keinsci.com/workshop/KEQC_content.html）和北京科音基礎（中級）量子化學培訓班（http://www.keinsci.com/workshop/KBQC_content.html）都專門有一節對基組做詳細介紹，后者講得明顯更深、更廣。 3 任務類型與理論方法的結合 第2節介紹的各種理論方法原則上都可以用于第1節介紹的各種任務。例如，分子力場、DFT、QM/MM等理論方法全都可以用于幾何優化、振動分析、產生反應路徑等任務。這些組合在程序實現上一般沒有什么需要特殊考慮的，一個程序支持什么理論方法，就能基于它們做什么依賴于勢能面的任務。 搞明白了前面介紹的名詞，從頭算動力學（ab initio molecular dynamics, AIMD）和第一性原理動力學（first-principle molecular dynamics, FPMD）顧名思義自然就知道是描述什么計算了。這兩個詞其實沒有嚴格的區分，實際中經常混淆使用，但從嚴格的習俗來說，AIMD和FPMD分別最適用于描述使用量子化學方法對孤立體系和周期性體系做分子動力學模擬。所以北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）里面講CP2K做基于GFN-xTB和DFT的分子動力學的地方標題用的是第一性原理動力學，但若說成AIMD也完全可以，不會造成什么誤會。 由于基于分子力場做分子動力學的文章數目和流行程度遠高于非常昂貴的AIMD/FPMD，因此平時說的分子動力學模擬默認是指的基于經典力場的分子動力學，這也往往被叫做經典分子動力學。但要注意“經典”這個詞在不同語境下有不同含義，有的場合下是特指把原子核當做經典粒子來考慮，此時基于BOMD形式做基于量子化學方法的分子動力學也是經典分子動力學。而與之相對的是“量子動力學”，強調把原子核以量子力學方式描述，原子核有對應的波函數。顯然，此處“量子動力學”和“基于量子化學做分子動力學”又不是同一個概念。 QM/MM MD指什么也不言而喻，顯然是用的QM/MM理論方法做分子動力學任務，這結合了QM/MM對體系描述的特點和分子動力學研究問題的能力。顯然這比所有原子都用QM描述（相當于AIMD/FPMD）要便宜得多得多，當然代價就是MM區域在動力學過程中描述得相對較糙，而且預期不會發生化學反應，同時還得想辦法去找力場參數。因此能否用QM/MM MD代替AIMD/FPMD當然得根據被研究的體系和問題判斷，盲目瞎用純粹是自討苦吃（吐槽一下，很多初學者沒有基本理論常識，一看文獻里有QM/MM就覺得好像高大上似的，就整天想著QM/MM而不知道其局限性，太naive了）。 文獻里的許多用詞并不標準，必須結合語境和常識理解說的是什么。比如可能有的文章沒明確用QM/MM MD這個詞而只說了QM/MM，若實際上在此基礎上跑了動力學，當然就是QM/MM MD。再比如有的文章表面上說是做的AIMD，但從計算細節描述來看作者把一部分當做MM來描述以節約時間，顯然實際上也是QM/MM MD。理解名詞必須有基本的應變能力。 如果用的理論方法描述的是激發態，那么前述各種任務就是對激發態做的。比如用TDDFT做激發態分子動力學，那就是動力學的每一步用的原子受力都是用TDDFT算的激發態的受力。順帶一提，有些文獻里聲稱做的是激發態動力學研究，但其實內容里根本沒做分子動力學，只是對激發態計算了勢壘，文中這么說只是由于勢壘和研究反應速率的反應動力學領域里的反應速率常數密切相關。另外值得注意的是研究激發態動力學往往需要考慮內轉換、系間竄越等勢能面切換的過程，考慮這些需要做非絕熱動力學，需要額外的算法，這方面的坑很深。所以并不是一個程序支持了計算激發態的理論方法和分子動力學模擬后就能完美、嚴格地做激發態動力學研究相關問題，情況遠沒那么簡單。 http://www.shanxitv.org/679 2023-08-03T02:19:00+08:00 江蘇科技大學的劉澤玉等人和北京科音自然科學研究中心的盧天等人近期發表了一篇非常系統性的通過理論計算研究碳鏈體系的電子光譜文章： Jiaojiao Wang, Zeyu Liu,* Qing Zhou, Tian Lu,* Xia Wang, Xiufen Yan,* Mengdi Zhao, Aihua Yuan, Size dependence of electronic spectrum for H-capped carbon chains, H-(C≡C)n-H (n = 3–9, 15): Analysis of its nature and prediction for carbyne, Comput. Theor. Chem., 1227, 114255 (2023) DOI: 10.1016/j.comptc.2023.114255 https://pan.baidu.com/s/1S8-wURaxp9_VEbUN1q-H_g?pwd=x2e1 下文是此文的中文版，歡迎閱讀和引用。本文研究的碳鏈體系和18碳環（cyclo[18]carbon）有非常緊密的聯系，18碳環的相關理論研究工作匯總見http://www.shanxitv.org/carbon_ring.html 文中使用的LOL-pi分析的介紹見《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432），文中用的空穴-電子分析介紹見《使用Multiwfn做空穴-電子分析全面考察電子激發特征》（http://www.shanxitv.org/434），文中用的BLA、BOA分析見《使用Multiwfn計算Bond length/order alternation (BLA/BOA)和考察鍵長、鍵級、鍵角、二面角隨鍵序號的變化》（http://www.shanxitv.org/501）。 氫封端碳鏈H-(C≡C)n-H (n = 3-9, 15)的電子光譜的尺寸依賴性：性質分析及對碳炔的預測 摘要 碳鏈的電子光譜的本質迄今還沒有得到很好的理解，值得進一步研究。在本工作中，采用（含時）密度泛函理論TDDFT研究了H-(C≡C)n-H (n = 3-9)的幾何結構、電子結構以及11Σu+ ← X1Σg+電子躍遷。計算的幾何參數和垂直激發能與已知的實驗觀測完全一致。通過對碳鏈電子激發的波函數分析，深入解釋了碳鏈的電子和光學特性的尺寸依賴性，并對它們進行了數據擬合和外推，通過n = 15的特定分子證實了擬合公式的準確性。根據提出的外推公式，可以估算出碳鏈分子電子和光學特性在鏈長增加時呈現出的特性，甚至獲得sp雜化碳同素異形體碳炔的極限值。 關鍵詞：碳鏈，聚炔，碳炔，電子光譜，尺寸依賴性 1. 前言 碳鏈是由sp雜化的碳原子連接的一維線性分子，通式為R-(C≡C)n-R。氫封端的H-(C≡C)n-H是各種末端取代的碳鏈的最簡單形式。在文獻中，它們有時被稱為聚乙炔，這是錯誤的，因為分子式為(H-C≡C-H)n的聚合物才是聚乙炔，而H-(C≡C)n-H的確切名稱應該是聚炔。 早期的探索發現，氫封端的碳鏈與密集星際云化學有關，這些高度不飽和的聚合物對星際介質的天體化學演化有顯著貢獻[1]。此外，還檢測到它們是形成多環芳烴和富勒烯的中間體[2]。最近通過凝聚相實驗[3，4]成功合成并表征了環[18]碳，這喚起了碳鏈的另一個長期隱藏的身份，即它們可能是此類化學性質不穩定的碳原子環的開環產物[5]。在對環[18]碳及其類似物和衍生物[6-18]進行理論研究的同時，我們也對不同長度碳鏈的性質感興趣。由于環狀碳簇的電子光譜已經得到了廣泛研究[6，8，10]，具有類似電子結構的鏈狀分子H-(C≡C)n-H的相關性質就成為一個值得研究的課題。 在氣相（n = 1-13）[19-21]、溶液（n = 2-10，12）[19，22]和氖（n = 6-12）[23]等不同介質中的實驗已經觀察到了幾種H-(C≡C)n-H鏈的吸收光譜，并證明了它們的電子和光學特性與尺寸有關。此外，還有對氫封端碳鏈的結構特征、電子激發能和振子強度的尺寸依賴性進行的一些理論研究[24-27]。但是，也許是受當時分析手段的限制，幾乎所有的報道都僅僅是關于吸收帶的位置和分布的，對碳鏈尺寸依賴的光學特性還缺乏全面的理解和分析。此外，通過外推有限長度碳鏈性質的尺寸依賴性來預測假想的碳同素異形體碳炔的相關性質具有重要的理論和實際意義。 延續我們對具有特殊拓撲結構的材料分子的研究興趣[6-18,28-32]，本文詳細考察了碳鏈H-(C≡C)n-H (n = 3-9, 15)在氣相中的光學激發性質，并在電子結構水平上揭示了其尺寸依賴性的本質。所提出的電子和光學特性的擬合公式可以推廣到更長的碳鏈上，并有望用于預測sp雜化碳同素異形體碳炔的性質。 2. 計算細節 我們最近的研究發現，ωB97XD/def2-TZVP水平可以很好地再現實驗觀測以及UCCSD/def2-TZVP水平上[7,9,16]計算得到的環[18]碳的幾何結構，這也保證了其對當前具有相同sp雜化構型的體系的優化具有準確性和穩健性。因此，本工作在氣相條件下，使用ωB97XD交換相關泛函[33]和def2-TZVP基組[34]，通過密度泛函理論（DFT）對碳鏈H-(C≡C)n-H (n = 3-9, 15)的幾何結構進行了優化。所有優化的幾何結構都是鍵長交替的線型結構，并且是沒有虛頻的勢能面極小點。碳鏈的激發能通過含時密度泛函理論（TD-DFT）計算得到，與幾何優化的計算水平ωB97XD/def2-TZVP相同。 上述所有量子化學計算均采用Gaussian 16 (A.03)程序包[35]完成。此外，我們還利用ORCA軟件[37]進行了高精度的STEOM-DLPNO-CCSD/def2-TZVP計算[36]以驗證TDDFT電子激發計算的可靠性。本文還使用Multiwfn 3.8(dev)進行了電子結構分析[38]。基于Multiwfn導出的文件，通過VMD軟件[39]繪制了分子軌道（MO）和實空間函數的等值面圖。 3. 結果與討論 3.1 幾何和電子結構 碳鏈H-(C≡C)n-H具有11Σg+基態，且具有閉殼層電子構型[20,25]。在ωB97XD/def2-TZVP水平上優化的碳鏈H-(C≡C)n-H (n = 3-9)基態幾何結構如圖1所示，H-(C≡C)n-H (n = 3-9)的笛卡爾坐標見表S1。 圖1. 在ωB97XD/def2-TZVP水平上優化的碳鏈H-(C≡C)n-H (n = 3-9)的基態幾何結構。圖中還顯示了計算的（括號外）和實驗的（括號內）[40]鍵長（單位：?）。 如圖1所示，基態碳鏈H-(C≡C)n-H (n = 3-9)是具有中心對稱性的長短C-C鍵交替的線型結構，屬于D∞h點群。這種鍵長交替（BLA）的結構被以前關于碳鏈體系的光學和理論研究所認可[21,26,27,40-44]，這種對稱性降低的變形可以理解為Peierls畸變的結果[26,45]。從數值上看，計算出的H-(C≡C)3-H的C-C鍵長為1.201/1.366/1.207 ?，與僅有的實驗觀察結果1.20/1.38/1.20 ?一致[40]。 完全基于π分子軌道計算的定域化軌道定位函數（LOL）[46]，即LOL-π，是一個用于揭示分子中π電子的離域情況的非常流行的實空間函數[6,11,13,47]。由于其獨特的sp雜化電子構型，碳鏈被認為有兩套類似于環碳的相互垂直的π共軛系統[6]。表S2列出了碳鏈H-(C≡C)n-H (n = 3-9)的兩組π分子軌道編號。圖2和圖S1分別顯示了在鏈上方0.5 ?的平面上的一組π系統（LOL-πx）和在包含鏈的平面上的另一組?系統（LOL-πy）的LOL-π的填色圖。盡管π電子的整體離域作用是顯著的，但由于長C-C鍵的性質更類似于單鍵，因此LOL-π在短C-C鍵附近的分布比在長C-C鍵附近的分布要充分得多，這很大程度上體現了長C-C鍵附近的電子離域作用的受阻。 圖2. 碳鏈H-(C≡C)n-H (n = 3-9)的LOL-πx在碳鏈上方0.5 ?平面上的填色圖，原子核位置用粉紅色字母標出，化學鍵用棕色直線標出。 Wiberg鍵級反映了兩個成鍵原子共享電子對的平均數目，通過它可以確定成鍵特性[48]。對于所研究的碳鏈H-(C≡C)n-H (n = 3-9)，Wiberg鍵級可以進一步分解為π分子軌道貢獻和σ分子軌道貢獻以研究它們各自的特性。圖3(a)和(b)分別繪制了計算得出的碳鏈π和σ分子軌道的Wiberg鍵級。可以看出，雖然長C-C鍵的π軌道的Wiberg鍵級明顯小于短C-C鍵，這表明長C-C鍵的π共軛作用較差，但長C-C鍵的π電子離域絕不能被輕易地忽略，因為相應的π-Wiberg鍵級達到了約0.48，是個不小的值。從鏈末端到鏈中心，π分子軌道的鍵級和σ分子軌道的鍵級都略有下降。尤其是最末端的C-C鍵（短鍵）的鍵級頗大，這顯然是參考文獻[26]中描述的末端效應造成的。碳鏈的鍵級分析結論與前面提到的通過BLA和LOL-π進行的鍵長和軌道離域性的分析是一致的。 圖3. 碳鏈H-(C≡C)n-H (n = 3-9, 15)的 (a) π分子軌道和 (b) σ分子軌道貢獻的Wiberg鍵級。 3.2 11Σu+ ← X1Σg+激發的電子光譜 除了紅外（IR）光譜，電子光譜是檢測和確定分子存在的另一種方法[49]。根據選擇規則，sp雜化碳鏈的允許躍遷被歸結為11Σu+ ← X1Σg+激發。本文采用兩種理論方法計算了碳鏈H-(C≡C)n-H (n = 3-9, 15)的50個最低激發態的激發能和振子強度。由TD-ωB97XD/def2-TZVP和非常穩健但昂貴得多的STEOM-DLPNO-CCSD/def2-TZVP計算獲得的相關數據以及11Σu+ ← X1Σg+躍遷的實驗值列于表1中。 表1. 計算的（括號外）和實驗的（括號內）吸收波長（λ）、碳鏈H-(C≡C)n-H (n = 3-9, 15) 11Σu+ ← X1Σg+躍遷所涉及的激發態（Sn）、主要躍遷軌道（MTOs）、躍遷能（Eg）和振子強度（f） a 本工作，TD-ωB97XD/def2-TZVP水平。 b 本工作，STEOM-DLPNO-CCSD/def2-TZVP水平。 c 參考文獻[19]。 d 參考文獻[20]。 e 參考文獻[21]。 從表1可以看出，對于所研究的碳鏈，兩種理論方法得到的11Σu+ ← X1Σg+波長的趨勢一致，兩者的系統偏差為14.3 ± 1.0 nm。在TD-ωB97XD/def2-TZVP水平上，11Σu+ ← X1Σg+波長的計算值與所有實驗觀測值驚人地一致。具體而言，相對于實驗值，模擬的波長總體小了5.6 ± 1.8 nm。這種一致性證實了我們計算方式的可靠性，并確保了所有后續特性分析的準確性。 雖然表1中的STEOM-DLPNO-CCSD結果也顯示出理想的精度，但其性能輕微劣于TD-ωB97XD。我們認為，ωB97XD在全局離域聚炔的TD-DFT計算中表現如此出色的關鍵原因之一在于它在長程極限（100%）時的Hartree-Fock成分非常高，這極大地避免了嚴重影響大范圍共軛體系描述精度的自相互作用誤差問題。本研究中采用的非常有效的TD-ωB97XD/def2-TZVP水平對于研究具有類似一維共軛特性的其它系統也很有適用性。 運用TD-DFT計算的結果模擬這些聚炔的吸收光譜如圖4(a)所示。結果表明，碳鏈H-(C≡C)n-H (n = 3-9)的11Σu+ ← X1Σg+吸收帶都是由基態S0到較高激發態Sn（n = 7、8或10）的躍遷產生的。隨著鏈長的增加，吸收帶明顯地向長波方向移動。碳鏈中11Σu+ ← X1Σg+躍遷的振子強度也隨著鏈長的增加而增大，這正如Vuitton和Scemama等人在他們的理論工作中所發現的[24,49]。振子強度之所以呈現出如此變化的趨勢，原因有以下兩點：(1)我們的計算結果清楚地表明，碳鏈的激發能隨著鏈長的增加而降低，而激發能是振子強度的分母；(2)由于參與躍遷的軌道在整個碳鏈上離域，如圖4(b)所示，碳鏈越長，躍遷偶極矩越大，振子強度與躍遷偶極矩的模的平方成正比。巧合的是，我們計算的H-(C≡C)n-H (n = 3-9)的振子強度值非常接近鏈中的n值，且整體略小于相應的半經驗結果[49]。值得一提的是，環[18]碳在相同計算水平下的振子強度為f = 3.02[6]，遠小于相近尺寸碳鏈的振子強度。 圖3. (a)碳鏈H-(C≡C)n-H (n = 3-9)的電子光譜（曲線）和振子強度（豎線）。采用半峰半寬（HWHM）為0.333 eV的高斯函數將理論計算的數據展寬為光譜曲線。(b)碳鏈H-(C≡C)n-H (n = 3-9)在各自關鍵激發態的主要分子軌道（等值面為0.006）。綠色和藍色區域分別表示正和負軌道相位。 如表1所示，所研究體系在紫外范圍內的強吸收幾乎都是由前線分子軌道HOMO-1/HOMO → LUMO/LUMO+1之間的電子激發產生的，它們都對應于典型的π-π*激發。由于碳鏈H-(C≡C)n-H (n = 3-9, 15)的對稱性，分子的HOMO-1/HOMO和LUMO/LUMO+1分別在能量上相互簡并。在圖4(b)中，占據的前線分子軌道（HOMO-1/HOMO）和未占據的前線分子軌道（LUMO/LUMO+1）的能級分別隨著鏈長的增加而增加和減小，導致分子的前線分子軌道能隙顯著減小。能隙的減小反映了較長碳鏈中電子相關性的增加，在其他π共軛低聚物中也觀察到了相同的行為[8,50-52]。可以預期其它分子軌道能級也會以與這些前線分子軌道相同的方式發生變化，因此隨著碳鏈的增長，吸收帶會發生明顯的紅移。值得注意的是，隨著鏈長的增加，HOMO-1/HOMO → LUMO/LUMO+1電子躍遷對最大吸收的貢獻逐漸減小，這表明在較長的碳鏈中有更多的其它軌道參與了11Σu+ ← X1Σg+躍遷，這是前線分子軌道的能級分布隨著碳鏈長度的增加而變得更密集的必然結果。順帶一提，在氖介質和溶液中通常可以觀察到相對于氣相中吸收波長的相當恒定的紅移，這與鏈的大小沒有明顯的關系[21-23]。關于非真空中電子躍遷的討論不屬于本工作的研究范圍。 非常流行的空穴-電子分析方法可以更具體地展示電子躍遷過程的特征[6,8,10,17,28,53]。關于空穴-電子分析方法的相關描述，參見參考文獻[6]。各種碳鏈的關鍵的電子激發的空穴-電子分布如圖5所示，由此圖展現的特征我們可以再次確定，所有碳鏈的最大吸收都是由π-π*躍遷引起的，因為涉及11Σu+ ← X1Σg+激發的相應空穴和電子都環繞分布在C-C鍵周圍，相關特性與環碳分子的電子躍遷特性相一致[6,8]。圖5還顯示，在11Σu+ ← X1Σg+激發過程中，空穴和電子的分布相對更加局域化，即集中在碳鏈的中心，而不是沿整個碳鏈均勻分布。這說明碳鏈中間的π電子共軛比兩端的強，這與前面討論的它們的幾何特征一致。此外，高度重疊的空穴和電子分布與環[18]碳[6]的相關特征極其相似，這解釋了碳鏈具有超強吸收，即巨大的振子強度的原因。 圖5. 碳鏈H-(C≡C)n-H (n = 3-9, 15)在11Σu+ ← X1Σg+激發（等值面為0.005 au）的空穴和電子分布的實空間表示。綠色和藍色區域分別表示空穴和電子分布。 3.3 預測碳炔的結構和激發特性 由于長鏈聚炔在一般條件下具有很高的反應活性，因此很難在實驗室中合成，更不用說難以捉摸的碳的同素異形體碳炔了。因此，它們的幾何和光學特征并不容易獲得。考慮到上文討論的碳鏈特性隨其長度發生規律性變化的特點，將相關性質準確外推到更長尺寸可以為未來可能會發現的更長的碳鏈甚至為預測sp雜化的碳炔的特性提供參考。對于一條特定的碳鏈，正如之前的DFT計算所揭示的那樣[25,43]，長C-C鍵和短C-C鍵的長度差從鏈的兩端向中心逐漸減小，圖S2中的BLA圖更清楚地表明了這一點。對于不同的碳鏈，隨著鏈長的增加，分子中相應位置的長C-C鍵和短C-C鍵分別有逐漸縮短和增長的趨勢。然而，C-H鍵并沒有隨著鏈的增長而發生明顯改變，這也與之前的理論研究報道一致[25,43]。雖然在目前研究的有限長度碳鏈中，上述兩點關于BLA的規律性似乎并不特別顯著，但可以想象，當鏈的長度無限增加時，這些規律的積累效應必然會非常明顯。可以預期無限長碳鏈上的所有長短C-C鍵會分別收斂到一個恒定值。根據Scemama等人給出的碳鏈H-(C≡C)n-H中C-C鍵長度的外推公式[49]，可以推斷出無限長碳鏈中長C-C鍵和短C-C鍵的長度分別約為1.329和1.229 ?。 以往的研究表明，11Σu+ ← X1Σg+躍遷波長與碳鏈中碳原子的數量之間存在線性或非線性的關系[21,24]。圖3(a)展示的不同碳鏈H-(C≡C)n-H (n = 3-9)的允許的電子躍遷的演變過程清楚地表明碳鏈的11Σu+ ← X1Σg+吸收波長并不隨碳鏈長度的增加而線性增加。圖6(a)顯示了11Σu+ ← X1Σg+躍遷的波長與C-C單元數量的關系。將ωB97XD/def2-TZVP水平得到的吸收波長用指數函數進行擬合可得到： λ(nm) = -305.807exp(-n/8.153) + 389.522               (1) 該方程與之前相關性質的理論研究[25,49]保持了相同的形式，R2 > 0.999，擬合質量相當完美。將n = 15代入擬合方程得到的數值結果為341.0 nm，這與DFT直接計算得到的結果（338.9 nm）也非常吻合。因此，此處擬合的方程是可靠的，并可用于預測任意長度碳鏈的11Σu+ ← X1Σg+吸收波長。根據公式(1)的推斷，無限長碳鏈的漸近值為389.5 nm。該值與之前半經驗計算的結論（384 nm）相當一致[49]。通過修正計算值和實驗值之間5.6 nm的系統誤差，我們可以推斷，碳炔-(C≡C∞)-的最大吸收將收斂到約395.1 nm。由于其絕大部分吸收帶可能位于可見光區域之外，可推斷碳鏈H-(C≡C)n-H應該像環[18]碳[6]一樣是無色的，至少對于那些不是很長的碳鏈來說是這樣。值得一提的是，我們之前對環[2N]碳（N = 3-15）光物理性質的研究表明，環碳分子中N的奇偶性不同，電子能譜的尺寸依賴性也遵循不同的規律[8]，而具有奇數和偶數C-C單元的碳鏈H-(C≡C)n-H (n = 3-9, 15)的變化特征則是一致的。 圖6. 碳鏈H-(C≡C)n-H中11Σu+ ← X1Σg+電子躍遷的最大波長(a)和振子強度(b)隨鏈中C-C單元的變化。 參照參考文獻[49]中的方程形式，我們提出了振子強度隨碳鏈H-(C≡C)n-H中C-C單元數變化的外推公式： f = -1.205exp(-n/1.823) + 0.931n + 0.480               (2) 如圖6(b)所示，擬合方程的相關系數為R2 > 0.999。通過n = 15得到的外推振子強度為14.445，與DFT計算出的結果（14.488）完全一致。 4. 結論 本文通過(TD-)DFT計算，研究了具有廣泛π電子離域的碳鏈H-(C≡C)n-H (n = 3-9, 15)的幾何和電子結構以及吸收光譜。在所有研究的碳鏈中都發現了鍵交替現象，但隨著碳鏈長度的增加，長C-C鍵和短C-C鍵之間的差異逐漸減小。所研究的每條碳鏈都有一個可觀測到的吸收帶，其振子強度非常大，與π-π*躍遷相對應。吸收波長和強度隨碳鏈長度的變化而有規律地變化。對電子結構的深入分析表明，碳鏈的強吸收歸因于電子激發過程中空穴和電子的高度重疊分布。對于較長的碳鏈，本文還得出了HOMO-LUMO gap、允許的11Σu+ ← X1Σg+吸收波長和振子強度的外推公式，準確性通過特定分子H-(C≡C)15-H進行了證實。這一工作比以往對碳鏈電子激發特性的研究更加系統和深入，有助于在實驗中針對特定性質的需要合理設計不同長度的碳鏈。 參考文獻 [1] J. Cernicharo, The polymerization of acetylene, hydrogen cyanide, and carbon chains in the neutral layers of carbon-rich proto-planetary nebulae, Astrophys. J., 608 (2004) L41-L44. [2] D.K. Bohme, PAH and fullerene ions and ion/molecule reactions in interstellar and circumstellar chemistry, Chem. Rev., 92 (1992) 1487-1508. [3] K. Kaiser, L.M. Scriven, F. Schulz, P. Gawel, L. Gross, H.L. Anderson, An sp-hybridized molecular carbon allotrope, cyclo[18]carbon, Science, 365 (2019) 1299-1301. [4] L.M. Scriven, K. Kaiser, F. Schulz, A.J. Sterling, S.L. Woltering, P. Gawel, et al., Synthesis of cyclo[18]carbon via debromination of C18Br6, J. Am. Chem. Soc., 142 (2020) 12921-12924. [5] D. Huang, S.L. Simon, G.B. McKenna, Chain length dependence of the thermodynamic properties of linear and cyclic alkanes and polymers, J. Chem. Phys., 122 (2005) 084907. [6] Z. Liu, T. Lu, Q. Chen, An sp-hybridized all-carboatomic ring, cyclo[18]carbon: Electronic structure, electronic spectrum, and optical nonlinearity, Carbon, 165 (2020) 461-467. [7] Z. Liu, T. Lu, Q. Chen, An sp-hybridized all-carboatomic ring, cyclo[18]carbon: Bonding character, electron delocalization, and aromaticity, Carbon, 165 (2020) 468-475. [8] Z. Liu, T. Lu, A. Yuan, X. Wang, Q. Chen, X. Yan, Remarkable size effect on photophysical and nonlinear optical properties of all-carboatomic rings, cyclo[18]carbon and its analogues, Chem. Asian. J., 16 (2021) 2267-2271. [9] Z. Liu, T. Lu, Q. Chen, Comment on "Theoretical investigation on bond and spectrum of cyclo[18]carbon (C18) with sp-hybridized", J. Mol. Model., 27 (2021) 42. [10] X. Wang, Z. Liu, X. Yan, T. Lu, H. Wang, W. Xiong, et al., Photophysical properties and optical nonlinearity of cyclo[18]carbon (C18) precursors, C18-(CO)n (n = 2, 4, and 6): Focusing on the effect of the carbonyl groups, Phys. Chem. Chem. Phys., 24 (2022) 7466-7473. [11] X. Wang, Z. Liu, X. Yan, T. Lu, W. Zheng, W. Xiong, Bonding character, electron delocalization, and aromaticity of cyclo[18]carbon (C18) precursors, C18-(CO)n (n = 6, 4, and 2): Focusing on the effect of carbonyl (-CO) groups, Chem. Eur. J., 28 (2022) e202103815. [12] T. Lu, Z. Liu, Q. Chen, Accurate theoretical evaluation of strain energy of all-carboatomic ring (cyclo[2n]carbon), boron nitride ring, and cyclic polyacetylene, Chin. Phys. B, 31 (2022) 126101. [13] Z. Liu, T. Lu, Q. Chen, Intermolecular interaction characteristics of the all-carboatomic ring, cyclo[18]carbon: Focusing on molecular adsorption and stacking, Carbon, 171 (2021) 514-523. [14] Z. Liu, T. Lu, Q. Chen, Vibrational spectra and molecular vibrational behaviors of all-carboatomic rings, cyclo[18]carbon and its analogues, Chem. Asian J., 16 (2021) 56-63. [15] T. Lu, Q. Chen, Ultrastrong regulation effect of the electric field on the all-Carboatomic ring cyclo[18]carbon, Chem. Phys. Chem., 22 (2021) 386-395. [16] T. Lu, Z. Liu, Q. Chen, Comment on “18 and 12 – Member carbon rings (cyclo[n]carbons) – A density functional study”, Mat. Sci. Eng. B, 273 (2021) 115425. [17] Z. Liu, X. Wang, T. Lu, A. Yuan, X. Yan, Potential optical molecular switch: Lithium@cyclo[18]carbon complex transforming between two stable configurations, Carbon, 187 (2022) 78-85. [18] X. Wang, Z. Liu, J. Wang, T. Lu, W. Xiong, X. Yan, M. Zhao, M. Orozco-Ic, Electronic structure and aromaticity of an unusual cyclo[18]carbon precursor, C18Br6, Chem. Eur. J., 29 (2023) e202300348. [19] E. Kloster-Jensen, H.-J. Haink, H. Christen, The electronic spectra of unsubstituted mono- to penta-acetylene in the gas phase and in solution in the range 1100 to 4000 ?, Helv. Chim. Acta, 57 (1974) 1731-1744. [20] C. Apetrei, R. Nagarajan, J.P. Maier, Gas phase 11Σu+ ← X1Σg+ electronic spectra of polyacetylenes HC2nH, n = 5-7, J. Phys. Chem. A, 113 (2009) 11099-11100. [21] T. Pino, H. Ding, F. Güthe, J.P. Maier, Electronic spectra of the chains HC2nH (n = 8–13) in the gas phase, J. Chem. Phys., 114 (2001) 2208-2212. [22] R. Eastmond, T.R. Johnson, D.R.M. Walton, A General Synthesis of the Parent Polyyens H(C=C)nH (n = 4-10,12), Tetrahedron, 28 (1972) 4601-4616. [23] M. Grutter, M. Wyss, J. Fulara, J.P. Maier, Electronic absorption spectra of the polyacetylene chains HC2nH, HC2nH-, and HC2n-1N- (n = 6-12) in neon matrixes, J. Phys. Chem. A, 102 (1998) 9785-9790. [24] V. Vuitton, A. Scemama, M.-C. Gazeau, P. Chaquin, Y. Benilan, IR and UV spectroscopic data for polyynes: predictions for long carbon chain compounds in titan's atmosphere, Adv. Space Res., 27 (2001) 283-288. [25] C. Zhang, Z. Cao, H. Wu, Q. Zhang, Linear and nonlinear feature of electronic excitation energy in carbon chains HC2n+1H and HC2nH, Int. J. Quantum Chem., 98 (2004) 299-308. [26] X. Fan, L. Liu, J. Lin, Z. Shen, J.-L. Kuo, Density functional theory study of finite carbon chains, ACS Nano, 3 (2009) 3788-3794. [27] C.D. Zeinalipour-Yazdi, D.P. Pullman, Quantitative structure-property relationships for longitudinal, transverse, and molecular static polarizabilities in polyynes, J. Phys. Chem. B, 112 (2008) 7377-7386. [28] Liu Z, Lu T. Optical properties of novel conjugated nanohoops: Revealing the effects of topology and size. J Phys Chem C,124 (2020) 7353-7360. [29] Liu Z, Lu T, Hua S, Yu Y. Aromaticity of Hu?ckel and Mo?bius topologies involved in conformation conversion of macrocyclic [32]octaphyrin(1.0.1.0.1.0.1.0): Refined evidence from multiple visual criteria. J Phys Chem C, 123 (2019) 18593-18599. [30] Liu Z, Lu T. Controllable photophysical and nonlinear properties in conformation isomerization of macrocyclic [32]octaphyrin(1.0.1.0.1.0.1.0) involving Hu?ckel and Mo?bius topologies. J Phys Chem C, 124 (2020) 845-853. [31] Liu Z, Yan X, Li L, Wu G. Theoretical investigation of the topology and metalation effects onthe first hyperpolarizability of rosarins. Chem Phys Lett., 641 (2015) 5-8. [32] Liu Z, Hua S, Wu G. Extended first hyperpolarizability of quasi-octupolar molecules by halogenated methylation: Whether the iodine atom is the best choice, J Phys Chem C, 122 (2018) 21548-21556. [33] J.D. Chai, M. Head-Gordon, Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections, Phys. Chem. Chem. Phys., 10 (2008) 6615-6620. [34] F. Weigend, R. Ahlrichs, Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy, Phys. Chem. Chem. Phys., 7 (2005) 3297-3305. [35] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, et al., Gaussian 16 Rev. A.03, (2016). [36] M. Nooijen, R. J. Bartlett, A new method for excited states: Similarity transformed equation-of-motion coupled-cluster theory, J. Chem. Phys., 106 (1997) 6441-6448. [37] F. Neese, Software update: The ORCA program system, version 4.0. WIREs: Comput. Mol. Sci., 8 (2018), e1327. [38] T. Lu, F. Chen, Multiwfn: A multifunctional wavefunction analyzer, J. Comput. Chem., 33 (2012) 580-592. [39] W. Humphrey, A. Dalke, K. Schulten, VMD: Visual molecular dynamics, J. Mol. Graph., 14 (1996) 33-38. [40] R. Hoffmann, Extended Hückel theory-v cumulenes, polyenes, polyacetylenes and Cn, Tetrahedron, 22 (1966) 521-538. [41] Q. Fan, G.V. Pfeiffer, Theoretical study of linear Cn (n=6-10) and HCnH (n=2-10) molecules, Chem. Phys. Lett., 162 (1989) 472-478. [42] T.D. Poulsen, K.V. Mikkelsen, J.G. Fripiat, D. Jacquemin, B. Champagne, MP2 correlation effects upon the electronic and vibrational properties of polyyne, J. Chem. Phys., 114 (2001) 5917-5922. [43] A. Scemama, P. Chaquin, M.-C. Gazeau, Y. Benilan, Theoretical study of the structure and properties of polyynes and monocyano- and dicyanopolyynes: Predictions for long chain compounds, J. Phys. Chem. A, 106 (2002) 3828-3837. [44] S. Szafert, J.A. Gladysz, Carbon in one dimension: Structural analysis of the higher conjugated polyynes, Chem. Rev., 106 (2006) PR1-PR33. [45] S, Suhai, Bond alternation in infinite polyene: Peierls distortion reduced by electron correlation. Chem. Phys. Lett., 96 (1983) 619-625. [46] H.L. Schmider, A.D. Becke, Chemical content of the kinetic energy density, J. Mol. Struct., 527 (2000) 51-61. [47] T. Lu, Q. Chen, A simple method of identifying π orbitals for non-planar systems and a protocol of studying p electronic structure, Theor. Chem. Accounts, 139 (2020) 25. [48] K.B. Wiberg, Application of the pople-santry-segal CNDO method to the cyclopropylcarbinyl and cyclobutyl cation and to bicyclobutane, Tetrahedron, 24 (1968) 1083-1096. [49] A. Scemama, P. Chaquin, M.-C. Gazeau, Y. Benilan, Semi-empirical calculation of electronic absorption wavelengths of polyynes, monocyano- and dicyanopolyynes. Predictions for long chain compounds and carbon allotrope carbyne, Chem. Phys. Lett., 361 (2002) 520-524. [50] S.S. Zade, N. Zamoshchik, M. Bendikov, From short conjugated oligomers to conjugated polymers. Lessons from studies on long conjugated oligomers, Acc. Chem. Res., 44 (2011) 14-24. [51] S. Hemmatiyan, D.A. Mazziotti, Unraveling the band gap trend in the narrowest graphene nanoribbons from the spin-adapted excited-spectra reduced density matrix method, J. Phys. Chem. C, 123 (2019) 14619-14624. [52] Z. Liu, X. Yan, L. Li, G. Wu, Modulation of the optical properties of D-π-A type azobenzene derivatives by changing the π-conjugated backbones: A theoretical study, J. Theor. Comput. Chem., 14 (2015) 1550041. [53] S. Zhang, Y. Wang, H. Xu, A new naphthalimide-picolinohydrazide derived fluorescent “turn-on” probe for hypersensitive detection of Al3+ ions and applications of real water analysis and bio-imaging, Spectrochim. Acta A Mol. Biomol. Spectrosc., 275 (2022) 121193. http://www.shanxitv.org/677 2023-07-29T18:32:00+08:00 下面是第12001~13000篇使用了筆者開發的功能非常全面而且在世界范圍流行的Multiwfn波函數分析程序（http://www.shanxitv.org/multiwfn）發表的文章的全文pdf合集（其中漏了幾篇文章對應的文件，不好找了）。 百度網盤下載地址：https://pan.baidu.com/s/1vbW8R5iSHkeEv2w_hANzaw?pwd=y8gf 這些文章的列表見文件包里的full list.txt，或者看http://www.shanxitv.org/multiwfn/all_citation3.html里面的第12001~13000號文章。 這些文章可以視為Multiwfn的例子庫，便于用戶了解Multiwfn在量子化學實際研究中能發揮的重要價值。如果想找某一類應用實例，可以利用Acrobat對指定文件夾里全部文檔進行搜索的功能搜索相關關鍵詞。 由于本文件包里的文章都是它們剛在期刊網站上刊登的時候我就下載的，所以文檔中大多沒有頁碼、卷號，且多數都是還沒經過proofing的原始狀態，大家只要到google學術搜索里搜索標題就能找到文章最終版鏈接。在Multiwfn主頁（http://www.shanxitv.org/multiwfn）的Cited by頁面里也有所有引用了Multiwfn的文章的匯總列表和文章鏈接。 之前的引用Multiwfn的第1~1010篇文章可以從這里下載：http://bbs.keinsci.com/thread-2850-1-1.html 之前的引用Multiwfn的第1011~2001篇文章可以從這里下載：http://bbs.keinsci.com/thread-6754-1-1.html 之前的引用Multiwfn的第2002~3000篇文章可以從這里下載：http://bbs.keinsci.com/thread-11023-1-1.html 之前的引用Multiwfn的第3001~4000篇文章可以從這里下載：http://bbs.keinsci.com/thread-14181-1-1.html 之前的引用Multiwfn的第4001~5000篇文章可以從這里下載：http://bbs.keinsci.com/thread-16840-1-1.html 之前的引用Multiwfn的第5001~6000篇文章可以從這里下載：http://bbs.keinsci.com/thread-19818-1-1.html 之前的引用Multiwfn的第6001~7000篇文章可以從這里下載：http://bbs.keinsci.com/thread-22597-1-1.html 之前的引用Multiwfn的第7001~8000篇文章可以從這里下載：http://bbs.keinsci.com/thread-25497-1-1.html 之前的引用Multiwfn的第8001~9000篇文章可以從這里下載：http://bbs.keinsci.com/thread-27978-1-1.html 之前的引用Multiwfn的第9001~10000篇文章可以從這里下載：http://bbs.keinsci.com/thread-31224-1-1.html 之前的引用Multiwfn的第10001~11000篇文章可以從這里下載：http://bbs.keinsci.com/thread-34124-1-1.html 之前的引用Multiwfn的第11001~12000篇文章可以從這里下載：http://bbs.keinsci.com/thread-36387-1-1.html 再次順帶強調關于恰當引用的問題。在Multiwfn程序主頁的首頁和Download頁面、程序手冊、程序啟動時屏幕上顯示的信息、程序包內的諸多文件（Multiwfn quick start.pdf文檔、How to cite Multiwfn.pdf文檔、LICENSE.txt文件）、Multiwfn入門貼（http://www.shanxitv.org/167）、Multiwfn FAQ（http://www.shanxitv.org/452）等盡可能多、用戶可以輕易看到的所有地方都非常非常非常明確強調如何對Multiwfn進行正確的引用，但使用Multiwfn的文章中不引用或者胡亂引用的現象依然十分嚴重（特別是中國的文章尤甚）！第12001~13000篇文章里沒按要求恰當引用Multiwfn的文章多達156篇。經常是只提及Multiwfn但不引用，甚至在引用Multiwfn的地方引用的卻是和Multiwfn根本沒有任何直接聯系的其他人的文章。還有不少文章里作者明明在研究中用了Multiwfn做分析，在論文里竟然連Multiwfn的名字都不提。不按明確聲明的方式恰當引用，對免費而且沒有任何經費支持的程序的發展十分不利，同時也是嚴重缺乏學術道德的行為。借此機會再次強烈呼吁用戶重視對Multiwfn的恰當引用。最恰當的引用方式見Multiwfn可執行文件包內的How to cite Multiwfn.pdf文檔的說明。沒有恰當引用Multiwfn的文章在本次的文件包里的full list.txt里都已經注上了Multiwfn was not properly cited!的字樣，在此予以強烈批評。 久久精品国产99久久香蕉