Effect of Cross Linking on Molecular Structure of Polydimethylsiloxane/Hydroxyapatite: Molecular Dynamics Simulation

1  Introduction

Thermoplastic polymers [1,2] have widely been used in various engineering and medical applications owing to their favorable electrical, thermal, and mechanical characteristics [3,4]. These polymeric materials typically exhibit high electrical resistance but low thermal conductivity. In materials science, molecular simulations are used to explore the potential of combining different material bonds and to predict their properties. Molecular simulation research has been conducting for several decades. During this period, researchers developed some reliable and specialized methods for modeling systems based on scale, state of matter, and other influencing factors. The two main simulation approaches used to study the material properties are continuum-based and atomistic-based methods [5–7]. Consequently, it is highly important to investigate how the temperature affects binding energy and mechanical behavior of materials. To study the properties and energy of polymer structures via molecular dynamics (MD) simulations, Materials Studio (MS) software is widely regarded as one of the most effective simulation tools [8–11].

MD is widely employed in advanced computational materials research to understand the material’s behavior at atomic and molecular scales. MD simulation is also a popular computational method for examining the intrinsic properties of both crystalline and amorphous materials [12,13]. Composites with their exceptional strength-to-weight ratio, corrosion resistance, and high-temperature performance compared to a monophase material (metal, alloy, polymer or ceramic) are therefore extensively applied in aerospace, biomedical, and defense fields. Specially in metals or alloys, the dislocations present in the crystals or lattices significantly influence their mechanical properties [14]. Recently, several MD techniques have increasingly been used to study high-energy-density materials, high-entropy materials, and various novel energetic materials. Many researchers focus on system minimization or optimization studies using MD simulations. For instance, structural optimization is commonly performed after constructing a molecular model. Modeling a structure often generates a high-energy molecule, and simulations are typically initiated from the optimized structure to identify accurate characteristics of its various parameters. BIOVIA Materials Studio provides multiple optimization methods, including conjugate gradient, steepest descent, and Newton-Raphson approaches. Computational simulations serve as valuable tools for gaining deeper insights into molecular dynamics and overcoming experimental limitations. In particular, mesoscopic dynamics (MesoDyn) simulations are employed as a powerful analytical tool for investigating polymer material properties and their dynamic evolution under realistic conditions [15,16].

The process of crystallization is largely influenced by its surrounding environment. Different factors such as confinement and solution composition are known to affect various aspects of crystal formation. Hydroxyapatite (HA), with the molecular formula [Ca10(PO4)6(OH)2] is a key biomineral component of bone [17–19]. Consequently, there is considerable interest in understanding the structure, dynamics, and reactivity of fluids around HA nanoparticles to gain insights into human bone. It is well established that water exchange around the constituent cations (Ca2+ in HA) is the rate-limiting step in the growth of ionic crystals from aqueous solutions. MD simulations of mineral–water interfaces provide a valuable tool for probing the mechanisms that govern crystallization [20].

Nontoxic elastomers such as polydimethylsiloxane (PDMS), with the chemical formula CH3[Si(CH3)2O]nSi(CH3)3, where ‘n’ represents the number of repeating units, are among the most widely used versatile polymers in the organosilicone group [21,22]. While siloxanes generally interact with their surroundings, PDMS is biocompatible, non-toxic, and relatively inert. Its mechanical properties are attracted significant attention in recent years, particularly for microfluidic devices, which can be easily and economically produced with varying degrees of crosslinking. The phase behavior, structural parameters, and temporal evolution of PDMS polymers have been studied by simulating their aggregation in aqueous solutions using the MesoDyn technique [23–25]. This polymer is widely utilized due to its excellent mechanical strength, high dielectric properties, thermal stability, and consistent electrical performance. Given the growing demand for this material, researchers focus on methods to enhance its overall performance, including incorporating nanoparticles to improve specific properties for targeted applications. With recent advancements in computational technology, MD simulation has become an essential tool for investigating material behavior at specific scales [26–29]. Consequently, MD simulation serves as a key technique for predicting the mechanical properties of polymeric materials. Materials Studio (MS), a widely used materials science simulation software developed by Accelrys (USA), has been applied extensively across pharmaceutical, petrochemical, automotive, aerospace, and educational research fields. The software provides a comprehensive implementation of advanced simulation techniques, including molecular dynamics (MD), Monte Carlo (MC), quantum mechanics (QM), dissipative particle dynamics (DPD), and others. It facilitates the realization of three-dimensional (3D) molecular modeling and structural optimization. The Microsoft-standard user interface of Materials Studio software offers an effective platform for analyzing results, allowing precise control through its parameter settings and control panel. MS software mainly allows multiscale modelling to customize material properties. VAMP can predict the electronic features including HOMO (i.e., highest occupied molecular orbital) and LUMO (i.e., lowest unoccupied molecular orbital) gaps, and chemical reactivity using semi-empirical quantum mechanics for organic light-emitting diode (OLED) and pharmaceutical design. On the other hand, Forcite analyses atomistic interactions using classical molecular dynamics to determine aircraft composite and lubricant mechanical strength and adhesion levels. In contrast, MesoDyn simulates mesoscale morphology using dynamic density functional theory (DFT) to forecast mixture stability or phase separation. This is essential for optimizing nanocomposites and consumer emulsions (viz., food, cosmetics, pharmaceuticals, and so on). These techniques link electrical structure to bulk performance in real-world production, allowing researchers to make virtual prototype materials and save experimental costs.

In this context, the crosslinking behavior in a composite system consisting of PDMS polymer and HA ceramics has not yet been modeled through simulation. Some MD studies have examined polymer-based composite interfacial interactions, but they focused on non-reactive systems, physical adsorption mechanisms, or composite systems without explicitly modelling the crosslinking process. The first simulation-based unit cell models of PDMS, HA, and PDMS-HA molecular structures were created using Materials Studio (MS) software in this study. The present study is the first to directly mimic molecular crosslinking between PDMS chains and HA surfaces. This study models the reactive crosslinking mechanism in the presence of filler material (HA), quantifies how the HA surface affects network formation, interfacial bonding, and structural evolution, and provides atomistic insight into how crosslink density and interactions affect composite mechanical and structural properties. This study simulates PDMS/HA composite interface chemistry more realistically than prior non-reactive polymer simulations by including reactive interactions and specifically tracking network development near the HA surface.

Therefore, the present study primarily focuses on understanding the structural changes, interactions, and energetics of molecules. Specifically, it investigates silicone derivatives with porous bioceramic powder substituents attached to the silicone chain. The main objective of this study is to determine their adsorption behavior on phosphorus surfaces, adsorption capacity, and active sites using MD simulations, while examining properties such as crosslinking/bonding effects, bond lengths, bond angles, and torsion angles, along with quantum chemical parameters including deformation energy, kinetic energy, potential energy, non-bond energy, and total energy at optimized geometries within the Materials Studio (MS, BIOVIA Software) environment. Additionally, the thermal properties and scattering behavior of the PDMS-HA composite molecular structure are analyzed. Molecular structural models were constructed in MS software, where all calculations were performed. The molecular simulations were carried out using the VAMP, Forcite, and MesoDyn modules, and these MD methods were applied to evaluate the thermal, physical, and chemical characteristics and energies of PDMS, HA, and PDMS-HA composite molecular models under different conditions using MS software.

2  Modeling of Molecular Structure

2.1 Molecular Structure of PDMS

PDMS elastomer was prepared as the polymer material. Its molecular structure, represented as a “ball-and-stick” model of the PDMS repeat unit, was constructed using Materials Studio software, as shown in Fig. 1 [3]. The PDMS molecular structure, represented by the empirical formula (C2H6OSi)n was developed.

Figure 1: Molecular model of pristine PDMS.

2.2 Molecular Structure of HA

To create hydroxyapatite as a porous bioceramic material, MS software was used to construct the molecular structure of the HA repeat unit in a “ball-and-stick” model, as shown in Fig. 2 [3]. The molecular structure of HA represented by the empirical formula Ca10P6O26H2 was developed.

Figure 2: Molecular model of HA.

2.3 Molecular Structure of PDMS-HA

To create the PDMS-HA crosslinked composite material, the “ball-and-stick” molecular model of the PDMS-HA repeat unit was built in MS software, as shown in Fig. 3 [3]. The molecular structure of PDMS-HA represented by the empirical formula (C2H6OSi)n–Ca10P6O26H2 was constructed.

Figure 3: Molecular model of PDMS-HA.

The PDMS-HA repeat unit (C2H6OSi)n–Ca10P6O26H2 forms a representative polymeric chain, constructed with an appropriate degree of polymerization. Previous studies showed that a polymer chain containing ten to fifteen repeating units meets accuracy requirements. After careful analysis and validation of the simulation, the researcher selected ten repeating units or monomers for the polymer chain. Therefore, the present PDMS-HA bonding structure was designed as a polymeric chain comprising a total of ten monomer units. The PDMS-HA molecular structure model was then initialized in the dynamic module using the Forcite module.

3  Simulation Steps of Materials Studio

The steps for constructing the molecular model and performing calculations [5,6] in molecular simulations using MS software are schematically illustrated in Fig. 4.

Figure 4: Steps for molecular simulation and calculation of PDMS-HA.

The molecular structure of the polymer and porous bioceramic material was determined using the MD simulation technique. Based on the input parameters listed in Table 1, the required model and particle sizes were obtained. It was necessary to adjust the simulation time according to the size of the molecular structure. The MD simulation results include progress reports, molecular size diagrams, and result files. The simulation methodology, techniques, and parameters used in MS software for the materials in this study are summarized in Table 1 [9]. Molecular dynamics simulations were performed using the Universal Force Field (UFF). PDMS chains consisting of repeat units were constructed and introduced into the simulation cell containing hydroxyapatite. After structural optimization, simulations were conducted under the NVE ensemble. The total simulation time step was 0–50 ns. These parameters were selected to ensure structural equilibration and stable energy convergence of the system.

4  Results and Discussion

4.1 Physical, Chemical and Thermal Molecular Properties via VAMP

In this study, convergence was evaluated by performing geometry optimization of PDMS, HA, and PDMS-HA using the VAMP method, as illustrated in Fig. 5. The molecular property results indicate that both PDMS and HA gradually approach convergence and achieve stability as the number of optimization steps increases. Accordingly, the PDMS-HA model was developed and subjected to geometry optimization simulation.

Figure 5: Convergence obtained from PDMS-HA after geometry optimization run for (a) PDMS, (b) HA, and (c) PDMS-HA composite molecules.

Fig. 6 illustrates the relative differences in enthalpy and entropy for PDMS, HA, and PDMS-HA as representative examples. The heat capacity, enthalpy, and entropy were determined using the VAMP method, and the results demonstrate a linear dependence of these properties on temperature (Fig. 6). The thermal properties of PDMS and HA were also calculated using the VAMP approach. At 700 K, the heat capacity of PDMS-HA (305 cal/mol·K) was compared with the corresponding values of PDMS (145 cal/mol·K) and HA (225 cal/mol·K). Likewise, the enthalpy of PDMS-HA (150 kcal/mol) was evaluated against PDMS (55 kcal/mol) and HA (110 kcal/mol) at the same temperature. In addition, the entropy of PDMS-HA (580 cal/mol·K) was compared with PDMS (300 cal/mol·K) and HA (440 cal/mol·K) at 700 K. Although biomedical applications operate near physiological temperatures (~310 K), the thermodynamic properties were evaluated up to 700 K to investigate the thermal stability and structural resilience of the PDMS–HA composite. Elevated temperatures in molecular dynamics simulations facilitate enhanced conformational sampling and provide insight into the robustness of interfacial interactions under extreme conditions. The results therefore reflect intrinsic thermodynamic trends rather than direct physiological behavior. The thermodynamic properties at elevated temperatures were intended to evaluate intrinsic thermal stability and temperature-dependent trends rather than direct biomedical operating conditions. The composite’s thermal ceiling and structural endurance at high temperatures depend on these conditions. HA particles “pin” PDMS chains against heat degradation because the PDMS-HA composite has higher entropy than its constituents. Curing energy when cooling to body temperature is estimated by enthalpy. Heat capacity measures energy dissipation during sterilization or high-energy laser protection, determining thermal shock resistance. Interfacial stability and manufacturing constraints are essential for biomedical coating reliability. These results indicate that, in terms of entropy, the thermal performance of PDMS-HA is 93.33% higher than that of PDMS and 31.82% higher than that of HA. In this context, recent breakthroughs in self-healing and high entropy materials coating include atomic-scale selective oxidation and crack-sealing [30]. It could useful for predicting the self-healing and mechanical bonding properties of the PDMS-HA composites, which were developed by previous studies.

Figure 6: Thermal properties including entropy, heat capacity and enthalpy obtained from VAMP MD simulation for (a) PDMS, (b) HA, and (c) PDMS-HA composite molecules.

The structural changes in the PDMS-HA molecules following the VAMP process indicate notable variations in bond length, bond angle, and torsional effects, as presented in Table 2.

The geometry optimization results, including energy and gradient norm (i.e., a potential energy gradient’s magnitude in relation to the atomic locations) values obtained from the VAMP process for PDMS, HA, and PDMS-HA molecules, are presented in Table 3. It can be observed that the energy of the PDMS molecular model is higher than that of HA for each molecule during geometry optimization.

4.2 Equilibrium of States via Forcite Method

To thoroughly analyze the static equilibrium state of the PDMS-HA molecular structure, the system must undergo dynamic simulation over time before reaching equilibrium. Both temperature and energy are used together to assess the equilibrium condition of the molecular structure. As the system approaches equilibrium, fluctuations in temperature and energy are observed. As shown in Fig. 7, the static potential energy and total energy initially increase, then fluctuate within a defined range, and eventually stabilize as time progresses for both PDMS and HA molecules. The system is considered to have reached equilibrium when the energy variation falls within ±5%–10%.

Figure 7: Energies (viz., total energy, potential energy, kinetic energy and nonbonded energy) curves for (a) PDMS, (b) HA and (c) PDMS-HA after dynamics run using Forcite method.

The energy profile obtained after the dynamic run indicates that PDMS and HA molecules produce lower kinetic energy compared to their potential and/or total energy. Furthermore, both kinetic and potential energies are significantly lower than the nonbond energy observed during the dynamic simulation. Therefore, when using the Forcite method, the kinetic and nonbond energies are more stable than the potential energy for both PDMS and HA molecules during dynamic runs [31]. The stability of the PDMS-HA system was evaluated based on interfacial interaction energy rather than total energy alone. The total energy profile primarily confirms structural equilibration of the system. Despite less fluctuations in total energy for PDMS-HA system confirm system equilibration (see Fig. 7c), but they do not directly quantify interfacial stability. Interfacial energy evaluates the specific “affinity” between the two PDMS and HA phases on the other hand, the crosslink density measures the network connectivity formed by water elimination. Interfacial stability cannot be inferred solely from total energy behavior and would require calculation of binding or interaction energies. The binding energy mainly depends on the potential energy of the material. In the PDMS-HA composite, the binding energy can be evaluated by the difference between potential energy of the composite and the combined potential energy of HA and PDMS. Therefore, binding energy can measure the strength of interaction between PDMS and HA. A more negative binding energy indicates stronger adhesion and enhanced interfacial stability.

The variation in the radius of gyration was analyzed during the simulation to better understand the molecular structural behavior of PDMS and HA beads at different concentrations. As illustrated in Fig. 8, the radius of gyration initially increases, followed by fluctuations within a certain range until it reaches a maximum value. Larger radius of gyration values indicates that the PDMS and HA molecules are oriented perpendicular to the interface. More pronounced fluctuations are observed at lower concentrations, as the interface has not yet reached saturation with PDMS and HA molecules. The overall increase in the radius of gyration suggests that the polymer chains become more extended. Once the interfacial region becomes saturated with PDMS and HA molecules, fewer variations are observed at higher concentrations, as shown in Fig. 8. The magnified image in Fig. 8, which depicts molecular orientation at the interface, indicates that HA molecules are dispersed within and interact with the PDMS matrix. Additionally, the radial distribution function, g(r), representing interactions between non-bonded atoms in PDMS-HA, describes the probability of finding neighboring particles at a distance r (radius of gyration) from a reference particle.

Figure 8: Distribution function vs. Radius of gyration (i.e., probability density) obtained from MD simulation for (a) PDMS, (b) HA, and (c) PDMS-HA composite molecules.

The dynamic energy values of PDMS and HA molecules are presented in Table 4. It is observed that, during dynamic simulation, the energy levels of HA molecules are higher than those of PDMS molecules in both the initial and final states.

The dynamic total energy values of the PDMS-HA composite molecules are presented in Table 5. It shows that during dynamic simulation, the energy of the PDMS-HA molecules is lower than that of both PDMS and HA molecules as obtained from geometry optimization runs.

4.3 Dynamic Density Function via MesoDyn Method

Density functional theory (DFT), based on the assumption that the free energy of an inhomogeneous liquid depends on the local density function, serves as the fundamental principle of the MesoDyn approach. All thermodynamic properties can be derived from the free energy. In this method, each bead corresponds to a specific component type and represents covalently bonded groups of atoms, such as one or several structural units of a polymer chain. The technique dynamically updates the component density fields according to Langevin noise and chemical potential gradients [32]. In MesoDyn simulations, two sets of parameters must be defined to describe the chemical nature of the system: the chain topology in terms of repeat segments (or beads) and the interaction energy between different components. For the first set, MesoDyn uses a Gaussian chain representation with springs and beads. Each bead in the Gaussian chain acts as a statistical unit representing several actual monomers, and different bead types correspond to distinct components. The chain topology depends on the degree of coarse-graining applied to the original system, while the springs model the stretching behavior of chain fragments. Since both parameters must be integers, the Gaussian chain A4B4 was chosen for the PDMS-HA system. Here, A4B4 denotes a coarse-grained model in which A represents 4 PDMS particles or molecules and B represents 4 HA particles or molecules within a composite polymer chain. This approach allows MesoDyn simulations to model large systems at a mesoscopic scale while maintaining computational efficiency.

Fig. 9 shows the fluctuation of the order parameters for beads A and B of the PDMS-HA composite molecule over time steps ranging from 0 to 50 ns. Consequently, the order parameter of the beads in the composite system increased as the aggregation process advanced. The figure clearly illustrates the two-stage bead formation process. In the first stage (bead A), the order parameters approach 0, indicating that the components of the system are widely dispersed. During the second stage (bead B), the order parameters of each component increase significantly, reflecting the formation of ordered aggregates. At this stage, the order parameter fluctuates within a narrow range rather than remaining constant. Fig. 9 shows that the initial stage for bead B is nearly equivalent to that of bead A [16]. It is to be noted that although MesoDyn simulations are less accurate and stable compared to atomistic MD simulations, the MesoDyn method is much faster and can efficiently cover longer timescales. The atomistic MD simulations provide detailed insight into the interfacial interactions between PDMS chains and hydroxyapatite, including binding energies, hydrogen bonding behaviour, and local structural organization. These molecular-level interactions determine the thermodynamic compatibility between the two components. Based on these interaction characteristics, MesoDyn simulations were employed to investigate the molecular-level interactions influence mesoscopic morphology, including phase distribution and domain evolution. In this way, the mesoscopic results do not stand independently but rather extend the atomistic findings to larger length and time scales, enabling prediction of dispersion quality and structural homogeneity in the composite. MesoDyn mainly works on the mesoscale, which ranges from nanometres to micrometres. It predicted the organization of the polymer chains and ceramic particles over a longer period and volume. In contrast, Forcite functions at the atomistic level (angstroms to nanometres). It computed certain interactions, such as the van der Waals forces and the covalent Si–O–P bonding.

Figure 9: Stages of the order parameters for PDMS-HA molecular structure.

The free energy density, defined as energy per unit volume (~RT/V, where R is the universal gas constant in kcal/mol·K, T is the absolute temperature in K, and V is the molar volume in Å3) for PDMS, HA, and the PDMS-HA composite molecules calculated during the MD simulations, is presented in Fig. 10. Compared to values obtained after 25 ns, the optimized conformations of HA adsorbed on PDMS surfaces show minimal conformational changes by the end of the MD run, resulting in a nearly linear curve. This is important because the adsorption process is influenced not only by the specific surface chemistry of the crystalline polymer exposed to the bioceramic powder but also by the relative size of the crystalline facets over time. Surface area plays a critical role in adhesion, particularly over extended timescales. Fig. 10a depicts the free energy density for PDMS. The curve begins slightly above 0.3727 kcal/Å3 and gradually decreases over time, reaching approximately 0.3715 kcal/Å3 at around 49 ns, indicating a reduction in free energy density with time. Fig. 10b shows the free energy density for HA. The curve starts at ~0.3719 kcal/Å3, rises slightly to ~0.3722 kcal/Å3 at 12 ns, peaks near 0.3725 kcal/Å3 at 23 ns, and then decreases to 0.3715 kcal/Å3 by ~49 ns. The fluctuations reflect the presence of stronger intermolecular interactions, such as hydrogen bonding or ionic forces. Fig. 10c presents the free energy density for the PDMS-HA composite. The curve starts at a high value of ~0.3727 kcal/Å3 and decreases similarly to PDMS, but the decline is slightly more gradual, reaching as 0.3718 kcal/Å3 at ~30 ns due to the incorporation of the bioceramic material. At equilibrium (~30 ns), the free energy densities for PDMS-HA, PDMS, and HA are 0.3718, 0.3720, and 0.3719 kcal/Å3, respectively, indicating that the equilibrium free energy density of PDMS-HA is 0.054% and 0.027% lower than that of PDMS and HA, respectively. These modest percentage shifts occur in this simulation since it uses coarse-grained beads and mean-field density functionals instead of discrete atoms. It is to be noted that although the change in equilibrium free energy density (0.054% and 0.027%) seems very small, but the thermodynamics indicates that a negative shift in free energy (ΔG < 0) implies that the “mixed” state is more stable than the “separated”. In the PDMS-HA composite, the PDMS polymer chains and HA ceramic molecules form a novel interaction due to the newly formed interfaces. The polymer chains with a high affinity for the ceramic surface (owing to Van der Waals forces, hydrogen bonding, or electrostatic interactions) would release energy resulting to a decreased system energy. The change in the equilibrium free energy density tends to more in PDMS since the polymer chains near the HA ceramic surface often transition from a random coil to a more “ordered” or adsorbed state. In comparison to the pristine polymer, this local densification or alignment may result in a lower energy state of greater magnitude (0.054%). In contrast, the surface energy of HA ceramic is often high and when the PDMS chains “wetted” the HA molecules (i.e., surface energy is partially met), they helped only in “relaxing” the HA particle’s surface-atoms, led to decrease in equilibrium free energy density at lower scale (0.027%).

Figure 10: Time evolution of the free energy density in kcal/Å3 for (a) PDMS (b) HA, and (c) PDMS-HA molecular structure.

The Langevin equations were applied to account for the frictional and random forces in both equilibrium and non-equilibrium dynamic systems, as summarized in Table 6. The step size factor is defined as 1 for a full step, 0.5 for a half step, and 0.25 for a quarter step. In this context, polymer and porous bioceramic particles undergo random motion due to the dynamics of the system. Based on the dynamic simulation results, the molecular models of the studied molecules are updated according to classical mechanics, following Newton’s or Hamilton’s equations (Eq. (1)) [33].

mdvdt=FTotal(1)

here, m represents the particle mass, dv/dt is the acceleration, v is the velocity, t is the time, and F is the force applied to the particle. By using Eq. (1), the final results of the MD simulation within the MesoDyn method are obtained. At each step of the random motion, the particle size factor may vary depending on the positions of the atoms.

4.4 Scattering Effect of Molecular Structure

The scattering behavior of PDMS, HA, and the PDMS-HA molecular structures is shown in Fig. 11. Simulated X-ray scattering patterns, obtained using the scattering effect module of Materials Studio, reveal the amorphous nature of PDMS through a broad, diffuse pattern, while HA displays weak oscillatory features due to short-range crystalline order. In the PDMS–HA composite, the suppression and broadening of HA diffraction features indicate strong interfacial interactions and disruption of long-range order, confirming a homogeneous dispersion of HA within the polymer matrix. The contributions to the X-ray scattering peaks are shown separately for PDMS (Fig. 11a), HA (Fig. 11b), and the PDMS-HA composite (Fig. 11c). The simulated X-ray scattering patterns were analysed quantitatively by extracting peak positions. The major diffraction peaks at 2θ = 25.9°, 31.7°, and 32.9° correspond to the (002), (211), and (300) planes of hydroxyapatite and are consistent with reported experimental data. Minor peak shifts observed in the PDMS–HA composite suggest interfacial interaction induced lattice distortion. Changes in relative peak intensities further indicate modifications in local structural ordering near the interface. The polymerization peaks observed for PDMS, HA, and PDMS-HA primarily arise from intermolecular backbone-backbone correlations. In the PDMS-HA composite, some sharp HA crystal peaks are still visible, while the lower-intensity peak corresponding to the amorphous PDMS indicates interchain packing. The simulation results using the united-atom model closely reproduce experimental X-ray scattering data, particularly capturing the polymerization peak and the amorphous characteristics [34]. The properties obtained from MD simulations closely matched with reported experimental data. For example, the peak positions of major X-ray scattering peaks of the present simulated PDMS/HA composites evidently resembled the experimental XRD results obtained in our reported study [34]. Small deviations are attributed to limitations of the force field, finite-size effects, and idealized simulation conditions, but the overall agreement supports the reliability of the molecular-level insights. This effect may result from the inclusion of explicit hydrogen atoms in the PDMS-HA model, which influences the scattering pattern.

Figure 11: Scattering effect of (a) PDMS (b) HA and (c) PDMS-HA molecular structure.

4.5 Cross Linking Simulation

To crosslink polymer chains, a chemical reaction must occur between them. Uncrosslinked polymer chains are more mobile within the polymer matrix and are connected through crosslinking. Crosslinked chains behave differently from uncrosslinked ones. Crosslinking enhances chemical and solvent resistance, as the chains resist flow under stress and swell less in solvents compared to loose-chains. When two or more molecules are covalently bonded, as illustrated in Fig. 12, a strong connection and durable structure are formed. Chemical modification involves attaching or removing functional groups to change solubility or other properties of the original molecule. The chemical reaction between PDMS and HA is shown in Eq. (2). In this reaction, water molecule can be eliminated after crosslinking occurs between the two molecular models. Reactive silanol groups (–Si–OH) are present on PDMS chains, and surface hydroxyl groups exist on HA, and a condensation reaction expressed in Eq. (3). Since the two oxygen-bearing groups unite to create a stable siloxane-type bridge (Si–O–P or Si–O–Ca complex), a water molecule (H2O) is eliminated.

(C2H6OSi)n+Ca10(PO4)6(OH)2=(C2H6OSi)n−(Ca10P6O26H2)(2)

Si−OH + HO−surface→Si−O−surface +H2O(3)

Figure 12: Cross linking effect of PDMS-HA composite using Forcite and VAMP methods.

In this case, water elimination is chemically reasonable because it corresponds to a dehydration condensation reaction, forming a Si–O–Ca or Si–O–P linkage. The interfacial interaction between PDMS and hydroxyapatite is carefully analyzed from a chemical perspective. In the absence of reactive silanol groups, the dominant interactions are hydrogen bonding and electrostatic attractions between the PDMS siloxane backbone and HA surface hydroxyl and calcium sites. Covalent bonding is only possible when silanol-functionalized PDMS undergoes condensation with surface hydroxyl groups of HA, leading to Si–O surface linkages accompanied by water elimination. Therefore, the nature of the interfacial interaction depends on the chemical functionality explicitly introduced in the model. The crosslinking mechanism was simulated by forming bonds between nearby reactive sites. After selecting a reactive OH group on HA, adjacent reactive hydrogen sites on PDMS are identified. Bonds are then formed between the closest available pairs, and new potential pairs are generated for subsequent reactions. At each step of the simulation, bonds are formed between the pair of reactive sites (one on PDMS and one on HA) that are closest to each other. Once a bond is formed, these two sites are removed from the pool of potential pairs, and the distances between the remaining unreacted sites are recalculated to identify the next closest pair. This process is repeated until the desired crosslinking density is reached. The system is relaxed using cycles of molecular mechanics (MM) and molecular dynamics (MD) to remove any unfavorable interactions caused by the formation of new bonds. Additionally, any potential chain spearing or ring catenation is detected and removed from the system.

The polymer system (PDMS) was modeled in Materials Studio software using the chain growth method, as illustrated in Fig. 13. The resulting composite system exhibits more natural characteristics with increased amorphousness while retaining a certain degree of molecular crosslinking. Crosslinked polymer molecules are packed into a periodic cell at a specified density, either manually or using an algorithm such as the Amorphous Cell module in Materials Studio. Residual unreacted functional groups are then utilized for final crosslinking, and molecular dynamics (MD) simulations are applied to relax and stabilize the molecular structure.

Figure 13: Cross linking effect of PDMS simulation.

The modelling method for porous bioceramic (HA) systems in Materials Studio software is illustrated in Fig. 14. To relax the system, energy minimization was performed after 1000 steps of MD simulation. The system was then analyzed to identify nearby reactive sites. Subsequently, residual unreacted functional groups were utilized to complete the final crosslinking, and MD simulations were applied to relax the molecular structure using the Forcite and VAMP methods.

Figure 14: Cross linking effect of HA simulation.

5  Conclusions

In this study, all molecular models of PDMS, HA, and PDMS-HA composites were constructed using Materials Studio software. The crosslinking behavior between PDMS and HA in a composite system is modeled through molecular dynamics (MD) simulations. The research primarily focused on investigating the equilibrium states and energy of PDMS, HA, and PDMS-HA under different conditions using MD simulations in Materials Studio. The results indicated that the total energy of PDMS-HA increased over time, confirming that the PDMS-HA composite exhibits greater stability and compatibility compared to pure PDMS and HA due to its consistent energy fluctuations. The intermolecular interaction energy between PDMS-HA crystal surfaces was also found to be strong.

Phase behavior and mesoscopic structure analysis using the MesoDyn simulation method revealed that the equilibrium energy density of the PDMS-HA composite is 0.054% lower than that of PDMS but 0.027% lower than HA. Apparently, this change in energy is very small, but the composite’s thermodynamic stability and interfacial bonding are enhanced by the interface of polymer chains and ceramic molecules. This change comes owing to the differing surfaces as the PDMS polymer chains and HA ceramic molecules interact differently in the composite. High-affinity polymer chains for ceramic surfaces lower the system energy compared to the polymer. The local densification or alignment saves 0.054% energy over pure polymers. In contrast, PDMS chains “wetted” and “relaxed” HA molecules or particles, lowering the equilibrium free energy density 0.027% compared to HA ceramic. Furthermore, the crosslinking in PDMS-HA achieved appropriate bond lengths, bond angles, and torsion angles compared to the original chemical structures of PDMS and HA. These results suggest that the crosslinked PDMS-HA composite can serve as a stronger and tougher coating or binder through chemically docking of two or more molecules via covalent bonds.

This research study provides molecular- and mesoscopic level insights into the interactions between PDMS and hydroxyapatite, highlighting an interfacial bonding, crosslinking behaviour, and chain mobility that govern composite structure and stability. While the simulations offer valuable guidance, certain limitations, including simplified reactive assumptions, finite chain lengths, and high-temperature thermodynamic evaluations, must be considered when extrapolating to experimental systems. Nevertheless, the findings have direct implications for material design: optimizing PDMS functionalization and HA dispersion can enhance adhesion, mechanical performance, and structural uniformity in biomedical devices, bone scaffolds, and bioactive coatings. In recent studies, the advanced self-healing high entropy coating are investigated by atomic-scale selective oxidation and crack-sealing methods. Furthermore, the multiscale simulation approach demonstrates a pathway to bridge molecular-level understanding with mesoscopic morphology, providing a framework for rational design of PDMS/HA-based self-healing composites.

Acknowledgement: The authors acknowledge the SRM Institute of Science and Technology, Kattankulathur for providing all the facilities utilized for this study.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: The authors confirm contribution to the paper as follows: study conception and design: Chellaiah Ayyanar, Sumit Pramanik; data collection: Chellaiah Ayyanar; analysis and interpretation of results: Chellaiah Ayyanar, Sumit Pramanik; draft manuscript preparation: Chellaiah Ayyanar, Sumit Pramanik. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Data available on request from the authors. The data that support the findings of this study are available from the corresponding author, [Sumit Pramanik], upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Yu KQ, Li ZS, Sun J. Polymer structures and glass transition: a molecular dynamics simulation study. Macromol Theory Simul. 2001;10(6):624–33. doi:10.1002/1521-3919(20010701)10:6624::AID-MATS624>3.0.CO;2-K. [Google Scholar] [CrossRef]

2. Gartner TE III, Jayaraman A. Modeling and simulations of polymers: a roadmap. Macromolecules. 2019;52(3):755–86. doi:10.1021/acs.macromol.8b01836. [Google Scholar] [CrossRef]

3. Nouri N, Ziaei-Rad S. A molecular dynamics investigation on mechanical properties of cross-linked polymer networks. Macromolecules. 2011;44(13):5481–9. doi:10.1021/ma2005519. [Google Scholar] [CrossRef]

4. Schneider F, Fellner T, Wilde J, Wallrabe U. Mechanical properties of silicones for MEMS. J Micromech Microeng. 2008;18(6):065008. doi:10.1088/0960-1317/18/6/065008. [Google Scholar] [CrossRef]

5. Sharma S, Kumar P, Chandra R. Molecular dynamics simulation of polymer-matrix composites using BIOVIA materials studio, LAMMPS, and GROMACS. In: Molecular dynamics simulation of nanocomposites using BIOVIA materials studio, lammps and gromacs. Amsterdam, The Netherlands: Elsevier; 2019. p. 141–225. [Google Scholar]

6. Sharma S, Kumar P, Chandra R. Molecular dynamics simulation of metal matrix composites using BIOVIA materials studio, LAMMPS, and GROMACS. In: Molecular dynamics simulation of nanocomposites using BIOVIA materials studio, lammps and gromacs. Amsterdam, The Netherlands: Elsevier; 2019. p. 101–40. [Google Scholar]

7. Sharma S, Kumar P, Chandra R. Molecular dynamics simulation of ceramic matrix composites using BIOVIA materials studio, LAMMPS, and GROMACS. In: Molecular dynamics simulation of nanocomposites using BIOVIA materials studio, lammps and gromacs. Amsterdam, The Netherlands: Elsevier; 2019. p. 227–58. [Google Scholar]

8. Alder BJ, Wainwright TE. Studies in molecular dynamics. I General Method J Chem Phys. 1959;31(2):459–66. doi:10.1063/1.1730376. [Google Scholar] [CrossRef]

9. Han Y, Elliott J. Molecular dynamics simulations of the elastic properties of polymer/carbon nanotube composites. Comput Mater Sci. 2007;39(2):315–23. doi:10.1016/j.commatsci.2006.06.011. [Google Scholar] [CrossRef]

10. Krishna S, Sreedhar I, Patel CM. Molecular dynamics simulation of polyamide-based materials—a review. Comput Mater Sci. 2021;200(5):110853. doi:10.1016/j.commatsci.2021.110853. [Google Scholar] [CrossRef]

11. Adnan A, Sun CT, Mahfuz H. A molecular dynamics simulation study to investigate the effect of filler size on elastic properties of polymer nanocomposites. Compos Sci Technol. 2007;67(3–4):348–56. doi:10.1016/j.compscitech.2006.09.015. [Google Scholar] [CrossRef]

12. Urata S, Li S. A multiscale model for amorphous materials. Comput Mater Sci. 2017;135(3–4):64–77. doi:10.1016/j.commatsci.2017.03.029. [Google Scholar] [CrossRef]

13. Massobrio C. The structure of amorphous materials using molecular dynamics. Bristol, UK: IOP Publishing; 2022. [Google Scholar]

14. Gowthaman S. A review on mechanical and material characterisation through molecular dynamics using large-scale atomic/molecular massively parallel simulator (LAMMPS). Funct Compos Struct. 2023;5(1):012005. doi:10.1088/2631-6331/acc3d5. [Google Scholar] [CrossRef]

15. Fraaije JGEM, van Vlimmeren BAC, Maurits NM, Postma M, Evers OA, Hoffmann C, et al. The dynamic mean-field density functional method and its application to the mesoscopic dynamics of quenched block copolymer melts. J Chem Phys. 1997;106(10):4260–9. doi:10.1063/1.473129. [Google Scholar] [CrossRef]

16. Li Y, Xu G, Zhu Y, Wang Y, Gong H. Aggregation behavior of Pluronic copolymer in the presence of surfactant: mesoscopic simulation. Colloids Surf A Physicochem Eng Aspects. 2009;334(1–3):124–30. doi:10.1016/j.colsurfa.2008.10.029. [Google Scholar] [CrossRef]

17. de Leeuw NH. Computer simulations of structures and properties of the biomaterial hydroxyapatite. J Mater Chem. 2010;20(26):5376. doi:10.1039/b921400c. [Google Scholar] [CrossRef]

18. Di Tommaso D, Prakash M, Lemaire T, Lewerenz M, De Leeuw N, Naili S. Molecular dynamics simulations of hydroxyapatite nanopores in contact with electrolyte solutions: the effect of nanoconfinement and solvated ions on the surface reactivity and the structural, dynamical, and vibrational properties of water. Crystals. 2017;7(2):57. doi:10.3390/cryst7020057. [Google Scholar] [CrossRef]

19. Mostafa NY, Brown PW. Computer simulation of stoichiometric hydroxyapatite: structure and substitutions. J Phys Chem Solids. 2007;68(3):431–7. doi:10.1016/j.jpcs.2006.12.011. [Google Scholar] [CrossRef]

20. Likhachev I, Balabaev N, Bystrov V, Paramonova E, Avakyan L, Bulina N. Molecular dynamics simulation of the thermal behavior of hydroxyapatite. Nanomaterials. 2022;12(23):4244. doi:10.3390/nano12234244. [Google Scholar] [PubMed] [CrossRef]

21. Fojtíková J, Kalvoda L. Molecular dynamics simulations of poly(dimethylsiloxane) elasticity. Acta Phys Pol A. 2018;134(3):857–8. doi:10.12693/aphyspola.134.857. [Google Scholar] [CrossRef]

22. Sharfeddin A, Volinsky AA, Mohan G, Gallant ND. Comparison of the macroscale and microscale tests for measuring elastic properties of polydimethylsiloxane. J Appl Polym Sci. 2015;132(42):42680. doi:10.1002/app.42680. [Google Scholar] [CrossRef]

23. Yang S, Yuan S, Zhang X, Yan Y. Phase behavior of tri-block copolymers in solution: mesoscopic simulation study. Colloids Surf A Physicochem Eng Aspects. 2008;322(1–3):87–96. doi:10.1016/j.colsurfa.2008.02.029. [Google Scholar] [CrossRef]

24. Weisgraber TH, Gee RH, Maiti A, Clague DS, Chinn S, Maxwell RS. A mesoscopic network model for permanent set in crosslinked elastomers. Polymer. 2009;50(23):5613–7. doi:10.1016/j.polymer.2009.09.046. [Google Scholar] [CrossRef]

25. Fojtíková J, Kalvoda L, Sedlák P. Molecular dynamics simulations of poly(dimethylsiloxane) properties. Acta Phys Pol A. 2015;128(4):637–40. doi:10.12693/aphyspola.128.637. [Google Scholar] [CrossRef]

26. Guo Y, Zhang D, Zhang X, Wu Y. Thermal properties and mechanical behavior of functionalized carbon nanotube-filled polypropylene composites using molecular dynamics simulation. Mater Today Commun. 2023;37(3):107510. doi:10.1016/j.mtcomm.2023.107510. [Google Scholar] [CrossRef]

27. Chua J, Tu Q. A molecular dynamics study of crosslinked phthalonitrile polymers: the effect of crosslink density on thermomechanical and dielectric properties. Polymers. 2018;10(1):64. doi:10.3390/polym10010064. [Google Scholar] [PubMed] [CrossRef]

28. Alder BJ, Wainwright TE. Phase transition for a hard sphere system. J Chem Phys. 1957;27(5):1208–9. doi:10.1063/1.1743957. [Google Scholar] [CrossRef]

29. Yu R, Wang Q, Wang W, Xiao Y, Wang Z, Zhou X, et al. Polyurethane/graphene oxide nanocomposite and its modified asphalt binder: preparation, properties and molecular dynamics simulation. Mater Des. 2021;209(5):109994. doi:10.1016/j.matdes.2021.109994. [Google Scholar] [CrossRef]

30. Luo H, Jin K, Tao J, Wang H. Properties prediction and design of self-healing epoxy resin combining molecular dynamics simulation and back propagation neural network. Mater Res Express. 2021;8(4):045308. doi:10.1088/2053-1591/abf66b. [Google Scholar] [CrossRef]

31. Kumar U, Rathi R, Sharma S. Carbon nano-tube reinforced nylon 6,6 composites: a molecular dynamics approach. Eng Solid Mech. 2020;8(4):389–96. doi:10.5267/j.esm.2020.2.002. [Google Scholar] [CrossRef]

32. Wu H, Xin Y. Molecular dynamics and MesoDyn simulations for the miscibility of polyvinyl alcohol/polyvinyl pyrrolidone blends. Plast Rubber Compos. 2017;46(2):69–76. doi:10.1080/14658011.2017.1280642. [Google Scholar] [CrossRef]

33. Essiz S, Coalson RD. A rigid-body Newtonian propagation scheme based on instantaneous decomposition into rotation and translation blocks. J Chem Phys. 2006;124(14):144116. doi:10.1063/1.2158996. [Google Scholar] [PubMed] [CrossRef]

34. Ayyanar C, Pramanik S. Effect of graphene on self-healing performance of hydroxyapatite/polydimethysiloxane composites. Mater Manuf Process. 2023;38(9):1068–80. doi:10.1080/10426914.2022.2146714. [Google Scholar] [CrossRef]