Interfacial Binding Energy between Calcium-Silicate-Hydrates and Epoxy Resin: A Molecular Dynamics Study

Microcapsules encapsulated within epoxy as a curing agent have been successfully applied in self-healing materials, in which the healing performance significantly depends on the binding behaviour of the epoxy curing agent with the cement matrix. In this paper, the binding energy was investigated by molecular dynamics simulation, which could overcome the shortcomings of traditional microscopic experimental methods. In addition to the construction of different molecular models of epoxy, curing agents, and dilutants, seven models were established to investigate the effects of chain length, curing agent, and epoxy resin chain direction on the interfacial binding energy. The results showed that an increase of chain length exhibited had limited effect on the binding energy, while the curing agent and the direction of the epoxy significantly affected the interfacial binding energy. Among different factors, the curing agent tetrethylenepentamine exhibited the highest value of interfacial binding energy by an increment of 31.03 kcal/mol, indicating a better binding ability of the microcapsule core and the cement matrix. This study provides a microscopic insight into the interface behaviour between the microcapsule core and the cement matrix.


Introduction
As the most widely used building material worldwide [1], the safety and service life of concrete is important. Unfortunately, during its service life, microcracks or defects inevitably occur in concrete. Without proper control, this may lead to a serious impact on the structures, and further reduce the durability and shorten the service life of the structures [2,3]. To solve this issue, many scholars have attempted to repair microcracks by employing various methods [4][5][6][7]. Among them, according to the principle of bionics, self-healing concrete with microcapsules have been developed, which can be used to automatically locate and repair the microcracks inside concrete and further improve the durability of concrete structures [8]. Moreover, compared with traditional repairing technology, the microcapsule self-healing system exhibits the advantages of automatic repair, low cost, and renewability. In general, this self-healing system is designed based on three steps to achieve the repairing effect: (1) embedding the encapsulated healing agent in the cement matrix, (2) releasing the healing agent after the breaking of the container when the matrix produces microcracks, and (3) bonding the cementitious matrix with the cured healing agent to block the microcracks [9].
Currently, microcapsules with a urea-formaldehyde (UF) shell have been proven to be one of the most effective microcapsules for repairing concrete cracks [10,11]. According to research, the core of the microcapsules was prepared by a mixture of bisphenol A type epoxy resin (E-51), and diluent (n-butyl glycidyl ether, BGE) was used as the healing agent, the binding behaviour between the organic curing agent and the inorganic cementitious matrix in a microcapsule-based concrete with the assistance of MD.
This paper aims to investigate the binding behaviour between tobermorite, the hydration product of cement and epoxy resin, from the perspective of the molecular level. To obtain a comprehensive understanding, the effects of epoxy chain length, curing agent type, and epoxy direction on the interfacial behaviour were investigated by determining the mean square displacement (MSD), the radial distribution function (RDF), and the interfacial binding energy of the established composite models between tobermorite and epoxy resin with different chemical structures.

Forcefield
The force field illustrates the interaction of an atom with its surrounding atoms, which is the core of mechanical simulation. In this study, the COMPASS force field is used to describe the interaction between atoms, which is an ab initio force field, and its algorithm is largely derived from the early force field, namely the CFF force field. The potential function of the COMPASS force field can be expressed in Equation (1) [22]: where b is the bond, θ is the angle, φ is the torsion angle, χ is the out-of-plane angle, q is the charge, and r the distance between two atoms. In this formula, the first seven terms are the bonding effect of atoms, which are: bond length, bond angle, dihedral angle, bond expansion-bond expansion coupling, bond expansion-bond angle bending coupling, bond expansion-dihedral angle twisting coupling, bond angle bending-dihedral angle torsion coupling, and bond angle bending-bond angle bending coupling, the latter two are non-bonding interactions of atoms, van der Waals and Coulomb interactions, respectively.
The COMPASS force field covers a wide range of covalent molecules, including polymers; organic molecules and small inorganic molecules; and non-covalent bond models covering a range of ionic materials including metals, metal oxides, and metal halides. It has been proven that the COMPASS force field shows good compatibility with cementitious materials [23][24][25][26].

Model Construction
The simulation in this paper was carried out using Biovia Materials Studio 2016 (Biovia Co., Shanghai, China). It has been widely agreed that tobermorite has a similar interlayer structure to C-S-H [27][28][29][30][31][32], therefore the model of tobermorite has been commonly applied to simulate C-S-H due to its complexity [33][34][35][36]. In this paper, based on the monoclinal ordered layer structure of the spatial group B11b and the chemical formula of Ca 5 Si 6 O 16 (OH) 2 •7H 2 O [37], the tobermorite 14 Å was modelled, which is shown in Figure 1. It should be noted that although the calcium atoms and water molecules randomly occupy the interlaminar space of tobermorite 14 Å, in this model, they were assumed to alternately reside in these positions in an ordered model, which is a common approach in the related study [37]. on the monoclinal ordered layer structure of the spatial group formula of Ca5Si6O16(OH)2•7H2O [37], the tobermorite 14 Å was m in Figure 1. It should be noted that although the calcium atom randomly occupy the interlaminar space of tobermorite 14 Å, i assumed to alternately reside in these positions in an ordered mo approach in the related study [37]. Since epoxy is commonly used as a healing agent in self-h was selected to simulate the binding behaviour with a cementitio in Figure 2, the molecular chain of epoxy resin contains unique groups, such as epoxy groups, hydroxyl groups, and ether gro advantages of bisphenol A type epoxy resin such as high bond surface, good stability, low shrinkage, and good mechanical prop as the core of the microcapsule [41]. According to the literature, or equal to 2 were a fluid phase [12], therefore, n was set as 0, 1, 2 i Moreover, in order to study the influence of the curing agent and resin models were established, namely, the epoxy resin model epoxy resin model after the reaction with diluent and curing agen according to the chemical structure formula of epoxy resin is sho temperature, E-51 epoxy resin exhibited a high viscosity, which epoxy's ability to flow into cracks after rupture, and requires viscosity to infiltrate the crack surface. For a self-healing proce flows into the crack, a curing reaction then occurs with the cur ( Figure 4) and TEPA ( Figure 5), while the diluent (BGE, as s participate in the curing reaction. According to the curing rea models of the cured epoxy resin with MC120D and TEPA wer shown in Figures 7 and 8, respectively.  Since epoxy is commonly used as a healing agent in self-healing systems [38,39], it was selected to simulate the binding behaviour with a cementitious matrix. As illustrated in Figure 2, the molecular chain of epoxy resin contains unique active groups and polar groups, such as epoxy groups, hydroxyl groups, and ether groups [40]. In view of the advantages of bisphenol A type epoxy resin such as high bond strength, wide bonding surface, good stability, low shrinkage, and good mechanical properties, E-51 was selected as the core of the microcapsule [41]. According to the literature, epoxies with n less than or equal to 2 were a fluid phase [12], therefore, n was set as 0, 1, 2 in the epoxy resin model. Moreover, in order to study the influence of the curing agent and chain length, six epoxy resin models were established, namely, the epoxy resin model without curing and the epoxy resin model after the reaction with diluent and curing agent. The model established according to the chemical structure formula of epoxy resin is shown in Figure 3. At room temperature, E-51 epoxy resin exhibited a high viscosity, which is unfavourable for the epoxy's ability to flow into cracks after rupture, and requires a diluent to reduce the viscosity to infiltrate the crack surface. For a self-healing process, after the epoxy resin flows into the crack, a curing reaction then occurs with the curing agent, i.e., MC120D ( Figure 4) and TEPA ( Figure 5), while the diluent (BGE, as shown in Figure 6) may participate in the curing reaction. According to the curing reaction principle [42], the models of the cured epoxy resin with MC120D and TEPA were established, which are shown in Figures 7 and 8, respectively.

Model Construction
The simulation in this paper was carried out using Biovia Materials Studio 2016 (Biovia Co., Shanghai, China). It has been widely agreed that tobermorite has a similar interlayer structure to C-S-H [27][28][29][30][31][32], therefore the model of tobermorite has been commonly applied to simulate C-S-H due to its complexity [33][34][35][36]. In this paper, based on the monoclinal ordered layer structure of the spatial group B11b and the chemical formula of Ca5Si6O16(OH)2•7H2O [37], the tobermorite 14 Å was modelled, which is shown in Figure 1. It should be noted that although the calcium atoms and water molecules randomly occupy the interlaminar space of tobermorite 14 Å, in this model, they were assumed to alternately reside in these positions in an ordered model, which is a common approach in the related study [37]. Since epoxy is commonly used as a healing agent in self-healing systems [38,39], it was selected to simulate the binding behaviour with a cementitious matrix. As illustrated in Figure 2, the molecular chain of epoxy resin contains unique active groups and polar groups, such as epoxy groups, hydroxyl groups, and ether groups [40]. In view of the advantages of bisphenol A type epoxy resin such as high bond strength, wide bonding surface, good stability, low shrinkage, and good mechanical properties, E-51 was selected as the core of the microcapsule [41]. According to the literature, epoxies with n less than or equal to 2 were a fluid phase [12], therefore, n was set as 0, 1, 2 in the epoxy resin model. Moreover, in order to study the influence of the curing agent and chain length, six epoxy resin models were established, namely, the epoxy resin model without curing and the epoxy resin model after the reaction with diluent and curing agent. The model established according to the chemical structure formula of epoxy resin is shown in Figure 3. At room temperature, E-51 epoxy resin exhibited a high viscosity, which is unfavourable for the epoxy's ability to flow into cracks after rupture, and requires a diluent to reduce the viscosity to infiltrate the crack surface. For a self-healing process, after the epoxy resin flows into the crack, a curing reaction then occurs with the curing agent, i.e., MC120D ( Figure 4) and TEPA ( Figure 5), while the diluent (BGE, as shown in Figure 6) may participate in the curing reaction. According to the curing reaction principle [42], the models of the cured epoxy resin with MC120D and TEPA were established, which are shown in Figures 7 and 8, respectively.

Simulation
In this study, a total of 7 models of composites were established, which are shown in Table 1. In order to optimize the tobermorite 14 Å structure, the smart algorithm was used for energy minimisation. To establish the interlayer structure, the cleave layer tool was used for the optimized tobermorite 14 Å along [0 0 1] (Z direction) and [0 1 0] (Y direction) to obtain the XY plane and XZ plane, respectively. Then, a three-dimensional space structure was created with vacuum space. In order to build the layered structure, the

Simulation
In this study, a total of 7 models of composites were established, which are shown in Table 1. In order to optimize the tobermorite 14 Å structure, the smart algorithm was used for energy minimisation. To establish the interlayer structure, the cleave layer tool was used for the optimized tobermorite 14 Å along [0 0 1] (Z direction) and [0 1 0] (Y direction) to obtain the XY plane and XZ plane, respectively. Then, a three-dimensional space structure was created with vacuum space. In order to build the layered structure, the epoxy models were built as crystals with length c (along the Z direction) of 10 Å. The layered structure was built with layer 1 and layer 3 of tobermorite and layer 2 of epoxy placed in between. Therefore, the width of the crack in the final model was obtained as 14 Å after matching the model structure, which was in the range of 0.5 nm-100 nm as reported in C-S-H by Ma et al. [43]. As shown in Figure 9, the interlayer structure, in which tobermorite 14 Å was on both sides and epoxy resin was in the middle, was constructed. After the whole composite model was established, the structure of the whole model was optimised again.  During the entire MD simulation process, the time step was set as 1 fs and the temperature was fixed at 298 K. The temperature and pressure were monitored and controlled by Nose Thermostat and Berendsen Barostat. Meanwhile, periodic boundary conditions were used [44]. The entire simulation process was divided into two steps. First, the models were relaxed for 1000 ps in the NPT (N-number of atoms, P-pressure, Ttemperature) isothermal-isobaric ensembles, in which the pressure was set to 0 GPa. Second, a further 2000 ps run was conducted in the NVT (N-number of atoms, Vvolume, T-temperature) ensembles to output the result of MD simulation. During the entire MD simulation process, the time step was set as 1 fs and the temperature was fixed at 298 K. The temperature and pressure were monitored and controlled by Nose Thermostat and Berendsen Barostat. Meanwhile, periodic boundary conditions were used [44]. The entire simulation process was divided into two steps. First, the models were relaxed for 1000 ps in the NPT (N-number of atoms, P-pressure, T-temperature) isothermal-isobaric ensembles, in which the pressure was set to 0 GPa. Second, a further 2000 ps run was conducted in the NVT (N-number of atoms, V-volume, T-temperature) ensembles to output the result of MD simulation.

Structural Analysis
The initial cell lengths in the Z direction before and after MD are listed in Table 2. It can be seen from the table that after MD simulation, the length of all models in the Z direction decreased, which could be due to the reduction in the layer spacing of each compound. Moreover, it should be noted that the cell length in the Z direction of the seven models all decreased after MD simulation, but the reduction in cell length may be attributed to the insertion of the epoxy resin chain, which was smaller than that of the model without epoxy resin. Moreover, as shown in Figure 10, the relative concentration distribution of the three layers of the overall model (the tobermorite layer on both sides and the epoxy resin layer in the middle, models 2 to 6) also reflects a similar phenomenon. In Figure 10a, the atoms of layer 2 were distributed between 51.74−61.52 Å, where there were no atoms of layer 1 and layer 3. On the contrary, the atoms of layer 2 were distributed between 25.80−32.32 Å and the atoms of layer 1 and layer 3 were distributed between 0−35.58 Å and 26.35−59.38 Å, from which we can see that the atoms of the three layers infiltrated each other after the MD simulation. These might be due to the molecular interaction between the different layers, resulting in the decrease in distance and the cracks between the models being reduced. interaction between the different layers, resulting in the decrease in distance and the cracks between the models being reduced.

Mean Square Displacement
The binding energy mainly depends on the interaction between O and N in cured epoxy resin and Ca in tobermorite, which can be reflected by the mean square displacement (MSD) to indicate the repairing efficiency of the cured epoxy resin. The MSD is often used to analyse the displacement of atoms with time, and the estimation of the motion parameters can be obtained in Equation (2). MSD (t) could be applied to describe atoms deviating from their initial position as the function of time.
where r i (0) is the original position for atom i and r i (t) represents the position for atom i at time t, N is the number of atoms in the system. In the polymer composite, the MSD for the oxygen atoms in polymers and the O w atoms in water are often used to indicate the faster movement of interlayer species, including the rotation and vibration of the polymer chains and their branch structures [16]. The MSD of oxygen atoms in each model epoxy resin chain are presented in Figure 11a. It showed that the displacements of the six models were 0.12 Å, 0.12 Å, 0.19 Å, 0.23 Å, 0.20 Å, and 0.22 Å, respectively. It can be clearly seen from the MSD curves of models 2, 3, and 4 that when the chain length was 0 or 1 (models 2 and 3), the displacement was still the same; when the chain length increased to 2 (model 4), a noticeable increment on the MSD curve could be observed, which might be attributed to the increasing chain length enhancing the atomic interaction between the epoxy resin and the tobermorite. Meanwhile, it was found that the MSD also increased when the epoxy resin was cured and the epoxy resin chain was parallel to the Z axis (models 5, 6, and 7). This may be because the O atoms in epoxy received more attraction from tobermorite than in model 2 and model 3. The deduction can be further proven with the increased binding energy between layers, which will be discussed in Section 3.3.
Similarly, the MSD of Si atoms in tobermorite is also shown in Figure 11b. It can be seen from the figure that the displacements of the seven models were 0.13 Å, 0.10 Å, 0.13 Å, 0.21 Å, 0.15 Å, 0.13 Å, and 0.20 Å, respectively. Compared with Figure 11a,b, it is obvious that the MSD of the Si atoms was smaller than that of the O atoms, which indicates that tobermorite showed a more stable skeleton function.
Since N-atoms only existed in the epoxy with a curing agent, the MSD of N atoms in the epoxy resin chain in models 5 and 6 are presented in Figure 11c. It is obvious that the MSD of N atoms in model 6 (0.18 Å) was greater than that in model 5 (0.12 Å), which indicates that the curing agent, TEPA, was more affected by the interaction of tobermorite than that by MC120D.

Binding Energy
The different configurations and dynamic behaviour of the composites could b attributed to the change in interaction between the polymers and the substance

Binding Energy
The different configurations and dynamic behaviour of the composites could be attributed to the change in interaction between the polymers and the substances. Therefore, in this study, the binding energy, which is the opposite of the interaction energy between the tobermorite and the epoxy, was calculated by following Equation (3) [45]: where E b is the binding energy between the tobermorite and the epoxy, E I is the interaction energy, E total is the total energy of the whole model, and E Tobermorite and E epoxy are the energy of the tobermorite and epoxy, respectively. The interfacial binding energy between the epoxy resin and the tobermorite on both sides is presented in Table 3. It can be seen that the interfacial binding energy of models 2 and 3 was close, while the interfacial binding energy of model 4 was larger than the two models, which is consistent with the previous investigations on MSD. Thus, it can be concluded that increasing the chain length can strengthen the interaction strength between epoxy resin and tobermorite. Moreover, in models 5, 6, and 7, the type of curing agent and the direction of the chain increased the interfacial binding energy, with significant effect from the curing agent. As presented in Table 3, after being cured with MC120D and TEPA, the interfacial binding energy increased by 21.33 kcal/mol (Model 5) and 31.03 kcal/mol (Model 6), respectively, while it only increased by 8.43 kcal/mol (Model 7) after changing direction. Therefore, it can be safely deducted that the use of a curing agent can strengthen the interfacial binding energy between epoxy resin and tobermorite, with a higher binding energy from TEPA, which is consistent with Zhang's experimental results [12].

Radial Distribution Function
The radial distribution function (RDF), which deals with spatial atomic correlations, can provide abundant structural information for the tobermorite/polymer composites [46]. Therefore, in addition to the analysis in Section 3.3, the RDF of model 5 and model 6 was investigated to analyse the spatial correlation of the different atoms between the two layers. According to the literature, in which Ca atoms are commonly selected to indicate the connection between polymers and C-S-H [47], the RDF of Ca atoms of tobermorite in the models and other atoms in the epoxy resin is shown in Figure 12a,b. Since it is widely accepted that the position of the first peak is related to the connection between the two atoms [48,49], i the lower value of r indicates a stronger connection. It can be seen from the figure that the first peak of Ca T -H EX , Ca T -O EX , Ca T -C EX , and Ca T -N EX of model 5 appeared at 2.35 Å, 2.35 Å, 2.65 Å, and 3.05 Å, respectively, while model 6 appeared at 2.45 Å, 2.25 Å, 2.65 Å, and 2.35 Å, respectively. Obviously, Ca-O appears earlier than others, indicating that there is a good spatial correlation between Ca and O. Therefore, the Ca-O bond played a major role in the interaction between the tobermorite and the epoxy resin layers. Moreover, by comparing the two figures, it can be clearly seen that the value of r of the Ca and O atoms of model 6 is smaller than that of model 5, which may suggest that TEPA offered a better interaction between O in the epoxy and Ca in the tobermorite. Furthermore, it can be clearly seen that the appearance of the first peaks of N and Ca in model 6 was earlier than that in model 5, which suggested that N in the TEPA played a better role in connecting with Ca. Therefore, the epoxy resin using TEPA acting as the curing agent had better interfacial binding energy with tobermorite compared with MC120D.
between O and H atoms of less than 2.45 Å may indicate the existence of a hydrogen bond [51]. Therefore, based on the current results, it can be deducted that the interfacial binding energy not only consists of the interaction between Ca and other atoms in the epoxy resin, but also of the hydrogen bonding between O atoms and H atoms. Similarly, the positions of first peak of model 5 and model 6 are the same, and the peak height of model 6 (1.31) is greater than that of model 5 (0.88), which further proves a better connection of epoxy with TEPA due to the hydration bond. By calculating the area of the first peak of the RDF, it can be observed that compared with model 5 (0.32), a larger value was obtained in model 6 (0.37), which therefore proved that the hydrogen bond of model 6 was stronger than that of model 5.  In addition, because the O and H atoms between the interface may form hydrogen bonds and increase the interfacial binding energy [50], the RDF of H in the tobermorite and O in the epoxy resin in model 5 and model 6 was simulated, and the results are shown in Figure 13. Obviously, the first peak of both models appeared at around 1.75 Å, indicating that the hydrogen bond existed, since, according to the literature, an RDF value between O and H atoms of less than 2.45 Å may indicate the existence of a hydrogen bond [51]. Therefore, based on the current results, it can be deducted that the interfacial binding energy not only consists of the interaction between Ca and other atoms in the epoxy resin, but also of the hydrogen bonding between O atoms and H atoms. Similarly, the positions of first peak of model 5 and model 6 are the same, and the peak height of model 6 (1.31) is greater than that of model 5 (0.88), which further proves a better connection of epoxy with TEPA due to the hydration bond. By calculating the area of the first peak of the RDF, it can be observed that compared with model 5 (0.32), a larger value was obtained in model 6 (0.37), which therefore proved that the hydrogen bond of model 6 was stronger than that of model 5.

Conclusions
In this work, seven composite models were established to investigate the in binding energy between C-S-H and epoxy resin via molecular dynamics sim Based on the simulation, the following conclusions can be summarized: (1) After the molecular dynamic simulation, the layers of the three-layer structu tobermorite and epoxy resin in the model permeate each other due to the int between atoms, and the cracks between the models were reduced. (2) The MSD results showed that there was a larger displacement of O in the cure which indicated the strong mutual attraction between epoxy resin and tobe Meantime, the N atom in the TEPA showed a larger displacement than MC120D, revealing that TEPA can offer better binding behaviour between th resin and the cementitious matrix. (3) Based on the MD simulation, the incorporation of a curing agent sign increased the binding energy, with a higher value of TEPA (50.53 k suggesting that the system of epoxy with a TEPA curing agent may provide repairing performance. (4) According to RDF analysis, the earliest appearance of the first peak o tobermorite and the O in epoxy was observed, indicating the strongest int between O atoms and Ca atoms. Additionally, the RDF of Ca in tobermorite a cured epoxy further confirmed that the epoxy cured by TEPA can offer binding effect on repairing cementitious materials. In addition, O and H between the interfaces will form hydrogen bonds, thus increasing the in binding energy. (5) The results revealed that the chain length and chain direction showed

Conclusions
In this work, seven composite models were established to investigate the interfacial binding energy between C-S-H and epoxy resin via molecular dynamics simulation. Based on the simulation, the following conclusions can be summarized: (1) After the molecular dynamic simulation, the layers of the three-layer structure of the tobermorite and epoxy resin in the model permeate each other due to the interaction between atoms, and the cracks between the models were reduced. (2) The MSD results showed that there was a larger displacement of O in the cured epoxy, which indicated the strong mutual attraction between epoxy resin and tobermorite. Meantime, the N atom in the TEPA showed a larger displacement than that in MC120D, revealing that TEPA can offer better binding behaviour between the epoxy resin and the cementitious matrix. (3) Based on the MD simulation, the incorporation of a curing agent significantly increased the binding energy, with a higher value of TEPA (50.53 kcal/mol) suggesting that the system of epoxy with a TEPA curing agent may provide a better repairing performance. (4) According to RDF analysis, the earliest appearance of the first peak of Ca in tobermorite and the O in epoxy was observed, indicating the strongest interaction between O atoms and Ca atoms. Additionally, the RDF of Ca in tobermorite and N in cured epoxy further confirmed that the epoxy cured by TEPA can offer a better binding effect on repairing cementitious materials. In addition, O and H atoms between the interfaces will form hydrogen bonds, thus increasing the interfacial binding energy. (5) The results revealed that the chain length and chain direction showed limited influence on the binding energy between the epoxy and the tobermorite, while the type of curing agent significantly boosted the binding energy. Based on the simulation, it exhibited a better ability to repair cracks in cementitious material.
Although the binding effect between epoxy resin and a cementitious matrix has been simulated, which may guide the selection of a curing agent to offer better repairing performance, an experimental investigation should be conducted to verify the simulation.