Cooperative electronic effects of Na⁺ and Ca²⁺ on an oxygenated aromatic model

Abstract

This study uses molecular modeling to investigate the individual and combined effects of sodium (Na⁺) and calcium (Ca²⁺) ions on a representative molecular model of Yallourn lignite. Lignite, a geologically young and chemically reactive coal rich in oxygen-containing functional groups, can be chemically modified through ion exchange with metal cations. Experimental studies have suggested that Na⁺ and Ca²⁺ ions alter the thermal behavior and gasification reactivity of lignite, particularly when co-added. To provide molecular-level insights into these effects, we performed density functional theory calculations at the B3LYP/6-31G* level, including full geometry optimization, electrostatic potential distributions, and Löwdin bond-order analysis. The results reveal that Na⁺ and Ca²⁺ exert distinct effects on molecular geometry and electron localization, while their co-addition induces cooperative stabilization and polarization. These findings deepen our understanding of ion-induced structural and electronic modifications and present a reproducible computational framework for designing catalytic upgrading strategies for low-grade carbonaceous materials.

---

Introduction

Alongside the expansion of renewable energy, improving the efficiency and environmental performance of existing carbon-based resources remains important. In this context, low-rank coals are widely distributed and attractive feedstocks, but their high moisture content and low energy density limit direct utilization. Lignite, a geologically young form of coal rich in oxygen-containing functional groups, can be chemically modified through ion exchange with metal cations. However, practical issues such as spontaneous heating and ignition during drying, pulverization, and transportation necessitate pretreatment and upgrading.

To enhance lignite utilization, various methods, such as drying, thermal upgrading, solvent extraction, catalytic treatment, and gasification, have been explored globally. Among these, chemical treatments that leverage selective interactions between functional groups (e.g., carboxyl (–COOH) and phenolic hydroxyl (–OH) groups) and metal ions have drawn attention for their potential to stabilize the structure of lignite and control its reactivity. A recent study by Wang et al.¹ demonstrated that cations such as Na⁺ and Ca²⁺ can interact with functional groups in lignite, improving thermal decomposition and promoting gasification reactions. Their results showed that treating Shengli lignite with Na₂CO₃ or NaOH generates –COONa structures, markedly enhancing ignitability and combustion performance. Similarly, experimental and computational studies by Zhang et al.² and Yin et al.³ revealed Na⁺- and Ca²⁺-induced structural changes in lignite. These findings highlight the need for theoretical investigations into ion-induced electronic effects.

Tsubouchi et al.⁴ demonstrated that adding nanoscale Ca/Na composite catalysts to lignite markedly promotes gasification compared with single-ion addition, providing experimental evidence of a synergistic effect. These findings suggest that similar mechanisms may apply to biomass-derived materials at intermediate stages of humification—between plant-based organics and coal—offering opportunities for early-stage resource conversion through metal-ion interactions.⁴ For example, Li et al.⁵ reported that chemical and thermal activation of coal gangue with CaO and CaSO₄ at 800 °C led to structural amorphization and improved pozzolanic reactivity. The catalytic effects of alkali and alkaline earth metals on biomass and coal pyrolysis have also been comprehensively reviewed.⁶ In addition, Yin et al.³ investigated how heating alters aromatic condensation and functional group distributions during the pyrolysis of lignite and bituminous coal. They showed that inert minerals, such as SiO₂ or potassium feldspar, can transform into reactive amorphous phases, implying potential relevance of lignite activation mechanisms.³ Thus, the metal co-addition effects elucidated in this study may offer promising strategies for the early utilization of plant-derived carbon resources, thereby contributing to a more sustainable carbon cycle.

Although previous studies have mainly focused on the effects of single-metal cations (e.g., Na⁺ or Ca²⁺), the combined or competing effects of simultaneous cation addition remain poorly understood. Given the differences in their valence, ionic radius, and electron affinity, Na⁺ and Ca²⁺ are expected to interact differently with carbon backbones. Specifically, Na⁺, being more mobile, tends to engage in π interactions with aromatic rings. In contrast, owing to its divalency, Ca²⁺ may form bridge-like structures through strong electrostatic interactions with oxygen atoms. Therefore, understanding the combined influence of these cations on charge distribution, molecular polarity, bond localization, and reactivity through theoretical analysis holds scientific and engineering significance. These theoretical approaches are increasingly being used to predict, design, and evaluate chemical modification processes for carbon-based materials.⁷ Shinohara and Tsubouchi^8,9 performed detailed density functional theory (DFT)-based studies demonstrating how ion exchange with Ca²⁺ or Na⁺ alters the electronic structure and reactivity of low-rank coals. Shi et al.¹⁰ constructed a molecular model of Xianfeng lignite based on experimental data and evaluated its structural features, supporting the validity of the Yallourn lignite model and the theoretical approach employed in this study. Bai et al.¹¹ combined experimental and DFT-based studies to investigate the effects of Na⁺ and Ca²⁺ on carbon structure and gasification reactivity, while Li et al.⁵ focused on experimental approaches to evaluate the reactivity enhancement by Ca-based additives. Bai et al.¹¹ used carbon black in DFT and experimental investigations to reveal the distinct catalytic behaviors of Na⁺ and Ca²⁺ ions during coal gasification, while He et al.⁷ demonstrated how Na⁺ modification influences bonding and electronic states. Additionally, Tsubouchi et al.⁴ reported that steam gasification with nanoscale Ca²⁺/Na⁺ composite catalysts achieved higher reactivity than individual cation addition, suggesting site-specific and complementary electronic and structural effects of the two ions in lignite. For instance, during steam gasification at 700 °C, the char conversion of raw Adaro coal was only approximately 18%, whereas Na⁺- and Ca²⁺-exchanged coals reached the values of ∼75% and ∼82%, respectively. Moreover, Ca²⁺/Na⁺ co-exchanged coal achieved nearly complete conversion (Table 1).⁴

Table 1 Char conversion of ion-exchanged Adaro coal samples during steam gasification at 700 °C (dry-ash-free, daf basis). Data adapted from the literature⁴

Sample code	Char conversion, mass% (daf^a)
a daf = dry ash-free basis. b AD = Adaro coal.
Raw AD^b coal	18
Na^+-AD coal	75
Ca^2+-AD coal	82
Ca^2+/Na^+-AD coal	100

---

To provide theoretical support for these experimental findings, this study employs quantum chemical calculations to analyze the structural stabilization and electronic modifications in a Yallourn lignite model following single- and dual-cation addition (Na⁺ and Ca²⁺). Recent advances in quantum chemical calculations allow direct visualization of molecular structural and electronic changes. In particular, DFT and molecular orbital methods serve as powerful tools to pre-evaluate carbon material reactivity and guide process design. Indicators such as optimized geometry, electrostatic potential (ESP) maps, and Löwdin bond orders (BOs) provide valuable insights into the reactivity of metal-modified lignite models. In particular, Ugwumadu et al.¹² employed molecular dynamics simulations with the Simulation of Thermal Emission of Atoms and Molecules method to elucidate the structural evolution of carbonized and graphitized materials under the influences of heteroatoms (O, N, S), an approach resonant with the DFT-based analysis in the present study.

Accordingly, this study focuses on a Yallourn lignite-based molecular model to investigate how the single and co-addition of Na⁺ and Ca²⁺ ions affect molecular structure, electronic distribution, and bonding characteristics. Specifically, we performed geometry optimizations using Becke's three-parameter Lee–Yang–Parr hybrid functional (B3LYP)/6-31G* level of theory, quantified local relaxation based on Löwdin BO analysis, and visualized charge redistribution through ESP mapping. These methods provide a theoretical framework for understanding the cooperative electronic effects and enhanced reactivity resulting from Na⁺/Ca²⁺ co-addition. The findings of this study are expected to contribute to advanced lignite utilization in gasification and pyrolysis processes by not only validating the effectiveness of additive design but also informing functional site targeting and electronic complementarity at the molecular level. From a theoretical chemistry perspective, this study offers new insights into the efficient use of carbon resources in addressing societal energy needs.

Computational details

Molecular model construction

To investigate the electronic effects of cation exchange, we constructed a model of Yallourn low-rank coal (H) and compared it to Na⁺-exchanged (Na), Ca²⁺-exchanged (Ca), and Na⁺/Ca²⁺ co-exchanged structures (Na/Ca-1 and Na/Ca-2, corresponding to two distinct cation arrangements).

The molecular model employed in this study was derived from lignite obtained from the Yallourn coal mine, Victoria, Australia. This low-rank coal is characterized by its relatively high-oxygen content, low coalification, and abundance of polar functional groups, such as –COOH and –OH. These functional groups strongly influence the chemical reactivity and metal ion affinity of the carbon framework, making them central to the design of the molecular model. Based on prior experimental findings, the present theoretical analysis focuses on Yallourn lignite because of its high reactivity and well-characterized functional group composition.

Experimental investigations have demonstrated that ion exchange with Na⁺ and Ca²⁺ enhances the gasification reactivity of low-rank coals, such as Yallourn, especially when both cations are introduced simultaneously.^4,8,9 These effects are attributed to interactions between the metal ions and oxygen-containing functional groups that play a crucial role in determining catalytic behavior and form the rationale for the present theoretical analysis. This interpretation is supported by classical studies that demonstrate a strong correlation between CaO dispersion and char gasification reactivity in lignite,¹³ highlighting the key role of calcium dispersion in catalytic enhancement.

In recent years, structural characterization of lignite has increasingly utilized representative structural units derived from molecular dynamics and quantum chemical modeling. These computational approaches are often combined with experimental techniques, such as nuclear magnetic resonance, Fourier transform infrared spectroscopy (FT-IR), and X-ray photoelectron spectroscopy (XPS). The structural model of lignite employed in this study is consistent with previous semiquantitative models constructed using FT-IR and XPS analyses,^14–16 in conjunction with simulation approaches such as those reported by Zhang et al.¹⁷ Additionally, this model builds on earlier computational studies on coal structure evolution during carbonization and graphitization.¹⁸ Following this approach, a single molecule was constructed with an aromatic ring as the core framework and carboxyl and phenolic –OH groups as functional side chains (Fig. 1). This design was inspired by the molecular orbital calculations of Isoda et al.,¹⁹ which provide a theoretical basis for evaluating reactivity features such as π-electron delocalization, decarboxylation of carboxylic acids, and the hydrogen bonding potential of hydroxyl groups. Moreover, the selected molecular model maintains near-electrical neutrality while allowing localized structural and electronic changes upon metal ion incorporation. In particular, the –COOH group serves as the primary site for electrostatic interactions with Na⁺ and Ca²⁺ ions, whereas the phenolic –OH group, which can conjugate with the aromatic ring, facilitates selective coordination behavior.

Fig. 1 Optimized molecular structures of the Yallourn low-rank coal model and its Na⁺/Ca²⁺ co-exchanged derivatives, computed at the Becke's three-parameter Lee–Yang–Parr hybrid functional (B3LYP)/6-31G* level of theory. (a) Yallourn unit structure; (b) unmodified model (H); (c) Na⁺/Ca²⁺ co-exchanged model with Ca²⁺ bridging carboxylate and phenolate sites (Na/Ca-1), and (d) Na⁺/Ca²⁺ co-exchanged model with Ca²⁺ chelating the carboxylate group (Na/Ca-2). Dashed lines indicate ionic coordination between Ca²⁺ and deprotonated functional groups (–COO⁻ and –O⁻), reflecting distinct cation-binding modes in co-exchanged structures.

The model was intentionally constructed without imposing symmetry constraints, enabling three-dimensional relaxation during geometry optimization. This approach provides a more realistic representation of post-optimization conformational changes and intramolecular interactions among functional groups. This molecular design not only mimics the original lignite structure but also enables the quantitative evaluation of phenomena such as structural relaxation, electronic redistribution, and reactivity variations induced by metal ion addition.

In summary, the molecular model adopted in this study captures the key chemical characteristics of Yallourn lignite and includes structural features well-suited for quantum chemical analysis. Therefore, it provides a valid foundation for subsequent geometry optimization, electronic state analysis, and bonding evaluation.

Software program and computational conditions

All quantum chemical calculations in this study were performed using Spartan’18 for Windows (Build 1.3.0; Wavefunction, Inc., Irvine, CA, USA).²⁰ Spartan is a commercially available computational chemistry software program that provides an integrated platform for molecular model construction, geometry optimization, energy evaluation, and visualization of electron distribution. It is extensively used in educational and research contexts, particularly for its balance of computational accuracy and flexibility, making it well-suited for modeling complex, heterogeneous systems such as lignite.

Before quantum chemical calculations, initial molecular structures were generated using molecular mechanics based on the Merck molecular force field (MMFF).²¹ This preliminary step aimed to identify low-energy conformers and generate a physically reasonable starting geometry for subsequent calculations. MMFF is a well-established classical force field optimized for organic molecules and has proven effective in predicting realistic three-dimensional molecular structures. A precedent for this approach in coal research can be found in the work of Carlson,²² who applied molecular dynamics simulations to a bituminous coal model, demonstrating the stabilizing roles of van der Waals interactions and hydrogen bonding.

Subsequently, the geometries were refined through quantum chemical geometry optimizations. Two principal computational methods were employed: The first is the restricted Hartree–Fock (HF) method,²³ which provides a fundamental ab initio framework for many-electron systems via a mean-field approximation. The second is the B3LYP method, an extensively adopted hybrid DFT approach that uses the B3 and LYP correlation functionals.^24,25 This method is recognized for its balance between accuracy and computational efficiency, making it particularly suitable for chemically diverse and structurally complex systems such as lignite.

Natural orbital and density matrix-based population analyses were employed to evaluate electronic redistribution and local variations in electron density quantitatively, caused by metal ion incorporation. Based on Löwdin's quantum theory of many-particle systems,²⁶ this approach utilizes symmetric orthogonalization to convert non-orthogonal basis functions into an orthogonal set, ensuring more stable and transferable BO values across different optimized geometries.

When selecting the basis set, we considered not only the balance between computational cost and accuracy but also the reproducibility of the infrared (IR) spectra. The lightweight 3-21G* basis set was used for the initial exploration of molecular structures. For selected candidates, IR spectra were computed using the 6-311G* and 6-31G* basis sets. Comparison with the experimental spectra reported by Wang et al.¹ revealed that the 6-31G* basis set provided better agreement in both vibrational peak positions and intensities. Accordingly, we selected the B3LYP/6-31G* level of theory for the main analyses in our study. This choice ensures a consistent and computationally tractable protocol for comparing relative structural and spectral trends across all models. The B3LYP functional with polarized basis sets has been widely validated for reproducing vibrational properties and metal–oxygen interactions in organic and coordination systems.^{24,25,27–29}

For each molecular configuration—unexchanged, Na⁺-exchanged, Ca²⁺-exchanged, and Na⁺/Ca²⁺ co-exchanged—geometry optimizations and vibrational calculations were performed using the B3LYP/6-31G* level of theory, iteratively refined until the local energy minima were reached. Stephens et al.³⁰ demonstrated that B3LYP combines the B3 exchange functional—which incorporates exact HF exchange and empirical parameters—with the LYP correlation functional, and is well known for its high predictive accuracy.^24,25

To evaluate the electronic characteristics of each structure, ESP maps and Löwdin BOs were calculated. ESP maps visually represent charge localization by mapping charge distributions onto isosurfaces of electron density, making them particularly useful for identifying changes in polarity and charge concentration induced by metal ion incorporation. Conversely, Löwdin BOs offer a quantitative measure of electronic bonding strength in each bond, yielding useful indicators of structural flexibility and electronic stabilization effects. Together, these analyses help visualize and quantify the electronic effects and bonding changes caused by metal ion addition, thus providing molecular-level insight into the electronic effects of metal-ion addition in the present model. We note that absolute charge and bond-order values can depend on the chosen population/partitioning scheme; therefore, our discussion focuses on consistent, relative trends obtained using the same settings across all structures.

Geometry optimizations were considered converged when the energy change was <1.0 × 10⁻⁶ Hartree, the maximum gradient was <3.0 × 10⁻⁴ Hartree Bohr⁻¹, and atomic displacements were <1.2 × 10⁻³ Bohr. All calculations were performed under gas-phase conditions, treating the systems as isolated molecules without explicit solvent effects. As a solvent-sensitivity check, representative single-point PCM (water) calculations were also performed (SI, Table S10). Harmonic vibrational frequencies and IR intensities were calculated for all the optimized structures. For comparison with the experimental spectra, a uniform scaling factor (0.95×) was applied during postprocessing; unless otherwise noted, all discussed peak positions correspond to the scaled frequencies. In response to the reviewers, supplementary gas-phase calculations at the M06-2X/def2-TZVP level were performed for representative structures (H, Na, Ca, and Na/Ca-1; labeled M, M–Na, M–Ca, and M–NaCa in the SI) as a sensitivity check for the IR assignments. The corresponding unscaled harmonic frequencies and IR intensities, together with values scaled by 0.95, are provided in the SI (Tables S2–S5), along with representative mode assignments confirmed by normal-mode inspection (Tables S6–S9).

These gas-phase calculations reflect standard computational practices in lignite modeling but may introduce limitations in simulating condensed-phase behavior.²⁰

Results and discussion

Optimized structures and total energies

In this study, we investigated the effects of Na⁺/Ca²⁺ co-addition on the structure and energetic stability of a simplified lignite model inspired by Yallourn coal. Geometry optimizations were performed at the B3LYP/6-31G* level, with particular focus on the structural diversity and energetic implications of the Na⁺/Ca²⁺ co-exchanged configurations. As shown in Fig. 1(a–d), two configurations were considered for the Na⁺/Ca²⁺ co-exchanged structure: Na/Ca-1 (Fig. 1c) and Na/Ca-2 (Fig. 1d). In the Na/Ca-1 structure (Fig. 1c), one of the –COOH groups undergoes deprotonation to form a carboxylate anion (–COO⁻), which then coordinates with Ca²⁺ to form a –COOCa complex, resulting in an overall charge-neutral coordination environment. Simultaneously, a phenolic –OH group is deprotonated to form a phenolate anion (–O⁻), allowing bridging coordination between the phenolate and carboxylate oxygens via Ca²⁺ (Ph–O⁻–Ca²⁺–OOC). Selected metal–oxygen distances are summarized in the SI (Table S1). In this structure, Na⁺ and Ca²⁺ preferentially bind to carboxyl and phenolic hydroxyl groups, respectively, resulting in spatially separated electrostatic and electronic influences. This arrangement enhances intramolecular structural stability and reduces electron density polarization, providing molecular-level insight into factors governing framework stability.

In contrast, in the Na/Ca-2 structure (Fig. 1d), Ca²⁺ simultaneously coordinates with the phenolate oxygen (from deprotonated Ph–OH) and oxygen atom of the carboxylate group, forming a stable bidentate bridging structure (Ph–O⁻⋯Ca²⁺⋯–OOC). Meanwhile, Na⁺ is electrostatically associated with a different phenolate oxygen, resulting in a Ph–ONa structure. In this configuration, Na⁺ and Ca²⁺ are in close proximity, potentially causing overlapping local electrostatic fields and competition for coordination sites. This may induce local structural strain within the molecule. Compared with the structure of Na/Ca-1 (Fig. 1c), the Na/Ca-2 structure (Fig. 1d) exhibits stronger interactions between the metal ions.

A comparison of the total energies of the two configurations (Fig. 1c and d) revealed that the Na/Ca-1 configuration (Fig. 1c) is more stable by approximately 110 kJ mol⁻¹ (Table 2). This substantial energy difference cannot be attributed solely to the presence of metal ions. Instead, it reflects the optimized binding positions and spatial arrangement, which allow each ion to exert its specific electronic influence. Specifically, the bridging effect of divalent Ca²⁺ together with the local relaxation associated with monovalent Na⁺ is consistent with a cooperative relief of structural stress throughout the molecule. These results indicate that the stabilizing effect of the Na⁺ and Ca²⁺ co-addition arises from the presence of these two ions and the optimization of their binding sites and spatial distribution. This suggests that, at the molecular level, the positional arrangement and functional group selectivity of the ions play a more critical role than the total quantity of added ions.

Table 2 Computed total energies of Na⁺/Ca²⁺ co-exchanged structures optimized at the Becke's three-parameter Lee–Yang–Parr hybrid functional (B3LYP)/6-31G* level of theory. The Na/Ca-1 configuration is thermodynamically more stable than that of Na/Ca-2, with an energy difference of 110.2 kJ mol⁻¹ (0.04196 Hartree^b)

Structure	Relative energy (Hartree)	Relative energy (kJ mol^−1)
a Energy of Na/Ca-1 is taken as the reference (0.0 kJ mol^−1). b Conversion factor used is 1 Hartree = 2625.5 kJ mol^−1.
Na/Ca-1^a	−2178.330695	0.0
Na/Ca-2	−2178.288731	110.2

---

Based on these findings, we selected the Na/Ca-1 structure (Fig. 1c) as the primary focus for further structural and electronic analyses.

Comparison of IR spectra

To evaluate structural changes in the lignite model, IR spectra were calculated for the optimized structures at the HF/6-311G* and B3LYP/6-31G* levels. Unless otherwise noted, the discussion below refers to scaled harmonic frequencies (0.95×), while Fig. 2 and 3 show the corresponding unscaled spectra. Full unscaled/scaled frequency–intensity lists and representative assignments are provided in the SI (Tables S2–S5 and S6–S9).

Fig. 2 Calculated infrared (IR) spectra of model compounds related to Yallourn low-rank coal, computed at the B3LYP/6-31G* level of theory: (a) reference molecule (1-hydroxy-2-naphthoic acid), (b) unmodified coal model (H), (c) Ca²⁺-exchanged model (Ca), and (d) Na⁺/Ca²⁺ co-exchanged model (Na/Ca-1). Spectra are plotted as unscaled harmonic frequencies and IR intensities. Peak labels in the figure indicate the unscaled harmonic frequencies. Peak positions used for comparison with experiment are based on the scaled frequencies (0.95×) and are summarized in Table 3; full scaled frequency/intensity lists and representative assignments confirmed by normal-mode inspection are provided in the SI (full scaled frequency/intensity lists: Tables S2–S5; representative assignments confirmed by normal-mode inspection: Tables S6–S9).

Fig. 3 Calculated IR spectra of model compounds related to Yallourn low-rank coal, computed at the RHF/6-311G* level of theory: (a) reference molecule (1-hydroxy-2-naphthoic acid), (b) unmodified coal model (H), (c) Ca²⁺-exchanged model (Ca), and (d) Na⁺/Ca²⁺ co-exchanged model (Na/Ca-1). Spectra are plotted as unscaled harmonic frequencies and IR intensities, and peak labels in the figure indicate the unscaled harmonic frequencies. Peak positions used for comparison with experiment are based on the scaled frequencies (0.95×) and are summarized in Table 3. Full scaled frequency/intensity lists and representative assignments confirmed by normal-mode inspection are provided in the SI (full scaled frequency/intensity lists: Tables S2–S5; representative assignments confirmed by normal-mode inspection: Tables S6–S9).

Wang et al.¹ reported a strong IR band at 1705 cm⁻¹ and discussed carboxylate-related features in the 1610–1550 and ∼1400 cm⁻¹ regions (often summarized as bands near 1570 and 1400 cm⁻¹ for Na⁺-exchanged lignite). Because the experimental spectra in Wang et al. were obtained for Na⁺-exchanged lignite without Ca²⁺, we used these bands only as qualitative reference points for the carboxylate-related region instead of assuming a strict peak-to-peak correspondence.

In the present study, B3LYP calculations predicted a peak at 1423 cm⁻¹ (scaled), which was assigned to the stretching vibration of the C=C–O bond in the carboxyl group coordinated with Na⁺ or Ca²⁺. Another band was observed at ∼1264 cm⁻¹ (scaled), corresponding to C–O bond stretching vibrations in the coordinated carboxyl group. Although the predicted peak positions deviate slightly from the experimental values, the overall results show reasonable agreement (Fig. 2 and Table 3).

Table 3 Comparison of calculated infrared absorption peaks with experimental values reported by Wang et al.¹ Vibrational modes were assigned based on molecular structures and characteristic wavenumber regions

Method^d	Structure	Peak position (scaled by 0.95 ×, cm^−1)	Vibrational mode^b	Experimental peak^c (Wang et al., cm^−1)	Agreement^a
Notes:a Agreement is defined as follows: ○ = reasonable match (within ±50 cm^−1), △ = tentative match (within ±100 cm^−1 or trend-based), × = poor or no match. Peak positions listed are scaled (0.95×); agreement is evaluated using the scaled values. When two experimental bands are listed (1570 and 1400 cm^−1), agreement is assessed based on the closer band.b Symbol ν denotes stretching vibration.c Vibrational mode assignments are based on theoretical predictions, as experimental assignments were not explicitly provided by Wang et al.^1d Frequencies calculated using Hartree–Fock methods are generally overestimated, owing to the absence of electron correlation.
HF/6-311G*	H	1910	νC=O	1705	×
		1464	νC–O	—	—
	Ca	1643	νC=O	1570, 1400 (Na)	△
	Na/Ca-1	1632	νC=O	1570, 1400 (Na)	△
B3LYP/6-31G*	H	1664	νC=O	1705	○
		1264	νC–O	—	—
	Ca	1435	νC=O	1570, 1400 (Na)	○
	Na/Ca-1	1423	νC=O	1570, 1400 (Na)	○

---

In contrast, the HF method predicted a major absorption peak at approximately 2010 cm⁻¹ (unscaled; 1910 cm⁻¹ scaled), which is substantially blue-shifted with respect to the experimental peak (1705 cm⁻¹). This discrepancy is attributed to the neglect of electron correlation effects in the HF method, which leads to overestimated vibrational frequencies. This blue shift may be further amplified by the gas-phase approximation used in the calculations, which excludes condensed-phase effects such as intermolecular hydrogen bonding and solvation. In addition, the simplified molecular model employed may not fully capture the complex structural and electronic environments present in real lignite samples, contributing to deviations from experimental spectra. Moreover, for the Ca²⁺ and Na +/Ca²⁺ structures, the HF spectra lacked distinct signals corresponding to either C=O or C–O vibrations, rendering mode assignments unreliable (Fig. 3).

Moreover, the B3LYP spectra exhibited an absorption band at 1664 cm⁻¹ (scaled; 1752 cm⁻¹ unscaled), attributable to the C=O stretching vibration of unreacted carboxylic acid (–COOH) groups. In contrast, the HF spectra yielded a corresponding peak at 1632 cm⁻¹ (scaled; 1718 cm⁻¹ unscaled), which is likely overestimated due to the inherent neglect of electron correlation effects in the HF method. Notably, these –COOH-related peaks were absent in the experimental spectra reported by Wang et al.,¹ suggesting that these functional groups were either not present or below the detection threshold under their experimental conditions.

In summary, the B3LYP/6-31G* method provided IR spectral predictions that were in better agreement with experimental observations, especially in reproducing vibrational features associated with carboxyl and hydroxyl functional groups. In comparison, the HF method resulted in large deviations and limited assignability due to its oversimplified treatment of electron correlation. Therefore, the B3LYP method was selected for further analyses of the electronic structure and bonding characteristics. The scaled IR frequencies (0.95×) and intensities are provided in the SI (Tables S2–S5), and the key band assignments discussed here are summarized in the SI (Tables S6–S9).

Visualization of electronic properties via ESP mapping

ESP maps were generated to elucidate the electron density distribution and local electronic characteristics of the optimized lignite model structures. ESP is a useful indicator for visualizing electron-rich and electron-deficient regions on the molecular surface, providing insights into electron localization/delocalization, molecular polarity, and potential reactive sites. Calculations were performed at the B3LYP/6-31G* level, and the ESP distributions (within ±200 kJ mol⁻¹) were mapped onto an isodensity surface of 0.002 e Bohr⁻³ using the Spartan‘18 software program (Fig. 4 and 5).²⁰

Fig. 4 Electrostatic potential (ESP) maps of Yallourn low-rank coal model structures calculated at the B3LYP/6-31G* level of theory. (a) Pristine model (H), exhibiting a uniform ESP distribution across the molecule; (b) Ca²⁺-exchanged model (Ca), displaying localized positive potential near the Ca²⁺ coordination site, accompanied by contraction of negative regions around the aromatic framework, and (c) Na⁺/Ca²⁺ co-exchanged model (Na/Ca-1), exhibiting enhanced charge polarization and more distinct negative potentials around carboxylate and phenolic sites, indicating cooperative electronic perturbation by both cations. These results suggest that Na⁺/Ca²⁺ co-addition induces the strongest electronic redistribution, reflecting intensified cation–framework interactions and potential modulation of reactivity. Atom colors: H (light gray), C (gray), O (dark red), Ca (red), Na (blue). Calcium (Ca²⁺) is distinguishable from oxygen by its larger sphere size and labeled position. The color scale represents electrostatic potential values ranging from −200 kJ mol⁻¹ (deep red) to +200 kJ mol⁻¹ (deep blue). Thus, red regions indicate negative electrostatic potential (electron-rich areas), while blue regions represent positive electrostatic potential (electron-deficient areas).

Fig. 5 ESP maps of Yallourn low-rank coal model structures calculated at the B3LYP/6-31G* level of theory. (a) Pristine structure (H), displaying a uniform electrostatic potential distribution; (b) Na⁺-exchanged structure (3Na), exhibiting intensified negative potentials near the carboxylate and phenolic groups, indicating localized charge accumulation induced by Na⁺ coordination, and (c) Na⁺/Ca²⁺ co-exchanged structure (Na/Ca-1), showing enhanced charge polarization and distinct spatial separation of positive and negative regions, reflecting synergistic effects of the two cations. These findings highlight the critical importance of cation identity and spatial arrangement in modulating the electronic environment of lignite. Atom colors: H (light gray), C (gray), O (dark red), Ca (red), and Na (blue). Ca²⁺ can be distinguished from O by its larger atomic size and labeled position. The color scale represents ESP values from −200 kJ mol⁻¹ (deep red) to +200 kJ mol⁻¹ (deep blue). Thus, red regions denote negative electrostatic potential (electron-rich areas), while blue regions indicate positive electrostatic potential (electron-deficient areas).

In the untreated Yallourn lignite structure, a uniform ESP was observed across the aromatic plane, with moderately negative regions around the carboxyl groups. This indicates an overall balanced charge distribution, without marked polarization across the molecule (Fig. 4a and 5a). In the Ca²⁺-exchanged structure, Ca²⁺ bridged two oxygen atoms, creating a concentrated region of positive potential near the ion (Fig. 4b). This led to a relative contraction of the negative potential region across the aromatic ring, reflecting a more localized electron density distribution. The strong electrostatic influence of the divalent Ca²⁺ ion induces pronounced local polarization of the electron cloud. In the Na⁺-exchanged structure, Na⁺ coordinated near the phenolic hydroxyl and carboxyl groups, producing a localized region of positive potential directly above its binding site. Simultaneously, the region around the aromatic ring distant from Na⁺ showed extended negative potential, indicating an intramolecular charge redistribution (Fig. 5b). In this work, the Na⁺–π interaction is inferred from the optimized geometries and associated electronic descriptors (e.g., ESP distributions and bond-order trends), rather than from isolated Na⁺–π complex models.³¹ As a result, a dipolar potential gradient can be observed across the molecule.

The most pronounced changes were observed in the Na⁺/Ca²⁺ co-exchanged structure (Fig. 4c and 5c). In this configuration, Na⁺ was positioned near the aromatic plane, while Ca²⁺ coordinated near the carboxyl group. Consequently, the positive charge distribution was spatially separated within the molecule. A strong negative potential remained on the aromatic plane, while a high positive potential emerged near the carboxyl region. This evident spatial separation of positive and negative regions in the ESP map suggests the development of a pronounced intramolecular electrostatic dipole.

This dipole formation arises from the site-specific binding of Na⁺ and Ca²⁺, with Na⁺ near the electron-rich aromatic domain and Ca²⁺ near the electron-deficient carboxyl region. This asymmetric coordination considerably enhances the molecular dipole moment, contributing not only to local polarization but also to long-range electrostatic interactions. These effects may influence intermolecular associations and could potentially impact reactivity trends, although condensed-phase effects are beyond the scope of the present model.

These results indicate that the introduction of cations markedly reconstructs the charge distribution within the molecule, potentially alleviating internal strain and modulating electronic interactions. These ESP maps reveal how Na⁺ and Ca²⁺ ions induce spatially distinct electrostatic regions, providing a qualitative basis for subsequent analyses, such as BO evaluation and reactivity descriptors.

Changes in BO

To gain deeper insights into the bonding characteristics and electronic structural changes within the molecule, Löwdin BO analyses were conducted for each optimized structure at the B3LYP/6-31G* level. The analyses focused on the aromatic C–C bonds, C=O and C–O bonds in the carboxyl group, and the C–O bonds in phenolic –OH groups. Variations in BOs across different ionic environments are illustrated in Fig. 6–8.

Fig. 6 Löwdin bond order analysis of Yallourn low-rank coal model structures and associated reference compounds, computed at the B3LYP/6-31G* level of theory. Bond orders of the unmodified coal model (a: H) and Na⁺/Ca²⁺ co-exchanged model (b: Na/Ca-1) are compared to evaluate electronic and structural perturbations induced by cation coordination. Reference compounds (c) benzene, (d) (1R)-1,2,3,4-tetrahydronaphthalene-1-carboxylic acid, and (e) 4-methylbenzene-1,3-diol serve as benchmarks for typical bond orders of aromatic and functional moieties. Notable variations in C–C and C–O bond orders highlight localized electron redistribution and structural relaxation effects, particularly within the aromatic and carboxyl/phenolic domains, upon metal ion exchange.

Fig. 7 Difference analysis of Löwdin bond orders between the Na⁺/Ca²⁺ co-exchanged coal model (Na/Ca–1) and the unmodified Yallourn low-rank coal structure (H), calculated at the B3LYP/6-31G* level of theory. Bond-order differences (ΔBO) from the pristine model are shown in the right panel, with only bonds exhibiting |ΔBO| ≥ 0.05 highlighted for clarity. Bonds with ΔBO ≥ +0.05 are indicated by bold underline symbols, while those with ΔBO ≤ −0.05 are indicated by dashed underline symbols. These changes reflect localized electronic redistribution induced by simultaneous Na⁺ and Ca²⁺ incorporation.

Fig. 8 Löwdin bond order differences (ΔBO) in Yallourn low-rank coal structures upon mono- and co-addition of Na⁺ and Ca²⁺ ions, calculated at the B3LYP/6-31G* level of theory. Bond order changes are presented with respect to the unmodified structure (a: H), with only magnitude values ≥±0.05 marked for emphasis. Bonds exhibiting ΔBO ≤ −0.05 are indicated by dashed underline symbols to denote weakening, while those with ΔBO ≥ +0.05 are indicated by bold underline symbols to denote strengthening. For sites where ΔBO could not be determined due to structural mismatch, the original BO values are shown. The structures include Na⁺-exchanged models (d: Na, e: 2Na, f: 3Na), a Ca²⁺-exchanged model (b: Ca), and a Na⁺/Ca²⁺ co-exchanged model (c: Na/Ca-1). This figure demonstrates progressive bond polarization and electronic redistribution, particularly at carboxyl and phenolic functional groups, in response to increasing levels and combinations of cation exchange.

In the unmodified (pristine) lignite model (Fig. 6a), the aromatic C–C bonds exhibited BO values of approximately 1.40, while the C=O bond in the carboxyl group had a value of 1.88, and the C–O bond was 1.05. These values are consistent with those found in resonance-stabilized systems, indicating partial aromaticity and delocalized bonding in the polar functional groups. To interpret these values, we compared them with representative BO values from structurally related reference compounds—benzene (Fig. 6c), (1R)-1,2,3,4-tetrahydronaphthalene-1-carboxylic acid (Fig. 6d), and 4-methylbenzene-1,3-diol (Fig. 6e)—which served as benchmarks for the lignite model.

In the Na⁺/Ca²⁺ co-doped structure (Fig. 6b), Na⁺ was bound near the phenolic -OH group, while Ca²⁺ coordinated with the carboxyl group. These cations induced localized, complementary electronic effects, resulting in site-specific bond relaxation and enhanced electron delocalization.

Fig. 7 compares the bond orders (BOs) of the pristine and Na⁺/Ca²⁺ co-doped (Na/Ca-1) structures, along with the corresponding BO differences (ΔBO).

A ±0.05 threshold was adopted to differentiate meaningful BO variations from minor numerical fluctuations arising from the inherent precision limits of quantum-chemical calculations. This criterion was based on our previous study,⁹ in which a ΔBO of 0.05 was shown to correspond to bond-strength changes on the order of tens of kcal mol⁻¹ (≈10² kJ mol⁻¹).

Fig. 8 displays the ΔBO between the pristine and ion-exchanged structures, clarifying how specific bonds are weakened or strengthened. Marked bond relaxation is evident near the carboxyl and phenolic –OH groups. For example, in the Na⁺/Ca²⁺ co-doped model (Na/Ca-1: Fig. 8c), ΔBO values of −0.57 (C=O) and −0.24 (C–O) reflect a notable increase in flexibility. Additionally, the figure compares the effect of increasing Na⁺ or Ca²⁺ ions on BO variation. Specifically, as the number of Na⁺ ions increases from one to three (Fig. 8d–f), both the C=O and C–O bonds progressively weaken over broader regions, suggesting increased electron delocalization. Conversely, Ca²⁺ induces localized BO enhancement near its coordination points (Fig. 8b). In the co-doped models (Fig. 8c), these effects are spatially distinct, with Na⁺ acting near phenolic –OH sites and Ca²⁺ targeting carboxyl sites, enabling selective bond modulation.

These ΔBO trends agree well with the ESP distributions in Fig. 4 and 5 (Section “Visualization of electronic properties via ESP mapping”), indicating that Na⁺/Ca²⁺ co-addition induces both charge redistribution and bond flexibility. This level of control shows great potential for tuning reactivity and guiding catalytic design in lignite upgrading systems.

Mechanistic insights into ion-induced stabilization and reactivity

The total energy comparisons (Section “Optimized structures and total energies”) reveal that Na/Ca-1 exhibits a relative energy that is 110.2 kJ mol⁻¹ (0.04196 Hartree) lower than that of Na/Ca-2 (Fig. 1b). This stabilization suggests that Na⁺ and Ca²⁺ respectively preferentially coordinate with distinct functional groups—aromatic hydroxyl and carboxyl groups, resulting in a more favorable electrostatic and electronic configuration.

As discussed in the Section “Comparison of IR spectra,” the B3LYP/6-31G*-calculated IR spectra revealed characteristic shifts in the C=O and C–O stretching bands upon Na⁺ and Ca²⁺ coordination, which are in reasonable agreement with experimental observations. These vibrational features further support the hypothesis that ionic interactions induce structural changes at polar functional groups, consistent with the observed BO relaxations and charge redistributions described below. ESP maps (Section entitled “Visualization of electronic properties via ESP mapping”) revealed a clear spatial separation of positive and negative charge regions in the co-doped structure. This electronic asymmetry originates from the site-selective binding of the two cations and contributes to localized polarization effects. Specifically, π-electron delocalization was promoted around the aromatic ring, while a strong positive potential formed near the carboxyl group. This charge separation is expected to modulate intramolecular interactions and enhance chemical reactivity.

BO analysis (Section “Changes in BO”) further supports this interpretation, showing that the carboxyl C=O and C–O bonds are markedly relaxed in the co-doped structure (Fig. 8c). This weakened bonding environment implies increased molecular flexibility, which is consistent with the charge redistribution observed in the ESP maps (Section entitled “Visualization of electronic properties via ESP mapping”). The relaxation of key polar bonds suggests that ion exchange facilitates structural destabilization at specific sites, potentially lowering the activation barriers for decomposition reactions, such as those related to CO₂ or H₂O elimination. In addition, energy calculations for these elimination reactions (Section “Changes in BO”) demonstrated that the co-doped model (Na/Ca-1) exhibited lower dissociation energies compared with the other structures. This implies that electron density redistribution and bond relaxation facilitate the initial bond cleavage steps required for gasification. These theoretical findings provide a molecular-level explanation for the enhanced gasification performance observed in our previous study, where Ca²⁺/Na⁺ co-exchanged lignite samples exhibited markedly faster steam gasification at lower temperatures than their singly exchanged counterparts.⁴ The selective bond weakening and electron redistribution reported herein suggest that the co-addition of these two ions effectively lowers the activation energy for thermal decomposition, thereby promoting reactivity under gasification conditions.

Taken together, these findings demonstrate that the co-addition of Na⁺ and Ca²⁺ induces complementary and site-selective electronic effects, leading to local structural relaxation and electron redistribution, which in turn contribute to enhanced chemical reactivity. Na⁺ and Ca²⁺ do not interfere with each other but instead act selectively on different regions of the molecule, enabling spatially controlled modulation of the electronic structure. This molecular-level insight provided a theoretical basis for the experimental observations reported in our previous work, where nanoscale Na/Ca co-catalysts markedly enhanced the steam gasification of lignite.⁴ The current quantum chemical results explain that this enhancement arises from selective coordination, bond relaxation, and charge redistribution at the atomic scale. In conclusion, the results presented herein offer valuable guidance for catalyst design in coal upgrading. The strategic co-addition of metal ions with complementary binding behaviors may improve lignite utilization and provide broader applicability to other low-grade carbonaceous resources, such as biomass.

However, it is important to note that the molecular model used in this study represents a simplified approximation of lignite structure and does not explicitly account for condensed-phase effects, such as intermolecular interactions or mineral matter. Additionally, the calculations were performed under gas-phase conditions, which may have not fully captured the behavior of lignite during actual gasification processes. These limitations should be considered when extrapolating the findings to real-world systems.

Conclusions

This study theoretically elucidates how the individual and combined addition of Na⁺ and Ca²⁺ ions modifies the structural and electronic properties of a Yallourn lignite model. The results demonstrate that ion exchange induces molecular relaxation, electron redistribution, and bond localization changes. Na⁺ preferentially interacts with π-systems, whereas Ca²⁺ exhibits strong electrostatic coordination with oxygen atoms. Notably, the co-addition of Na⁺ and Ca²⁺ produces complementary effects, stabilizing the structure and enhancing molecular polarity to a greater extent than either ion alone. These findings are consistent with previous experimental reports showing enhanced gasification reactivity upon dual-ion treatment. The combination of structural analysis, electrostatic mapping, and BO quantification provides a solid theoretical foundation for understanding the electronic basis of lignite activation by metal ions. These insights contribute to the rational design of catalytic systems for the efficient utilization of lignite and similar carbonaceous resources. Future studies should validate the predicted electronic effects experimentally and explore other metal ions with distinct valence states and coordination chemistries.

Author contributions

Yuji Shinohara: conceptualization, methodology, quantum chemical calculations, data curation, formal analysis, visualization, writing – original draft, writing – review & editing. Naoto Tsubouchi: supervision, project administration, scientific advice, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

No new experimental data were generated or analyzed in this study. All results are based on density functional theory calculations.

Data supporting the findings of this study are provided within the article and the supplementary information (SI). Supplementary information: including representative optimized Cartesian coordinates (Section S1), selected metal–oxygen distances (Table S1), full vibrational frequency/intensity datasets (unscaled and scaled; Tables S2–S5), representative IR mode assignment tables (Tables S6–S9), single-point PCM (water) total energies (Table S10), and a limited quantitative deviation analysis for selected experimental reference bands (Table S11). See DOI: https://doi.org/10.1039/d5nj04547a.

Acknowledgements

This research study did not receive any financial support from public, commercial, or not-for-profit funding agencies. The authors thank the Steel Foundation for Environmental Protection Technology for technical support and helpful discussions.

References

1. Y. Wang, Y. Zhao, R. He, Z. Yan, X. Li, H. Zhou, N. Li, K. Zhi, Y. Song, Y. Teng and Q. Liu, ACS Omega, 2021, 6(38), 24848–24858,  DOI:10.1021/acsomega.1c03695.
2. B. Zhang, Z. Tian, Q. Wang, R. Jia, D. Ma and X. Guilin, Energy, 2025, 327, 136352,  DOI:10.1016/j.energy.2025.136352.
3. Y. Yin, Z. Wu, J. Tao, C. Qi, W. Zhang and S. Cheng, Investigation of the evolution of the chemical structure of bituminous coals and lignite during pyrolysis, Crystals, 2022, 12(4), 444,  DOI:10.3390/cryst12040444.
4. N. Tsubouchi, Y. Mochizuki, Y. Shinohara, Y. Hanaoka, T. Kikuchi and Y. Ohtsuka, Energy Fuels, 2018, 32(1), 226–232,  DOI:10.1021/acs.energyfuels.7b02964.
5. Y. Li, Y. Yao, X. Liu, H. Sun and W. Ni, Fuel, 2013, 109, 527–533,  DOI:10.1016/j.fuel.2013.03.010.
6. W. Wang, R. Lemaire, A. Bensakhria and D. Luart, J. Anal. Appl. Pyrolysis, 2022, 163, 105479,  DOI:10.1016/j.jaap.2022.105479.
7. W. He, C. Hao, C. Deng, J. Gao, Y. Wang and Y. Qin, Fuel, 2026, 404, 136341,  DOI:10.1016/j.fuel.2025.136341.
8. Y. Shinohara and N. Tsubouchi, J. Mol. Struct., 2020, 1218, 128544,  DOI:10.1016/j.molstruc.2020.128544.
9. Y. Shinohara and N. Tsubouchi, ACS Omega, 2020, 5(3), 1688–1697,  DOI:10.1021/acsomega.9b03780.
10. Y. Shi, Y. Zhu, S. Chen, Y. Wang and Y. Song, Energies, 2024, 17(5), 1049,  DOI:10.3390/en17051049.
11. R. Bai, L. Liu, N. Li, Q. Liu, Y. Meng, Y. Teng, H. Zhou and X. Zhou, Fuel, 2023, 349, 128654,  DOI:10.1016/j.fuel.2023.128654.
12. C. Ugwumadu, R. Olson III, N. L. Smith, K. Nepal, Y. Al-Majali, J. Trembly and D. A. Drabold, Nanotechnology, 2023, 35(9), 095703,  DOI:10.1088/1361-6528/ad1058.
13. L. R. Radovic, P. L. Walker Jr. and R. G. Jenkins, Fuel, 1983, 62(2), 209–212,  DOI:10.1016/0016-2361(83)90200-4.
14. Y. Zhao, C. Xing, C. Shao, G. Chen, S. Sun, G. Chen, L. Zhang, J. Pei, P. Qiu and S. Guo, Fuel, 2020, 278, 118229,  DOI:10.1016/j.fuel.2020.118229.
15. C. Wang, Y. Xing, K. Shi, S. Wang, Y. Xia, J. Li and X. Gui, ACS Omega, 2024, 9(1), 1881–1893,  DOI:10.1021/acsomega.3c08574.
16. X. Wang, S. Wang, Y. Zhao, R. Mu and Y. Shao, ACS Omega, 2023, 8(49), 46528–46539,  DOI:10.1021/acsomega.3c04985.
17. D. Zhang, Y. Li, C. Zi, Y. Zhang, X. Hu, G. Tian and W. Zhao, ACS Omega, 2021, 6(15), 10281–10287,  DOI:10.1021/acsomega.1c00627.
18. K. Iwata, H. Itoh and K. Ouchi, Fuel Process. Technol., 1989, 3(3–4), 221–229,  DOI:10.1016/0378-3820(80)90021-1.
19. T. Isoda, H. Takagi, H. Saiki, K. Kusakabe and S. Morooka, J. Energy Soc. Jpn, 2000, 79(6), 511–521,  DOI:10.3775/jie.79.511.
20. SPARTAN’18 for Windows (Revision 1.3.0), Wavefunction Inc., Irvine, CA, 2018, https://www.wavefun.com/.
21. T. A. Halgren, J. Comput. Chem., 1996, 17(5–6), 490–519 , https://onlinelibrary.wiley.com/doi/10.1002/(SICI)1096-987X(199604)17:5/6%3C490::AID-JCC1%3E3.0.CO;2-P.
22. G. A. Carlson, Energy Fuels, 1992, 6(6), 771–778,  DOI:10.1021/ef00036a012.
23. A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, Courier Corporation, 1996, ISBN 0-486-69186-1.
24. A. D. Becke, J. Chem. Phys., 1993, 98(7), 5648–5652,  DOI:10.1063/1.464913.
25. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37(2), 785–789,  DOI:10.1103/PhysRevB.37.785.
26. P.-O. Löwdin, Phys. Rev., 1955, 97(6), 1474–1489,  DOI:10.1103/PhysRev.97.1474.
27. M. U. Khan, M. R. S. A. Janjua, J. Yaqoob, R. Hussain, M. Khalid, A. Syed, A. M. Elgorban and N. S. S. Zaghloul, J. Photochem. Photobiol., A, 2023, 440, 114667,  DOI:10.1016/j.jphotochem.2023.114667.
28. C. Lei, Z. Yang, B. Zhang, M.-H. Lee, Q. Jing, Z. Chen, X.-C. Huang, Y. Wang, S. Pan and M. R. S. A. Janjua, Phys. Chem. Chem. Phys., 2014, 16, 20089–20096,  10.1039/C4CP02539C.
29. A. Umar, J. Yaqoob, M. U. Khan, R. Hussain, A. Alhadhrami, A. S. A. Almalki and M. R. S. A. Janjua, J. Phys. Chem. Solids, 2022, 169, 110859,  DOI:10.1016/j.jpcs.2022.110859.
30. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627,  DOI:10.1021/j100096a001.
31. P. K. Ramya, A. Shivhare, M. M. Deshmukh and C. H. Suresh, Phys. Chem. Chem. Phys., 2024, 26(42), 27054–27065,  10.1039/D4CP03371J.

---