基于色散校正的密度泛函理論的食品添加劑太赫茲光譜研究

李天瑩，李 春，章 龍，馬卿效，蔣 玲

南京林業大學信息科學技術學院，江蘇 南京 210037

引 言

食品添加劑的違規使用或濫用會對身體健康產生嚴重的危害，傳統的高效液相色譜法(HPLC)等檢測方法存在耗時長、操作難等缺點，太赫茲光譜技術具有方便快捷、準確率高等優勢，被廣泛應用于食品添加劑檢測領域。目前，食品添加劑的太赫茲光譜研究多側重于實驗測量及結合化學計量學方法的定性定量分析[1-3]，缺乏深入的理論研究。

梁成森等[4]利用太赫茲時域光譜儀測試了甜味劑D-木糖及其加氫后得到的木糖醇，理論計算采用B3LYP泛函，與實驗結果吻合較好且能夠準確預測木糖醇在低頻處相對于D-木糖的紅移現象。Ge等[5]利用GAUSSIAN09軟件結合B3LYP/6-311++G(d, p)算法對苯甲酸等六種抗氧化劑的單分子進行計算，部分吸收峰與實驗結果吻合較好，但無法預測分子間作用力或晶格整體振動模式。以往的研究多基于B3LYP泛函，沒有根據研究體系的結構特點對比不同泛函的計算性能[6]。近年來，許多研究表明廣義梯度近似(GGA)泛函在太赫茲光譜計算中表現出良好的計算性能[7-8]，尤其是對于富含氫鍵的體系。King等[9]利用GGA泛函中的PBE泛函對腺嘌呤和胸腺嘧啶的共晶體進行理論模擬，由于共晶體中存在弱色散作用力，在PBE泛函的基礎上添加了色散校正系數，結果表明，色散修正后的PBE泛函對于得到收斂的晶體結構和高質量的振動頻譜有顯著影響。然而，對于不同計算方法添加色散修正后的計算精度、計算時長，缺乏深入的對比分析。

本文采用太赫茲時域光譜儀測試了三種富含氫鍵的食品添加劑，即苯甲酸、山梨酸和木糖醇，理論計算采用CRYSTAL14軟件包構建周期性晶體結構?？紤]到描述分子間相互作用的色散作用力對氫鍵系統至關重要，在傳統的B3LYP和PBE泛函基礎上添加了色散校正系數，建立了B3LYP-D和PBE-D模型，對比分析了色散校正前后四種模型的預測精度。從結構變化、計算精度、計算時長等多個角度評估了B3LYP-D和PBE-D模型的計算性能。

1 實驗部分

1.1 儀器

采用日本Advantest公司TAS7500SP型太赫茲時域光譜儀(THz-TDS)，儀器包含透射、反射、衰減全反射三個模塊，本文采用透射模塊。系統分辨率7.6 GHz, 動態范圍高于60 dB，理想狀態下可測得0.1～4 THz范圍內的高質量光譜。儀器自帶樣品腔，測量時通過充入氮氣減少空氣中水汽的影響，所有測試在室溫下進行。

1.2 樣品

所用樣品(純度≥99%)均購于Sigma-Aldrich公司，使用前未經進一步純化，三種食品添加劑的結構信息如表1所示。為了將吸收峰強度降低到合適的水平，將苯甲酸、山梨酸和木糖醇分別與聚乙烯以1∶10的質量均勻混合，混合后的粉末放置于直徑為13 mm的真空模具中壓制成厚度約為1 mm的薄片。同時將100 mg的聚乙烯粉末按照相同的方法壓制成相同規格的薄片，作為參考背景。

表1 三種食品添加劑的結構信息Table 1 The structures of three food additives

2 理論計算

2.1 模型的建立

采用北京并行科技股份有限公司提供的HPC系統(基于中國國家網格的超算云平臺)進行計算，利用ParaCloud(并行云桌面)平臺運行CRYSTAL14軟件[12]。

理論計算初始結構來源于劍橋晶體數據庫(CCDC)，在室溫(283～303 K)下獲得。CRYSTAL14軟件依據平移對稱性構建周期性晶胞構型，相比于單分子、多聚體及單晶胞構型，更能還原物質的實際形態。在傳統的B3LYP、PBE泛函基礎上添加了GRIMME倫敦色散校正系數，建立了B3LYP-D和PBE-D模型，基組函數均采用6-311G(d, p)。

2.2 模型評價方法

從結構優化、計算精度、計算時長等多角度分析不同計算模型的性能。其中結構優化部分采用均方根偏差(root-mean-square deviation, RMSD)來評估計算模型，計算精度部分采用平均絕對誤差(average absolute error, AAE)作為評價標準，見式(1)

(1)

其中，pcal_i和pexp_i分別代表幾何結構參數的計算值和對應的實驗值，RMSD的值越小，證明匹配效果越好，得到的幾何構型越穩定。

(2)

其中，pcal_j和pexp_j分別代表吸收峰位置(THz)的計算值和對應的實驗值，AAE越小，證明計算模型的模擬精度越高，計算結果越可靠。

3 結果與討論

3.1 收斂標準

CRYSTAL14計算過程包含兩個重要的收斂標準： 自洽場收斂閾值(TOLDEE)和雙電子積分截斷標準(TOLINTEG)，對計算時長影響較大。表2為不同收斂標準下基于B3LYP算法的苯甲酸模擬結果。當TOLDEE為默認值10-6，TOLINTEG為默認值(6, 6, 6, 6, 12)時，計算結果出現虛頻(-0.25 THz)。當TOLDEE=10-9，TOLINTEG=(8, 8, 8, 8, 16)時，吸收峰位的平均絕對誤差(AAE)為0.103，可以滿足基本的計算需求。進一步提高收斂標準至TOLDEE=10-12、TOLINTEG=(10, 10, 10, 10, 20)，模擬精度沒有明顯提高，且計算時間為原來的二倍以上。綜合考慮，后續計算收斂標準設置為TOLDEE=10-9、TOLINTEG=(8, 8, 8, 8, 16)。

表2 不同收斂標準下苯甲酸的模擬結果Table 2 Simulation results of benzoic acid under different convergence criteria

3.2 計算結果

3.2.1 結構優化

穩定的幾何結構是頻率計算的前提，優化后的結構與實驗數據越接近，證明得到的結構越穩定。為便于舉例說明，圖1給出了苯甲酸的晶胞結構、氫鍵二聚體及原子標號，其中虛線部分表示分子間氫鍵。苯甲酸晶體為單斜晶體，空間組為P21/c，單晶胞包含四個分子。CRYSTAL14程序通過識別空間組號自動生成對稱算子，結合分子坐標、晶格參數等信息還原晶體的實際結構。氫鍵是指氫原子介于兩個電負性較大的原子之間，形成的一種X—H…Y型分子間或分子內相互作用。在苯甲酸晶體中，存在大量的分子間氫鍵，即位于二聚體羧基之間的O—H…O型氫鍵，與色散作用力密切相關，考慮色散作用有利于得到收斂的結構以及高質量的THz光譜，因此在傳統泛函的基礎上添加了色散校正系數[11]。

表3給出了苯甲酸基于B3LYP-D和PBE-D泛函優化后的結構與實驗值的對比，包含鍵長、鍵角及晶胞體積變化。采用同樣的計算方法可以得到木糖醇優化后鍵長值RMSD為0.031 9(B3LYP-D)和0.005 5(PBE-D)，鍵角值RMSD為0.696(B3LYP-D)和0.749(PBE-D)，體積變化為-6.53%(B3LYP-D)和-8.80%(PBE-D)。結果表明，兩種經色散校正后的泛函均能得到較為穩定的結構，其中PBE-D泛函計算得到的結構體積更小。

需要注意的是，CRYSTAL14支持兩種結構優化方式： 完全優化方式更為精確，能夠同時優化晶格參數和原子坐標，但對于初始結構和基態相距較遠的體系，會造成長時間無法收斂且浪費大量機時的情況； 固定晶格參數法只優化原子坐標，精度略低但能夠保證結構快速收斂。根據三種食品添加劑各自的體系特點，苯甲酸和木糖醇采用完全優化方式，而山梨酸采取固定晶格參數法，(CCDC中山梨酸初始體積為1 191.175 ?3，優化后參考基態體積為595.588 ?3，體積預計減小50%)。山梨酸優化后鍵長值RMSD為0.0111(B3LYP-D)和0.0147(PBE-D)，鍵角值RMSD為0.860(B3LYP-D)和0.749(PBE-D)，與實驗結果較為吻合，能夠滿足計算需求。

圖1 苯甲酸的晶胞結構(a)及其氫鍵二聚體(b)Fig.1 Molecular packing in unit cell (a) and hydrogen bonding dimer of benzoic acid (b)

表3 苯甲酸優化后結構參數與實驗值對比Table 3 Optimized structural parameters of benzoic acid compared with experimental values

3.2.2 頻率計算

頻率計算可以得到研究對象在THz波段的吸收峰，從而進行深入的振動來源分析，本文采用MOLDRAW軟件分析吸收峰的振動來源。采用上述收斂標準TOLDEE=10-9、TOLINTEG=(8, 8, 8, 8, 16)，計算結果均無虛頻。

理論計算得到的吸收峰位與實驗結果的匹配程度能夠最直觀的反映計算模型的精度，表4對比了三種食品添加劑基于四種模型的計算結果及對應的振動模式。結果表明，色散校正模型能夠提高吸收峰預測精度，其中B3LYP-D泛函得到的結果更優。振動來源分析表明，三種食品添加劑的吸收峰振動來源以晶格整體的平移或旋轉為主，分子內振動模式貢獻較小。需要說明的是，振動模式往往是耦合的，表4列出的是吸收峰主要的振動模式，而省略了振動強度較弱的模式，例如山梨酸在1.79 THz處的主要振動模式為晶格沿a軸旋轉，但該頻點處伴隨著—CH3強度較弱的扭轉。

圖2給出了基于不同模擬方法的吸收峰位平均絕對誤差(AAE)，相比于色散校正前的模型，色散校正后的模型AAE降低了0.03～0.1。結果表明，提出的色散校正模型能夠提高模擬光譜預測精度，其中，B3LYP-D模型精度更高。然而，理論計算仍存在誤差，一方面源于振動頻率計算基于絕對零度，而實驗數據則在室溫下獲得； 另一方面，對于高估或低估振動強度的情況，可能是由于優化前后分子結構的細微改變，同時實驗測試中食品添加劑與聚乙烯混合比例的不同也會對吸收峰強度造成影響。

表4 吸收峰位及對應的振動模式Table 4 The position of peaks (THz) and vibration modes

圖2 基于不同模擬方法的吸收峰位平均絕對誤差

采用的服務器CPU型號為Intel Xeon E5-2692(2.20 GHz), 單節點可完成24核運算，所有計算均基于4節點96核計算，基于兩種色散校正后的泛函的CPU總耗時如表5所示。三種食品添加劑基于B3LYP-D泛函的耗時分別為PBE-D泛函的3.56倍、2.89倍、2.07倍，PBE-D泛函所用的計算機時明顯更少，對于較大的研究體系，可以節省大量的計算時間。因此，在體系較大、精度要求不高的情況下，可以優先考慮PBE-D泛函。

表5 兩種經色散校正后泛函的計算時長Table 5 The CPU time consumed by two dispersion- correction functionals (s)

4 結 論

利用THz-TDS系統測試了苯甲酸、山梨酸、木糖醇的吸收光譜，理論計算采用CRYSTAL14軟件構建了基于周期性邊界條件的晶體結構，討論了兩個重要的收斂標準。根據三種食品添加劑富含氫鍵的結構特點，在傳統的B3LYP和PBE泛函的基礎上添加了色散校正系數，構建了B3LYP-D和PBE-D模型。結果表明，色散校正模型能夠提高吸收峰指認精度。B3LYP-D模型的精度更高，但較為耗時，PBE-D模型的精度略低于B3LYP-D模型，但可以節省計算機時，需要根據具體的體系情況靈活運用。振動來源分析表明，三種食品添加劑的振動來源均以晶格整體的平移和旋轉為主，是以往的單分子或多聚體模型無法預測到的。本文提出的CRYSTAL14部分參數設置以及適用于氫鍵系統的固態色散校正模型，對于其他類似體系的理論研究具有重要的參考價值。