Polymer-metal oxide interfaces in XHNBR/PA6 blends: computational insights toward sustainable crosslinking

Abstract

Polyamide 6 (PA6) is a high-performance thermoplastic widely used in engineering applications, while carboxylated hydrogenated nitrile rubber (XHNBR) provides viscoelastic damping and reactive carboxyl groups for efficient ionic crosslinking with metal oxides. Their combination enables thermoplastic elastomer systems with improved toughness and vibration attenuation. In this work, the interfacial interactions between PA6, XHNBR, and representative metal oxides (ZnO, MgO, CaO, and MgO₂) were investigated using density functional theory (DFT) combined with a conformational clustering approach. Cluster models of the oxides, (MO)₁₂ and (MO₂)₆, were employed to reproduce local coordination environments and enable detailed electronic characterization through natural bond orbital (NBO) and quantum theory of atoms in molecules (QTAIM) analyses. The results reveal distinct binding preferences among polymer functional groups: carboxyl moieties exhibit the strongest and most complex interactions, involving proton transfer and metal–oxygen coordination, while amine and amide groups form weaker, primarily electrostatic contacts. Among the oxides, CaO produces the most exothermic and predominantly ionic interactions, making it a promising, less toxic alternative to ZnO for crosslinking applications. To account for conformational flexibility, a clustering-based sampling strategy was applied to PA6-XHNBR dimers, allowing exploration of the configurational landscape and evaluation of Boltzmann-weighted interaction energies. The analysis demonstrates that conformers with significant population weights govern the effective interfacial stabilization, underscoring the importance of conformational diversity in accurately describing polymer–oxide and polymer–polymer interfaces. These findings provide molecular-level guidelines for designing greener damping materials for automotive and industrial applications.

---

1 Introduction

The reduction of structural vibrations remains a key performance requirement across industries such as automotive, railway, and aeronautics. A widely adopted strategy to address this challenge is the incorporation of viscoelastic materials, which can dissipate vibrational energy and thereby improve damping efficiency.^1–3 Within this context, thermoplastic elastomers (TPEs) have gained significant attention. These materials combine the elasticity of rubbers with the processability of thermoplastics, offering a unique balance between mechanical performance and manufacturing versatility.⁴ A related class of materials, thermoplastic vulcanizates (TPVs), are obtained through dynamic vulcanization of rubber dispersed in a thermoplastic matrix, which further enhances their damping behavior and broadens their application window.⁵

Among the practical routes to TPEs and TPVs, melt mixing of rubbers and thermoplastics is particularly attractive for its simplicity, scalability, and compatibility with standard compounding operations. In this work, we focus on a hybrid system formed by polyamide 6 (PA6) as the thermoplastic matrix and carboxylated hydrogenated nitrile butadiene rubber (XHNBR) as the elastomeric component. XHNBR is a strong candidate for damping applications owing to its viscoelastic response and the presence of carboxyl functionalities that enable efficient chemical or ionic crosslinking.⁶ A conventional approach to enhance crosslink density and damping in carboxylic rubbers is the addition of metal oxide (MO) nanoparticles as crosslinking agents. Besides, when blended with PA6, XHNBR can impart improved damping while contributing to the overall toughness of the composite, making the blend relevant for demanding engineering uses.

The processing method strongly influences the microstructure and, therefore, the properties of the material. If metal oxide nanoparticles are first dispersed within the rubber and only then blended with PA6, the system tends to form MO-XHNBR domains embedded in a continuous PA6 matrix (Fig. 1a). Conversely, when the three components are combined simultaneously, configurations that promote direct contact among rubber, PA6, and MO at a shared boundary become more likely (Fig. 1b).

Fig. 1 Schematic representation of the two composite architectures considered in this study. (a) MO nanoparticles coated with XHNBR and dispersed within a PA6 matrix. (b) Fully mixed system in which MO, XHNBR and PA6 are homogeneously blended without a defined interphase.

These two limiting scenarios, schematically depicted in Fig. 1, naturally lead to three distinct interphases: metal nanoparticle-rubber (MO-XHNBR), metal nanoparticle-thermoplastic matrix (MO-PA6), and rubber-thermoplastic matrix (XHNBR-PA6). The MO-PA6 and XHNBR-PA6 interphases might influence morphology, load transfer, and energy dissipation and must be considered in a comprehensive description of the material. In contrast, the MO-XHNBR interphase is expected to dominate the crosslinking chemistry and, in turn, the damping response, because the oxide surface sites interact strongly with the carboxyl groups of XHNBR.⁷ Metal oxide nanoparticles, such as ZnO and Al₂O₃, provide mechanical reinforcement, attributed to enhanced MO-PA6 interfacial adhesion.⁸ Meanwhile, the XHNBR-PA6 interfacial interaction may induce changes in morphology and cristallinity due to hydrogen bonding and dipolar interactions between the carboxyl and nitrile groups of XHNBR and the amine or amide end groups of PA6.⁹

Despite the technological importance of such multiphase systems, the molecular mechanisms that control polymer–metal oxide interactions and their translation into macroscopic damping remain insufficiently understood. Much of the existing literature emphasizes macroscopic performance metrics, offering limited insight into how oxide surface chemistry, dispersion state, and polymer conformations couple at the interface. Computational studies have provided valuable clues, but often at the cost of oversimplification—e.g., using small organic fragments or idealized oxide clusters that do not fully capture the structural complexity of the polymer matrix or the extended nature of realistic oxide surfaces.^10–12

To fill these gaps, we perform interphase-resolved density functional theory (DFT) calculations to study the atomic-scale interactions between XHNBR/PA6 blends and different metal oxides. Our objective is to characterize the interactions at the MO-XHNBR, MO-PA6, and XHNBR-PA6 interphases, and to relate these molecular signatures to processing routes and targeted damping performance. Historically, zinc oxide (ZnO) has been widely employed as a crosslinking agent because of its strong reactivity toward carboxyl groups; however, it presents drawbacks including premature reactivity, dispersion challenges, and environmental toxicity concerns.¹³ These limitations motivate the search for safer and more sustainable alternatives. Therefore, rather than facusing only on ZnO nanoclusters, we also perform a screening of different MOs to investigate their crosslinking performance and identify greener substitutes.¹⁴

To address the complexity of the XHNBR/PA6 blend and the variety of metal oxides considered, we develop a systematic workflow that integrates a clustering procedure to identify representative configurations, reducing the effective configurational space while preserving chemically relevant diversity. This methodology enables a consistent comparison across oxides and interphases, with particular emphasis on the MO-XHNBR interactions that primarily govern crosslink formation. By establishing molecular-level guidelines for interfacial design, our results provide insights for the development of safer, greener, and more effective damping materials for advanced engineering applications.

2 Computational procedure and models

2.1 Quantum chemical calculations

All quantum chemical calculations were performed within the framework of density functional theory (DFT) using the Gaussian 16 program package.¹⁵ Geometry optimizations and harmonic frequency calculations were carried out at the ωB97XD/6-31+G(d,p) level of theory. Zero-point vibrational energy (ZPVE) and thermal corrections to the enthalpy and Gibbs free energy at 298 K were obtained from the frequency analyses. Single-point energy refinements were performed on the optimized geometries at the ωB97XD/6-311++G(2df,2p) level to improve the electronic energy description. The total enthalpy (H) of each species was evaluated as H = E_elec + H^corr, where E_elec is the single-point electronic energy and H^corr represents the vibrational and thermal corrections. Interaction enthalpies (ΔH) were calculated using the supermolecular approach.

Periodic DFT calculations were conducted using the VASP package^16–19 with the PBE0 functional and projector-augmented wave (PAW) potentials.

2.1.1 Analysis of bonding interactions. The nature of the interactions was analyzed by means of the natural bond orbital (NBO) method,^20–23 using the Gaussian 16 suite of programs,¹⁵ and the quantum theory of atoms in molecules (QTAIM),^24–27 using the AIMAll program.²⁸ In QTAIM, the bond critical point (BCP) defines the bond path connecting two atoms within a molecule, and its properties serve to characterize the nature of their interaction.²⁴ The BCP parameters analyzed include the electron density (ρ), the Laplacian of the electron density, (∇²ρ), the kinetic (G), the potential (V), and the total (H) electron energy densities. These parameters provide information on the density of electrons in a given bond path. The negative value of ∇²ρ denotes electron density concentration in the interatomic region of a bond, corresponding to covalent bonds, whereas a positive value indicates depletion of electron density usually attributed to non-covalent interactions. If the ∇²ρ is positive but H is negative, the interaction is classified as partially covalent. If ρ > 0.03 a.u. and H is negative, the interaction is electrostatic with partially covalent characteristics.²⁹

2.2 Clustering

The conformational complexity of the XHNBR-PA6 interphase is very high due to the large number of internal degrees of freedom, especially in dimer models. To efficiently explore this configurational space, we first generated a conformer ensemble using the CREST program,³⁰ employing an energy window of 9 kcal mol⁻¹. As CREST retains only conformers within this window relative to the lowest-energy structure, the number of sampled conformers reflects not only the size of the system but also the strength and directionality of the intermolecular interactions present. In particular, for the combined PA6-XHNBR system, strong stabilizing interactions restrict the accessible low-energy configuration space, resulting in fewer retained conformers compared to the more flexible isolated systems. The resulting conformers were then subjected to a clustering analysis in order to identify representative structures for further refinement at a higher level of theory. This procedure allows a systematic reduction of redundancy while retaining chemically meaningful diversity.

2.2.1 K-means clustering. Given a set of n observations {x₁, x₂, …, x_n}, where each x_i is a d-dimensional real vector, the K-Means algorithm seeks to partition the data into k(≤n) disjoint subsets S = {S₁, S₂, …, S_k}. The objective is to minimize the within-cluster sum of squares (WCSS), i.e., the total variance within each cluster. Formally, the optimization problem iswhere µ_i denotes the centroid of cluster S_i,

|S_i| is the number of elements in S_i, and ‖·‖ is the Euclidean (L²) norm.

This criterion is equivalent to minimizing the average pairwise squared distance between points within the same cluster:

The equivalence follows from the identity

Since the total variance of the dataset is fixed, minimizing the WCSS is equivalent to maximizing the variance between clusters, also known as the between-cluster sum of squares (BCSS).³¹

In the present work, the number of clusters was fixed to k = 10 as a practical compromise between structural diversity and computational feasibility. Selecting too few clusters would not adequately capture the geometric variability of the conformational ensemble generated by CREST, whereas a larger number of clusters would substantially increase the number of subsequent high-level DFT optimizations required for representative structures, rendering the analysis computationally prohibitive. This choice allows the main conformational basins to be represented while keeping the electronic-structure calculations tractable and consistently applied across all systems.

2.2.2 Dimensionality reduction. To visualize the clustering and highlight structural trends, the descriptor space was reduced to two principal components using Principal Component Analysis (PCA). PCA projects the data into an orthogonal basis that maximizes variance, obtained by solving the eigenvalue problem

Cv_j = λ_jv_j, (5)

where C is the covariance matrix of the descriptor set, λ_j are the eigenvalues, and v_j the corresponding eigenvectors. The first two components (v₁, v₂) capture the largest variance and were used to generate 2D projections of the conformational space.^32,33

This clustering and selection protocol ensures that the chosen structures are both computationally tractable and structurally representative of the configurational diversity of the XHNBR-PA6 interphase.

2.2.3 Boltzmann-weighted population. After identifying the representative conformers through clustering, their relative thermodynamic stability was assessed using Boltzmann statistics. This analysis allows quantifying the population distribution of each conformer at thermal equilibrium and provides a more realistic weighting for subsequent energy-averaged properties.³⁴

Given a set of conformers with electronic energies E_i (in kcal mol⁻¹), the Boltzmann population p_i of each conformer is expressed as

where E_min is the energy of the most stable conformer, k_B is the Boltzmann constant (1.9872 × 10⁻³ kcal mol⁻¹ K⁻¹), and T is the temperature (298.15 K).

The Boltzmann-weighted average energy 〈E〉 of the ensemble is then given by

This quantity represents the thermodynamic mean energy at equilibrium, accounting for the statistical contribution of all conformers.

2.3 Molecular models

To capture the variety of interactions present in the multiphase system, we employed a set of representative molecular models of the metal oxides (MO), XHNBR, and PA6, as illustrated in Fig. 2.

Fig. 2 Molecular models used to represent the different components of the XHNBR/PA6/metal oxide system. (a–c) Metal oxide clusters with increasing size, (MO)_x = 6, 9, and 12, used to approximate small nanoparticles (metal atoms shown in purple and oxygen atoms in red). (d) Periodic slab model of the metal oxide, representative of extended surfaces and larger particles. (e–h) XHNBR models: (e) carboxyl fragment, (f) nitrile fragment, (g) XHNBR monomer, and (h) XHNBR dimer. (i and j) PA6 models: (i) amide fragment, (j) PA6 dimer.

Different levels of structural complexity were considered depending on the specific goal of each simulation. For metal oxides, small cluster models were chosen to approximate the reactivity of nanoparticles of reduced size and quantum effects, while periodic surface models were used to describe interactions with larger particles. For the polymers, simplified molecular fragments were introduced to reduce computational cost while retaining chemically relevant functionality: characteristic functional groups of XHNBR and PA6 were studied individually, and monomeric and dimeric units were also included to capture larger-scale conformational effects.

2.3.1 MO models. Four different representations of the metal oxides were considered (Fig. 2a–d). To approximate small nanoparticles, cluster models of (MO)_x with x = 6, 9, and 12 were constructed. These clusters provide a simple yet chemically meaningful way to probe the size dependence of local binding motifs while keeping the calculations computationally affordable. In addition to these finite clusters, a periodic slab model was included to describe the crystalline surfaces of the oxides (Fig. 2d). This representation captures long-range periodicity and extended coordination environments, thus allowing the investigation of interfacial interactions more relevant to larger nanoparticles. For MgO, we selected the (100) surface, which is both the thermodynamically most stable facet and the most widely studied in the literature.³⁵ The lattice parameter was fixed at 4.190 Å, in agreement with experimental data. In the case of ZnO, the wurtzite crystal structure was modeled using the (0001) polar surface, which is characteristic of hexagonal ZnO. Lattice parameters of a = 3.249 Å and c = 5.206 Å were adopted, consistent with reported crystallographic values.³⁶ These two oxides were prioritized because ZnO has been the conventional choice in rubber formulations, while MgO is the most frequently proposed greener alternative. Together, the cluster and slab models provide a systematic framework for analyzing both local and extended interaction environments.

2.3.2 XHNBR models. The elastomeric component was represented through four molecular models (Fig. 2e–h), chosen to balance chemical detail and computational tractability. First, individual functional groups were isolated to probe their specific roles in oxide interactions: a carboxyl fragment, which is the primary site for crosslinking reactions, and a nitrile fragment, which may contribute to secondary interactions or influence polarity. Beyond these fragments, a single monomer unit of XHNBR was modeled to study the combined influence of multiple functional groups within a chemically realistic repeating unit. Finally, a dimer model was introduced, which incorporates conformational flexibility and begins to reflect the structural complexity of the polymer chain. Together, these representations allow us to systematically dissect the contribution of functional chemistry, local environment, and conformational freedom to the overall behavior of XHNBR at oxide interfaces.

While these models do not capture long-range chain flexibility on the scale of multiple persistence lengths, they are designed to isolate local chemical interactions at the interface, which are the primary focus of the present electronic-structure analysis.

2.3.3 PA6 models. To represent the thermoplastic PA6 matrix, a simplified set of models was employed (Fig. 2i and j). Since the polymer backbone is overwhelmingly composed of amide repeat units, the central model used in this work is an isolated amide fragment. This model was analyzed through two complementary interaction motifs: one in which the metal oxide engages primarily with the carbonyl oxygen (–O), and another in which the interaction occurs via the amide nitrogen (–N). Treating these two binding modes separately allows a clear distinction between the oxygen-driven and nitrogen-driven coordination patterns that may arise along the polymer backbone. Finally, a dimeric unit of PA6 was considered, which introduces backbone flexibility and enables the study of local conformational effects beyond single functional groups. This hierarchical modeling strategy ensures that both reactive end groups and backbone-specific interactions are explicitly accounted for in the simulations, while maintaining a level of simplicity suitable for atomistic studies.

We note that the use of fragment and dimer models does not aim to represent the full polymer-chain flexibility expected at larger length scales, but rather to capture local binding motifs and short-range conformational effects, which can be reliably addressed within a quantum-chemical framework.

3 Results and discussion

To systematically unravel the interfacial phenomena governing the multiphase behavior of the studied system, the results are organized in a progressive manner from the simplest to the most complex interactions. Section 3.1 presents the initial screening of MO/MO₂ clusters to identify representative candidates for subsequent analysis. Section 3.2 focuses on the interactions between XHNBR and the selected MO/MO₂, including both cluster and periodic surface models, and is subdivided into three parts covering molecular fragments, monomeric and dimeric units, and extended surfaces. Section 3.3 then addresses the PA6-MO interface following a similar analytical framework, encompassing geometric analysis, interaction energetics, and electronic characterization through NBO and QTAIM analyses, and the PA6-XHNBR interface.

3.1 Screening of (MO)₆ and (MO₂)₆ clusters as an alternative to ZnO

A preliminary computational screening was conducted on alkaline earth oxide clusters (MO)₆ and their corresponding peroxide analogues (MO₂)₆ (M = Mg, Ca, Sr) to evaluate their potential as substitutes for ZnO in rubber crosslinking applications. The interaction of these clusters with the XHNBR carboxyl model (Fig. 2e) was used as a representative test of their ability to form stable interfacial bonds with the polar functional groups of the polymer.

The computed interaction enthalpies (ΔH), displayed in Fig. 3, reveal a consistent trend in which alkaline earth oxides exhibit stronger interactions with the carboxyl group than their peroxide counterparts. This behavior reflects the higher basicity and stronger Lewis acidity of the metal cations in the oxides, which favor more robust ionic and coordinative bonding with the carboxylate oxygen atoms.

Fig. 3 Interaction enthalpies (ΔH), in kcal mol⁻¹, between the XHNBR carboxyl model and different oxide and peroxide clusters, (MO)₆ and (MO₂)₆ (M = Mg, Ca, Sr). The red bar corresponds to the reference (ZnO)₆, while the blue bars indicate the candidates selected for further analysis.

Based on the comparison with (ZnO)₆, three systems were selected for detailed analysis in the subsequent sections: MgO, CaO, and MgO₂ (blue bars in Fig. 3). The selection criterion was primarily the similarity of their interaction energies to that of ZnO, ensuring comparable reactivity and interfacial stability. In particular, MgO and MgO₂ display interaction enthalpies within a few kcal mol⁻¹ of the reference, suggesting that they could reproduce ZnO-like behavior while offering improved environmental compatibility.

The choice of MgO₂ (instead of CaO₂, which exhibits a slightly closer ΔH value) is further justified by its strong industrial relevance. Magnesium-based additives are widely used in rubber and polymer formulations, where MgO acts as an acid acceptor, cure activator, and stabilizer in halogenated elastomers.³⁷ This well-established role of magnesium compounds in elastomer chemistry makes MgO₂ a more pragmatic candidate for translation from computational models to experimental applications, aligning the theoretical screening with realistic material considerations.

3.2 Interaction of XHNBR with MO/MO₂

Following the preliminary screening, the interaction of XHNBR with the selected MO/MO₂ was investigated in detail. From this point onward, only the (MO)₁₂/(MO₂)₆ cluster models were considered, as they provide a suitable compromise between computational cost and the accurate representation of the local coordination environment in oxide nanoparticles. This section is organized as follows: first, we analyze the interactions between the individual XHNBR functional models (Fig. 2e and f) and the MO/MO₂ clusters, followed by the monomeric and dimeric XHNBR units (Fig. 2g and h) to account for cooperative effects and intramolecular constraints. Subsequently, (MO)_x/(MO₂)_x periodic surface models (Fig. 2d) are examined to capture long-range interfacial effects.

Representative interaction geometries between (ZnO)₁₂ and the four XHNBR models are shown in Fig. 4 (the coordinates of all optimized structures can be seen in Section 3 in SI). Each configuration illustrates the preferred binding arrangement for the carboxyl and nitrile functional groups, both individually and in combination within the monomer and dimer models. These structures serve as reference geometries for the subsequent comparison with the alternative oxides, providing a baseline for assessing how the oxide composition influences the interfacial binding motifs.

Fig. 4 Optimized structures of (ZnO)₁₂ and the different models of XHNBR, obtained at the ωB97XD/6-31+G(d,p) level of theory using Gaussian 16 (see Section 2.1). (a) Carboxyl group, (b) nitrile group, (c) monomer and (d) dimer.

3.2.1 Interactions of XHNBR functional models with MO/MO₂ cluster models. To gain deeper insight into the nature of the interactions between XHNBR and the metal clusters, a detailed structural analysis was performed. The following section presents an evaluation of the interaction focusing on key interatomic distances and interaction energies, and their implications for binding strength and stability.

The initial analysis of the XHNBR-MO/MO₂ interface focused on the fundamental interaction motifs between the isolated functional groups of XHNBR and the metal clusters. Two fragment models were employed to represent the reactive sites of the polymer: the carboxyl and nitrile groups (Fig. 2e and f). These simplified systems enable a direct evaluation of the intrinsic binding behavior of each functionality, isolating their individual contributions to the overall interfacial interaction. Selected bonds and interactions have been used for the discussion and labeled according to Fig. 5.

Fig. 5 Schematic representation of the interaction and the relevant bonds between the XHNBR functional groups, (a) carboxyl and (b) nitrile, and the (MO)₁₂ clusters.

The optimized geometries and interaction enthalpies (ΔH) are summarized in Table 1. For the nitrile group, the R₁ distance remains nearly constant at approximately 1.15 Å, indicating preservation of the C≡N triple bond and negligible perturbation of its electronic structure upon adsorption. The R₂ distance, describing the coordination between the nitrile nitrogen and the cluster metal atom, exhibits moderate variations among the different MO/MO₂, primarily reflecting the intrinsic size and Lewis acidity of the cation. In particular, CaO presents a notably longer R₂ distance and a smaller interaction enthalpy compared with MgO and ZnO. This behavior arises from the lower charge density and larger ionic radius of Ca²⁺ (r_Ca²⁺ ≈ 0.99 Å vs. r_Mg²⁺ ≈ 0.72 Å and r_Zn⁺² ≈ 0.74 Å), which reduce its ability to polarize and attract the nitrile lone pair, leading to a weaker coordination. As a result, the nitrile group interacts more weakly and at a longer distance from the Ca centers.

Table 1 Bond distances (R_i, in Å) of the optimized structures, following the nomenclature defined in Fig. 5, and corresponding interaction enthalpies (ΔH, in kcal mol⁻¹). Cluster-based results were obtained from Gaussian 16 calculations at the ωB97XD/6-311++G(2df,2p)//ωB97XD/6-31G+(d,p) level of theory, while periodic surface results were obtained from VASP calculations (see Section 2.1). The upper section of the table corresponds to cluster models ((MO)₁₂ and (MO₂)₆), while the lower section corresponds to periodic surface models ((MO)_x and (MO₂)_x). Cluster models were employed to capture local interaction characteristics at the oxide surface, whereas periodic models were used to assess the influence of extended surface periodicity

	Nitrile			Carboxyl
	R1	R2	ΔH	R1	R2	R3	R4	R5	ΔH
(ZnO)12	1.151	2.093	−15.61	1.515	1.021	1.929	1.252	1.279	−37.92
(MgO)12	1.152	2.209	−14.73	1.577	1.010	1.982	1.255	1.273	−48.20
(CaO)12	1.154	2.595	−9.13	2.413	0.965	2.392	1.263	1.263	−68.44
(MgO2)6	1.151	2.147	−21.45	1.064	1.408	2.006	1.292	1.249	−40.52
(ZnO)x	1.16	2.05	−50.88	1.32	1.12	1.97	1.27	1.28	−32.84
(MgO)x	1.16	2.25	−42.55	1.42	1.08	2.06	1.29	1.27	−27.17
(CaO)x	1.16	2.62	−40.89	1.82	0.99	2.39	1.28	1.27	−46.20
(MgO2)x	1.16	2.19	−47.54	1.08	1.44	2.05	1.32	1.25	−25.60

---

In contrast, the carboxyl group exhibits more pronounced structural rearrangements upon interaction with the oxide clusters. Proton transfer from the carboxylic OH group to the cluster occurs in all cases except for MgO₂, as indicated by the short O–H distance to the cluster (R₂ ≈ 1.0 Å) in MgO, ZnO, and CaO. This transfer leads to the formation of a surface hydroxyl species and a deprotonated carboxylate moiety. In these systems, the internal C–O bond lengths (R₄ and R₅) become nearly equivalent, consistent with increased electronic delocalization within the carboxylate group. In MgO₂, however, the proton remains largely associated with the carboxyl group (R₂ = 1.408 Å), and the two C–O bond distances (R₄ and R₅) are more different, reflecting the absence of full deprotonation and limited delocalization. These differences highlight the lower proton affinity and distinct reactivity of the peroxide surface, which also contributes to its smaller interaction enthalpy relative to the oxide clusters.

Thus, carboxyl-functionalized fragments exhibit markedly stronger binding to the oxide clusters than their nitrile counterparts. This difference arises from the carboxyl group's capacity to form multiple stabilizing interactions, electrostatic, hydrogen bonding, and, in most cases, proton transfer, whereas the nitrile primarily coordinates to surface metal atoms through its nitrogen lone pair. Consequently, the carboxyl group is expected to play a dominant role in the interfacial adhesion and crosslinking activity of XHNBR on metal oxide surfaces.

To further characterize the nature of the interactions at the metal oxide cluster/XHNBR interface, NBO (Table 2) and QTAIM (Table 3) analyses were performed. For the nitrile group, the interaction corresponding to the coordination between the nitrile nitrogen and the metal atom, is predominantly electrostatic in all cases except for ZnO, where it acquires a partially covalent character. This stronger polarization in ZnO originates from the higher polarizability and moderate Lewis acidity of Zn²⁺, which facilitate enhanced charge delocalization from the nitrile lone pair into the metal orbitals. In contrast, Mg²⁺ and especially Ca²⁺, characterized by lower charge density and larger ionic radii, form weaker and more distant electrostatic contacts, consistent with the longer R₂ distances and smaller interaction enthalpies reported in Table 1. This trend reflects the decreasing ability of the metal center to stabilize donor–acceptor interactions as the cation becomes more ionic and less polarizing. For the carboxyl group, the nature of the interactions is influenced by the occurrence of proton transfer. In MgO, ZnO, and CaO, where proton transfer occurs, the O⋯H interaction in the carboxyl group (bond 1) is now electrostatic, while the O⋯H interaction with the cluster (bond 2) exhibits covalent character. The M⋯O interactions (bond 3), representing interactions between the carboxylate oxygen and nearby metal atoms, are predominantly electrostatic. In MgO₂, proton transfer is absent, leading to a different bonding pattern: bond 1 shows covalent character corresponding to the carboxylic O–H bond, while bond 2 remains electrostatic.

Table 2 NBO results for the interactions between characteristic XHNBR model molecules and (MO)₁₂ clusters (M = Zn, Mg, Ca). Covalent bond occupancy (occ), atomic contribution of the covalent bond (in %) and second-order interaction energies (ΔE⁽²⁾), in kcal mol⁻¹. Bond nomenclature from Fig. 5

Group	Bond		(ZnO)12	(MgO)12	(CaO)12	(MgO2)6
Nitrile	2	occ	1.972	—	—	—
		% M – % N	6.5–93.5%	—	—	—
		ΔE^(2)	—	18.9 (Nsp → Mg4s)	8.13 (Nsp → Ca5s)	23.83 (Nsp → Mg4s)
Carboxyl	1	occ	—	—	—	1.979
		% O – % H	—	—	—	83.24–16.76%
		ΔE^(2)	10.46 (Osp → σ*(O–H))	42.66 (O2p → σ*(O–H))	1.02 (O2p → σ*(O–H)	—
	2	occ	1.993	1.998	1.999	—
		% O – % H	82.2–17.8%	80.7–19.3%	75.8–24.16%	—
		ΔE^(2)	—	—	—	92.47 (Osp^2 → σ*(O–H)
	3	occ	—	—	—	—
		% O – % M	—	—	—	—
		ΔE^(2)	17.1 (Osp → Zn4p)	12.6 (Osp^2 → Mg4s)	6.19 (Osp → Ca5s)	14.23 (Osp → Mg4s)

---

Table 3 QTAIM results, in a.u., for the interaction between characteristic XHNBR model molecules and (MO)₁₂ clusters (M = Zn, Mg, Ca). Electron density (ρ), its Laplacian (∇²ρ) and total electron energy density (H). Bond nomenclature from Fig. 5

Group	Bond		(ZnO)12	(MgO)12	(CaO)12	(MgO2)6
Nitrile	2	ρ	0.066	0.025	0.021	0.030
		∇^2ρ	0.250	0.173	0.104	0.214
		H	−0.021	0.008	0.004	0.010
Carboxyl	1	ρ	0.072	0.062	—	0.260
		∇^2ρ	0.163	0.159	—	−1.196
		H	−0.014	−0.007	—	−0.380
	2	ρ	0.291	0.302	0.357	0.098
		∇^2ρ	−1.538	−1.640	−2.006	0.110
		H	−0.464	−0.489	−0.584	−0.044
	3	ρ	0.093	0.040	0.032	0.041
		∇^2ρ	0.392	0.355	0.168	0.324
		H	−0.037	0.015	0.004	0.014

---

QTAIM analysis corroborates the trends observed in the geometry and NBO results, providing a clear picture of the balance between electrostatic and covalent contributions. For the nitrile group, ZnO exhibits negative total energy density (H) at the bond critical point (BCP), confirming partial covalency, whereas MgO, CaO, and MgO₂ show positive H and ∇²ρ, characteristic of closed-shell electrostatic interactions. The particularly low electron density (ρ) and large Laplacian (∇²ρ) values for CaO reflect its weaker and more long-range coordination to the nitrile nitrogen, consistent with the larger cation radius and lower polarizing power compared to Mg and Zn. For the carboxyl group, QTAIM confirms the correlation with proton transfer observed in the geometries: when proton transfer occurs (MgO, ZnO, CaO), bond 1 is predominantly electrostatic and bond 2 covalent, reflecting the formation of a stable M–O bond with the deprotonated carboxylate. In MgO₂, the absence of proton transfer leads to bond 1 adopting covalent character and bond 2 remaining electrostatic, in agreement with the retained proton on the carboxyl group and the unequal C–O bond lengths (R₄ and R₅). Bond 3 shows partial covalency only for ZnO, while in the remaining systems it is purely electrostatic, consistent with the lower polarizability of Mg, Ca, and the peroxide surface.

3.2.2 Interactions of XHNBR monomers and dimers with MO/MO₂ cluster models. General interaction distance trends remain consistent with those observed for the fragment models (See Tables 1 and 2 in the SI). In several optimized structures, some interatomic distances become significantly elongated as a consequence of the molecular flexibility, particularly in the dimeric systems, where only a portion of the molecule effectively interacts with the cluster surface.

When comparing different molecular models, monomeric systems consistently exhibit greater stabilization than their respective fragments, with the largest increase observed for the oxides and a more modest effect for the peroxide system. This suggests that the additional degrees of freedom in the monomer allow for more favorable binding conformations. In the dimer systems, interaction energies remain slightly more stabilizing than in the monomers, but the effect is less pronounced. In the case of the peroxide system, the dimer interaction energy is even less stabilizing than that of the monomer, indicating that conformational constraints limit the extent of interaction with the cluster. Given the small size of the nanoparticles, only a portion of the dimer is able to effectively interact with the surface, reducing the overall stabilization. These findings highlight the influence of molecular flexibility and steric effects on binding strength, reinforcing the role of system conformation in determining the observed interaction energies.

3.2.3 Interactions of XHNBR functional models with MO/MO₂ surface models. The analysis of the optimized geometries at the XHNBR/MO interface reveals fundamental differences in interaction behavior compared to the cluster models. While the same bonding motifs observed in the finite systems are generally retained, the crystalline periodicity and rigidity of the surfaces substantially limit their structural flexibility, preventing local reorganization upon adsorption (Fig. 6). As a result, the metal oxide surfaces exhibit more localized and less adaptable interactions with the functional groups.

Fig. 6 Optimized structures of (ZnO)_x surface and different models of XHNBR, obtained from VASP calculations (see Section 2.1). (a) Carboxyl group and (b) nitrile group.

In the nitrile group, the C≡N triple bond remains structurally intact (R₁ = 1.16 Å, see Table 1 lower section), consistent with the cluster-level observations. The metal–nitrogen distances (R₂) fall within a narrow range of 2.1–2.6 Å across all surfaces, confirming the formation of a stable coordination bond. This uniformity indicates that the nitrile group can effectively adapt to the more rigid environment of the crystalline oxide surfaces, establishing localized and well-defined interactions even under restricted relaxation conditions. This structural trend correlates with the gradual weakening of the interaction enthalpy, which becomes less exothermic along the same sequence (from −50.88 kcal mol⁻¹ to −40.89 kcal mol⁻¹).

By contrast, the carboxyl group exhibits a reduced capacity for multidentate coordination and proton transfer. While such processes were energetically favorable in the cluster models, they are largely suppressed on the periodic surfaces due to the limited structural flexibility of the lattice. For all analyzed periodic surfaces except (MgO₂)_x, a clear proton transfer to the surface oxygen atom is observed, consistent with the short R₂ distances (typically <1.2 Å) characteristic of newly formed O–H bonds at the oxide surface. The resulting R₄ and R₅ distances for MgO, ZnO, and CaO reflect a partially delocalized carboxylate motif. The broader range of M–O distances (R₃–R₅ = 1.25–2.39 Å), reflects the different coordination environments at the surface. The CaO surface, despite having the longest metal–oxygen bonds (up to 2.39 Å), achieves one of the strongest interaction enthalpies (−46.20 kcal mol⁻¹), suggesting that electrostatic stabilization compensates for the weaker geometric overlap. Conversely, MgO₂ exhibits the weakest binding (−25.60 kcal mol⁻¹), in line with its lower proton affinity and reduced capacity to form stable surface interactions.

Unlike the cluster models, the electronic structure analyses based on NBO and QTAIM could not be applied to the periodic surfaces due to methodological limitations. In periodic calculations, NBO analysis is not compatible with PAW pseudopotentials and therefore cannot be performed. QTAIM can, in principle, be applied to periodic systems, but the topology analysis differs substantially from the molecular case and does not provide the same type or quality of information we obtain from cluster models. For these reasons, NBO and QTAIM analyses were not carried out for the periodic surfaces.

3.3 Interaction of PA6 with MO/MO₂ cluster models and XHNBR dimers

Following the same methodological approach applied to the XHNBR-MO systems, the interaction between PA6 and the selected metal oxides (ZnO, MgO, CaO) and peroxide (MgO₂) clusters was investigated. Periodic models were not further explored, since the cluster-based approach was previously shown to provide an accurate local description of the metal–oxygen coordination environment while enabling detailed electronic characterization through NBO and QTAIM analyses.

This section is organized as follows: in Section 3.3.1, the interaction between PA6 molecular fragments and the oxide clusters is examined, focusing on the optimized geometries, interaction enthalpies, and the nature of the interfacial interactions. In Section 3.3.2, the interaction energy of the PA6-XHNBR interface is discussed, emphasizing the influence of the clustering methodology in accurately describing the conformational and energetic diversity of the interface.

3.3.1 Interactions of PA6 functional model with MO/MO₂ cluster models. To elucidate the interfacial behavior of polyamide-6 (PA6) with metal oxide clusters, model interactions were investigated using the amide functional group as the representative structural unit of the polymer backbone. Two distinct conformations were considered: amide(–N), where the metal interacts with the amide nitrogen lone pair, and amide(–O), where coordination occurs through the carbonyl oxygen. Representative optimized geometries of the PA6 amide models interacting with (ZnO)₁₂ are shown in Fig. 7. These structures illustrate the preferred binding modes of the amide functional group with the metal cluster, serving as reference for comparison among different metal oxides and peroxides. The analysis of these models enables the evaluation of how the oxide composition influences both the interaction strength and the electronic character of the interface. The optimized geometries and interaction enthalpies (ΔH) are summarized in Table 4. Selected bonds and interactions have been used for the discussion and labeled according to Fig. 8.

Fig. 7 Optimized structures of (ZnO)₁₂ and the two coordination modes of the amide group in PA, obtained at the ωB97XD/6-31+G(d,p) level of theory using Gaussian 16 (see Section 2.1). (a) Amide(–N) and (b) amide(–O).

Table 4 Bond distances (R₁ and R₂, in Å) of the optimized structures, following the nomenclature defined in Fig. 8, and corresponding interaction enthalpies (ΔH, in kcal mol⁻¹) between the two conformations of PA6, amide(–N) and amide(–O), with (MO)₁₂ (M = Zn, Mg, Ca) and (MgO₂)₆ clusters. Results obtained at the ωB97XD/6-311++G(2df,2p)//ωB97XD/6-31+G(d,p) level of theory, using Gaussian 16 (see Section 2.1)

	Amide(–N)		Amide(–O)
	R1	ΔH	R2	ΔH
(ZnO)12	2.182	−17.07	2.034	−23.72
(MgO)12	2.305	−16.64	2.064	−23.83
(CaO)12	2.461	−27.12	2.418	−20.77
(MgO2)6	2.199	−23.74	1.986	−33.98

---

Fig. 8 Schematic representation of the interaction and the relevant bonds between the PA6 amide model and the (MO)₁₂ clusters.

The interaction of PA6 with metal oxide clusters reveals two distinct coordination modes through the amide functional group: amide(–N) and amide(–O). The results summarized in Table 4 show that both coordination strength and geometry depend on the nature of the metal center and the coordination site involved. In all cases except for (CaO)₁₂, the amide(–O) configuration is more stable than the amide(–N) one, with interaction enthalpies between −20.77 and −33.98 kcal mol⁻¹. This trend reflects the intrinsic electronic structure of the amide group: the carbonyl oxygen carries a significant partial negative charge and possesses a localized lone pair, making it a stronger Lewis base than the partially delocalized nitrogen lone pair. Consequently, metal cations such as Zn²⁺ and Mg²⁺ preferentially coordinate to the carbonyl oxygen, forming M–O bonds with shorter distances (≈2.0 Å) and larger stabilization energies. The resulting interactions are largely electrostatic in nature, with minor charge-transfer contributions, consistent with the moderate Lewis acidity of ZnO and MgO surfaces.

In contrast, the (CaO)₁₂ cluster exhibits an inverted trend, showing a notably higher stabilization for the amide(–N) configuration (ΔH = −27.12 kcal mol⁻¹) compared to amide(–O) (ΔH = −20.77 kcal mol⁻¹). Structural analysis reveals that this stronger interaction arises from proton transfer from the amide N–H group to a basic surface oxygen atom (R(O–H) = 0.975 Å). The resulting partial deprotonation of the amide generates an anionic nitrogen that coordinates strongly to Ca²⁺, forming an ionic Ca–N bond. This process is driven by the high basicity and proton affinity of (CaO)₁₂ and is consistent with the previous observations for amine and carboxylate models of XHNBR, where CaO surfaces favor proton abstraction and strong ionic coordination.

For the (MgO₂)₆ cluster, both orientations form stable complexes, but the interaction strengths follow the opposite trend. The amide(–O) configuration displays the strongest stabilization (ΔH = −33.98 kcal mol⁻¹), whereas the amide(–N) interaction is significantly weaker (ΔH = −23.74 kcal mol⁻¹). In contrast to the behavior observed for CaO, no proton transfer from the amide N–H group is detected in either configuration.

To further characterize the nature of the interactions at the metal oxide cluster/XHNBR interface, NBO (Table 5) and QTAIM (Table 6) analyses were performed. The NBO results indicate that the main donor orbitals correspond to the amide N_2p and O_sp lone pairs, which interact with empty metal orbitals (Mg_4s, Ca_5s and Zn_4p). The second-order stabilization energies, ΔE⁽²⁾, reveal stronger donor–acceptor interactions for (ZnO)₁₂ (20.89 kcal mol⁻¹ for amide(–N)) and (MgO₂)₆ (17.64 kcal mol⁻¹ for amide(–O)). (CaO)₁₂ and (MgO)₁₂ exhibit weaker ΔE⁽²⁾ values (5.71–15.39 kcal mol⁻¹), reflecting interactions dominated by ionic rather than covalent contributions.

Table 5 NBO results for the interactions between characteristic PA6 model molecules and (MO)₁₂ clusters. Covalent bond occupancy (occ), atomic contribution of the covalent bond (in %) and second-order interaction energies (ΔE⁽²⁾), in kcal mol⁻¹. Bond nomenclature from Fig. 5

Group	Bond		(ZnO)12	(MgO)12	(CaO)12	(MgO2)6
Amide(–N)	1	occ	—	—	—	—
		% M – % N	—	—	—	—
		ΔE^(2)	20.89 (N2p → Zn4p)	15.39 (N2p → Mg4s)	13.23 (Nsp^2 → Ca5s)	17.46 (N2p → Mg4s)
Amide(–O)	2	occ	1.841	—	—	—
		% M – % O	1.2–98.8%	—	—	—
		ΔE^(2)	—	13.05 (Osp → Mg4s)	5.71 (Osp → Ca5s)	17.64 (Osp → Mg4s)

---

Table 6 QTAIM results, in a.u., for the interaction between characteristic PA6 model molecules and (MO)₁₂ clusters (M = Zn, Mg, Ca). Electron density (ρ), its Laplacian (∇²ρ) and total electron energy density (H). Bond nomenclature from Fig. 5

Group	Bond		(ZnO)12	(MgO)12	(CaO)12	(MgO2)6
Amide(–N)	1	ρ	0.060	0.024	0.034	0.032
		∇^2ρ	0.200	0.135	0.158	0.194
		H	−0.017	0.005	0.003	0.006
Amide(–O)	2	ρ	0.068	0.033	0.028	0.040
		∇^2ρ	0.293	0.251	0.153	0.337
		H	−0.020	0.011	0.005	0.015

---

The QTAIM analysis further confirms the trends observed in the NBO results. For the amide(–N) interaction mode, (ZnO)₁₂ displays the highest electron density at the bond critical point (ρ = 0.060 a.u.) and the only negative total energy density (H = −0.017 a.u), which indicates a partially covalent character. In contrast, (MgO)₁₂, (CaO)₁₂, and (MgO₂)₆ present lower electron densities and positive H values, reflecting closed-shell, predominantly ionic interactions. In the amide(–O), (ZnO)₁₂ again shows the strongest interaction, with the highest electron density (ρ = 0.068 a.u.) and a negative total energy density (−0.020 a.u). (MgO₂)₆ also shows a relatively high electron density (ρ = 0.040 a.u), although its positive H value indicates an interaction that remains ionic. (MgO)₁₂ and (CaO)₁₂ exhibit smaller ρ values and positive H, consistent with weaker, electrostatic interactions.

Overall, the NBO and QTAIM analyses reveal that the stronger interactions in the amide(–O) coordination arise from enhanced orbital overlap and larger electron density at the BCP, whereas amide(–N) coordination in (CaO)₁₂ remains weak and predominantly ionic, even after proton transfer. Thus, the bonding nature across these systems spans from partially covalent only in the case of (ZnO)₁₂, while (MgO)₁₂, (MgO₂)₆, and (CaO)₁₂ display fully ionic interactions, in direct correspondence with the oxides' increasing basicity and decreasing polarizing power, and consistent with the trends observed previously for XHNBR systems.

3.3.2 Interactions of PA6 dimer with XHNBR dimer. The conformational flexibility of both PA6 and XHNBR leads to a highly multidimensional configurational landscape, making the evaluation of interfacial interactions particularly sensitive to the choice of representative structures. To obtain statistically meaningful interaction energetics, we applied the conformational sampling and clustering protocol detailed in Section 2.2. This strategy combines extensive MD and CREST-based sampling with K-means clustering, ensuring that the most populated regions of the conformational space are adequately represented while avoiding redundancy. For each subsystem (PA6 dimer, XHNBR dimer, and the combined PA6-XHNBR complex), ten representative conformers were selected and subsequently optimized at the DFT level to compute accurate Boltzmann-weighted energetics (See Fig. S1–S3 in SI).

Fig. 9a and b illustrate the distribution of the sampled XHNBR dimer conformers projected onto the first two principal components. In Fig. 9a, each point is color-coded according to its cluster assignment, revealing the structural diversity captured by the descriptor space. Fig. 9b shows the corresponding Gaussian kernel density map, which highlights the most populated regions of the conformational ensemble. The cluster centroids marked with crosses correspond to the conformers selected for DFT refinement.

Fig. 9 Clustering and Boltzmann-weighted analysis of XHNBR dimer conformers. (a) PCA projection of all sampled conformers colored by cluster. (b) Kernel density map highlighting the most populated regions. (c) Relative DFT energies and cluster centroids. (d) Boltzmann populations of the ten DFT-refined conformers.

From a chemical standpoint, the PA6-XHNBR interface is expected to be strongly influenced by H-bonding interactions involving the PA6 amide groups and the polar functionalities of XHNBR. While no explicit hydrogen-bond analysis was performed for this system, such interactions are consistent with the chemical nature of the interacting groups and are expected to contribute significantly to the stabilization of the most populated conformers identified in the clustering analysis.

The relative DFT energies of the optimized XHNBR conformers and their corresponding Boltzmann populations are presented in Fig. 9c. These results clearly show that only a small subset of conformers contributes significantly to the thermodynamic ensemble. The Boltzmann distribution plotted in Fig. 9d reinforces this observation, demonstrating that a few low-energy minima dominate the equilibrium population. This confirms the need to account for both conformational diversity and statistical weighting when evaluating interfacial adhesion.

Using the Boltzmann-weighted energies of the isolated subsystems and the combined PA6-XHNBR complex, the interaction enthalpy was computed as: H =H(PA6-XHNBR) – H(PA6)-H(XHNBR). Using the weighted energies obtained from the ten DFT-refined conformers are, the Boltzmann-weighted interaction enthalpy is ΔH = −23.97 kcal mol⁻¹. Importantly, this value arises from a conformational ensemble rather than a single structure, highlighting the necessity of statistically weighted approaches for polymer–polymer interfaces. Analyses based solely on a single optimized structure would significantly underestimate or overestimate the adhesion energy, depending on the chosen geometry.

4 Conclusions

For XHNBR, the carboxyl moieties are identified as the dominant contributors to adhesion, exhibiting the strongest and most complex interactions with the oxide surfaces. These involve proton transfer, metal–oxygen coordination, and charge delocalization, leading to stabilized carboxylate-metal motifs characterized by short O–H contacts and significant π-electron delocalization. In contrast, nitrile groups form weaker, predominantly electrostatic interactions through coordination of their nitrogen lone pair to surface metal atoms. The interaction strength follows the order (CaO)₁₂ > (MgO)₁₂ > (ZnO)₁₂ > (MgO₂)₆, reflecting the increasing basicity and ionic character of the oxides. While ZnO promotes partial covalency, CaO forms largely ionic but stronger bonds, making it a promising, less toxic alternative to ZnO for crosslinking in XHNBR-based systems. At the periodic surface level, the same general binding motifs are maintained; however, surface rigidity limits interfacial reorganization and reduces overall stabilization. This restriction weakens the ability of carboxyl groups to form multiple contacts, shifting the dominant contribution toward nitrile interactions on crystalline surfaces.

The interaction of PA6 with MO clusters is governed by the dual coordination capability of the amide group, whose binding preference depends strongly on the electronic structure and basicity of the oxide surface. Across ZnO, MgO and MgO₂, the amide(–O) orientation generally provides stronger stabilization, driven by the localized carbonyl lone pair and its ability to form stronger contacts with metal centers. Only in CaO does this trend invert: the high basicity of the CaO surface promotes proton transfer from the amide N–H group, generating an anionic nitrogen that binds ionically to Ca²⁺ and shifting the preferred mode to amide(–N).

NBO and QTAIM analyses reveal a clear progression from partially covalent interactions only in ZnO to fully ionic bonding in MgO, MgO₂ and CaO, in line with the increasing basicity and decreasing polarizing power of the corresponding oxides. These results show that, although PA6 contains fewer functional groups than XHNBR, its amide units still establish meaningful and oxide-dependent interfacial interactions, with oxygen coordination dominating except in surfaces capable of driving proton transfer.

The conformational clustering and Boltzmann-weighted sampling methodology used here enables a realistic description of these flexible interfaces, capturing the contributions of multiple low-energy configurations that govern the overall interfacial stabilization.

Overall, these results provide molecular-level guidelines for tailoring interfacial chemistry in hybrid thermoplastic elastomers. By selecting oxides with appropriate basicity (e.g., CaO or MgO) and optimizing the carboxyl content of XHNBR, interfacial adhesion, crosslink density, and damping performance can be enhanced while improving environmental compatibility relative to conventional ZnO-based systems.

Author contributions

Conceptualization: JJE, AB, JMM and FR. Methodology: UC, JMM and FR. Investigation: UC. Formal analysis: UC. Visualization: UC. Supervision: JJE, AB, JMM and FR. Funding acquisition: JMM and FR. Writing – original draft: UC and FR. Writing – review and editing: All authors.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: Fig. S1–S3; Tables S1 and S2, and cartesian coordinates. See DOI: https://doi.org/10.1039/d5ra09279e.

Acknowledgements

Technical and human support provided by IZO-SGI, SGIker (UPV/EHU, MICINN, GV/EJ, ERDF and ESF) is gratefully acknowledged for assistance and generous allocation of computational resources. We acknowledge support from grant PID2023-148587NB-I00, funded by MICIU/AEI/10.13039/501100011033 and by ERDF/EU, from the María de Maeztu Excellence Unit CEX2023-001303-M funded by MCIN/AEI/10.13039/501100011033, and by the European Union NextGenerationEU/PRTR-C17.I1. This work was also supported by the IKUR HPC&AI Strategy, an initiative of the Basque Government aimed at developing high-performance computing and artificial intelligence research. UC acknowledge a PhD fellow under this strategy through the collaboration agreement with Ikerbasque.

Notes and references

1. Y. H. Elmoghazy, B. Safaei, M. Asmael, S. Sahmani, Q. Zeeshan and Z. Qin, Arch. Comput. Methods Eng., 2024, 31, 2611–2662.
2. A. Burgoa, B. Lekube and K. Zulueta, J. Appl. Polym. Sci., 2022, 139, e51748.
3. K. Zulueta, A. Burgoa, B. Lekube, J. L. Vilas and A. Arrillaga, J. Appl. Polym. Sci., 2022, 139, e53026.
4. G. Holden, in Thermoplastic Elastomers, Elsevier, 2024, pp. 97–113.
5. A. Reghunadhan, H. Akhina, A. Ajitha, N. Chandran, S. T. Nair, H. J. Maria and S. Thomas, in Thermoplastic Elastomers (TPEs) from Rubber-Plastic Blends, Elsevier, 2024, pp. 291–314.
6. H. Dong and Y. Zhang, Compos. Appl. Sci. Manuf., 2023, 175, 107783.
7. H. P. Brown, Rubber Chem. Technol., 1963, 36, 931–962.
8. X. H. Li, B. Kretzschmar, A. Janke, L. Häussler, K. Schneider and M. Stamm, Advanced Materials and Processes II, 2012, pp. 272–276.
9. Z. Starý, T. Pemsel, J. Baldrian and H. Münstedt, Polymer, 2012, 53, 1881–1889.
10. I. Y. Gotlib, D. V. Filyukov and S. W. de Leeuw, Polym. Sci., 2009, 51, 583–592.
11. E. Kubyshkina, M. Unge and B. L. G. Jonsson, J. Chem. Phys., 2017, 146, 051101.
12. E. Perrin, M. Schoen, F.-X. Coudert and A. Boutin, J. Phys. Chem. B, 2018, 122, 4573–4582.
13. L. Ibarra and M. Alzorriz, J. Appl. Polym. Sci., 2002, 84, 605–615.
14. A. Burgoa, R. Hernandez and J. L. Vilas, Polym. Int., 2020, 69, 467–475.
15. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, Wallingford CT, 2016, vol. 2.
16. G. Kresse and J. Hafner, Phys. Rev. B, 1993, 47, 558–561.
17. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
18. G. Kresse and J. Furthmüller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
19. G. Kresse and D. Joubert, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
20. A. E. Reed, R. B. Weinstock and F. Weinhold, J. Chem. Phys., 1985, 83, 735–746.
21. A. E. Reed and F. Weinhold, J. Chem. Phys., 1985, 83, 1736.
22. F. Weinhold, J. Chem. Educ., 1999, 76, 1141.
23. E. D. Glendening, C. R. Landis and F. Weinhold, WIREs Comput. Mol. Sci., 2012, 2, 1–42.
24. R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press: New York, USA, 1994.
25. R. F. W. Bader, J. Phys. Chem. A, 1998, 102, 7314–7323.
26. S. Jenkins and I. Morrison, Chem. Phys. Lett., 2000, 317, 97–102.
27. W. D. Arnold and E. Oldfield, J. Am. Chem. Soc., 2000, 122, 12835–12841.
28. T. A. Keith, AIMAll (Version 19.10.12), TK Gristmill Software, Overland Park KS, USA, https://aim.tkgristmill.com, 2019.
29. P. Popelier, Atoms in Molecules: an Introduction, Prentice Hall, 2000.
30. P. Pracht, S. Grimme, C. Bannwarth, F. Bohle, S. Ehlert, G. Feldmann, J. Gorges, M. Müller, T. Neudecker, C. Plett, S. Spicher, P. Steinbach, P. A. Wesołowski and F. Zeller, J. Chem. Phys., 2024, 160, 114110.
31. P. J. Lisboa, T. A. Etchells, I. H. Jarman and S. J. Chambers, BMC Bioinf., 2013, 14, year.
32. F. L. Gewers, G. R. Ferreira, H. F. D. Arruda, F. N. Silva, C. H. Comin, D. R. Amancio and L. D. F. Costa, ACM Comput. Surv., 2021, 54, 1–34.
33. I. T. Jolliffe and J. Cadima, Philos. Trans. R. Soc. A Math. Phys. Eng. Sci., 2016, 374, 20150202.
34. J. L. Phillips, M. E. Colvin and S. Newsam, BMC Bioinf., 2011, 12, 445.
35. H. Zhou, L. Zhang, S. Li and G. Lu, Phys. Chem. Chem. Phys., 2024, 26, 2269–2276.
36. A. Baillard, Y. Wang, L. Zhang and G. Lu, Surf. Sci. Rep., 2024, 79, 100540.
37. Akrochem Technical Paper, Magnesium Oxide Solutions for Rubber Compounding, https://www.akrochem.com/pdf/technical_papers/magnesium_oxide_solutions1.pdf, 2020, Accessed October 2025.

---