基于色散校正的密度泛函理论的食品添加剂太赫兹光谱研究
李天莹, 李春, 章龙, 马卿效, 蒋玲*
南京林业大学信息科学技术学院, 江苏 南京 210037
*通讯作者 e-mail: jiangling@njfu.edu.cn

作者简介: 李天莹, 1994年生, 南京林业大学信息科学技术学院硕士研究生 e-mail: m18362985762@163.com

摘要

食品添加剂违规使用或滥用会对身体健康产生严重危害, 相比于传统检测方法, 太赫兹光谱技术具有检测周期短、 准确率高等优势, 被广泛应用于食品添加剂检测领域, 但对其缺乏深入的理论分析。 主要致力于根据晶体结构特征改进传统计算方法, 提高食品添加剂模拟光谱预测精度并评估不同模型的计算性能。 选取了三种富含氢键的食品添加剂, 即苯甲酸、 山梨酸和木糖醇, 实验采用太赫兹时域光谱仪, 理论计算采用CRYSTAL14软件包计算周期性晶体结构。 鉴于富含氢键的体系中伦敦色散作用力不可忽略, 在传统的Becke-3-Lee-Yang-Par(B3LYP)泛函和Perdew-Burke-Ernzerhof(PBE)泛函的基础上添加了伦敦色散校正系数, 构建了B3LYP-D和PBE-D模型。 计算得到吸收峰位平均绝对误差(AAE)分别为苯甲酸0.073(B3LYP-D)、 0.096 (PBE-D), 山梨酸0.039(B3LYP-D)、 0.047(PBE-D), 木糖醇0.023(B3LYP-D)、 0.087(PBE-D), 相比于色散校正前的模型, AAE降低了0.03~0.1。 对比了B3LYP-D和PBE-D模型的计算时长, B3LYP-D泛函耗时均为PBE-D泛函的二倍以上。 结果表明, 经色散校正后的模型适用于富含氢键的系统, 能够提高模拟光谱预测精度, 其中B3LYP-D模型预测精度更高但较为耗时, PBE-D泛函预测精度略低于B3LYP-D模型, 但计算速度更快。 振动来源分析表明, 三种食品添加剂的振动模式均以晶格整体的平移和旋转为主, 分子内的振动贡献较小。 研究提出的色散校正模型对于其他类似体系的理论研究具有重要的参考价值。

关键词: 太赫兹光谱; 食品添加剂; 周期性晶体结构; 氢键系统; 色散校正
中图分类号:O433.5 文献标志码:A
Investigation on Terahertz Spectroscopy of Food Additives Based on Dispersion-Correction Functional Theory
LI Tian-ying, LI Chun, ZHANG Long, MA Qing-xiao, JIANG Ling*
College of Information Science and Technology, Nanjing Forestry University, Nanjing 210037, China
*Corresponding author
Abstract

Illegal use or abuse of food additives candamage human health. Terahertz (THz) spectroscopy was widely used for food additives detection with fast and accurate etc. advantages, but lacked in-depth theoretical analysis. This paperexplores to improve traditional methods to produce higher quality simulation spectra, and evaluate different models. THz spectra of sorbic acid, benzoic acid and xylitol were obtained by terahertz time-domain spectroscopy. Simulations were performed using CRYSTAL14 software to calculate the periodic crystals. Considering that dispersion force cannot be ignored in hydrogen-bonded systems, the dispersion correction term was used to augment the traditional functionals to construct the B3LYP-D and PBE-D models. The average absolute error (AAE) of benzoic acid absorption peak positions was 0.073 for B3LYP-D model, 0.096 for PBE-D model. For sorbic acid, AAE was 0.039 and 0.047, xylitol 0.023 and 0.087. The values of AAE were decreased by 0.03~0.1 compared with primitive functionals. The computation time of B3LYP-D model was more than twice of PBE-D model. Results showed that dispersion-corrected models can produce higher quality simulation spectra for hydrogen-bonded systems. The B3LYP-D model holds higher accuracy but is more time-consuming; the PBE-D model provides comparable accuracy to B3LYP-D model with much higher simulation speed. The spectral features were as signed as primarily lattice translations and rotations with lesser intramolecular torsions. The dispersion-corrected model proposed in this paper has important reference value for the theoretical research of other similar systems.

Keyword: Terahertz spectra; Food additive; Periodic crystal structure; Hydrogen-bonded system; Dispersion-correction
引言

食品添加剂的违规使用或滥用会对身体健康产生严重的危害, 传统的高效液相色谱法(HPLC)等检测方法存在耗时长、 操作难等缺点, 太赫兹光谱技术具有方便快捷、 准确率高等优势, 被广泛应用于食品添加剂检测领域。 目前, 食品添加剂的太赫兹光谱研究多侧重于实验测量及结合化学计量学方法的定性定量分析[1, 2, 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)

RMSD=i=1n(pcal_i-pexp_i)2n(1)

其中, pcal_ipexp_i分别代表几何结构参数的计算值和对应的实验值, RMSD的值越小, 证明匹配效果越好, 得到的几何构型越稳定。

AAE=j=1n|pcal_j-pexp_j|n(2)

其中, pcal_jpexp_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]

图1 苯甲酸的晶胞结构(a)及其氢键二聚体(b)Fig.1 Molecular packing in unit cell (a) and hydrogen bonding dimer of benzoic acid (b)

表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泛函计算得到的结构体积更小。

表3 苯甲酸优化后结构参数与实验值对比 Table 3 Optimized structural parameters of benzoic acid compared with experimental values

需要注意的是, 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), 与实验结果较为吻合, 能够满足计算需求。

3.2.2 频率计算

频率计算可以得到研究对象在THz波段的吸收峰, 从而进行深入的振动来源分析, 本文采用MOLDRAW软件分析吸收峰的振动来源。 采用上述收敛标准TOLDEE=10-9、 TOLINTEG=(8, 8, 8, 8, 16), 计算结果均无虚频。

理论计算得到的吸收峰位与实验结果的匹配程度能够最直观的反映计算模型的精度, 表4对比了三种食品添加剂基于四种模型的计算结果及对应的振动模式。 结果表明, 色散校正模型能够提高吸收峰预测精度, 其中B3LYP-D泛函得到的结果更优。 振动来源分析表明, 三种食品添加剂的吸收峰振动来源以晶格整体的平移或旋转为主, 分子内振动模式贡献较小。 需要说明的是, 振动模式往往是耦合的, 表4列出的是吸收峰主要的振动模式, 而省略了振动强度较弱的模式, 例如山梨酸在1.79 THz处的主要振动模式为晶格沿a轴旋转, 但该频点处伴随着— CH3强度较弱的扭转。

表4 吸收峰位及对应的振动模式 Table 4 The position of peaks (THz) and vibration modes

图2给出了基于不同模拟方法的吸收峰位平均绝对误差(AAE), 相比于色散校正前的模型, 色散校正后的模型AAE降低了0.03~0.1。 结果表明, 提出的色散校正模型能够提高模拟光谱预测精度, 其中, B3LYP-D模型精度更高。 然而, 理论计算仍存在误差, 一方面源于振动频率计算基于绝对零度, 而实验数据则在室温下获得; 另一方面, 对于高估或低估振动强度的情况, 可能是由于优化前后分子结构的细微改变, 同时实验测试中食品添加剂与聚乙烯混合比例的不同也会对吸收峰强度造成影响。

图2 基于不同模拟方法的吸收峰位平均绝对误差Fig.2 Average absolute error of absorption peaks based on different simulation methods

3.2.3 计算时长

采用的服务器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部分参数设置以及适用于氢键系统的固态色散校正模型, 对于其他类似体系的理论研究具有重要的参考价值。

参考文献
[1] Haddad J E, Miollis F D, Sleiman J B, et al. Analytical Chemistry, 2014, 86: 4927. [本文引用:1]
[2] Sun X D, Zhu K, Liu J B, et al. Journal of Infrared, Millimeter, and Terahertz Waves, 2019, 40: 466. [本文引用:1]
[3] Du C M, Zhang X, Zhang Z Y. Vibrational Spectroscopy, 2019, 100: 64. [本文引用:1]
[4] LIANG Cheng-sen, ZHAO Guo-zhong(梁成森, 赵国忠). Spectroscopy and Spectral Analysis(光谱学与光谱分析), 2011, 31(2): 323. [本文引用:1]
[5] Ge H Y, Jiang Y Y, Lian F Y, et al. Optik, 2016, 127: 4954. [本文引用:1]
[6] Su T F, Huang R, Su Y Q, et al. Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy, 2015, 139: 456. [本文引用:1]
[7] Yan H, Fan W H, Zheng Z P. Optics Communications, 2012, 285: 1593. [本文引用:1]
[8] Zhang H, Zhang Z H, Zhao X Y, et al. Chinese Physics B, 2015, 24(7): 073301. [本文引用:1]
[9] King M D, Ouellette W, Korter T M. J. Phys. Chem. A, 2011, 115(34): 9467. [本文引用:1]
[10] Dovesi R, Orland o R, Erba A, et al. Int. J. Quantum Chem. , 2014, 114(19): 1287. [本文引用:1]
[11] Grimme S, Antony J, Ehrlich S, et al. J. Chem. Phys. , 2010, 132: 154104. [本文引用:1]