实验用样品苯甲酸、 山梨酸均采购自阿拉丁试剂网, 聚乙烯(Polyethylene, PE)采购于西格玛奥德里奇(Sigma-Aldridge)公司。 首先将所用样品(苯甲酸与山梨酸)与聚乙烯按实验设定的质量比称量好并倒入玛瑙研钵中充分研磨, 然后用压片机以5 MPa的压力压制5 min, 制为厚度1.35~1.40 mm左右, 直径为13 mm的圆形薄片。 分别制作浓度梯度为10%的苯甲酸与山梨酸单质样片各3个, 以及苯甲酸和山梨酸的混合样品3个, 样片制备信息表如表1。

表1

Table 1

表1(Table 1)

表1 样品制备信息表 Table 1 Sample preparation information sheet Sample Thickness /mm Proportion /% 苯甲酸 (Benzoic acid) 1.40 10 1.38 10 1.40 10 山梨酸 (Sorbic acid) 1.38 10 1.38 10 1.42 10 混合样片 (Benzoic acid+Sorbic acid) 1.40 10、10 1.40 10、10 1.40 10、10

表1 样品制备信息表 Table 1 Sample preparation information sheet

实验时, 为降低样片颗粒度及均匀度所带来的误差, 对每个样片的同一点位均进行三次测量, 并取其平均值保存为最终结果, 3组9个样品共获得27组光谱数据。

经过太赫兹时域光谱系统获得样片的时域信号后, 根据Dorney等和Duvillaret等提出的光学参数提取模型, 结合快速傅里叶变换(fast Fourier transform, FFT)计算样片的折射率与吸收系数, 如式(1)和式(2)所示^{[8, 9]}。

ns(ω)=1-cωdφ(ω)(1)

式(1)中, c为光速, d为样片厚度, ω 为时域信号经傅里叶变换(FFT)后的角频率, φ (ω )为样品信号与参考信号之间的相位差。

A(ω)=-2dln|Esam(ω)|[ns(ω)+1]2|Eref(ω)|4ns(ω)(2)

式(2)中, A(ω )为吸收系数, E_sam(ω )为样品信号的幅频特性, E_ref(ω )为参考信号的幅频特性, d为样片厚度, n_s(ω )为试样的折射系数。

图1为苯甲酸和山梨酸在0.5~2.2 THz波段的特征吸收谱。 首先, 如图可知苯甲酸共有4个特征吸收峰, 分别位于0.619、 0.916、 1.531和2.069 THz, 而山梨酸则仅在1.603 THz有1个较明显的特征吸收峰, 故可根据特征吸收峰的峰位及峰数实现上述两种添加剂的定性鉴别。 其次, 图1中混合物(苯甲酸和山梨酸)样片的吸收光谱中属于苯甲酸的1.531 THz处的单质吸收峰在吸收谱中缺失了, 属于特殊吸收现象, 这将影响混合样片通过吸收峰进行定性识别的准确性, 进一步对后续的精确定量分析形成干扰, 故本文结合量子化学分析方法对上述特殊吸收现象进行了进一步的分析。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 苯甲酸、 山梨酸、 苯甲酸和山梨酸混合物的太赫兹实验谱Fig.1 Benzoic acid, sorbic acid, and benzoic acid-sorbic acid mixture terahertz experimental spectra

首先, 从剑桥晶体数据库中心(Cambridge Crystallographic Data Centre, CCDC)提取苯甲酸和山梨酸的单分子及晶胞构型, 并利用Gaussian 09W展开理论模拟谱的计算。 结果表明, 二者单分子构型的吸收模拟谱与实验谱的特征吸收峰不符, 这是由于单分子构型的结构过于单一, 无法充分体现其在自然条件下的微观状态^[10]。 相比之下, 晶胞构型中包含多个单分子构型的排列, 更全方位地分析思考了分子内化学键、 分子间的氢键及范德华作用的影响, 能更准确地模拟分子间和分子内的相互作用, 所以采用晶胞构型进行模拟谱的计算获取。

其次, 采用苯甲酸与山梨酸的晶胞构型, 基于杂化泛函的B3LYP(Becke, 3-parameter, Lee-Yang-Parr)^[11]结合Grimme的DFT-D3(BJ)^[12]色散矫正和基组6-311g对两种物质结构进行优化和频率的计算, 获得了它们的太赫兹理论模拟吸收谱, 见图2。 由图2(a)中可以看出苯甲酸的模拟谱和实验谱中吸收峰个数相同, 峰位相似。 由图2(b)可知, 山梨酸的模拟谱与实验谱基本吻合, 均有1个极为显著的吸收峰, 且峰位相似。 这表明本文选择的晶胞构型和计算基组是合适的。 而峰位上存在的偏差是由于理论模拟谱是在绝对零度(0 K)下计算获取的, 但实验谱是在常温(298 K)下测量得到的, 随着温度的升高, 热膨胀效应使两个原子之间的键长增加, 分子间距离也随之增大, 导致晶胞的晶格参数(如晶胞体积)发生改变。 随着晶胞体积的增加, 分子间相互作用的力度会降低, 导致实验光谱与模拟光谱在吸收峰峰位上出现偏移^[13]。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 苯甲酸(a)和山梨酸(b)实验谱、 晶胞构型模拟谱对比Fig.2 Comparison of experimental spectra and simulated spectra for benzoic acid (a) and sorbic acid (b)

为进一步研究苯甲酸与山梨酸混合物的吸收峰形成机理, 本文进行了混合物的模拟谱计算。 由于CCDC数据库中没有苯甲酸和山梨酸混合物的晶胞构型, 故本文利用Molclus中的Genmer^[14]组件构建了500个期望初始团簇构型。 以苯甲酸为团簇中心, 加入山梨酸单质, 通过计算获得了含有500个初始构型的能量值和分子结构坐标信息的xyz文件, 但在这些初始团簇构型中存在一些结构相似, 且不合理的构型, 故先用专门针对原子之间的弱相互作用的半经验方法PM6-DH+在MOPAC^[15]下进行了预优化, 优化完成后按结构能量由低到高排序, 筛选出124个团簇构型。 在此基础上再用密度泛函理论(density functional theory, DFT)-B3LYP-6-311g^{* *} 结合D3矫正的方法进行结构优化计算, 筛选得到50个苯甲酸和山梨酸混合物的团簇构型。 最后基于结构的稳定性要求, 按照结构能量从低到高排序选取了15个团簇构型, 排除掉其中出现虚频的构型, 最终获得图3所示的13个团簇构型作为模拟谱计算的输入构型。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 苯甲酸和山梨酸混合物13个团簇构型Fig.3 Thirteen cluster configurations of benzoic acid and sorbic acid mixture

利用Gaussian 09对上述的13个团簇构型分别在DFT-B3lyp-6-311g^{* *} 水平下进行频率计算, 获取团簇构型在0.4~2.5 THz的模拟谱, 表2为13个团簇构型吸收峰的汇总信息。

表2

Table 2

表2(Table 2)

表2 团簇构型Cluster1-13(sim)的吸收峰统计表 Table 2 Absorption peak statistics for cluster structures Cluster1-13(sim) Number Absorption peak Peak position of absorption/THz Number Absorption peak Peak position of absorption/THz Cluster1(sim) 1 0.543 Cluster2(sim) 4 0.929, 1.261, 1.899, 2.353 Cluster3(sim) 2 1.531, 2.252 Cluster4(sim) 2 1.710, 1.916 Cluster5(sim) 4 0.620, 0.921, 1.710, 2.381 Cluster6(sim) 3 0.820, 1.673, 2.303 Cluster7(sim) 5 0.744, 1.030, 1.576, 1.769, 2.193 Cluster8(sim) 3 1.068, 1.479, 1.975 Cluster9(sim) 3 0.984, 1.685, 2.118 Cluster10(sim) 4 0.623, 1.013, 1.618, 2.282 Cluster11(sim) 3 1.026, 1.505, 2.101 Cluster12(sim) 2 0.673, 1.778 Cluster13(sim) 2 0.917, 2.303

表2 团簇构型Cluster1-13(sim)的吸收峰统计表 Table 2 Absorption peak statistics for cluster structures Cluster1-13(sim)

由图1可知, 苯甲酸和山梨酸混合物实验谱中包含4个吸收峰, 由表2可知, Cluster2(sim)、 Cluster5(sim)、 Cluster10(sim)构型的模拟谱吸收峰个数与之前图1二者混合物实验谱中吸收峰数量吻合, 但根据二者分别的模拟谱可知, 二者混合物的吸收谱应包含5个特征吸收峰, 表2中只有Cluster7(sim)构型的模拟谱具有5个吸收峰, 且峰位分别为0.744、 1.030、 1.576、 1.769和2.193 THz, 与实验谱峰位基本吻合, 如图4所示。 故因此推断Cluster7(sim)团簇构型更接近实际混合物样品的晶体结构。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 基于混合物7号团簇构型的模拟谱与实验谱对比Fig.4 Comparison of simulated and experimental spectra based on the Cluster7(sim) configuration of the mixture

有研究指出, 太赫兹范围内的特征吸收峰是由单个分子或晶体单元的低频振动引起的, 这些振动是由分子间或分子内部的键长、 键角和扭转角的变化所致。 且振动模式和旋转方式的不同也会影响吸收峰的位置和强度。 势能分布(PED)是一种分析方法, 用于分解分子的振动模式, 能够确定各基团在振动模式中的贡献值, 从而加深对振动模式特征的理解, 并帮助识别出主导的振动模式。 本文使用VEDA4对以上研究的晶胞构型的模拟谱吸收峰的振动模式进行PED分析, 完成模拟谱特征吸收峰的振动模式匹配^[16]。 计算所得的苯甲酸、 山梨酸、 与Cluster7 (sim)团簇构型贡献率和振动模式见表3和表4。 表中, STRE表示键长伸缩, BEND表示键角弯曲, TORS表示二面角扭转, 括号内的数字为该振动模式占有的贡献百分比。 通过利用VMD^[17]软件显示分子的振动分量, 可更直观地显示苯甲酸和山梨酸晶胞内各个原子的振动方向, 如图5所示。 图中黄色箭头表示苯甲酸和山梨酸分子分别在特征吸收峰处的振转模式。 由图5(a)可知, 苯甲酸在0.627 THz处的振动模式经分析可指认为贡献率为14%的46O45H所在平面的键长伸缩。

表3

Table 3

表3(Table 3)

表3 苯甲酸、 山梨酸及其团簇构型振动模式PED分析 Table 3 PED analysis of benzoic acid and sorbate vibration modes Sample Atomic numbering bitmap Experimental frequency/THz Calculation frequency/THz PED analysis benzoic acid 0.619 0.627 STRE: 46O45H(14) 0.916 0.980 TORS: 16O18C19C20C(34) 1.531 1.320 STRE: 46O45H(21) BEND: 32O45H46O(-11) TORS: 35C34C33C32O(-16) 2.069 1.967 STRE: 46O25H(-35) TORS: 33C32O45H46O(-25) sorbic acid 1.603 1.446 BEND: 22H35O48C(10) TORS: 42C44C46C48C(-26)

表3 苯甲酸、 山梨酸及其团簇构型振动模式PED分析 Table 3 PED analysis of benzoic acid and sorbate vibration modes

表4

Table 4

表4(Table 4)

表4 苯甲酸和山梨酸Cluster7(sim)团簇构型振动模式PED分析 Table 4 PED analysis of Benzoic acid and sorbate Cluster7(sim) cluster configuration vibration mode Sample Atomic numbering bitmap The corresponding experimental spectra of benzoic acid and sorbic acid/THz Calculation frequency/THz PED analysis Mixture (Benzoic acid sorbic acid) 0.619 0.744 BEND: 25C27C29C(16) TORS: 17H16O15H2O(36) 0.916 1.030 BEND: 25C27C29C(-22) TORS: 15H2O3C4C(-16) 1.531 1.576 STRE: 16O15H(15) BEND: 25C27C29C(-12) TORS: 23C25C27C29C(-16) 1.603 1.769 STRE: 16O15H(-16) TORS: 17H16O15H2O(49) 2.069 2.193 STRE: 16O15H(21) BEND: 37C16O15H(10) TORS: 5C4C3C2O(-14) TORS: 23C25C27C29C(16)

表4 苯甲酸和山梨酸Cluster7(sim)团簇构型振动模式PED分析 Table 4 PED analysis of Benzoic acid and sorbate Cluster7(sim) cluster configuration vibration mode

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 振动模式 (a)、 (b)、 (c)、 (d)分别为苯甲酸在0.627、 0.980、 1.320和1.967 THz处的振动模式; (e)为山梨酸在1.446 THz处的振动模式Fig.5 Vibration modes (a), (b), (c) and (d) of benzoic acid at 0.627, 0.980, 1.320, 1.967 THz, respectively; Vibration modes (e) of sorbic acid at 1.446 THz

由表3、 表4和图5发现苯甲酸和山梨酸的太赫兹吸收峰是复杂的分子内和分子间相互作用和各自分子的多种振动模式共同作用的结果。 这些模式包括羧基之间的氢键伸缩和弯曲振动, 羧基的C=O和C— O键伸缩振动, 苯环的呼吸和扭转振动, 以及C=C双键的伸缩振动、 C— H键的伸缩与弯曲振动, 还有O— H…O氢键间的键角伸缩和弯曲。 这些分子的振动模式可以通过吸收峰的位置信息进行区分。 通过对比两者的吸收峰位置和强度, 可以对吸收峰形成机理实现有效解析。

为了进一步考察分子内外相互作用, 解析特征吸收峰形成机理, 本文利用基于Hirshfeld原子空间划分的独立梯度模型(IGMH)对苯甲酸和山梨酸分子之间和内部的相互作用类型、 位置等进行了可视化分析。 相比于约化密度函数(reduced density gradient, RDG), IGMH更加直观且可以自由选取片段, 分析特定片段间或片段内的相互作用, 避免无关区域的等值面妨碍, 容易获得高效果的等值面, 耗时也更低^{[18, 19]}。 IGMH函数表达式如式(3)和式(4)

δg(r)=∑i|∇ρi(r)|-∑i∇ρi(r)(3)

δg(r)=δginter(r)+δgintra(r)(4)

式中, r为坐标矢量, i为原子序号, ∇为梯度算子, ρ _i代表i原子的电子密度, |∇ρ _i(r)|为取矢量的模。 δ g可以对结构中各原子间特定片段之间(interfragment)和内部(intrafragment)的相互作用进行直观地展现, 并将其分割为δ g_inter和δ g_intra。

通过将ρ (r)与电子密度Hessian矩阵中的第二大本征值λ ₂的符号函数sign(λ ₂)相结合, 所得的乘积sign(λ ₂)ρ (r)可以映射到具有不同颜色的IGMH等值面上, 以此直观地揭示出共价作用、 弱相互作用的具体方位、 类型等^[20]。 色彩梯度设定为蓝色过渡到绿色, 再过渡到红色, 其中蓝色代表ρ (r)> 0, λ ₂< 0的区域, 这通常代表存在吸引作用; 绿色部分代表ρ (r)≈ 0, λ ₂≈ 0的位置, 一般为弱氢键作用或范德华作用; 红色代表λ ₂大于零的区域, 这通常关联于笼状结构或环状结构中较强的立体阻碍效应, 即位阻作用。 图6和图7展示的用波函数分析软件Multiwfn^[21]和VMD计算得到的结果。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 (a)苯甲酸填色IGMH散点图; (b)苯甲酸填色IGMH等值面图Fig.6 (a) Benzoic acid coloring IGMH scatter diagram; (b) benzoic acid coloring IGMH contour map

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 (a)山梨酸填色IGMH散点图; (b)山梨酸填色IGMH等值面图Fig.7 (a) Sorbic acid coloring IGMH scatter diagram; (b) Sorbic acid coloring IGMH contour map

图6(a)显示了苯甲酸晶胞的IGMH散点图, 图6(b)为其IGMH等值面图。 通过与PED分析结合, 针对苯甲酸吸收峰成因的特定原子进行了IGMH分析。 可以看出, 图中出现了多个不同位置的尖峰, 蓝色尖峰位于-0.60~-0.20 a.u, 体现出成键区域的电子密度较大, 这是共价键作用区域的典型特征。 结合图6, 可以推测该现象主要由C— C键(如19C— 24C、 20C— 21C、 34C— 39C、 35C— 36C等)、 C— H键(如20C— 25H、 34C— 40H等)及羧基中的C=O双键(48C=46O、33C=31O等)和其他共价键作用(如33C…32O— 45H、 18C…17O— 30H)导致。

在-0.04~0.04 a.u区域, 出现了较小的绿色尖峰, 这是弱相互作用区域, 主要归因于羧基中的O— H键和C=O双键与周围分子中的C— H键形成的弱氢键作用(如48C=46O…25H、18C=16O…55H、32O— 45H…44H等)。 此外, 苯环之间还存在π — π 堆积相互作用, 这是分子间芳香环之间的π 电子云相互吸引所导致的。

同时, 绿色区域也反映了范德华力的贡献, 尤其是在氢键和π — π 堆积较弱的区域, 范德华色散力通过分子间的瞬时偶极相互作用进一步增强了分子间的吸引力。 IGMH等值面图中对应化学键的区域周围有红色部分, 但这些区域没有明确的物理意义。

图7(a)为山梨酸晶胞IGMH散点图, (b)为山梨酸晶胞IGMH等值面图。 结合PED分析可知, 其山梨酸位于1.603 THz的太赫兹吸收峰成因可归结为: 羧基的C=O双键作为氢键受体与周围分子的C— H形成较弱的氢键(48C=35O…25H 等)作用, 这些氢键较弱, 但仍对分子间的空间排布和晶体的稳定性产生了一定影响。 其次, 山梨酸的分子骨架中存在多个C=C双键(如0C=42C), 这些作为共价键不仅维持了分子本身的稳定结构, 还通过分子振动参与了特定频率的太赫兹吸收。 与此同时, 山梨酸晶胞中的C— H键(42C— 43H、 44C— 45H、 46C— 47H)和羧基(48C…33O— 34H)也起到了共价键的作用, 进一步增强了分子结构的刚性。 分子间的范德华作用也广泛存在, 增强了晶体结构的稳定性。 IGMH图中的绿色区域代表了这些弱氢键和范德华力的作用, 它们共同维持了山梨酸晶体的有序性。

虽然IGMH对上述苯甲酸和山梨酸的晶胞构型进行了相互作用的可视化分析, 但在相同的等值面数值设置下难以清晰展示各种强度的相互作用, 而且需要高时长的计算。 故本文引入相互作用区域指示函数(IRI), 其既可以反映分子间的相互作用, 又可以反映出原子间的化学键^[22]。 IRI 的函数定义如式(5)所示

IRI=|∇ρ(r)|[ρ(r)]α(5)

式(5)中; ρ (r)是电子密度, ∇ρ (r)是梯度矢量, α 是可变参数。 与IGMH分析相似, IRI方法也是通过将sign(λ ₂)ρ (r)函数在IRI等值面用不同颜色映射, 从而实现对不同区域间相互作用的识别。 因此, 色彩的参考标准与上述IGMH 是相同的。

图8是苯甲酸和山梨酸混合物Cluster7(sim)团簇构型的IRI填色散点图和等值面图。 从图8可以看出, 两个单分子羧基之间交互相错, 在此区域中C、 H、 O之间形成了不一的弱相互作用, 参见图8(b)中大面积的黄绿色等值面。 苯甲酸单分子羧基上的15H与山梨酸单分子羧基的16O间等值面也为绿色, 说明存在弱氢键作用, 31C=18O…28H、3C=1O…10H、2O— 15H…14H等之间呈现绿色和红色并存的区域, 说明同时存在吸引的弱氢键作用和一定程度的位阻效应。 IRI等值面中位于C— C键、 C— H键和C=O键之间的蓝色圆柱状物, 表示此处存在强的相互作用, 即化学键作用。 苯环中的红色区域说明存在较强的位阻效应, 且两个单分子间出现的大片微红绿色, 说明此处存在范德华作用。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 (a)Cluster7(sim)团簇构型填色IRI散点图; (b)Cluster7(sim)团簇构型填色IRI等值面图Fig.8 (a) Cluster7(sim) configuration coloring IRI scatter diagram; (b) Cluster7(sim) configuration coloring IRI contour map

结合团簇振动模式分析可知, 混合物团簇构型在1.576 THz处出现的吸收峰, 可归因于分子之间的氢键振动引起的, 其大部分振动出现在两个分子羧基之间的16O和15H的氢键上, 由于团簇构型改变了羧基的键长, 导致该频段出现了吸收峰。