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

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

南京林業(yè)大學(xué)信息科學(xué)技術(shù)學(xué)院，江蘇 南京 210037

引 言

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

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

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

1 實驗部分

1.1 儀器

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

1.2 樣品

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

表1 三種食品添加劑的結(jié)構(gòu)信息Table 1 The structures of three food additives

2 理論計算

2.1 模型的建立

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

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

2.2 模型評價方法

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

(1)

其中，pcal_i和pexp_i分別代表幾何結(jié)構(gòu)參數(shù)的計算值和對應(yīng)的實驗值，RMSD的值越小，證明匹配效果越好，得到的幾何構(gòu)型越穩(wěn)定。

(2)

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

3 結(jié)果與討論

3.1 收斂標準

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

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

3.2 計算結(jié)果

3.2.1 結(jié)構(gòu)優(yōu)化

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

表3給出了苯甲酸基于B3LYP-D和PBE-D泛函優(yōu)化后的結(jié)構(gòu)與實驗值的對比，包含鍵長、鍵角及晶胞體積變化。采用同樣的計算方法可以得到木糖醇優(yōu)化后鍵長值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)。結(jié)果表明，兩種經(jīng)色散校正后的泛函均能得到較為穩(wěn)定的結(jié)構(gòu)，其中PBE-D泛函計算得到的結(jié)構(gòu)體積更小。

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

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

表3 苯甲酸優(yōu)化后結(jié)構(gòu)參數(shù)與實驗值對比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)，計算結(jié)果均無虛頻。

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

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

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

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

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

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

4 結(jié) 論

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