通過量子化學計算和Multiwfn程序預測化學物質的顏色
通過量子化學計算和Multiwfn程序預測化學物質的顏色
文/Sobereva@北京科音 2023-Mar-26
1 前言
顏色無疑是化學物質極其重要的特征，理論預測化學物質的顏色對于從頭設計和改良化學物質有重要意義。通過量子化學計算可以模擬分子的UV-Vis光譜、通過第一性原理計算可以模擬固體和液體的UV-Vis譜，在原理上可以進而預測物質對應的顏色。一個很常見、很粗糙的預測顏色的方式是：在可見光波長范圍內（一般取380-760 nm）考察吸收峰位置，根據其波長和對應的顏色的關系判斷物質吸收的光的顏色，然后根據補色圖得到它的補色，這就是物質實際顯現的顏色，即反射或者透射光對應的顏色。每次我在北京科音初級量子化學培訓班（http://www.keinsci.com/workshop/KEQC_content.html）里講電子激發計算的時候也都舉例介紹這種做法。這種做法雖然簡單，但是準確度很有限，而且僅對于在一個波長處有簡單的吸收峰的時候才適合。如果可見光區有不止一個吸收峰，或者需要考慮譜帶形狀對吸收的影響，這種做法無疑就不適用了。所以這種做法適合作為教學目的，而用于實際研究目的、試圖準確預測化合物的顏色的時候，明顯應當用更嚴格的做法。
本文將介紹Multiwfn支持的將物質的UV-Vis吸收光譜轉化為物質顯現的顏色的方法，具有極高的實用性，解決了很多人研究的實際困難。此方法原理嚴格，在Multiwfn中使用超級簡單方便，只需要敲幾下鍵盤，顏色就能立刻呈現在屏幕上！用戶可以直接將諸如Gaussian、ORCA、CP2K等程序做TDDFT等電子激發計算得到的輸出文件作為輸入文件，從而“從頭”預測物質的顏色，也可以將實驗或者其它方式得到的UV-Vis光譜曲線數據作為輸入文件來預測顏色。
本文在第2節將介紹光譜曲線轉化為對應的顏色的原理。第3節將以預測靛青染料水溶液的顏色演示如何在Multiwfn中操作，輸入文件是Gaussian程序在TDDFT下對其做電子激發計算的輸出文件。第4節將示例如何用通過Multiwfn將網絡上搜索到的誘惑紅的實驗UV-Vis譜轉化為它對應的顏色，
讀者必須使用在2023-Mar-26及以后更新的Multiwfn版本，否則沒有本文說的功能。Multiwfn可以在官網http://www.shanxitv.org/multiwfn免費下載，若用于發表文章必須按程序網站、手冊里的明確要求以正確方式引用Multiwfn。不了解Multiwfn者可參看《Multiwfn FAQ》（http://www.shanxitv.org/452）。
本文的例子涉及到的文件都可以在這里下載：http://www.shanxitv.org/attach/662/file.rar。
2 由吸收光譜轉化為顏色的方法
本節內容較長，如果只希望無腦地得到結果，直接看后面的例子就完了。
一種顏色對應于色彩空間里的一個點。色彩空間（color space）有不同的定義，常見的有CIE1931 XYZ、CIE xyY、CIELAB、CIE RGB、sRGB（standard RGB）、Adobe RGB、CMYK等。在這些色彩空間中，每種顏色通過特定的坐標描述，例如sRGB和Adobe RGB都用R（紅）、G（綠）、B（藍）作為分量表示，是范圍為(0-255,0-255,0-255)的整數，但每個分量也經常用0到1范圍的分數表示。CIE1931 XYZ通過X、Y、Z作為坐標表示顏色，CIE xyY通過x、y、Y表示顏色，CIELAB通過a*、b*、L*表示顏色。CIE各種色彩空間的分量的定義見https://en.wikipedia.org/wiki/CIE_1931_color_space和https://en.wikipedia.org/wiki/CIELAB_color_space。
下圖是CIE色度圖，里面的完整區域是人眼可識別的各種顏色。弧線輪廓上不同位置的顏色對應的是不同波長的單色光，而圖內側區域的顏色是不同波長的單色光混合所產生的顏色。圖中橫、縱坐標分別對應的是CIE xyY的x和y分量，二者一起體現的是色度（chromaticity），而Y分量對應的是亮度（luminance）。如果亮度為0，那么不管x和y是什么值顏色都為黑。亮度足夠大時，圖中(1/3,1/3)位置對應白色，而如果亮度不夠，則這個點對應的是灰色。
并不是所有人眼可以感受的顏色都是可以在屏幕上顯示出來的，也沒有任何屏幕能顯示出來上圖里的完整色域里的所有顏色。上圖中三角形包圍的區域是sRGB色域，也是絕大多數消費級顯示器能顯示出的顏色范圍。Adobe RGB色域也是上圖中的一個三角形的子集，但比sRGB包圍的面積更大。現在有越來越多的中檔顯示器能顯示出Adobe RGB色域的大部分顏色，因而比只支持sRGB色域的顯示器能顯示出更豐富的顏色。考慮到大多數人用的顯示器還是sRGB色域的，而且絕大多數圖片里的顏色也是基于sRGB色域記錄的，因此預測化學物質的顏色并在屏幕上呈現出來，重點就是怎么得到物質吸收的光所對應的sRGB色域里的RGB值，之后再取其補色是很簡單的事。
肯定有人會問，為什么我的顯示器是sRGB色域的，但上圖里在三角形以外還是顯示了顏色？這是因為此圖是把CIE色度圖轉化到了sRGB空間下，在三角形以內的顏色是你的屏幕上能如實顯示的，因而其中不同位置的顏色在你的屏幕上能夠區分。而在三角形以外的區域（sRGB色域無法表現的區域），上圖在作圖時做了特殊處理。具體來說，這些區域的x,y坐標（結合特定的Y）直接轉化成sRGB色域的(R,G,B)時，會出現數值為負或者超過255的分量值，將負值取0處理、超過255的當255處理之后就得到了上圖。因此上圖的三角形以外區域的顏色在你的屏幕上呈現的并不是其實際的顏色，而只是對應了sRGB色域三角形區域里邊界位置的顏色。
下面來說一下有了UV-Vis吸收光譜曲線后，怎么得到吸收的光對應的RGB值（下文一律說的是sRGB色彩空間的情況）。更多細節參看：
https://www.oceanopticsbook.info/view/photometry-and-visibility/chromaticity
https://www.oceanopticsbook.info/view/photometry-and-visibility/from-xyz-to-rgb
https://support.hunterlab.com/hc/en-us/article_attachments/201533555/an-1002b.pdf
下面的公式用于把依賴于波長λ的光譜曲線Λ(λ)轉化為CIE1931 XYZ的X、Y、Z值（三刺激值，tristimulus values），可見算每個值就是求一維積分而已。式中Km是發光效率，是個常數。式中帶上橫線的x、y、z叫做三刺激值函數（tristimulus function），或稱配色函數（color-matching function）。
CIE1931 2o三刺激值函數如下所示，列表數據可以在http://cvrl.ucl.ac.uk/cmfs.htm下載。顯然上式就是求光譜在這仨函數上的投影來得到X、Y、Z分量。
通過紅、綠、藍三種單色光并不能組合出所有人眼可以識別的顏色，但是通過以上三種三刺激值函數所分別對應的非單色光的混合則可以產生，因此任何人眼可識別的顏色都可以用X、Y、Z坐標來表示。
有了X、Y、Z后就可以按下式計算色度坐標x、y、z。這相當于對坐標做了歸一化，從而去除了對光譜強度的依賴性。由于x+y+z=1，因此只有x和y是獨立的，這兩個值就是前述的CIE1931 xyY色彩空間的色度值，用于描述獨立于亮度的顏色。CIE1931 xyY里的Y直接對應上面說的三刺激值里的Y，靠它的大小來體現亮度。
現在有了x和y，對照CIE色度圖中的相應位置，就知道光譜曲線對應的什么顏色了。當然，若要與人眼感受到的顏色相對應，還需要再額外指定亮度值Y。
光譜曲線對應的x、y色度值經常超過sRGB色域而無法在sRGB色域的顯示器上顯示出來，因此實際中得把X、Y、Z轉化為sRGB的R、G、B坐標，這只需要按下式做一個線性變換即可。各種色域顏色之間的變換關系見http://brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html。注意此處的X、Y、Z坐標被scale到了0至1之間。
注：人眼對亮度的響應不是線性的，對深色的差異比對淺色的差異更敏感。因此很多程序將XYZ轉換為RGB時還自動做Gamma校正使得RGB值被拉伸，使得0-255數值范圍內能更充分地展現出人眼能分辨的顏色。具體做法是R<=0.0031308時將R替換為12.92*R，R大于0.0031308時將R替換為1.055*R^(1/2.4)-0.055，其中2.4是gamma值。對G、B分量也都這么處理。對于本文的目的這個轉換沒必要，反倒有礙取補色，Multiwfn也不考慮這個校正。
這么轉化之后的R、G、B如果都在0到1之間，那么再乘上255并取最近的整數，就直接可以輸入到比如Windows的畫圖、Photoshop等程序的調色板里，在顯示器上顯示光譜對應的顏色了。但如上轉換后大概率會發現RGB的一個或多個分量小于0或大于1，也即超過了sRGB的色域范圍，這時候只能做近似處理得到sRGB色域能表達的顏色。比如乘了255并取整數后，發現數值是(-358,59,266)，就得把負值當成0，超過255的當成255，即取為(0,59,255)。或者在把B調回255的同時讓G也按照特定規則減小，比如等比例減小，沒有唯一的做法。
如上將物質的UV-Vis吸收光譜轉化為RGB坐標后，再取補色（complementary color）就是此物質展現出的顏色了。下面的圖叫補色圖，互相處在對面位置的兩種顏色彼此構成補色。一個顏色的補色是與之對比最強烈的顏色，也可以叫相反色。圖中標注的度數稱為色調（hue），紅、綠、藍分別是圖中0、120、240度的三個對稱位置。典型的顏色對應的RGB值在圖中也標注了。例如圖中的海藍色是210度的位置，在其對面位置（30度）的地方就是其補色，即橙色。
從算法上計算某個(R,G,B)值對應的顏色的補色很簡單，補色顏色為(255-R,255-G,255-B)。例如紅色是(255,0,0)，其補色就是(0,255,255)，即黃色。另外還要注意亮度問題。例如(255,255,155)是淺黃色，直接算出來的補色是(0,0,100)，對應的是深藍色，二者不僅在色度上有差別，在亮度上也有極大的差別。如果要得到亮度相同而僅在色調上有差異的補色，需要平移補色的RGB值使其最大分量與原顏色相同，對此例也即給(0,0,100)各分量都增加155成為(155,155,255)使得最大分量和之前一樣同為255，此顏色對應的是淡藍色。上面的補色圖中的所有顏色的亮度都是相同的，因為最大分量都是255。注意亮度值有不同計算方式，這里取最大分量值作為亮度的衡量標準是Hue Saturation and Brightness (HSB)色彩模型所使用的，其亮度（B）參數定義為max(R,G,B)/255*100%。
一種化學物質對光的吸收整體越強，由于反射和折射光越弱，無疑看到的顏色就越暗，反之亦然。因此考察其展現的顏色看似應當按(255-R,255-G,255-B)直接取補色。但實際中并無法得到物質透射和反射光的絕對的亮度，故考察物質顯現的顏色時在取補色時只需要體現色調的差異就夠了，可以把RGB最大分量都統一平移到255。
至此，從物質的吸收光譜轉化到其展現的顏色涉及的知識講述完畢，下面就是具體例子了。
3 實例：通過理論計算從頭預測靛青染料的顏色
下圖給出了靛青（indigo）的結構、水溶液的實際顏色、水溶液的實驗光譜（下圖中虛線）。下面就要基于量子化學計算和上一節介紹的原理在Multiwfn中從頭預測它的顏色，看看和實際是否相符。
此例的量子化學計算使用非常常用的Gaussian 16，通過最主流的TDDFT方法做電子激發計算，其輸出文件用于給Multiwfn繪制光譜并預測顏色。如果讀者對量子化學計算和光譜的模擬缺乏常識知識，請先參看《Gaussian中用TDDFT計算激發態和吸收、熒光、磷光光譜的方法》（http://www.shanxitv.org/314）和《使用Multiwfn繪制紅外、拉曼、UV-Vis、ECD、VCD和ROA光譜圖》（http://www.shanxitv.org/224）。如文中所述，用Gaussian不是必須的，還有很多其它選擇，比如ORCA、CP2K等。用戶還可以把任意程序計算的激發能和振子強度寫成標準的文本文件作為Multiwfn的輸入。凡是能用于Multiwfn繪制UV-Vis光譜的輸入文件一律都可以用來預測顏色。
本文的文件包里indigo_S0opt.gjf是靛青分子在IEFPCM模型表現的水環境下在常用的B3LYP/6-311G*級別下做基態的優化和振動分析的輸入文件，從相應的out文件可見沒有虛頻，體系是完全平面的。基于優化完的結構創建TDDFT的輸入文件indigo_TDDFT.gjf，使用B3LYP/def2-TZVP級別做TDDFT算最低20個激發態，依然用IEFPCM表現水環境。其輸出文件indigo_TDDFT.out在本文的文件包里提供了。
啟動Multiwfn，載入indigo_TDDFT.out，然后輸入
11 //繪制光譜
3 //UV-Vis
25 //基于光譜計算顏色
現在在屏幕上會看到在CIE1931 2o三刺激值函數有定義的360到830 nm區間內的UV-Vis吸收光譜，如下所示，在可見光范圍內吸收峰約在590 nm
在圖上點右鍵關閉光譜圖，然后馬上就蹦出來下圖。color是光譜對應的顏色，complementary color是其補色。下面的Maximum brightness of above colors是上面兩種顏色的RGB最大值都平移到255的情況，使得光譜顏色和補色只有色度的差異而沒有亮度的差異。對于基于吸收光譜預測物質顏色，一般就看上圖的右下角的顏色即可。上圖說明，水中的靛青吸收的是黃色光，顯現出的顏色是藍色，這和實際中靛青的顏色完全一致，靠理論計算預測靛青的顏色大成功！
此時在Multiwfn的文本窗口還顯示了顏色相關的具體數值，如下所示。可見由XYZ直接轉化成的RGB值超過了sRGB色彩空間，因此Multiwfn自動做了處理。
CIE1931 XYZ: 1077992.509 1123653.902 188673.482
Fractional CIE1931 XYZ: 0.9593634718 1.0000000000 0.1679106720
CIE1931 xy: 0.4509825283 0.4700851571
Note the R,G,B values show below correspond to standard RGB (sRGB) color space
RGB (0-1): 1.487926 0.953110 0.026876
RGB (0-255): 379 243 7
Note: The color exceeds sRGB color space! Now the R,G,B values are scaled into
valid range:
RGB (0-1): 1.000000 0.953110 0.026876
RGB (0-255): 255 243 7
RGB of complementary color: 0 12 248
RGB of original color (maximum brightness): 255 243 7
RGB of complementary color (maximum brightness): 7 19 255
在顯示顏色的窗口上點右鍵關閉后，Multiwfn會問你是否把這些顏色保存成圖像文件，如果選y，則當前目錄下就會出現以DISLIN為開頭的圖像文件，內容和窗口里顯示的一致。
如果在此之后再選0讓Multiwfn繪制光譜，看到的光譜將和預測顏色功能順帶給出的光譜相一致，因為使用選項25時Multiwfn自動會將橫坐標設為360到830 nm（與此同時繪制的點數設為了471，因此點的間隔為1 nm）。
4 使用Multiwfn對任意方式得到的UV-Vis光譜計算對應的顏色
為了讓Multiwfn基于光譜曲線預測顏色的功能有最大程度的靈活性和普適性，我令Multiwfn也能夠基于任意方式獲得的光譜曲線預測顏色，用戶只需要提供記錄光譜曲線的文本文件作為輸入文件即可，第一列為是以nm為單位的橫坐標，第二列是吸收曲線值。作為例子，此例用Multiwfn基于誘惑紅（Allura red）的實驗UV-Vis曲線預測它的顏色。
首先要獲得誘惑紅的UV-Vis光譜數據文件。用Google的圖片搜索功能搜Allura red UV-Vis，馬上就能找到好多此體系的UV-Vis譜，隨便點開一個，比如下圖
然后用WebPlotDigitizer（https://apps.automeris.io/wpd/）在線工具獲得此光譜曲線的XY數據點。做法是：
(1)先把搜到的光譜圖保存下來，然后拖到WebPlotDigitizer的窗口里，依次點Align Axes、Proceed按鈕
(2)依次點擊坐標軸這些位置：(190,0)、(890,0)、(190,0)、(190,4)，然后點擊網頁右邊的Complete按鈕，把X-axis的Point 1和Point 2分別設為190和890，Y-axis的Point 1和Point 2分別設為0和4，然后點OK
(3)在網頁右邊的顏色方塊上點左鍵，選擇黑色（即光譜圖中曲線的顏色），然后點Done。網頁右邊點Box，然后在圖片上拖拽一個方框以把圖片中曲線部分完全框住
(4)在網頁右下側的Algorithm里選擇X Step w/ Interpolation，ΔX Step設2，Smoothing設100，然后點Run。此時會看到下圖，可見光譜曲線被充分采樣了
點擊網頁左邊的View Data，在新窗口里點擊Copy to Clipboard，然后把數據粘貼到文本文件里，保存為光譜曲線數據文件spectrum.txt。此文件已經提供在了本文的文件包里。當前這個曲線數據對應190-890 nm，數據點間隔為2 nm。Multiwfn對用戶提供的曲線數據的數據范圍和數據點間隔沒有強制要求，在預測顏色的過程中會自動插值得到所需要的位置的光譜值。
啟動Multiwfn，載入spectrum.txt，然后依次輸入11、0，馬上預測出的光譜就顯示出來了，如下所示。與此同時具體顏色參數值也在文本窗口里顯示了。
上圖顯示出誘惑紅的吸收光譜曲線對應的補色是紅色，和誘惑紅的實際顏色完全一致！
5 總結&其它
本文介紹了基于物質的UV-Vis吸收光譜預測其顏色的原理，并介紹如何通過Multiwfn極其簡單方便地實現這個目的。由本文的例子可見整個過程容易極了，還不需要任何額外耗時。對于經常做電子激發計算的理論化學研究者，如果發現自己研究的物質在可見光區域內有吸收，都不妨用此文的方法預測顏色并且在文中提及，可以令文章的信息更加充實。本文的方法也可以指導研究者理論設計染料、色素類物質。
本文介紹的預測顏色的方法也非常有教學意義。在講量子化學計算時，如果能讓學生體會到靠簡單的步驟就能得到物質的準確的顏色，絕對會令他們很有快感和成就感。
本文雖然以分子體系為例，但同樣的過程也完全可以用于預測晶體等周期性體系的顏色，比如可以用CP2K做周期性體系的TDDFT計算的輸出文件作為Multiwfn的輸入文件然后按前例操作。CP2K做TDDFT的計算介紹見《使用CP2K結合Multiwfn對周期性體系模擬UV-Vis光譜和考察電子激發態》（http://www.shanxitv.org/634），筆者在北京科音CP2K第一性原理培訓班（http://www.keinsci.com/workshop/KFP_content.html）里也都有詳細講授。
從本文的例子可見用Multiwfn預測顏色的原理嚴格，但得到真實的物質顏色的前提之一是做電子激發用的體系的幾何和計算模型和實際有可比性。如果體系同時存在多種比例不小的構象，而且不同構象的UV-Vis光譜在可見光區域差異不小，就必須考慮構象權重平均才可能得到靠譜的顏色，參看《使用Multiwfn繪制構象權重平均的光譜》（http://www.shanxitv.org/383）。和Multiwfn常規方式繪制構象權重平均光譜一樣，用multiple.txt作為輸入文件，然后依次輸入11、3、25讓Multiwfn預測顏色即可。
考慮振動耦合有助于得到更準確的光譜曲線，某些情況下對于預測的顏色有不可忽視的改進。但由于考慮振動耦合又麻煩又耗時，故本文第3節的例子沒考慮此問題，感興趣的讀者可以看《振動分辨的電子光譜的計算》（http://www.shanxitv.org/223）。Gaussian等程序做這種計算給出的考慮振動耦合的光譜曲線數據改寫成兩列的txt數據后就可以直接作為Multiwfn的輸入文件按照第4節的例子預測顏色了。
本文的例子說的都是基于吸收光譜預測物質的顏色，取的是Multiwfn窗口中給出的complementary color下面的顏色。顯然也可以用Multiwfn預測物質發射的光的顏色。比如要預測熒光的顏色，可以按照《Gaussian中用TDDFT計算激發態和吸收、熒光、磷光光譜的方法》（http://www.shanxitv.org/314）里面說的繪制熒光的方法操作，只不過在繪制光譜界面里最后不是選0來繪制光譜圖，而是選25得到光譜曲線對應的顏色，屆時直接取Multiwfn給出的圖中的color下面的顏色即可（而此時一起給出的補色就不用管了）。
最后提醒一下，使用Multiwfn做任何事情，包括按本文預測顏色，請不要忘了恰當引用Multiwfn的原文。