Terahertz fingerprint spectrum and weak interaction analysis of nicotinamide–sorbic acid cocrystal

Abstract

Sorbic acid (SA) is a naturally occurring food additive that plays an important role in the food industry. The study of nicotinamide–sorbic acid (NIA–SA) cocrystal is significant for improving its dissolution rate, bioavailability, and physical stability. In this study, a liquid-assisted grinding technique was used to prepare the NIA–SA binary cocrystal compound, and powder X-ray diffraction and terahertz time-domain spectroscopy (THz-TDS) confirmed that this method can successfully obtain the cocrystal. Solid-state density functional theory (ss-DFT) calculations were utilized to obtain theoretical spectra, revealing the structural characteristics of the NIA–SA cocrystal and the relationship between terahertz absorption peaks and specific collective vibration modes. Subsequently, the reduced density gradient (RDG) analysis was employed to analyze energy interactions within the molecular system, identifying the characteristic absorption peaks associated with cocrystal formation and explaining their origins. The results demonstrate that combining THz-TDS with ss-DFT can comprehensively and dynamically reveal intermolecular interactions during the formation of cocrystals. The unique sensitivity of THz-TDS to low-frequency vibrations provides robust theoretical and experimental support for the design and development of novel cocrystal materials, highlighting the significant advantages and broad application prospects of this technique in cocrystal characterization.

---

1. Introduction

Preservatives extend the shelf life of food by inhibiting microbial activity, playing a crucial role in the food industry. Common preservatives include benzoic acid, sorbic acid, and its potassium salt (potassium sorbate).¹ Sorbic acid (C₆H₈O₂, CAS: 110-44-1, abbreviated as SA; Fig. 1) is a short-chain unsaturated fatty organic acid that participates in normal human metabolism and is a naturally occurring food additive. SA and its salts exhibit significant inhibitory effects against most bacteria, molds, and yeasts, and are commonly used as preservatives in food, animal feed, cosmetics, and pharmaceuticals. The antimicrobial activity of SA primarily depends on its undissociated carboxylic acid form, so increasing the concentration of the acid molecules can enhance its antimicrobial properties. However, under alkaline conditions, the reduced solubility of the carboxylic acid in water diminishes the antimicrobial effectiveness of SA. Therefore, increasing the water solubility of SA can effectively broaden its range of applications.² Nicotinamide (C₆H₆N₂O, CAS: 98-92-0, abbreviated as NIA; Fig. 1), also known as Vitamin B3, is highly soluble in water and is a commonly used dietary supplement. NIA is widely found in both plant and animal foods and has anti-inflammatory and antioxidant properties. NIA and SA can form a cocrystal through hydrogen bonding, thereby improving the physical properties of SA.³ In 2023, Chen et al. synthesized cocrystals and their salts of SA, NIA, piperazine (PIP), and 2,6-diaminopyridine (2,6-DAP) for the first time, and characterized them using XRD, FT-IR, DSC, and TG techniques.^4–6 They determined the solubility of SA and its multicomponent crystals in aqueous solution at 37 °C. The results showed that the solubility of the SA multicomponent crystals in water was significantly enhanced.⁷ Although this cocrystal structure has been well characterized, the intermolecular interactions responsible for cocrystal formation have not yet been fully explained.

Fig. 1 Chemical structure of nicotinamide and sorbic acid.

To comprehensively characterize the structural features and formation mechanisms of cocrystals, various characterization techniques have been widely applied by researchers. Among them, terahertz time-domain spectroscopy (THz-TDS) is particularly effective for probing low-frequency intermolecular vibration modes (such as intermolecular vibrations and lattice vibrations) by measuring the absorption and transmission properties of samples in the 0.1–10 THz frequency range. This technique allows for the effective analysis of changes in molecular structures, crystal structures, and intermolecular interactions in solids, demonstrating significant advantages over other characterization methods.⁸ In contrast to X-ray diffraction (XRD), which is limited to the analysis of static crystal structures, terahertz spectroscopy not only provides crystal structural information but also reveals the dynamic changes in weak intermolecular interactions (such as hydrogen bonds and van der Waals forces) during cocrystal formation.⁴ While Fourier transform infrared (FTIR) and Raman spectroscopy can detect the vibrational characteristics of functional groups, they are primarily sensitive to high-frequency vibration modes (such as C–H and O–H stretching vibrations).^9–11 Terahertz spectroscopy, on the other hand, offers higher resolution for low-frequency vibration modes, providing a more direct reflection of the changes in intermolecular interactions during cocrystal formation. Furthermore, although thermogravimetric analysis (TGA) and differential scanning calorimetry (DSC) can assess the thermal stability and phase transition behavior of cocrystals, they do not provide direct information on intermolecular interactions or crystal structures.^6,12 In contrast, terahertz spectroscopy offers a nondestructive means to simultaneously provide dielectric response and molecular vibration information of the cocrystals.

Based on the above advantages, THz-TDS has gradually become a novel and sensitive tool for studying molecular structures and understanding collective vibrational modes between molecules, particularly in investigating cocrystal polymorphism behaviors. It has proven useful in areas such as distinguishing drug polymorphs, characterizing pharmaceutical cocrystals, and differentiating chiral compounds from their racemic products.^13,14 In 2018, Wu et al. used THz-TDS to measure the characteristic absorption spectra of two structural isomers of dimethylurea (DMU) in the 0.6–1.8 THz range at room temperature. They also employed reduced density gradient (RDG) analysis to investigate the positions and types of intermolecular hydrogen bond interactions in the crystals of 1,1-DMU and 1,3-DMU. The results confirmed that THz-TDS can serve as an effective method for identifying crystal structural isomers and detecting intermolecular hydrogen bond interactions.¹⁵ In 2020, Li et al. investigated the cocrystals of indomethacin–nicotinamide (IND–NIC) and indomethacin–saccharin (IND–SAC), demonstrating that terahertz vibrational spectroscopy is a reliable tool for cocrystal detection.¹⁶ In 2023, Zhao et al. used a broadband air plasma THz-TDS system to obtain the terahertz fingerprint spectrum of adenosine in the range of 0.5 to 13 THz. Additionally, the density functional theory (DFT) and RDG analysis were employed to study the vibrational modes and weak intermolecular interactions of adenosine molecules. The results indicated that different absorption peaks in the THz spectrum correspond to specific vibrations of the molecules, and hydrogen bonds play a significant role in the vibrations in the terahertz range.¹⁷ Between 2022 and 2024, Du's research group carried out a series of systematic investigations into the construction and characterization of ethenzamide (ETZ) pharmaceutical cocrystals.¹¹ Wan et al. employed THz-TDS and Raman spectroscopy in combination with quantum chemical calculations to elucidate the structural features of hydrogen-bonding synthons in cocrystals formed between ETZ and various polyhydroxy acids, highlighting the crucial role of hydroxyl group number and position in supramolecular structure construction.¹⁸ Zhang et al. synthesized two polymorphic forms of the ETZ–ethylmalonic acid cocrystal and, through THz and Raman spectroscopy along with thermal phase transition and mechanical grinding studies, revealed the regulatory mechanisms of solvent polarity and mechanical force on polymorph formation pathways.¹⁹ Jing et al. compared two polymorphic forms of the ETZ–saccharin cocrystal, finding that although the primary hydrogen-bonding motifs were identical, differences existed in the extended hydrogen-bonding networks, and that THz spectroscopy outperformed Raman spectroscopy in distinguishing between the polymorphs.²⁰ In 2024, Su et al. measured XRPD, THz-TDS, and Raman spectra of salicylic acid, aspirin, their physical mixtures, and cocrystals. They also used DFT to simulate the structures and spectroscopic properties of salicylic acid and aspirin cocrystal molecules, which led to the assignment of molecular vibration modes.²¹ Wang et al. calculated the potential energy to distinguish the stability of different polymorphs of the pyrazinamide–malic acid (PZA–MA) cocrystal and matched the corresponding vibrational modes to the terahertz fingerprint spectrum. They combined terahertz spectroscopy with quantum chemical calculations to explore the relationship between THz absorption peaks and specific collective vibrational modes.²² The above works combine spectroscopy with quantum chemical calculations based on crystal structure, laying a theoretical foundation for the establishment of a terahertz spectroscopy database and the identification of biochemical molecules. Analyzing terahertz fingerprint spectra can provide information on lattice structures and internal interactions of substances. However, further in-depth research is needed to address questions about cocrystal formation and the detailed correspondence between terahertz fingerprint spectra and molecular vibration features.

This work combines terahertz spectroscopy experiments with first-principles calculations to systematically reveal the newly emerging low-frequency vibration modes and their intermolecular interaction mechanisms in the NIA–SA cocrystal system. In section 3.1 and 3.2, the cocrystals were obtained through liquid-assisted grinding, and their successful synthesis was confirmed by powder X-ray diffraction (PXRD) and terahertz spectroscopy. NIA, SA, and the NIA–SA cocrystal exhibit unique spectral features, allowing for easy identification of the NIA–SA cocrystal and its parent substances based on their spectral fingerprint characteristics. In section 3.3 and 3.4, solid-state DFT (ss-DFT) calculations were first employed to obtain the theoretical spectra, providing deep insights into the structural characteristics of the NIA–SA cocrystal and the relationship between the terahertz absorption peaks and specific collective vibration modes. Subsequently, the RDG method was used to further analyze the energy interactions within the NIA–SA cocrystal molecular system, identifying the most characteristic absorption peaks in the cocrystal formation and offering explanations for their origins. THz-TDS combined with ss-DFT provides reliable theoretical support and experimental validation for studying cocrystal formation and weak interactions.

2. Material and method

2.1 Material preparation

SA and NIA were purchased from J&K Scientific with purities not less than 98%. To minimize the effect of particle scattering on the experimental results, the samples were ground into fine powders in a mortar to meet the tablet compression requirements. At room temperature, 30 mg of each of NIA, SA, a physical mixture of NIA and SA, and a binary cocrystal of NIA–SA were mixed thoroughly with 120 mg of polyethylene powder using a pestle and mortar. The resulting mixture was then compressed into tablets at 10 MPa pressure. The mixture was pressed at 10 MPa for 5 minutes to form tablets with a diameter of 13 mm and thicknesses of 1.311 mm for NIA, 1.308 mm for SA, 1.353 mm for the physical mixture, and 1.540 mm for the cocrystal compound. Polyethylene was used as a substrate since it has a low absorption coefficient in the terahertz range. To ensure accuracy and reproducibility of the results, the thickness of the samples was averaged from three measurements taken at different positions. Each spectrum was obtained from 1200 scans. The reference signal (i.e., dry air as the background signal) and the sample signal were measured three times, and the average of these measurements was used to obtain the THz spectra for data analysis.

2.2 Cocrystal synthesis method

In this experiment, a liquid-assisted grinding technique was used to prepare the NIA–SA binary cocrystal compound. The molar mass of NIA is 122.125 g mol⁻¹, while that of SA is 112.13 g mol⁻¹. NIA and SA were added to 0.4 mL of ethanol solution in a 1 : 1 molar ratio and then ball-milled for 2 hours. The resulting mixture was placed in a drying oven for 24 hours to obtain the powdered cocrystal compound.

2.3 Terahertz time-domain spectroscopy (THz-TDS)

The samples were tested using the TAS7400 SU terahertz time-domain spectrometer from Advantest, Japan. This system supports ultra-broadband spectral analysis, with a dynamic frequency range of 0.5 to 7 THz. It offers various spectral analysis modes, including transmission, reflection, and attenuated total reflection (ATR), with transmission mode being used in this experiment. The spectrometer is equipped with an external dry air unit to eliminate water vapor absorption interference. The sample tablet was placed in the sample chamber of the spectrometer and continuously purged with dry air during detection to maintain relative humidity below 1.5%.

2.4 X-ray diffraction (XRD)

The samples were tested using a MiniFlex600 X-ray diffractometer. The scanning angle range was set from 5° to 40°, with a scanning speed of 4° per minute. Subsequently, the XRD patterns were simulated using Mercury software in conjunction with crystal structures from the Cambridge Crystallographic Data Centre (CCDC).^23,24

2.5 DFT-based calculations and analysis of weak interactions

The crystal structures of NIA, SA, and their cocrystal NIA–SA were obtained from the Cambridge Crystallographic Data Centre (CCDC) (Structural identification numbers: 131757,³ 1206459,²⁵ 2046374 (ref. 7)). In the asymmetric unit, SA molecules and NIA molecules form a cocrystal via N–O⋯H interactions. The NIA–SA cocrystal belongs to the triclinic system, with space group P1̄(2), as shown in Fig. 2. The calculations for terahertz spectroscopy were performed using the ss-DFT software CP2K.²⁶ The specific computational method selected was the Perdew–Burke–Ernzerhof (PBE) functional²⁷ within the generalized gradient approximation (GGA), combined with the pob-DZVP-rev2 basis set,²⁸ and supplemented with the D3 dispersion correction with BJ damping to describe the weak interactions within the crystal structure.^29,30 To balance computational accuracy and cost, we constructed supercells of 3 × 1 × 2, 1 × 3 × 1, and 3 × 2 × 1 for the three crystals and performed convergence tests on their energies, determining that the cut-off should not be less than 600 Ry and the rel-cut-off should not be less than 55 Ry. We first optimized the supercells to find the local minima on the potential energy surface. The geometry optimization was performed using the linear search method from the BFGS (Broyden–Fletcher–Goldfarb–Shanno) quasi-Newton algorithm, with a maximum displacement set to 0.003 Bohr, maximum force set to 1.0 × 10⁻⁵ Hartree/Bohr, root mean square displacement set to 1.5 × 10⁻³ Bohr, and root mean square force set to 3.0 × 10⁻⁴ Hartree/Bohr. To ensure the accuracy of low-frequency mode calculations, the self-consistent field convergence threshold was set to 1 × 10⁻⁸ Hartree. After geometry optimization was completed, the optimized structure was subjected to vibrational analysis using the same computational parameters.

Fig. 2 Structural formulas of NIA (a), SA (b), and their cocrystal NIA–SA (c).

Weak interactions are ubiquitous in various chemical systems and exist in many specific forms, such as hydrogen bonds (H-bonds), halogen bonds (X-bonds), van der Waals interaction, π–π stacking, steric effects, and more. Weak interactions are currently a hot topic in theoretical and computational chemistry and are key factors that need to be carefully considered in the fields of crystal engineering and drug design.³¹ The non-covalent interaction (NCI) method, also known as the RDG method, is a commonly used approach to study weak interactions. It was proposed by the research group of Johnson and Yang and is developed based on the atoms in molecules (AIM) method to describe deviations from uniform electron distribution.^32–34 This method uses RDG isosurface to reflect the types and strengths of non-covalent interactions, with RDG being a dimensionless form of the electron density gradient, expressed as follows:³⁵

Here, ∇ is the gradient operator, ρ(r) is the electron density, and |∇ρ(r)| is the norm of the electron density gradient. The value of ρ(r) can only reflect the strength of weak interactions, not their types. Therefore, the sign(λ₂) function which represents the second eigenvalue of the electron density Hessian matrix is introduced. By multiplying ρ(r) by sign(λ₂) and projecting the result onto the RDG isosurface, one can reflect both the types and strengths of non-covalent interactions.

To analyze the nature of the intermolecular weak interactions in NIA, SA, and their cocrystal NIA–SA, this study employed the wave function analysis software Multiwfn3.8(dev) (ref. 2023-Oct-8) for RDG visualization.^36,37 The input files were created using the wave function files after geometric optimization, with the same basis set and functional as before. In the calculations, the isosurface was set to 0.5, the grid spacing was set to 0.15, and RDG_maxrho was set to 0.05. Both the RDG scatter plots and isosurfaces were color-coded using RGB, with the color scale calibrated to blue-green-red, and the color scale range defined as (−0.035, −0.0075, 0.02). Additionally, VMD 1.9.3 (ref. 2016-Nov-30) software³⁸ was used to obtain scatter plots and isosurfaces of the non-covalent interactions, while gnuplot 6.0 software was employed to create the colored scatter plots. In this paper, we ignore the regions close to the atomic nuclei and the molecular edges, focusing only on the areas of weak interactions, specifically those near chemical bonds and non-covalent bonds. The values of |∇ρ(r)| and RDG in these two regions are relatively small, but ρ(r) differs. By combining the RDG values and ρ(r), we can identify the regions with weak interactions within the molecular system.

3. Results and discussions

3.1 XRD spectral characterization and analysis

Fig. 3 shows the diffraction peaks of NIA, SA and the NIA–SA cocrystal measured using an X-ray diffractometer and simulated with Mercury. It can be seen that the positions of the diffraction peaks in the experimental and simulation results correspond well, and they are consistent with the data previously published by Chen et al.⁷ For comparison, Table 1 lists the specific positions (2θ values) of the XRD diffraction peaks for the three substances obtained from both experimental and simulated results. As shown in Table 1, compared to the single crystals of NIA and SA, the NIA–SA cocrystal presents new diffraction peaks at 8.798, 14.040, 17.600, 26.260, and 27.021, which clearly confirms the formation of the cocrystal. The use of XRD allows for the characterization and analysis of the structures of NIA, SA and the NIA–SA cocrystal, laying a foundation for further research in THz-TDS.

Fig. 3 The XRD diffraction peaks of the NIA (red line), SA (green line), and NIA–SA cocrystal (blue line) cocrystal measured and calculated at room temperature.

Table 1 The comparison of XRD diffraction peak positions of NIA, SA, and NIA–SA cocrystals between experimental and simulated results

Material		2θ (deg)
NIA	Cal.	11.317	14.738	19.503	22.219	23.341	25.363	25.837	27.302	36.961
	Exp.	11.321	14.841	19.540	22.280	23.402	25.381	25.819	27.221	36.961
SA	Cal.	11.539	13.200	23.161	23.54	24.00	25.417	28.402
	Exp.	11.479	13.019	22.82	23.44	—	24.861	29.581
NIA–SA	Cal.	8.868	—	12.461	13.501	14.182	15.981	17.798	19.201	20.323
		22.876	24.1	25.138	25.459	26.223	26.499	26.937	27.96
	Exp.	8.798	11.477	12.32	13.441	14.04	15.94	17.6	19.119	20.081
		22.8	23.839	24.78	25.381	26.26	—	27.021	27.723

---

3.2 Terahertz spectroscopy measurement

Fig. 4 shows the absorption coefficient and refractive index of NIA, SA, their 1 : 1 mixture, and the NIA–SA cocrystal measured using a terahertz time-domain spectrometer in the range of 0.5 to 4 THz. To improve the accuracy of the analysis, we filtered out weaker absorption peaks caused by impurities and noise. As shown in the figure, NIA has absorption peaks at 0.595, 1.061, 1.923, 2.396, 2.922, 3.235 and 3.677 THz, while SA has absorption peaks at 1.669, 2.164, 2.538, 2.873 and 3.386 THz, which are consistent with previous reports.^39,40 The absorption peaks of the 1 : 1 mixture of NIA and SA occur at 0.602, 1.088, 1.945, 2.13, 2.436, 2.925, and 3.699 THz, indicating a linear superposition of the THz spectra from the NIA and SA components, which suggests that a simple physical mixture does not form new substances or chemical bonds. In contrast, the NIA–SA cocrystal exhibits new absorption peaks at 1.197 THz and 2.44 THz, indicating that NIA and SA interact through non-covalent bonds to form a new phase that is distinct from the parent constituents.

Fig. 4 Absorption coefficients and refractive index of NIA, SA, their 1 : 1 mixture, and the NIA–SA cocrystal in the range of 0.5 to 4 THz. ab: absorption coefficient; re: refractive index.

Additionally, the presence of impurities such as NIA and SA in the crystalline materials may have contributed to the observation of absorption peaks in the spectra of both the physical mixture and the cocrystal that resemble those of the parent constituents. Combined with the quantum chemistry theory calculation, we believe that the absorption peak of 1.045 THz in the NIA–SA cocrystal may be the residue of NIA monomer. Since terahertz spectroscopy can sensitively detect low-frequency molecular lattice vibrations, the experimental results can preliminarily confirm that NIA and SA interact through non-covalent bonds, particularly intermolecular hydrogen bonds. This interaction results in the formation of new absorption peaks in the THz fingerprint spectrum of the cocrystal, distinguishing it from both the parent materials and the physical mixture. The combination of experimental results and theoretical analysis indicates that the formation of hydrogen bonds in the cocrystal can alter the molecular structure and vibrational modes of certain functional groups, thereby changing the distribution of characteristic peaks in the vibrational spectrum.

3.3 Quantum chemical calculations

Quantum chemical calculations have proven to be an effective tool for interpreting experimental THz absorption spectra. To investigate how intermolecular forces function in cocrystals, we selected the supercell rather than isolated molecules as the basic structure for theoretical simulations. Fig. 5 compares the simulated THz spectra of NIA, SA, and the NIA–SA cocrystal with their experimental THz spectra. In the simulated spectra, each black vertical line represents the frequency of a characteristic vibrational mode, and its length is proportional to the vibrational intensity. As shown in Fig. 5, there is a certain discrepancy in peak frequencies between the experimental spectra and the theoretical calculations. This is because the spectral features primarily arise from the collective internal vibrations of the molecules in the low-frequency terahertz range of 0–4 THz, which are typically caused by weak bonds or non-bonding interactions such as van der Waals forces, dispersion forces, and hydrogen bonds. Due to the relatively weak nature of these intermolecular interactions, they are easily influenced by external environmental factors such as temperature. The simulations are typically performed at absolute zero, while the experimental THz spectra are obtained at room temperature. Therefore, the temperature difference may cause a shift in the position of the characteristic peaks to some extent. In addition, the differences between the molecular conformations in the computational model and the actual molecular structure may also contribute to the peak shift.

Fig. 5 Comparison of the simulated THz spectra and experimental THz spectra of NIA, SA, and the NIA–SA cocrystal.

The vibrational modes in the terahertz region are governed by intermolecular interactions, which often result in complex intramolecular and intermolecular vibrations. Fig. 6–8 show the vibrational mode diagrams of NIA, SA, and the NIA–SA cocrystal at different vibrational frequencies, calculated using CP2K. The corresponding coordinate axes are indicated in the bottom-right corner of each diagram, with the a-axis in red, the b-axis in green, and the c-axis in blue. Table 2 specifies the assignment of each vibrational mode in NIA, SA, and the NIA–SA cocrystal.

Fig. 6 The vibrational mode diagrams of the NIA supercell calculated using CP2K at the vibrational frequencies of (a) 0.843 THz, (b) 1.59 THz, (c) 2.415 THz, (d) 2.619 THz, (e) 3.031 THz, (f) 3.156 THz, (g) 3.345 THz, and (h) 3.915 THz.

Fig. 7 The vibrational mode diagrams of the SA supercell calculated using CP2K at the vibrational frequencies of (a) 0.987 THz, (b) 1.482 THz, (c) 1.674 THz, (d) 2.01 THz, (e) 2.268 THz, (f) 2.835 THz, (g) 3.219 THz, (h) 3.444 THz, (i) 4.17 THz.

Fig. 8 The vibrational mode diagrams of the NIA–SA cocrystal calculated using CP2K at the vibrational frequencies of (a) 1.485 THz, (b) 1.89 THz, (c) 2.361 THz, (d) 2.811 THz, (e) 3.198 THz.

Table 2 Comparison of the absorption peak positions between the experimental and calculated THz spectra and vibrational modes for NIA, SA and NIA–SA cocrystal

	Exp. (THz)	Cal. (THz)	Type
s: stretching, δtw: twisting, w: wagging, v: rotation, r: rocking.
NIA	0.595	0.843	Partial vibration: mainly from s (pyridine ring) (lattice along c-axis)
	1.061	1.59	Collective: r (NIA) (in-plane a–b)
	1.923	2.415	Collective: v (pyridine ring) + s (–C=O–NH2) (in-plane b–c)
	2.396	2.619	Collective: w (pyridine ring) + w (–C=O–NH2) (lattice along c-axis)
	—	3.031	Collective: δtw (pyridine ring) + δtw (–C=O–NH2) (out-plane)
	2.922	3.156	Collective: v (pyridine ring) + v (–C=O–NH2) (in-plane b–c)
	3.235	3.345	Partial vibration: v (pyridine ring) + v (–C=O–NH2)
	3.677	3.915	Collective: w (pyridine ring) + w (–C=O–NH2) (lattice along a-axis)
SA	—	0.987	Collective: δtw (SA) (lattice along c-axis)
	—	1.482	Collective: s (SA) (lattice along positive a-axis)
	1.669	1.674	Collective: s (SA) (lattice along negative a-axis)
	2.164	2.01	Partial vibration: mainly from δtw (–CH3)
	2.538	2.268	Collective: r (SA) (in-plane a–c)
	2.873	2.835	Collective: s (SA) (lattice along c-axis)
	3.386	3.219	Collective: v (SA) (in-plane b–c)
	—	3.444	Collective: v (SA) (in-plane a–c)
	—	4.17	Partial vibration: mainly from w (–CH3)
NIA–SA	1.197	1.485	Collective: v (NIA) (in-plane b–c) + w (SA) (out-plane)
	1.719	1.89	Collective: δtw (NIA) (out-plane) + δtw (SA) (out-plane)
	1.928	2.361	Collective: w (NIA) + mainly from δtw (SA) (out-plane)
	2.442	2.811	Collective: v (NIA) (in-plane b–c) + s (SA) (in-plane b–c)
	2.824	3.198	Collective: δtw (NIA) + w (SA) (partial vibration)

---

Fig. 6 presents the vibrational modes of NIA at eight different frequencies, calculated using CP2K. The experimental peak at 0.595 THz corresponds to the calculated mode at 0.843 THz, attributed to the collective stretching vibration of the NIA pyridine ring along the c-axis of the lattice (Fig. 6a). The experimental peak at 1.061 THz matching the calculated mode at 1.59 THz, arising from the collective rocking vibration of NIA in the a–b plane (Fig. 6b). The experimental peak at 1.923 THz aligns with the calculated mode at 2.415 THz, caused by pyridine ring rotation and in-plane stretching of the amide group (–C=O–NH₂) in the b–c plane (Fig. 6c). The experimental peak at 2.396 THz corresponds to the calculated mode at 2.619 THz, resulting from the collective wagging of two adjacent pyridine rings and amide groups along the c-axis (Fig. 6d). The calculated mode at 3.031 THz is not prominent in experiments, possibly due to insufficient resolution of the terahertz spectrometer; Fig. 6e reveals it stems from out-of-plane collective twisting of the amide group and C–H twisting in the pyridine ring. The experimental peak at 2.922 THz matching the calculated mode at 3.156 THz, driven by co-rotational vibration of the pyridine ring and amide group in the b–c plane (Fig. 6f). The experimental peak at 3.235 THz corresponds to the calculated mode at 3.345 THz, due to partial rotation of the pyridine ring and amide group (Fig. 6g). Finally, the experimental peak at 3.677 THz aligns with the calculated mode at 3.915 THz, caused by wagging vibrations of the amide group and C–H in the pyridine ring along the a-axis (Fig. 6h).

Fig. 7 illustrates the vibrational modes of sorbic acid at nine different frequencies, calculated using CP2K. The experimental peak at 1.669 THz corresponds to the calculated mode at 1.674 THz, caused by lattice stretching along the a-axis (Fig. 7c). The peak at 2.164 THz aligns with the calculated mode at 2.01 THz, primarily due to partial twisting of the –CH₃ group (Fig. 7d). The peak at 2.538 THz matches the calculated mode at 2.268 THz, resulting from collective rocking of SA in the a–c plane (Fig. 7e). The peak at 2.873 THz corresponds to the calculated mode at 2.835 THz, driven by collective stretching of SA along the c-axis (Fig. 7f). The peak at 3.386 THz aligns with the calculated mode at 3.219 THz, attributed to collective rotation of SA in the b–c plane (Fig. 7g). Combined with Fig. 5, theoretical peaks at 0.987 THz, 1.482 THz, 3.444 THz, and 4.17 THz (corresponding to sequentially Fig. 7a, b, h and i) are absent in experimental results, likely due to insufficient resolution of the terahertz spectrometer or noise obscuring weaker absorption peaks. We observed that SA exhibits numerous absorption peaks, some of which are relatively weak, though most can be accurately predicted. Vibration source analysis indicates that these peaks predominantly arise from lattice translations, rotations, and rocking, with minor contributions from intramolecular modes.

Next, we focus on analyzing the assignments of the vibrational modes corresponding to the characteristic peaks at various frequencies in the NIA–SA cocrystal. The experimental spectral feature peak at 1.197 THz corresponds to the calculated spectral feature peak at 1.485 THz, which is caused by the collective rotational vibration of NIA molecules in the b–c plane and the out-of-plane collective wagging vibration of SA molecules (Fig. 8a). The experimental spectral feature peak at 1.719 THz corresponds to the calculated spectral feature peak at 1.89 THz, caused by the out-of-plane collective twisting vibration of NIA and SA molecules (Fig. 8b). The experimental spectral feature peak at 1.928 THz corresponds to the calculated spectral feature peak at 2.361 THz, caused by the out-of-plane collective wagging vibration of NIA molecules and the collective twisting vibration of SA molecules (Fig. 8c). The experimental spectral feature peak at 2.442 THz corresponds to the calculated spectral feature peak at 2.811 THz, caused by the collective rotational vibration of NIA molecules in the b–c plane and the collective stretching vibration of SA molecules (Fig. 8d). These vibrations involve the motion of the hydrogen bond between the O⋯H bond in the carboxyl group of SA and the nitrogen atom in the pyridine ring of NIA, which is a key factor in the formation of the NIA–SA cocrystal. Therefore, this absorption peak can serve as a distinctive signature for the formation of the NIA–SA cocrystal in the 0–4 THz range. The experimental spectral feature peak at 2.824 THz corresponds to the calculated spectral feature peak at 3.198 THz, caused by the collective twisting vibration of NIA molecules and the out-of-plane collective wagging vibration of some SA molecules (Fig. 8e). Based on the experimental results and theoretical analysis, we can conclude that in the NIA–SA cocrystal, the vibrational modes are mainly due to collective vibrations involving all the molecules. The formation of hydrogen bonds in the cocrystal alters the molecular structure and the vibrational modes of certain functional groups, thereby affecting the distribution of the characteristic peaks in the vibrational spectrum.

3.4 Visualization analysis of weak interactions

In the terahertz frequency range, the low-frequency vibrational and rotational modes of molecules are significantly influenced by weak interactions such as van der Waals forces and hydrogen bonds. These weak interactions can affect the potential energy surface of the molecules, thereby causing changes in the position, intensity, and width of absorption peaks in the spectrum, as already shown in Fig. 5. To further analyze the weak interactions in the NIA–SA cocrystal, as well as the mechanisms of hydrogen bonds and van der Waals forces during the cocrystal formation process, this study employed the RDG method and analyzed the energy interactions in the molecular system using Multiwfn and VMD software.

Fig. 9 presents scatter plots and RDG isosurface maps for NIA, SA, and their cocrystal NIA–SA, visually revealing the regions involved in weak interactions within the system and illustrating the relationship between molecular structures and these weak interactions. The types of non-covalent interactions can be determined by the peak positions in the scatter plots, while their strengths are represented by different colors shown on the RDG isosurfaces. In this context, blue regions indicate a larger ρ(r) with sign(λ₂) = −1, signifying stronger attractive weak interactions such as hydrogen bonds and halogen bonds. Green regions exhibit a small ρ(r), with sign(λ₂) possibly being positive or negative, indicating areas of weaker van der Waals forces. Red regions have a larger ρ(r) with sign(λ₂) = +1 indicating areas of stronger steric hindrance.

Fig. 9 The scatter diagram sign λ₂·ρ and RDG isosurface maps for NIA, SA, and NIA–SA. Fig. 9a and a′ are the scatter plots and RDG isosurface maps for NIA. Fig. 9b and b′ are the scatter plots and RDG isosurface maps for SA. Fig. 9c and c′ are the scatter plots and RDG isosurface maps for NIA–SA cocrystal.

The scatter diagram of NIA crystals is shown in Fig. 9a, while Fig. 9a′ presents the isosurface map. In Fig. 9a, multiple peaks can be observed within the range of sign(λ₂)·ρ from −0.05 to 0.05 a.u. Based on the color range of the scatter plot, the isosurfaces (0 < RDG < 0.5, i.e. the area below the black dashed line) are classified into three types: The first type is the blue region located between −0.03 and −0.025 a.u., indicating the presence of strong hydrogen bond interactions in the system. Referring to Fig. 9a′, we can see that these interactions arise from the formation of N–H⋯O bonds between the O⋯H bonds of adjacent amide groups and the nitrogen atom on the pyridine ring, as well as the hydrogen atoms in the adjacent amide groups. By combining the terahertz vibration mode diagrams of NIA molecule shown in Fig. 6, we can find that the weak interaction at the theoretical absorption peak of 3.156 and 3.915 THz mainly originates from hydrogen bond interactions. The second type of isosurfaces consists of green regions located between −0.015 and −0.010 a.u., and between −0.010 and 0.000 a.u., as well as a yellow-green region between 0.00 and 0.010 a.u. The isosurfaces in the range of −0.015 to −0.010 a.u. and −0.010 to 0.000 a.u. exhibit van der Waals weak interactions, primarily derived from π–π stacking interactions between the opposing pyridine rings (offset face-to-face stacking) and the opposing amide groups (–CONH₂). In Fig. 9a′, it can be observed that significant cumulative effects occur within the multiphase system, affecting the physicochemical properties of the material. The yellow-green isosurface in the 0.00 to 0.010 a.u. range indicates the presence of weaker van der Waals interactions, mainly arising from the dispersion forces between the hydrogen atoms on the pyridine rings and the amide groups, as well as dipole–dipole interactions between amide groups. By combining Fig. 6, it can be found that the theoretical absorption peaks of the NIA molecule at 0.843, 1.59, 2.415, 2.619, 3.031, and 3.345 THz mainly originate from van der Waals interactions. The third type of isosurfaces is represented by the red-green region located between 0.015 and 0.020 a.u. and the red region between 0.020 and 0.028 a.u. In the 0.015 to 0.020 a.u. range, both hydrogen bond interactions and steric effects are present. Combining with Fig. 9a′, we can see that the oxygen atom in the amide group forms hydrogen bonds with the hydrogen atoms on the pyridine ring, while the spatial arrangement between the amide group and the pyridine ring may be constrained, leading to steric effects. The red isosurface in the 0.020 to 0.028 a.u. range mainly arises from the tension generated by the pyridine ring, resulting in steric hindrance.

Fig. 9b represents the scatter diagram of SA, while 9b′ shows the isosurface map. Using the same method as above, the isosurfaces (0 < RDG < 0.5, i.e. the area below the black dashed line) are classified into three types. The first type is the blue region located between −0.070 and −0.045 a.u., indicating the presence of strong hydrogen bond interactions in the system. Combining with the isosurface map in 9b′, we can see that this interaction mainly arises from the hydrogen atoms in the hydroxyl groups (–OH) of the carboxylic acid groups acting as hydrogen bond donors, forming hydrogen bonds with the carbonyl groups (C=O) in adjacent molecules. By combining the terahertz vibration mode diagrams of SA molecule in Fig. 7, it can be found that the generation of the theoretical absorption peaks at 1.482 THz and 1.674 THz mainly originates from hydrogen bond interactions. The second type of isosurfaces consists of the green region located between −0.012 and −0.008 a.u. and the yellow-green region between 0.000 and 0.010 a.u., corresponding to van der Waals interactions in the system. The isosurfaces in the −0.012 to −0.008 a.u. range correspond to π–π stacking interactions between olefinic groups and the interactions through dispersion forces among the alkyl chains (–CH₃, –CH₂) and the carbon–hydrogen backbone of the olefin part. The isosurfaces in the 0.000 to 0.010 a.u. range correspond to dipole–dipole interactions between carboxylic acid groups (–COOH), indicating van der Waals interactions. It can be found that the theoretical absorption peaks of SA at 0.987, 2.01, 2.268, 2.835, 3.219, 3.444, and 4.17 THz in Fig. 7 mainly originate from van der Waals interactions. The third type of isosurface is the red-green region located between 0.011 and 0.015 a.u., where both hydrogen bond interactions and steric effects are present. The hydrogen bonds are formed by the interaction between hydrogen atoms on the double bond (–CH=CH–) connecting the carboxylic acid groups within the SA molecule and the oxygen atoms on the carbonyl double bonds in the carboxylic acid groups. The steric effects arise from the tension generated by the double bond within the SA molecule.

Fig. 9c shows the scatter plot of the NIA–SA cocrystal, while 9c′ presents the isosurface map. Using the same method as above, the isosurfaces (0 < RDG < 0.5, i.e. the area below the black dashed line) are classified into three types. The first type is the dark blue region located between −0.100 and −0.070 a.u., indicating the presence of strong hydrogen bond interactions in the system. Combining with the isosurface map in 9c′, we can see that this interaction mainly arises from the hydrogen bond interactions between the carboxylic acid groups of the SA molecules and the pyridine rings in the NIA molecules, forming a stable cocrystal structure. Secondly, the blue region located between −0.045 and −0.035 a.u. also indicates hydrogen bond interactions in the system, primarily arising from the O⋯H bonds on the amide groups of the NIA molecules, which is consistent with the NIA monomer type. Additionally, the light blue region located between −0.03 and −0.02 a.u. indicates relatively weaker hydrogen bond interactions in the system. These interactions mainly arise from the following situations: the amino groups (–NH₂) in the amide groups of NIA molecules acting as hydrogen bond donors, forming weak hydrogen bonds with the carbonyl groups in the adjacent SA carboxylic acid groups; weak hydrogen bonds formed between the C–H bonds on the pyridine rings of NIA molecules and the carbonyl groups on adjacent SA molecules; and weak hydrogen bonds formed between the carbonyl groups on adjacent NIA groups and the C–H bonds on the pyridine rings. By combining the terahertz vibration mode diagrams of the NIA–SA binary cocrystal shown in Fig. 8, the weak interactions corresponding to the theoretical absorption peaks at 2.811 and 3.198 THz primarily arise from hydrogen bonding interactions. The hydrogen bond interaction between the O⋯H bond of the carboxylic acid group in the SA molecule and the nitrogen atom on the pyridine ring in the NIA molecule is the most crucial factor in the formation of the NIA–SA cocrystal. Therefore, the absorption peak at 2.811 THz (corresponding to the experimental absorption peak at 2.442 THz), which represents this interaction, is considered the main characteristic peak of the NIA–SA cocrystal. The second type of isosurfaces consists of the green region located between −0.015 and 0.00 a.u. and the yellow-green region between 0.00 and 0.010 a.u., indicating the presence of van der Waals weak interactions in the system. The green isosurface in the −0.015 to 0.00 a.u. range primarily reflects van der Waals forces arising from π–π stacking interactions between the olefinic groups opposite to SA molecules, and interactions involving the oxygen atoms on the amide groups between adjacent NIA molecules. In the yellow-green region between 0.00 and 0.010 a.u., the van der Waals forces arise from π–π stacking interactions between the pyridine rings opposite to NIA molecules. As shown in Fig. 8, the weak interactions of the theoretical absorption peaks located at 1.485, 1.89, and 2.361 THz, mainly arise from van der Waals forces. The third type of isosurfaces is represented by the red-green region located between 0.010 and 0.015 a.u. and the red region between 0.020 and 0.028 a.u. The red-green regions located between 0.010 and 0.015 a.u. correspond to the O⋯H bonds formed between the hydrogen atoms of the double bonds connecting the carboxyl groups in the SA molecules and the oxygen atoms of the carbonyl double bonds in the carboxyl groups, as well as the interactions between the C=O of the amide groups and the C–H bonds on the pyridine rings. These interactions exhibit both attractive and repulsive forces. The red isosurface in the 0.020 to 0.028 a.u. range corresponds to steric hindrance interactions, primarily arising from the tension generated by the pyridine rings. The RDG analysis allows for a clear observation and comparison of the types and positions of intermolecular interactions in the NIA, SA, and NIA–SA binary cocrystal, as well as their correlation with the vibrational spectra.

4. Conclusions

This study reports the terahertz vibrational spectra of NIA, SA, their physical mixture, and the NIA–SA binary cocrystal. The results reveal the emergence of new absorption peaks in the terahertz spectrum of the NIA–SA cocrystal, which differ from those observed in the individual crystals and the physical mixture, indicating the formation of a new phase distinct from the original single crystals. Vibrational analysis was conducted using THz-TDS experiments in conjunction with ss-DFT calculations. The results demonstrate that the characteristic peaks of the NIA–SA cocrystal in the 0–4 THz range primarily originate from collective intermolecular motions. The hydrogen bonding interaction between the O⋯H bond in the carboxyl group of SA and the nitrogen atom in the pyridine ring of NIA plays a crucial role in the formation of the NIA–SA cocrystal. These hydrogen bonds not only directly drive the formation of the cocrystal but also lead to the absorption peak at 2.811 THz (corresponding to the experimentally observed peak at 2.442 THz) through the collective rotational vibration of NIA molecules in the b–c plane and the collective stretching vibration of SA molecules. This absorption peak serves as a distinctive signature of NIA–SA cocrystal formation within the 0–4 THz range. In summary, the combination of THz-TDS and ss-DFT provides an effective approach for cocrystal identification. The analysis of vibrational modes and weak interactions facilitates the elucidation of the origin of absorption peaks in the spectrum, offering valuable insights into molecular structure exploration and material property understanding.

Data availability

Crystallographic data for Sorbic acid has been deposited at the CCDC under 1206459 and can be obtained from https://doi.org/10.1107/S0108270194002891.

Crystallographic data for Nicotinamide has been deposited at the CCDC under 131757 and can be obtained from https://doi.org/10.1107/S0108768198007848.

Crystallographic data for Nicotinamide-Sorbic acid has been deposited at the CCDC under 2046374 and can be obtained from https://doi.org/10.1039/D3CE00392B.

The code for CP2K can be found at https://www.cp2k.org. The version of the code employed for this study is version 2023.1.

The code for Multiwfn can be found at http://sobereva.com/multiwfn. The version of the code employed for this study is version Multiwfn 3.8(dev).

The code for Mercury can be found at https://www.ccdc.cam.ac.uk/solutions/software/mercury/. The version of the code employed for this study is version Mercury 3.7.

The code for VMD can be found at http://www.ks.uiuc.edu/Research/vmd/. The version of the code employed for this study is version VMD 1.9.3.

The code for gnuplot can be found at http://www.gnuplot.info. The version of the code employed for this study is version gnuplot 6.0.

Author contributions

Di Zhu: conceptualization, funding acquisition, writing – review & editing. Yali Liu: software, writing – original draft. Lei Wang: methodology. Hanming Zhang: investigation. Qiuhong Qu: project administration. Xiaodong Sun: supervision. Yizhu Zhang: funding acquisition.

Conflicts of interest

The authors declare no competing financial interest.

Acknowledgements

This paper is supported by National Natural Science Foundation of China (12174284) and Tianjin Metrology Science and Technology Project (2024TJMT044).

References

1. K. Danchana, P. Jitthiang, K. Uraisin and V. Cerdà, Food Chem., 2021, 361, 130086.
2. M. B. Liewen and E. H. Marth, J. Food Prot., 1985, 48, 364–375.
3. Y. Miwa, T. Mizuno, K. Tsuchida, T. Taga and Y. Iwata, Acta Crystallogr., Sect. B:Struct. Sci., 1999, 55, 78–84.
4. S. Zhang, H. Chen and Åke. C. Rasmuson, CrystEngComm, 2015, 17, 4125–4135.
5. H. Zhang, Y. Zhu, N. Qiao, Y. Chen and L. Gao, Pharmaceutics, 2017, 9, 54.
6. J. E. K. Schawe, S. Pogatscher and J. F. Löffler, Thermochim. Acta, 2020, 685, 178518.
7. C. Li, D. Wu, Z. Gao and W. Chen, CrystEngComm, 2023, 25, 4889–4901.
8. M. Koch, D. M. Mittleman, J. Ornik and E. Castro-Camus, Nat. Rev. Methods Primers, 2023, 3, 48.
9. E. Fardelli, A. D'Arco, S. Lupi, D. Billi, R. Moeller and M. C. Guidi, Spectrochim. Acta, Part A, 2023, 288, 122148.
10. Q. Zhao, J. Zhang, Y. Jing, J. Xue, J. Liu, J. Qin, Z. Hong and Y. Du, Chem. Phys., 2025, 591, 112584.
11. J. Zhang, M. Wan, J. Fang, Z. Hong, J. Liu, J. Qin, J. Xue and Y. Du, Spectrochim. Acta, Part A, 2023, 295, 122623.
12. S. Koranne, A. Sahoo, J. F. Krzyzaniak, S. Luthra, K. K. Arora and R. Suryanarayanan, Mol. Pharmaceutics, 2018, 15, 3297–3307.
13. H. A. Hafez, X. Chai, A. Ibrahim, S. Mondal, D. Férachou, X. Ropagnol and T. Ozaki, J. Opt., 2016, 18, 093004.
14. P. Wang, J. Zhao, Y. Zhang, Z. Zhu, L. Liu, H. Zhao, X. Yang, X. Yang, X. Sun and M. He, Int. J. Pharm., 2022, 620, 121759.
15. Y. Zhao, Z. Li, J. Liu, T. Chen, H. Zhang, B. Qin and Y. Wu, Spectrochim. Acta, Part A, 2018, 192, 336–342.
16. L. Xu, Y. Li, P. Jing, G. Xu, Q. Zhou, Y. Cai and X. Deng, Spectrochim. Acta, Part A, 2021, 249, 119309.
17. J. Zhang, Z. Zhu, Y. Wu, T. Ji, J. Wang, H. Zhu, W. Peng, M. Chen, S. Li and H. Zhao, J. Infrared, Millimeter, Terahertz Waves, 2023, 44, 814–829.
18. M. Wan, J. Fang, J. Xue, J. Liu, J. Qin, Z. Hong, J. Li and Y. Du, Int. J. Mol. Sci., 2022, 23, 8550.
19. J. Zhang, Y. Jing, M. Wan, J. Xue, J. Liu, J. Li and Y. Du, Spectrochim. Acta, Part A, 2024, 305, 123478.
20. Y. Jing, J. Zhang, M. Wan, J. Xue, J. Liu, J. Qin, Z. Hong and Y. Du, IEEE Trans. Terahertz Sci. Technol., 2024, 14, 152–161.
21. B. Peng, J. Shu, Z. Hou, S. Qian, B. Yan, B. Zhang, Y. Zhao, B. Su and C. Zhang, Int. J. Pharm., 2024, 651, 123767.
22. P. Wang, Y. Li, W. Han, Y. Yan, C. Zhang, Q. Qu, X. Zhang, L. Liu, X. Sun, X. Yang and M. He, Talanta, 2024, 278, 126489.
23. C. F. Macrae, I. Sovago, S. J. Cottrell, P. T. A. Galek, P. McCabe, E. Pidcock, M. Platings, G. P. Shields, J. S. Stevens, M. Towler and P. A. Wood, J. Appl. Crystallogr., 2020, 53, 226–235.
24. C. F. Macrae, P. R. Edgington, P. McCabe, E. Pidcock, G. P. Shields, R. Taylor, M. Towler and J. Van De Streek, J. Appl. Crystallogr., 2006, 39, 453–457.
25. P. J. Cox, Acta Crystallogr., Sect. C:Cryst. Struct. Commun., 1994, 50, 1620–1622.
26. T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöβ, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack and J. Hutter, J. Chem. Phys., 2020, 152, 194103.
27. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
28. R. Krishnan, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 1980, 72, 650–654.
29. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
30. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
31. T. Lu and Q. Chen, in Comprehensive Computational Chemistry, Elsevier, 2024, pp. 240–264.
32. E. R. Johnson, S. Keinan, P. Mori-Sánchez, J. Contreras-García, A. J. Cohen and W. Yang, J. Am. Chem. Soc., 2010, 132, 6498–6506.
33. P. Hohenberg and W. Kohn, Phys. Rev. A, 1964, 136, B864–B871.
34. A. J. Cohen, P. Mori-Sánchez and W. Yang, Science, 2008, 321, 792–794.
35. J. P. Perdew and W. Yue, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 8800–8802.
36. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
37. T. Lu, J. Chem. Phys., 2024, 161, 082503.
38. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
39. T. Y. Li, C. Li, L. Zhang, Q. Xiao and L. Jiang, Guangpuxue Yu Guangpu Fenxi, 2021, 41, 100–104.
40. Q. Zhou, Y. Shen, Y. Li, L. Xu, Y. Cai and X. Deng, Spectrochim. Acta, Part A, 2020, 236, 118346.

---