Combined DFT and MD Simulation Protocol to Characterize Self-Healing Properties in Disulfide-Containing Materials: Polyurethanes and Polymethacrylates as Case Studies

The introduction of dynamic bonds in polymeric materials facilitates the emergence of new functionalities, such as self-healing capacity. Understanding the role of the molecular structure in the efficiency of the self-healing process is fundamental to design new materials with improved features. Computational chemistry has emerged as a valuable tool for the characterization of polymeric materials. In this work, computational chemistry is used to analyze the observed self-healing capacity of a set of disulfide-containing polyurethanes and polymethacrylates, including different hard segments and dynamic bonds. For this purpose, a recently developed theoretical protocol has been used. This protocol is based on three parameters: the probability of generating radicals by cleavage of the disulfide bond, the energetic barrier of the exchange reaction among disulfides and the dynamics of the polymeric chains. This protocol is able to qualitatively explain the experimental self-healing properties of these materials. In particular, it explains both the great performance of two materials and the lack of self-healing capacity of another two. Besides, it can also describe the improvement of the self-healing capacity with increasing temperature. These results demonstrate the robustness and usefulness of this approach for the analysis and prediction of self-healing properties in polymeric materials. Therefore, this protocol allows to predict new materials with improved properties and will help the experimental community in the development of these improved materials.


INTRODUCTION
The use of dynamic covalent chemistry, introduced by Rowan et al. (2002), has emerged as a powerful tool in the design of new materials in different research fields, ranging from biochemistry to polymer chemistry. It is based on a reversible cleavage and reformation of certain covalent bonds under thermodynamic control, and its use in the design of new materials may introduce new functionalities which could assist creating materials with improved features. This improved features may vary from solid state recycling, property regulation, controllable degradation and reprocessing to self-healing capacity (Zhang et al., 2018). The dynamic character can be introduced by means of non-covalent interactions like hydrogen bonds (Cordier et al., 2008), π−π stacking (Burattini et al., 2009) or metal-ion interactions (Burnworth et al., 2011), but in the last few years, the use of dynamic covalent bonds has increased since, in addition to ensuring dynamic behavior, it also provides robustness to materials (Hager et al., 2010;Billiet et al., 2013;Hillewaere and Du Prez, 2015). Polymeric materials including dynamic bonds as crosslinks in their network have been labelled as covalent adaptive networks (CANs) and the activation of the bond-exchange reactions, which could be triggered by any external stimuli, can induce changes in their topology (Kloxin et al., 2010;Bowman and Kloxin, 2012;Kloxin and Bowman, 2012). Examples of CANs are the Diels-Alder reaction (Chen et al., 2002(Chen et al., , 2003Liu and Chen, 2007), boronic esters (Cash et al., 2015;Cromwell et al., 2015), olefin metathesis (Lu and Guan, 2012), polysiloxanes (Cho et al., 2012;Schmolke et al., 2015) or disulfides (Sarma et al., 2007;Belenguer et al., 2011;Matxain et al., 2016;Nevejans et al., 2016;Rekondo et al., 2014;Martin et al., 2014;Aguirresarobe et al., 2017). The use of disulfides as crosslinks, especially the aromatic ones Martin et al., 2014;Azcune and Odriozola, 2016;Takahashi et al., 2016), has been revealed as a useful technique to generate responsive materials due to the dynamics and reversibility of this bond, which can be activated without external stimuli (Canadell et al., 2011;Lafont et al., 2012;Amamoto et al., 2012).
One of the most studied polymers in this field are the polyurethanes (PU), because of its great versatility and the ease of functionalization with many chemical moieties (Engels et al., 2013). As a consequence, PUs are very promising class of materials to be used in self-healing applications due to the range of dynamic bonds that can be incorporated (Aguirresarobe et al., 2021). The dynamic character can be introduced in different ways, for example, by urethane exchange (transcarbamoylation) (Fortman et al., 2015;Zheng et al., 2016Zheng et al., , 2017, urea exchange (Ying et al., 2014;Zhang et al., 2016;Erice et al., 2018), or by using dichalcogenide bonds as crosslinkers, which can be both aliphatic (Otsuka et al., 2010;Xu et al., 2016;Pan et al., 2017) and aromatic Martin et al., 2016;An X. et al., 2017;Shin et al., 2020;Eom et al., 2021). The mechanical properties of PUs are greatly affected by supramolecular interactions like hydrogen bonds and by the packaging of the hard segments (Brunette et al., 1982;Mattia and Painter, 2007;Zhang et al., 2014;Erice et al., 2019). The presence of dynamic bonds may also lead to the decrease of the mechanical strength and, therefore, a balance has to be found when designing new self-healing materials (Nevejans et al., 2019).
In this way, Kim et al. studied a set of PUs with different hard segments containing disulfide bonds, in the pursuit of a robust and easily processable polymer that shows self-healing features at room temperature (Kim et al., 2017). The hard segment in PUs is generated by the reaction of isocyanates with short chain diols, and provides strength, hardness and high temperature performance. In particular, these authors synthesized three PUs with different hard segments including aromatic disulfides (Figure 1), providing varying rigidity and labeled after the isocyanate used: IP (isophorone diisocyanate), an asymmetric alicyclic structure, HM [4,4′-methylenebis (cyclohexyl isocyanate)], a symmetric alicyclic structure and M [4,4′methylenebis (phenyl isocyanate)], an aromatic structure. The IP material was the only one to show self-healing at room temperature with mechanical properties such as robustness, stretchability and durability, surpassing those of roomtemperature self-healable materials to date. The HM material showed improved self-healing efficiency with increasing temperature, while the M material showed no self-healing ability even at high temperatures. This lack of self-healing capacity has been explained due to the packaging and the restricted mobility of the chains in this material.
FIGURE 1 | Hard segments of self-healing PUs including disulfide bonds proposed in Ref. Kim et al. (2017). The nomenclature is related to the isocyanate used to generate the hard segment: IP = isophorone diisocyanate; HM = 4,4′-methylenebis (cyclohexyl isocyanate) and M = 4,4′methylenebis (phenyl isocyanate). Another approach to modify or to provide with self-healing features is changing the dynamic bond. This approach has been used by Takahashi et al. (2017a), whose work focused on polymethacrylate networks including different disulfide bonds (An S. Y. et al., 2017;Liu et al., 2020;Rusayyis and Torkelson, 2020). In particular, these authors synthesized two materials: BTA, including a thermally exchangeable bis(2,2,6,6tetramethylpiperidin-1-yl)disulfide (BiTEMPS) crosslink (Takahashi et al., 2017b;Tsurumi et al., 2020;Yokochi et al., 2021;Takashima et al., 2020Takashima et al., , 2021Kataoka et al., 2021), and ADSA, including and alkyllic disulfide, Figure 2. The authors observed that BTA, with N atoms linked to the disulfide, presented healing capacity both at room and at high temperature, while ADSA, with aliphatic carbon chains linked to the disulfides, did not show healing capacity even at high temperature, even though the hardener structure was the same. This is an expected result, since it has been previously demonstrated that the presence of nitrogen atoms linked to the sulfur weakens the S-S bond, favoring the disulfide exchange reaction and promoting the self-healing behavior .
The simultaneous design of dynamic networks combined with reasonable mechanical properties is becoming the landmark of this field. Nevertheless, exploring the multiple polymer/dynamic chemistry combinations to obtain materials with the targeted properties is not trivial and only achievable by trial and error. Thus, computational chemistry appears as a very powerful predictive tool to guide experimental research. It was by means of theoretical chemistry that the good performance of the disulfide bond in the self-healing process was elucidated. Matxain et al. (2016) explored the reaction mechanism involved in aromatic disulfide-based self-healing materials, such as those developed by Rekondo et al. (2014). A (2 + 1) radical-mediated reaction mechanism was found to be responsible for the selfhealing reaction, instead of the expected (2 + 2) metathesis mechanism, and it was later verified experimentally in the absence of catalyst (Nevejans et al., 2016). The first step in the (2 + 1) radical-mediated mechanism is the cleavage of the S-S bond, creating sulfenyl radicals, which may further attack neighbouring disulfide bonds. This attack takes place via a three-membered transition state, which eventually leads to a new sulfenyl radical in a single-displacement reaction and another disulfide compound, and thus, to exchanging sulfur atoms in the process. Following this work, the authors were able to establish a theoretical protocol to assess the theoretical selfhealing capacity in these kind of materials, consisting on three different parameters (Formoso et al., 2017;Irigoyen et al., 2019b): 1) the probability to generate sulfenyl radicals, which can be controlled by tuning the strength of the S-S bond, 2) the energetic barrier of the exchange process, and 3) the dynamics of the polymeric chains, governing the ability to flow and reshuffle. In previous works (Matxain et al., 2016;Irigoyen et al., 2019a;Ruipérez et al., 2018;Irigoyen et al., 2019b), different aliphatic and aromatic disulfide-and diselenidebased materials have been studied by using this theoretical protocol. These works showed that the proper choice of the dynamic bond, the non-covalent interactions and the packaging are key features to enhance the self-healing capacity of a polymeric material.
Up to now, efforts have been made concerning the understanding of the dynamic character, at the molecular level, of dichalcogenide-containing materials. The main goal of this work is to test the validity of the proposed theoretical protocol to qualitatively predict the self-healing capacity of real systems. In order to do that, we have performed computational simulations of the materials previously mentioned, IP, HM, M, BTA and ADSA, using the theoretical protocol developed in our group (Formoso et al., 2017;Irigoyen et al., 2019b).

Quantum Chemical Calculations
All geometry optimizations were carried out in gas phase within density functional theory (DFT) (Hohenberg and Kohn, 1964;Kohn and Sham, 1965), by using the long-range corrected ωB97XD functional (Chai and Head-Gordon, 2008) combined with the 6-31+G(d,p) basis set (Hehre et al., 1972). Harmonic vibrational frequencies were obtained by analytical differentiation of gradients in order to determine whether the structures found were minima or transition states. The frequencies were then used to evaluate the zero-point vibrational energy (ZPVE) and the thermal (T = 298 K) vibrational corrections to the enthalpy (H) and Gibbs free energy (G) in the harmonic oscillator approximation. Single point calculations using the 6-311++G(2df,2p) basis set (Krishnan et al., 1980) were performed on the optimized structures to refine the electronic energy. All DFT calculations were carried out using the Gaussian 16 package (Frisch et al., 2016).
Ab initio Born-Oppenheimer molecular dynamics (QMD) calculations were carried out for the optimized species in order to determine their thermal stability, using the PBE functional 1997) combined with a DZP quality basis set and the RI formalism, with the corresponding auxiliary basis sets (Deglmann et al., 2004;Eichkorn et al., 1995;Sierka et al., 2003). The calculations were carried out at 298 K by means of the Nose-Hoover thermostat. All these simulations were as long as 40,000 a.u. (9.651 ps), with a time step of 40 a.u. (1.93 ps), and were performed using the TURBOMOLE package (TURBOMOLE Gmbh, 2014).

Molecular Dynamics Simulations
The molecular dynamics (MD) simulations have been performed using GAFF ff14SB (Maier et al., 2015) amber force-field in the AMBER 14 molecular dynamics simulation package (Case et al., 2014). All structures were built via the LEaP module of Ambertools and the charges were computed using the restrained electrostatic potential (RESP) fitting procedure (Bayly et al., 1993). First, the ESP was calculated with the Gaussian 16 package using the 6-31G* basis set, at the Hartree-Fock level of theory, and then the RESP charges were obtained. All the simulations were carried out in vacuum in a canonical ensemble (NVT) with a 2 fs timestep. 100,000 steps of minimization (50,000 steps steep descent minimization plus 50,000 steps of conjugate gradient minimization) were followed by heating from 0 to 300 K over 200 ps. Covalent bond lengths involving hydrogen were constrained using the SHAKE algorithm.

Theoretical Self-Healing Capacity
The theoretical protocol to estimate the self-healing capacity of disulfide-based materials is based on three parameters (Formoso et al., 2017;Irigoyen et al., 2019b). The first parameter is the probability to generate sulfenyl radicals (ρ), which is governed by the dissociation of the S-S bond. This parameter is calculated from quantum molecular dynamics, and is defined as the number of simulation steps where the S-S bond distance is longer than 2.3 Å (N S−S ), a value for which the bond may be considered dissociated, divided by the total number of simulation steps, N tot .
The second parameter is the energetic barrier of the exchange reaction, which is used to estimate the rate constant (k), according to the Wigner, Eyring, Polanyi and Evans formulation of the Transition State Theory (Evans and Polanyi, 1935;Eyring, 1935): where R, k B and h are the ideal gas, Boltzmann and Planck constants, respectively, and ΔG is the activation energy corresponding to the attack of a sulfenyl radical to a disulfide bond. Finally, the mobility of the polymeric chains is estimated by analyzing the non-covalent interactions among chains, mainly hydrogen bonds. In order to do that, the parameter ω is defined, which is based on the distance between disulfides in the material, that will be governed by the interactions among chains: Using the radial distribution function of the sulfur atoms, three regions may be defined in order to calculate ω: the reacting region (I i ), where disulfides are close enough to undergo the exchange reaction (R ≤ 4.5 Å), the neighboring region (I ii ), where disulfides are far to react but with a non-negligible probability to approach the reacting region (4.5 < R < 20 Å), and the external region (I iii ), where the interaction among disulfides is neglected (R > 20 Å). Since the size of the external region is dependent on the size of the system considered, in order to have a sizeindependent parameter, this third region is not considered in the calculation of ω. The amount of sulfur atoms located in each region is calculated by integration of the radial distribution function within the limits defined above. These limits may change depending on the system analyzed.

System Models
Different model structures were used to carry out both quantum and molecular dynamics simulations. In quantum simulations small chain models were used, including the disulfide unit and the surrounding chemical structure. In particular, the model used for the PUs correspond to the hard segments depicted in Figure 1, while the corresponding to the polymethacrylate materials are represented in Figure 3. In the quantum chemical calculations, both one-and two-chain models were used.
The models used for the molecular dynamics simulations are built with two disulfide units including a small polymeric chain between them (Figure 4, where the model corresponding to the IP material is depicted). 16 chains are put randomly in a simulation box, following the procedure previously validated in our group (Irigoyen et al., 2019b), reproducing the experimental density of PUs and polymethacrylates after equilibration.

RESULTS AND DISCUSSION
In this section, the obtained results will be presented and discussed, in order to rationalize the experimental data obtained in the works of Kim et al. (2017) and Takahashi et al. (2017a), by using the theoretical protocol designed in our group. Concretely, this section will be divided into three subsections. First, we will summarize the experimental selfhealing properties of these materials. Then, we will make use of quantum chemical calculations to analyze the formation of sulfenyl radicals by cleavage of the disulfide bond. In order to do that, the bond dissociation energy (BDE) and the probability of radicals to be formed (ρ) will be evaluated. Also, in this subsection the exchange reaction barriers (ΔG) will be calculated. In the third subsection, by means of molecular dynamics simulations, the mobility of the polymeric chains will be analyzed by calculating the probability of sulfenyl radicals to find a neighboring disulfide bond (ω), which will be crucial to trigger the exchange reaction.

Self-Healing Properties of the Studied Materials
The polyurethanes synthesized by Kim et al. showed healing efficiency in the following order: IP > HM > M. Thus, IP fully recovered from surface scratch in 2 h, that is, this material showed 100% healing efficiency. The time needed to recover decreased with temperature and times of 30, 5, and 3 min were required at 40, 60, and 80°C, respectively. Besides, a sample of this material was cut into two pieces and then joined again at room temperature. After 2 h, the film resisted manual pulling and, after 4 h, a 5 kg load was lifted without tearing. M material showed no self-healing ability even at 80°C. Further evaluation of the self-healing properties was performed by tensile tests. Again, IP material showed the best efficiency, while HM material showed moderate mechanical recoveries at 25°C. The self-healing properties of the polymethacrylate networks studied by Takahashi et al. were tested by using strength and strain at break. First, a sample of BTA network was cut and the surfaces were put together at 120°C under mild pressure. The scar became invisible after 8 h and healing ratios of up to 85 and 92% for the strength and strain at break were recorded, respectively. Inferior recovery performance was observed for ADSA, and visible scars appeared. The difference in self-healing properties was also supported from fracture behavior in tensile tests. BTA fractured at different positions, similar to the pristine material, while ADSA fractured in the cut line.

Radical Formation and Disulfide Exchange
Several energetic and structural parameters have been considered in order to study the sulfenyl radicals formation and their influence in the disulfide exchange. The disulfide bond strength has been analyzed by calculating the BDE considering two different chain conformations, labeled as open and closed, in a single-chain model   ( Figure 5A). In addition to this, the interaction between two chains has been also calculated, to take into account the interchain noncovalent interactions, such as hydrogen bonding and π − π stacking among phenyl rings ( Figure 5B).

Bond Dissociation Energies
In Table 1 are collected the geometrical parameters obtained for the open and closed conformations, together with the BDEs. In the open conformation, the R-SS-R dihedral angles are large, more than 100°, which means that both chains linked to the sulfur atoms are rather far from each other, hindering the interaction among them, as depicted in Figure 5A. In the closed conformation, otherwise, the dihedral angle is small, of around 60°or less, which allows the intrachain interactions. See Figure 5B, where two chains interacting in a closed conformation are represented. In principle, a combination of different conformations would be present in the polymeric material due to the presence of multiple chains. However, in those materials including diphenyl disulfides (IP, HM, and M), the closed conformation seems to be the most favored one, since the interaction between the phenyl rings adjacent to the sulfur atoms helps reducing the steric repulsion with the rest of the chain. In the M and ADSA models both conformations can be obtained, while the steric repulsion in BTA prevents the formation of the closed conformation.
The calculated values in this work for BTA, a sulfenamide (33.77 kcal/mol), ADSA, an aliphatic disulfide (58.19 kcal/mol, open conformation, and 63.13 kcal/mol, closed conformation), and HM, an aromatic disulfide (47.34 kcal/mol), are in agreement with those values. Nevertheless, the BDEs calculated for IP and M (closed conformation) are notably larger, 59.31 and 65.65 kcal/mol, respectively. This suggests the presence of stronger intrachain interactions in this conformation (hydrogen bonding and π−π stacking). Note that the BDE is calculated as the energy difference between the whole molecule (disulfide) and the two separated sulfenyl fragments. This means that these non-covalent interactions are also included in the calculation of the BDE. Therefore, in the open conformation, where the non-covalent interactions are absent, the BDE values are expected to be lower and closer to the calculated in previous works. This effect is observed in the M model, for which the BDE in the open conformation is smaller (51.20 kcal/mol) than in the closed one (65.65 kcal/mol). In Table 1 is also observed that lower BDEs correspond to longer S-S bond distances.
Since intermolecular interactions have been found to be important in the self-healing capacity of materials, in this work the interchain interactions have been studied using a two-chain model ( Figure 5B). A deeper analysis will be carried out later by means of molecular dynamics simulations, including more and larger chains. In Table 2 are collected the geometrical parameters of each chain and the interaction energy between them. Inspecting the calculated interaction energies (ΔH), it is observed that in all cases interchain interactions are rather strong, around 30 kcal/mol. This suggests that the mobility of the chains may be limited.

Radical Formation
The BDE gives a measure of the strength of the bond and is relevant to estimate the probability to generate sulfenyl radicals (ρ, Eq. 1). Since the BDE is affected by chain conformation and interchain interactions, in order to estimate this probability, it is necessary to take into account the dynamics of many chains interacting at the same time in the material. Thus, ρ has been calculated for the five model systems by means of quantum molecular dynamics (QMD) simulations at different temperatures, those in which the experiments have been carried out, to see the influence of temperature in the formation of radicals. So that, four different temperatures have been used in the simulations (25, 40, 60, and 80°C) for the IP, HM and M systems, while two temperatures (25 and 120°C) have been used for the BTA and ADSA systems. In Table 3 are collected the results of these simulations.
Starting with the PUs, it can be seen that, at room temperature, IP is the only material presenting a value of ρ different from 0, which means that, theoretically, this would be the only one that could trigger the exchange reaction at that temperature due to the presence of radicals, in agreement with the experimental findings. At increasing temperatures, non-zero probabilities are calculated for the three PUs, but it is remarkable the almost negligible values calculated for the M material at any temperature, again consistent with the experimental results. This supports the assumption that the probability to break the S-S bond is not only dependent on the electronic features of the bond and its surroundings, but is also largely determined by the mobility of the chains, which in the case of the M material is notably reduced due to the packaging developed by the π-π interactions among phenyl rings. This has as a consequence a lower or even absent self-healing efficiency.
Inspecting the polymethacrylates, it can be seen that ρ for ADSA is zero at both room temperature and 120°C, meanwhile BTA presents high values at both temperatures, meaning that sulfenyl radicals are easily created and the exchange reaction promoted, in agreement with the self-healing efficiency observed experimentally in this material.
In Table 3, it can also be observed a non-regular trend for the variation of ρ with temperature. One could expect that increasing temperature would favor the dissociation of the bond. However, in some cases the opposite trend is observed. In Figure 6, the evolution of the disulfide bond distances during the simulation, used for the 2 | S-S bond length (R e ), in Å, R-SS-R dihedral angle (α), in degrees, and interaction energy (ΔH), in kcal/mol, between two polymeric chains. calculation of ρ, are depicted for both chains in a two-chain model of IP material, for which ρ shows a decreasing trend with increasing temperature. It is clearly seen that the disulfide of one of the chains shows a much larger variation of the S-S bond distance along the simulation (green line), while the other one (red line) remains with bond distances very similar to those calculated for the single-chain system. This can be related to the conformation adopted by the chain (open and closed), which, as it has been analyzed before, affects the bond distance (and the BDE). Thus, the red line would correspond to a chain with a closed conformation, where the π-π interactions between phenyl rings keep the sulfur atoms closer during the simulation, while the green line would correspond to a chain with an open conformation, providing more mobility to the sulfur atoms. This suggests that the bond distance and, therefore, ρ, may be affected by the initial conformation of the chains, that would be basically determined by both the starting conformation and the simulation length. In Table 3, ρ is extracted only from the closed conformation in order to relieve the computational effort. Nonetheless, several calculations were performed using both conformations as starting point, and it was observed that the value obtained for the closed one always corresponded to the smallest, therefore, we use it as the lower limit of this parameter. This feature hinders the observation of a regular trend in ρ, and could be partially solved by exploring a large number of conformations and using longer simulation times, which would be very costly computationally using quantum molecular dynamics simulations.

Barrier of the Exchange Reaction
Finally, the other parameter obtained by means of quantum chemical calculations is the exchange reaction barrier, calculated as the energy of the transition state, ΔG, which gives information about the kinetics of the exchange reaction (Eq. 2). In previous works (Matxain et al., 2016;Ruipérez et al., 2018;Irigoyen et al., 2019a), it was concluded that the chemical structure around the dynamic bond was a key feature in the calculation of the barrier. Since the calculation of the transition states is very demanding computationally and the systems studied in this work are rather large, we decided to compute only the reaction barrier for the IP material, obtaining a value of 12.02 kcal/mol, close to those previously calculated with a similar chemical structure around the disulfide bond, that is, an aromatic disulfide with a urethane unit attached to the ring in para position (12.21 kcal/mol) (Irigoyen et al., 2019a). So that, we can assume that the reaction barrier for all of IP, M, and HM will be similar, since all of them have the same chemical structure around the dynamic bond. Similarly, we can use our previous values calculated for sulfenamides (16-18 kcal/mol)  to estimate the approximate reaction barrier of BTA material. These results suggest that the exchange reaction barrier is not a relevant factor in the different self-healing performances observed in these materials.

Influence of the Hard Segment and Crosslinks in the Chain Mobility
The mobility of the chains is a key factor in the self-healing process, as it determines the ability of chains, and sulfur atoms in particular, to move and rearrange, which is essential for sulfur radicals to attack disulfide bonds and trigger the exchange reaction. This mobility is determined by the non-covalent interactions among chains which, as it has been demonstrated in the previous section, also affects the S-S bond dissociation process, the first step in the exchange mechanism.
Therefore, in this section, by means of classical molecular dynamics simulations, we will not only analyze the chain mobility by computing the ω parameter (Eq. 3), but also ρ will be reevaluated and compared with that obtained by QMD simulations. In classical simulations, many more and longer chains are included and, since the starting point is a random configuration, lots of different conformations are taken into account, removing the initial bias observed in the quantum simulations.  Inspecting the radial distribution functions of the sulfur atoms obtained from the molecular dynamics simulations (Figure 7), two different regions can be observed: one corresponding to the sulfur atoms within a disulfide bond, with a high, narrow peak at around the S-S equilibrium bond distance (~2.10 Å) and a large region, ranging from 4 to 25 Å, approximately (enlarged in the small box of Figure 7), which corresponds to the distances between sulfur atoms of different disulfide bonds. The parameter ω is obtained by integration of different regions of the radial distribution functions to estimate the number of disulfides in each region. In the same manner, a value for ρ can be defined by integrating the radial distribution function in the region of the S-S bond equilibrium distance (Figure 8), with the aim of obtaining a ratio of the sulfur atoms located away from the equilibrium distance, which could be considered as dissociated.

Reevaluation of ρ by Classical Molecular Dynamics
First, we will focus on the region of the radial distribution function corresponding to the sulfur atoms within each disulfide bond. In Figure 8A, the radial distribution functions of all the studied systems at 25°C are represented, while the region corresponding to larger bond distances are magnified in the right panel. It can be seen that both M (blue line) and BTA (black line) materials show the highest population of sulfur atoms at larger distances, while ADSA (orange line) shows the lowest one, with IP   Figure 8B). These results are a consequence of the chemical structure around the disulfide bond. BTA is a sulfenamide and, as it has been shown in the previous section, the N atoms linked to the sulfur produce weaker S-S bonds and larger bond distances ( Table 1).
The ADSA system consists on an aliphatic disulfide with stronger and, thus, shorter bonds. IP and HM have the same structure around the disulfide, that is, an aromatic ring attached to each sulfur atom, producing intermediate bond strengths and similar profiles. The exception to this trend corresponds to M, also an aromatic disulfide, but it presents the largest bond distances. In Table 1 can be observed that the M system, unlike IP and HM, may be found in both the open and closed conformation. The absence of the intrachain interactions in the open conformation makes the S-S bond dissociation easier and, thus, longer bond distances are observed.
In Figure 9A, the radial distribution functions for the IP system at different temperatures are represented. It can be observed that the higher the temperature, the lower the peak of the curve, as expected, since increasing the temperature leads to longer S-S bonds and the distribution is broader. In Figure 9B, an enlargement of the region corresponding to longer bond distances is represented and the same behavior can be observed, where higher temperatures show a larger amount of sulfur atoms at longer distances. This would correspond to a large probability for the disulfide bonds to be dissociated. This is the region of the radial distribution function that will be integrated in order to estimate the ratio of sulfur atoms laying at longer distances and, thus, with a higher probability to be dissociated. For that, a similar expression to that of ω parameter has been used, leading to the so-called ρ MD parameter (Eq. 4).
ρ MD I I I I + I II The regions selected for the integration are the following: I I ranges from 2.20 to 3.00 Å, where the sulfur atoms with higher probability to be dissociated are included, and I II , ranging from 1.80 to 3.00 Å. In this way, we obtain a proportion of the sulfur atoms that could be dissociated, similar to that obtained by using QMD simulations. The results are collected in Table 4. Now, it can be observed a more regular trend, where increasing temperatures corresponds to larger probabilities of the bond to be dissociated. This confirms that classical molecular dynamics simulations remove almost completely the limitations found in the quantum molecular dynamics calculations. Comparing Table 3 and Table 4, differences arise because ρ is calculated differently. In both cases, ρ is defined similarly in QMD and MD calculations, i.e., as the probability of S-S bond distance to be larger than a given threshold. In this way, using QMD we are able to study the formation and breaking of bonds, while using MD it is not possible. Thus, in the MD calculations, the S-S distance is defined by the contributions included in the force field, while in the QMD calculations, the S-S distances change according to quantum effects. Hence, for considering the formation of radicals, QMD methods are necessary. According to these results, the values for M and ADSA are 0 since in the QMD simulations no S-S distances are found to be longer than the threshold, due to electronic quantum effects. In both cases, the 0 value arise directly from the calculation. MD results are used to study the temperature effect, but it is not possible to perform a direct comparison with ρ from QMD calculations, since MD is not able to describe formation and breaking of the bond.

Mobility of the Polymeric Chains
Next, an analysis of the region of the radial distribution functions corresponding to longer distances, those between different disulfide units, will be performed (see the enlarged region in Figure 7). This will allow us to understand the mobility of the sulfur atoms and their ability to reach neighboring disulfides in order to trigger the exchange reaction. To do that, the following  regions are defined and integrated to evaluate the ω parameter (Eq. 3): I I , the reacting region, ranging from 3 Å to 4.5Å, where sulfurs are close enough to react, and I II , the neighboring region, from 4.5 Å to 20 Å, where the disulfides are further away, but with a non-zero probability to enter in the reacting region.
The results, at different temperatures, are collected in Table 5. It can be observed that the temperature has small influence in ω, as it fluctuates around a certain value. Since the materials are solids, the variations in the molecular structure caused by temperature are small in this range of temperatures, as they are well below the T g of polyurethanes (120-130°C). So that, it makes sense to define an average value (ω ave ), which should be characteristic of each material at temperatures below the T g .
The differences observed in the PUs can be explained by inspecting their chemical structures. The highest ω ave is obtained for HM, which includes aliphatic cyclohexanes, while the lowest one is obtained for M, which contains phenyl groups (Figure 1). The aromatic groups favor more rigid structures due to their planarity and the π-π stacking interactions, while the cyclohexanes possess more degrees of freedom and, therefore, more mobility. The intermediate value for the IP system can be explained by a more limited mobility of the substituted cyclohexanes, still larger than that of the phenyl rings.
Regarding the polymethacrylates, a large difference is observed between BTA and ADSA, which can also be ascribed to their molecular structure. ADSA contains aliphatic hydrocarbon chains attached to the disulfide bond, which can be effectively packed and favoring the hydrogen bonding between urethane moieties. This hinders the mobility of the chains and a very low value of ω ave is found (0.0006). In BTA, instead, the presence of a cyclohexane next to the disulfide leads to a much less packed molecular structure and prevents the formation of hydrogen bonds between urethanes. Thus, a much larger value of ω ave is calculated (0.0506).

CONCLUSION
In this work, we have analyzed the theoretical self-healing capacity of several materials and compared it to the experimental results. In general, a good agreement is found and the theoretical calculations are able to explain qualitatively the experimental self-healing capacities.
In particular, it explains the great performance of IP and BTA, as well as the lack of self-healing capacity of M and ADSA. More interestingly, it can also describe the improvement of the self-healing capacity of HM with increasing temperatures. However, an accurate representation of the temperature effects in the dissociation of the disulfide bond is problematic with the present protocol, since the initial conformation in the quantum molecular dynamics used to calculate ρ may introduce a conformational bias. This issue may be overcome by using classical molecular dynamics, that allow including many more conformations, reducing at a great extent the initial bias. Nevertheless, it is not possible to describe the cleavage of the bond and, therefore, the formation of radicals, by classical dynamics. This means that this technique can only be used to estimate the influence of the temperature by inspecting the variation of the bond distance during the simulation.
The mobility of the polymeric chains has been estimated using the ω parameter. Large values of ω correspond to greater mobility of the chains. In PUs, the following trend has been calculated for the three hardeners: ω HM > ω IP > ω M . This suggests that, for HM and IP, the sulfenyl radicals will have more mobility to approach another disulfide and favoring the self-healing process, while for M the opposite will take place, in agreement with the experimental findings. For the polymethacrylates, BTA shows a large ω value, while that of ADSA is almost zero, again in agreement with the differences observed in the self-healing behaviour of the two materials. Besides, we have observed that the increment of the temperature does not really have an impact in ω, as the properties of the bulky materials are not affected that much when working at temperatures well below their T g . We can conclude, then, that the effect of temperature in the self-healing capacity is mainly related to a greater generation of radicals in the material. Finally, the estimation of the kinetic barriers of the exchange reaction suggests that this is not a key factor in the self-healing capacity of these materials.
As a summary, we believe that this theoretical protocol will assist in the development of new materials with improved selfhealing properties.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

ACKNOWLEDGMENTS
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.