How to evaluate the interaction between epoxy and a free surface via metadynamics: A coarse-grained example

Abstract

Molecular investigations into polymer-surface interactions are essential for understanding and optimizing the mechanical behavior of, e.g., adhesive joints or polymer nanocomposites. In this context, metadynamics provides a powerful tool for determining the free energy of a molecular system through Gaussian sampling. While the literature demonstrates that metadynamics is well suited for adherend–adhesive systems, there is little discussion on the choice of the metadynamics parameters and their impact on the resulting free energy. This contribution provides a step-by-step methodology for determining the parameters of metadynamics, specifically the Gaussian’s depth, width, and deposition time, as well as the convergence. To illustrate this process, we utilize a simple coarse-grained system made up of a rigid substrate and either an epoxy or hardener molecule. Once we establish the optimal metadynamics parameters, we examine the effects of non-bonded polymer-surface interactions, which are modeled using 12-6 Lennard–Jones potentials, on the resulting free energy surface. As expected, we find that the distance parameter $σ_{ps}$ influences both the position and depth of the free energy well, while the energy parameter $ε_{ps}$ only affects the well depth. The methodology outlined here can be applied to more complex systems, making this article a useful guide for advanced studies on adhesive joints and polymer surface interactions in general.

Keywords

Epoxy coarse-graining metadynamics interface molecular dynamics

1. Introduction

Structural adhesives are high-performance bonding agents that provide strong, durable connections between various materials and are thus a key technology for cost-efficient and sustainable engineering solutions. Compared to conventional fastening methods, such as riveting, bolting, or welding, they offer a better strength-to-weight ratio and allow for bonding dissimilar materials without requiring holes [1,2]. This enables the use of high-performing fiber-reinforced polymers [3], making structural adhesive bonds widely used today in automotive [4], aerospace [5], and civil engineering [6]. The adhesive is often based on epoxy resins [7], as these are highly durable and cost-efficient [8]. Consequently, a curing process is required in which the epoxy and a hardening agent form a polymer network, which enables the load transfer. As a side effect of the curing, an interphase region forms in the adhesive near the adherend surfaces. This interphase differs from the bulk adhesive in terms of curing degree [9], glass transition [10], and mechanical properties [11,12]. Consequently, a comprehensive understanding of the interphase, particularly of its mechanical properties, is key to unlock the full potential of adhesive joints. However, this interphase, depending on material and process parameters, exhibits a thickness of only a few nanometers to several micrometers [10,13] and is therefore difficult to detect and characterize via experiments [11]. Similar interphase effects and challenges arise for other polymeric interphases, e.g., in polymer nanocomposites [14 –16].

As a remedy, particle-based numerical methods, such as molecular dynamics (MD), are commonly employed since these offer the required spatial resolution to accurately model the adhesive’s microstructure [17]. All-atom MD simulates every individual atom and thus offers detailed insights into the interphase formation [18,19], the impact of different adherend surfaces [20,21], or the detrimental effect of moisture [22,23]. However, all-atom MD is too computationally expensive to simulate samples large enough to evaluate the mechanical properties. The much more efficient coarse-grained MD (CGMD) solves this issue by substituting groups of atoms by superatoms, which considerably reduces the degrees of freedom [24]. Thus, CGMD facilitates the investigation of the mechanical behavior of adhesive joints [25,26] and their interphases [27 –29]. One crucial issue is the correct coarse-grained representation of non-bonded interactions, i.e., the Van der Waals and Coulomb interactions between the substrate surface and the polymer. One practical approach is to calibrate the coarse-grained potentials so that CGMD replicates the all-atom free energy surface, evaluated with metadynamics. Metadynamics [30] is a simulation method that introduces a time-dependent bias potential to a molecular system, facilitating the exploration of its free energy landscape. For example, Lau et al. [31] employ metadynamics to construct the free energy of a single DGEBA molecule detaching from a silica surface. They upscale their results by estimating the intrinsic strength at the interface via a worm-like chain-based fracture model [32] and calibrate a traction-separation law describing the epoxy-silica interface at the continuum scale. In a similar vein, Yang and Qu [33] use metadynamics to obtain the coarse-grained interactions at an epoxy-copper interface. They observe a mechanical interphase region where particle displacements differ from the bulk polymer and derive a traction-separation law from their results.

Coarse-graining is commonly used to model the mechanics of adhesive joints due to its computational efficiency [17]. However, most existing models do not feature the critical surface interactions. The literature indicates that metadynamics is an effective tool for integrating adherend–adhesive interactions into established models without a complete reparametrization. While transferability errors may arise, this metadynamics-based extension offers an efficient exploration of surface-chemistry effects, yielding valuable qualitative insights into adhesive joint behavior. However, the identification of the metadynamics parameters and their impact on the determined free energy surface often remains ambiguous in the literature. Therefore, this article uses the example of a simple coarse-grained adherend–adhesive system to illustrate the calibration of a metadynamics analysis in detail. Our research questions are: 1) How are the metadynamics parameters calibrated for an adherend–adhesive system? 2) How do the non-bonded interactions between the adherend and the adhesive influence their free energy surface? Hence, this contribution aims to deliver a reproducible and instructive protocol for conducting coarse-grained metadynamics investigations of polymer–surface interactions. By employing a minimalistic coarse-grained model, we ensure the primary focus remains on the methodology of free-energy sampling, rather than the physical interpretation of the results. In Section 2, we introduce the coarse-grained model, the basics of metadynamics, and the investigated sample setup. We then determine the Gaussian width, the time to reach convergence, and the Gaussian height in Section 3. With these optimal metadynamics parameters, we determine the free energy of the system and discuss the impact of the non-bonded substrate-polymer interactions. The methodology presented here relies on open-access tools and can be easily adapted to more complex systems, making this article a valuable starting point for future research.

2. Methods

2.1. Coarse-grained model

2.1.1. Epoxy

We perform the following investigations in coarse-grained resolution. To this end, we choose a common thermoset comprising the resin Diglycidyl ether of bisphenol A (DGEBA) and the hardener 4,4’-Diaminodicyclohexylmethane (PACM). Giuntoli et al. [34] calibrated the corresponding coarse-grained interaction potentials via the Energy Renormalization approach [5,35] with Gaussian process surrogate models [36]. They focused on reproducing the mechanical properties and obtained good agreement of the Debye-Waller factor, density, Young’s modulus, and yield stress with respect to their all-atom reference simulations. In addition, the model is computationally efficient and has thus been used to study the fracture behavior [37], impact of non-stoichiometric mixing ratios [38], and the interphase formation of adhesive joints [29]. The level of coarse-graining is visualized in Figure 1, and we provide the force fields in 4. In this study, we consider single DGEBA and PACM molecules and thus do not need to simulate the curing reactions.

Figure 1.

Simulation setup: rigid substrate with either PACM (top) or DGEBA (bottom) placed in simulation box, periodic boundary conditions in the x and z-directions, repulsive Lennard–Jones wall in the upper y-direction; insets visualize coarse-graining for PACM and DGEBA by Giuntoli et al. [34]; collective variable s measures the shortest distance of DGEBA’s and PACM’s center of mass from the substrate’s surface.

2.1.2. Substrate

We model the substrate, exemplarily, as a single crystal with a body-centered cubic lattice and lattice spacing $ℓ = 5 Å$ , representing a generic substrate material. Since we perform all simulations under the canonical ensemble (NVT), i.e., constant volume and temperature, the substrate beads would only fluctuate about their equilibrium position, which merely adds noise to the free energy evaluation. Hence, we consider the substrate as rigid by excluding it from the time integration, which additionally reduces the numerical cost. We model the non-bonded interactions between substrate and polymer via a shifted 12–6 Lennard–Jones potential with energy and distance parameters $ε_{ps}$ and $σ_{ps}$ , respectively. We choose a very large cutoff radius $r_{c} = 15 Å$ to ensure that we incorporate all relevant interactions. In order to study the impact of these interactions on the free energy, we evaluate all 30 combinations of $ε_{ps} \in {0.5, 1.0, 1.5, 2.0, 2.5, 3.0} kcal {mol}^{- 1}$ and $σ_{ps} \in {3.0, 3.5, 4.0, 4.5, 5.0} Å$ . The energy parameter $ε_{ps}$ sets the depth of the potential well, while $σ_{ps}$ determines the potential’s zero crossing and thus specifies the equilibrium distance at $r_{0} = \sqrt[6]{2} σ_{ps}$ . Thus, with a larger $σ_{ps}$ , beads located deeper in the substrate also influence the DGEBA or PACM molecule. For our system consisting of a rigid substrate and a mobile resin or hardener molecule, both larger $ε_{ps}$ and larger $σ_{ps}$ lead to a larger potential energy.

2.1.3. Transferability

Bottom-up CG models are usually calibrated based on the potential of mean force (PMF) obtained from a fine-scale, e.g., atomistic, reference simulation. The PMF is typically dependent on material composition and thermodynamic state, rendering bottom-up CG models similarly susceptible to changes in these conditions. Consequently, CG potentials are usually not transferable across different thermodynamic states (e.g., temperatures, densities, or pressures) or other local environments (bulk vs. interfacial regions) [39]. Hence, one must be extremely careful when extending an existing CG bulk model to incorporate an interface, as outlined here. The present CG model [34] employs the energy renormalization approach [35] to explicitly incorporate a wide range of curing degrees for the bulk epoxy. However, to show that the model can be extended to include an interface, it is necessary to carefully validate it using atomistic reference simulations. Because atomistic molecular dynamics (MD) simulations are usually too computationally expensive to directly analyze mechanical behavior, validation instead focuses on the structural features of the adhesive near the substrate. For example, Yang and Qu [33] derived the coarse-grained (CG) epoxy-copper interaction by detaching a single epoxy molecule in both CG and all-atom (AA) resolutions using metadynamics. They then validated their CG model by comparing the microstructural characteristics observed in densely packed AA and CG samples after curing [40]. Consequently, while calibration of adhesive-adherend interaction potentials can be performed with a single molecule, validation requires large-scale atomistic reference simulations to identify key microstructural features—an effort that is beyond the scope of this work. In general, the bi-material setup has to be considered already during the CG potential parametrization [39]. Addressing the limited transferability of coarse-grained models, novel strategies such as local density potentials [41] as well as dual [42] and microcanonical [43] approaches aim to enhance density and temperature transferability, respectively [44]. In this context, Jin et al. [45] provide a detailed review of enhanced bottom-up coarse-graining techniques and discuss transferability issues and possible remedies in depth.

2.2. Metadynamics

Metadynamics [30,46,47] identifies the free energy F of a molecular system with respect to a (small) number of collective variables s by adding external Gaussian potentials every $τ_{G}$ time steps along its trajectory [48]. By summing these Gaussians over time t and the number of collective variables d, we obtain the potential

\begin{matrix} V (s, t) = \sum_{k τ_{G} < t} ω \exp (- \sum_{i = 1}^{d} \frac{{[s_{i} - s_{i} (k τ_{G})]}^{2}}{2 δ s_{i}^{2}}) \end{matrix}

(1)

with Gaussian height ω and width $δ s_{i}$ of the ith collective variable. Practically, the repulsive Gaussian potentials gradually fill in the free energy surface, pushing the system out of local minima and allowing us to explore rare events. Hence, after sufficient time, we have completely mapped the collective variable space such that the potential V converges to the free energy F

\begin{matrix} \lim_{t \to \infty} V (s, t) = - F (s) + C_{fill} (t) \end{matrix}

(2)

where $C_{fill} (t)$ denotes a constant accounting for overfilling [48]. In the enhanced well-tempered metadynamics [49], the heights of the Gaussian hills are scaled at each step to improve convergence and avoid overfilling. To this end, a fictitious collective variable temperature $Δ T$ is introduced as an additional variable, that controls the extent of exploration. As a consequence, well-tempered metadynamics do no converge to the free energy $- F$ but to $- \frac{Δ T}{T + Δ T} F$ . Given the simplicity of our systems, the additional advantages offered by well-tempered metadynamics are not necessary, allowing us to avoid the complexity of selecting and calibrating the additional parameter $Δ T$ . The convergence time $t_{C}$ defines the period $t < t_{C}$ during which the potential well is not fully filled, while for $t > t_{C}$ we observe overfilling, i.e., convergence. Assuming a system with a single collective variable s whose free energy should be zero for $s_{C}^{\min} \leq s \leq s_{C}^{\max}$ , we can estimate $C_{fill} (t)$ as

\begin{matrix} C_{fill} (t) = \frac{1}{s_{C}^{\max} - s_{C}^{\min}} \int_{s_{C}^{\min}}^{s_{C}^{\max}} V (t) d s \end{matrix}

(3)

and compute the free energy via (2). In order to improve the estimation of the free energy, we compute the arithmetic average [50]

\begin{matrix} \bar{F} (s, t_{tot}) = - \frac{1}{t_{tot} - t_{C}} \int_{t_{C}}^{t_{tot}} V (s, t) - C_{fill} (t) d t \end{matrix}

(4)

for $t_{C} \leq t \leq t_{tot}$ , with total simulation time $t_{tot}$ . Analogously, we evaluate the associated standard deviation, to assess the statistical uncertainty of the averaged $\bar{F}$ .

The sampling of $\bar{F}$ becomes more accurate for smaller δs and is optimal at approximately half the fluctuations of the unbiased system, which can be determined in a conventional MD simulation [46,47]. Furthermore, the accuracy depends on the ratio $ω ∕ τ_{G}$ , meaning that the error can be reduced by either a smaller height ω or a larger Gaussian deposition time $τ_{G}$ [46]. Finally, we ensure convergence by time-averaging the estimated free energy $\bar{F}$ for $t > t_{C}$ [47], where $t_{C}$ marks the time when the free energy surface is completely mapped and overfilling occurs, which we discuss in Subsection 3.2. Furthermore, we compare the metadynamics solution to Umbrella sampling to validate our results. Alternatively, we could apply error estimators [47], which, however, exceeds the scope of this work.

2.3. Sample configurations

Similar to other contributions [31,33], we evaluate the free energy of a single resin or hardener molecule in the vicinity of a generic substrate. To this end, we initialize a simulation box $40 Å \times 100 Å \times 40 Å$ , and insert a single pre-equilibrated DGEBA or PACM molecule, referred to as mobile beads since only these are time-integrated. On the lower y-side of the simulation box, we place the space-centered cubic substrate comprising 576 beads, spanning the whole simulation box in the x- and z-direction. The substrate features a thickness of four crystal layers, i.e., $20 Å > r_{c}$ in the y-direction to ensure the lower y-end of the simulation box is beyond the cutoff radius of the mobile beads. To confine the mobile resin or hardener molecule to the simulation box, we apply periodic boundary conditions in the x- and z-directions and employ a repulsive Lennard–Jones wall to the upper y-boundary of the simulation box. To avoid interference with the repulsive wall, we choose $s_{C}^{\min} = 22 Å$ and $s_{C}^{\max} = 26 Å$ for the estimation of $C_{fill}$ according to (3). The sample setup is visualized in Figure 1 for resin (bottom) and hardener (top), respectively.

2.4. Simulation parameters

For the MD simulations, we employ LAMMPS [51,52] with the metadynamics plugin PLUMED [48,53]. We perform all simulations at 300 K under NVT ensemble with time step size $Δ t = 4 fs$ and temperature coupling parameter $t_{T} = 400 fs$ for the Nosè-Hoover thermostat [54], as suggested by [34].

According to [47], a successful metadynamics run requires (1) equilibrating the sample at the desired temperature, (2) choosing the collective variables, (3) setting the metadynamics parameters, and (4) ensuring convergence. Since we insert the mobile beads pre-equilibrated and the substrates are rigid, we already fulfill step 1. The collective variables have to describe the state of the system accurately, which makes their selection one of the most challenging aspects of metadynamics [47]. In our case, however, using the shortest distance between the mobile beads’ center of mass and the substrate’s surface beads, as depicted in Figure 1, appears as an obvious choice [31,33]. Due to the simple structure of the resin and hardener molecules in coarse-grained resolution, we neglect their rotation, leaving us with only one collective variable s. Using a single collective variable can lead to projection artifacts and may overlook important orientation information, potentially oversimplifying the free-energy representation of complex systems. While a two-dimensional surface could offer more detail, as demonstrated by [55], our emphasis on presenting a methodology for metadynamics sampling makes this simplification acceptable for this study. Hence, we only need to set the parameters for the actual metadynamics run, i.e., the Gaussian deposition time $τ_{G}$ , as well as the Gaussian height ω and width δs and ensure convergence by identifying the convergence time $t_{C}$ .

By inserting pre-equilibrated DGEBA and PACM molecules and due to the straightforward selection of collective variables, we already fulfill steps 1 and 2 [47]. We address steps 3 and 4 in Section 3 by determining δs, $t_{C}$ , and ω for different non-bonded interactions defined by $ε_{ps}$ and $σ_{ps}$ in Subsection 3.1, Subsection 3.2, and Subsection 3.3, respectively.

3. Results and discussions

3.1. Gaussian width δs

First, we determine the Gaussian width δs, which defines the resolution of the free energy surface in the collective variable space. To this end, we perform a standard (unconstrained) molecular dynamics run for 4 ns for all 30 $σ_{ps}$ - $ε_{ps}$ combinations. We monitor the standard deviation of our collective variable $std (s)$ , i.e., the fluctuations of the minimal distance of the resin or hardener molecule’s center of mass from the substrate’s surface. Figure 2 depicts $std (s)$ over the non-bonded well-depth $ε_{ps}$ for different non-bonded distance parameters $σ_{ps}$ . Evidently, these fluctuations decrease as $ε_{ps}$ and $σ_{ps}$ increase, since a higher potential energy, achieved through a deeper well $ε_{ps}$ or longer range $σ_{ps}$ , constrains the molecules more tightly. Resin and hardener molecules behave very similarly, since we obtain maxima of 3.75 Å and 3.43 Å, and minima of 0.08 Å and 0.15 Å for DGEBA and PACM, respectively. Since a smaller Gaussian width δs yields higher accuracy, we choose the overall minimum to set $δs = 0.5 std (s) = 0.04 Å$ as recommended by [47]. For this study, we use a uniform δs across all parameter combinations for simplicity and consistency. Please note that using such an extremely small Gaussian width is computationally inefficient and may lead to rugged free-energy estimates. While the former is not a problem due to our minimalist systems, the latter is partly mitigated by the time-averaging described by equation (4).

Figure 2.

Indentification of Gaussian width δs for (a) DGEBA and (b) PACM with standard deviation $std (s)$ of collective variable s, obtained in conventional molecular dynamics run for all $ε_{ps}$ - $σ_{ps}$ combinations; insets show s over t for low and high $ε_{ps}$ .

3.2. Convergence time $t_{C}$

Second, we ensure convergence by identifying the filling time $t_{C}$ the metadynamics run requires to completely fill the free energy landscape, which leads to diffusive, i.e., unconstrained, behavior of the collective variable. Figure 3(a) exemplifies the convergence behavior of our metadynamics run for $σ_{ps} = 5 Å$ and $ε_{ps} = 3 kcal {mol}^{- 1}$ . Initially, the free energy well fills rapidly, but converges for $t > 120 ns$ resulting in the converged $\bar{F}$ , which is averaged for $128 ns \leq t \leq 160 ns$ , and features negligible standard deviation for our simple setup. To validate this result, we perform Umbrella sampling with the same sample setup and following the LAMMPS tutorial by [56]. We evaluate the free energy in the interval $0 \leq s \leq 20 Å$ and move the harmonic potential in increments of $0.2 Å$ . In each increment, we first equilibrate the sample for 0.4 ns and then sample the collective variable for 1.6 ns. We employ the weighted histogram analysis method (WHAM) [57,58] implemented by [59] to derive the free energy profiles. As the Umbrella sampling depends on the strength of the harmonic potential k, we average for $k \in {7.5, 10, 12.5, 15, 17.5, 20} kcal mo l^{- 1} Å^{- 2}$ , cf. Figure 7 in the Appendix. Smaller k do not affect the molecules position, while larger k lead to inaccurate estimations of the free energy [56]. Figure 3(a) shows that the resulting free energy matches well with the one obtained via metadynamics and thus validates our approach. Figure 3(b) depicts the time evolution of the free energy minimum $\min (\bar{F})$ for different interaction strengths $ε_{ps}$ . Evidently, larger $ε_{ps}$ lead to a larger convergence time $t_{C}$ , as a deeper free energy well needs to be mapped. Larger interaction range $σ_{ps}$ has a similar effect, which renders $σ_{ps} = 5 Å$ and $ε_{ps} = 3 kcal {mol}^{- 1}$ the most demanding combination. Consequently, we perform all simulations for $t_{tot} = 160 ns$ , i.e., 40,000,000 time steps, to reach $t_{C}$ for all considered $ε_{ps}$ - $σ_{ps}$ combinations. Due to the simplicity of the considered molecules this can be easily achieved in our case, but we recommend individual determination of $t_{C}$ for more complex systems, to reduce the computational effort.

Figure 3.

Convergence for DGEBA molecule: (a) Free energy $\bar{F}$ over collective variable s for different times t for $σ_{ps} = 5 Å$ and $ε_{ps} = 3 kcal {mol}^{- 1}$ , converged $\bar{F}$ averaged over $128 ns \leq t \leq 160 ns$ , band indicates negligible standard deviation, reference solution obtained with umbrella sampling averaged over different spring constants, band indicates standard deviation. (b) Evolution of minimum energy $\min (\bar{F})$ over time t for $σ_{ps} = 5 Å$ and different $ε_{ps}$ .

3.3. Gaussian height ω

According to [46], metadynamics’ accuracy is proportional to the ratio $ω ∕ τ_{G}$ of Gaussian height ω and Gaussian deposition time $τ_{G}$ . Since the two parameters are not independent, we vary ω and set $τ_{G} = 4000 fs$ to avoid computational bottlenecks due to too frequently placed Gaussians [47]. We evaluate the Gaussian height from $0.05 k_{B} T$ to $0.5 k_{B} T$ resulting in the range $0.03 kcal {mol}^{- 1} \leq ω \leq 0.3 kcal {mol}^{- 1}$ for $T = 300 K$ . The effect of ω on the resulting averaged free energies $\bar{F} (s)$ is shown in Figure 4(a) and (b) for DGEBA and PACM, respectively. Considering DGEBA in Figure 4(a), ω does not influence the shape of $\bar{F} (s)$ but its level. If ω is too small, we underestimate the well-depths as we do not reach complete filling. Once we apply a minimum Gaussian height $ω \approx 0.1 kcal {mol}^{- 1}$ , the resulting $\bar{F}$ merely fluctuates around the converged solution, as evident from $\min (\bar{F})$ in Figure 4(a) inset. Larger ω sample the collective variable space quicker, but might also impair the accuracy [47]. Consequently, the choice of the Gaussian height ω represents a compromise between high efficiency (large ω) and high accuracy (low ω). Since we make the same observations for PACM in Figure 4(b), we choose $ω = 0.18 kcal {mol}^{- 1}$ , corresponding to $0.11 k_{B} T$ , and $τ_{G} = 4000 fs$ for all further simulations. The chosen small Gaussian width $δs = 0.04 Å$ resolves $\bar{F}$ in high detail, but also causes the large fluctuations obtained for $10 Å \leq s \leq 20 Å$ . However, these numerical fluctuations can be easily mitigated by applying, e.g., a Savitzky–Golay filter [60].

Figure 4.

Impact of Gaussian height ω on the converged free energy landscape $\bar{F}$ for (a) DGEBA and (b) PACM with color scheme given in (b); insets visualize the minimum free energy $\min (\bar{F})$ over ω.

3.4. Impact of non-bonded interactions

Having determined the metadynamics parameters $δs = 0.04 Å$ , $ω = 0.18 kcal {mol}^{- 1}$ , $τ_{G} = 4000 fs$ , and the filling time $t_{C}$ , we can now study the impact of the non-bonded interactions between substrate and resin or hardener on their free energy surface. To this end, we perform the metadynamics simulations for all 30 combinations of $ε_{ps} \in {0.5, 1.0, 1.5, 2.0, 2.5, 3.0} kcal {mol}^{- 1}$ and $σ_{ps} \in {3.0, 3.5, 4.0, 4.5, 5.0} Å$ with simulated time $t_{tot} = 160 ns$ .

3.4.1. DGEBA

Figure 5(a) (left) displays the effect of the non-bonded energy parameter $ε_{ps}$ on the averaged free energy $\bar{F} (s)$ . As expected, a larger $ε_{ps}$ leads to a lower free energy well $\min (\bar{F} (s))$ , while the minimum’s position $s (\min (\bar{F}))$ is barely affected, cf. Figure 5(b). In contrast, Figure 5(a) (right) reveals that larger $σ_{ps}$ decrease both $\min (\bar{F} (s))$ and its position. This is because a larger $σ_{ps}$ increases the effective range of non-bonded interactions, so that even more distant substrate beads have an impact on the DGEBA molecule. Figure 5(b) summarizes these trends and shows that $\min (\bar{F} (s))$ decreases proportional to $ε_{ps}$ and $σ_{ps}$ . The position of the minimum increases linearly with $σ_{ps}$ and is almost independent of $ε_{ps}$ .

Figure 5.

Impact of the DGEBA-substrate interactions on the free energy landscape: (a) Free energy $\bar{F}$ over collective variable s for varying $ε_{ps}$ (left) and $σ_{ps}$ (right). (b) Minima of free energy $\min (\bar{F})$ (top) and their location $s (\min (\bar{F}))$ over energy coefficient $ε_{ps}$ , for all $ε_{ps}$ - $σ_{ps}$ combinations.

3.4.2. PACM

In contrast to DGEBA, the free energy of the PACM molecule in Figure 6(a) (left) features two local minima, ${\bar{F}}_{A}^{\min}$ and ${\bar{F}}_{B}^{\min}$ , located at $s_{A}^{\min} = 0.98 Å$ and $s_{B}^{\min} = 1.99 Å$ , respectively. The insets in Figure 6(a) reveal that PACM is oriented convexly to the substrate surface in case A and can therefore reach a smaller distance $s_{A}^{\min}$ than in case B with a concave arrangement. The double wells form for $ε_{ps} \geq 1.0 kcal {mol}^{- 1}$ and deepen with increasing $ε_{ps}$ . In contrast, Figure 6(a) (right) demonstrates that with increasing $σ_{ps}$ , the two minima merge into one. We again attribute this to the increased range of the non-bonded interactions at larger $σ_{ps}$ , causing several substrate beads to exert an influence and their effect to overlap to only one minimum. Figure 6(b) shows that the general influence of $σ_{ps}$ and $ε_{ps}$ on $\min (\bar{F})$ and their position $s (\min (\bar{F}))$ is similar to that already discussed for DGEBA, cf. Figure 5(b). The jumps in the position of the minimum can be explained by the double well and the change from ${\bar{F}}_{A}^{\min} < {\bar{F}}_{B}^{\min}$ to ${\bar{F}}_{A}^{\min} > {\bar{F}}_{B}^{\min}$ .

Figure 6.

Impact of the PACM-substrate interactions on the free energy landscape: (a) Free energy $\bar{F}$ over collective variable s for varying $ε_{ps}$ (left) and $σ_{ps}$ (right); insets visualize the PACM configurations of the two local minima A and B. (b) Minima of free energy $\min (\bar{F})$ (top) and their location $s (\min (\bar{F}))$ over energy coefficient $ε_{ps}$ , for all $ε_{ps}$ - $σ_{ps}$ combinations.

Overall, the non-bonded interactions between the substrate and resin or hardener, parametrized in the Lennard–Jones parameter $ε_{ps}$ and $σ_{ps}$ , have a considerable impact on the resulting free energy $\bar{F} (s)$ . Depending on the examined molecule, $\bar{F} (s)$ features one or more minima, the depths of which depend on $ε_{ps}$ and $σ_{ps}$ , while their position depends only on $σ_{ps}$ .

4. Conclusion and outlook

In this contribution, we demonstrate how to determine the free energy landscape of an adherend–adhesive system via metadynamics and how the free energy is affected by the non-bonded adherend–adhesive interactions. To this end, we choose a simple system consisting of one resin or hardener molecule and a substrate surface in coarse-grained resolution. We show how the metadynamics parameters, i.e., the Gaussian width, filling time, Gaussian height, and Gaussian deposition time, are calibrated step by step. The optimal Gaussian width δs results from the fluctuations of the collective variable in an unbiased MD system. We estimate the filling time $t_{C}$ and thus the system’s convergence either from the fluctuations of the collective variable or from the filling constant, and identify a clear dependence of $t_{C}$ on the non-bonded interactions. The optimal Gaussian height ω represents a compromise between accuracy and efficiency, and results for our system in $ω = 0.18 kcal {mol}^{- 1}$ at which the lowest free energy occurs. With the identified optimal metadynamics parameters, we investigate the influence of the non-bonded interactions between substrate and resin or hardener, i.e., the Lennard–Jones parameters $ε_{ps}$ and $σ_{ps}$ on the free energy landscape $\bar{F}$ of the system. As expected, $\bar{F}$ strongly depends on the molecule type, $ε_{ps}$ and $σ_{ps}$ . While $σ_{ps}$ impacts both the position and depth of the free energy well, $ε_{ps}$ only affects its depth.

This article serves as a guide for advanced studies, and the methodology presented can be easily applied to more complex systems. We exclusively use open-source software and provide all the scripts utilized in this research via Zenodo [61], making this study a straightforward introduction to metadynamics in the context of adhesives. In addition, we intend to utilize the determined free energy dependencies of $ε_{ps}$ and $σ_{ps}$ to calibrate the coarse-grained non-bonded interactions between an Al₂O₃ substrate and PACM-DGEBA adhesive in a future contribution.

Simulation software and data

We employ LAMMPS (version 12 June 2025) [51,52] and the open-source, community-developed PLUMED library (v2.9.4) [48,53] for the MD and metadynamics simulations performed on the High-Performance-Computing (HPC) parallel CPU cluster Fritz of the Erlangen National High Performance Computing Center (NHR@FAU) of “Regionales Rechenzentrum Erlangen (RRZE).” We used for the evaluation, data processing and visualization of the simulations mainly MathWorks MATLAB R2023b and OVITO [62], as well as the WHAM algorithm [59]. We provide the data published in this contribution via Zenodo [61].

Footnotes

Appendix 1 Acknowledgements

We would like to express our gratitude to [] for publishing and documenting their research data, which served as the foundation for this contribution. In addition, we extend our thanks to the PLUMED consortium for their excellent documentation and tutorials on their metadynamics package.

ORCID iD

Maximilian Ries

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: M. Ries is funded by the DFG - 377472739 (Research Training Group GRK 2423 ‘Fracture across Scales - FRASCAL’). The author gratefully acknowledges the scientific support and HPC resources provided by the Erlangen National High Performance Computing Center (NHR@FAU) of the Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU) under the NHR project b275dc “Molecular Modeling of Adhesive Joints”. NHR funding is provided by federal and Bavarian state authorities. NHR@FAU hardware is partially funded by the German Research Foundation (DFG) - 440719683.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

Campilho

Banea

Chaves

, et al. eXtended finite element method for fracture characterization of adhesive joints in pure mode I. Nato Sc S Ss Iii C S 2011; 50(4): 1543–1549.

Pramanik

Basak

Dong

, et al. Joining of carbon fibre reinforced polymer (CFRP) composites and aluminium alloys – a review. Compos Part A Appl Sci Manuf 2017; 101: 1–29.

Taib

Boukhili

Achiou

, et al. Bonded joints with composite adherends. part I. effect of specimen configuration, adhesive thickness, spew fillet and adherend stiffness on fracture. Int J Adhes Adhes 2006; 26(4): 226–236.

Scattina

Peroni

, et al. Numerical analysis of hybrid joining in car body applications. J Adhes Sci Technol 2011; 25(18): 2409–2433.

Song

Hsu

Shull

, et al. Energy renormalization method for the coarse-graining of polymer viscoelasticity. Macromolecules 2018; 51(10): 3818–3827.

Zhang

Chen

. Reinforced concrete beams strengthened in flexure with near-surface mounted (NSM) CFRP strips: current status and research needs. Compos Part B Eng 2017; 131: 30–42.

Sancaktar

. Classification of adhesive and sealant materials. In: da Silva

Öchsner

Adams

(eds) Handbook of adhesion technology. Berlin and Heidelberg: Springer, 2018, p. 283–317.

Kim

Lim

Hwang

, et al. Composition of adhesives. In: da Silva

Öchsner

Adams

(eds) Handbook of adhesion technology. Berlin and Heidelberg: Springer, 2018, p. 319–343.

Wehlack

Possart

. Chemical structure formation and morphology in ultrathin polyurethane films on metals. In: Possart

(ed.) Adhesion current research and applications. Weinheim: Wiley, 2005, p. 71–88.

10.

Meiser

Kübel

Schäfer

, et al. Electron microscopic studies on the diffusion of metal ions in epoxy–metal interphases. Int J Adhes Adhes 2010; 30(3): 170–177.

11.

Possart

. Mechanical interphases in adhesives: experiments, modelling and simulation. Doctoral dissertation, LTM, 2014.

12.

Roche

Bouchet

Bentadjine

. Formation of epoxy-diamine/metal interphases. Int J Adhes Adhes 2002; 22(6): 431–441.

13.

Aufray

André Roche

. Epoxy–amine/metal interphases: influences from sharp needle-like crystal formation. Int J Adhes Adhes 2007; 27(5): 387–393.

14.

Mészáros

Horváth

Vas

, et al. Investigation of the correlations between the microstructure and the tensile properties multi-scale composites with a polylactic acid matrix, reinforced with carbon nanotubes and carbon fibers, with the use of the fiber bundle cell theory. Compos Sci Technol 2023; 242: 110154.

15.

Zhang

Kolluru

, et al. Determination of mechanical properties of polymer interphase using combined Atomic Force Microscope (AFM) experiments and finite element simulations. Macromolecules 2018; 51(20): 8229–8240.

16.

Ries

Possart

Steinmann

, et al. A coupled MD-FE methodology to characterize mechanical interphases in polymeric nanocomposites. Int J Mech Sci 2021; 204: 106564.

17.

Ries

. Mechanical behavior of adhesive joints: a review on modeling techniques. Comput Methods Mater Sci 2024; 24(4): 5–35.

18.

Lian

, et al. Influence of cross-linking density on the structure and properties of the interphase within supported ultrathin epoxy films. J Mater Sci 2016; 51(19): 9019–9030.

19.

Kanamori

Kimoto

Toriumi

, et al. On the cyclic fatigue of adhesively bonded aluminium: experiments and molecular dynamics simulation. Int J Adhes Adhes 2021; 107: 102848.

20.

Deng

Song

Lan

, et al. The surface modification effect on the interfacial properties of glass fiber-reinforced epoxy: a molecular dynamics study. Nanotechnol Rev 2022; 11(1): 1143–1157.

21.

Konrad

Zahn

. Interfaces in reinforced epoxy resins: from molecular scale understanding towards mechanical properties. J Mol Model 2023; 29(8): 243.

22.

Büyüköztürk

Buehler

Lau

, et al. Structural solution using molecular dynamics: fundamentals and a case study of epoxy-silica interface. Int J Solids Struct 2011; 48(14–15): 2131–2140.

23.

Tam

Lau

. Moisture effect on the mechanical and interfacial properties of epoxy-bonded material system: an atomistic and experimental investigation. Polymer 2015; 57: 132–142.

24.

Müller-Plathe

. Coarse-graining in polymer simulation: from the atomistic to the mesoscopic scale and back. ChemPhysChem 2002; 3(9): 754–769.

25.

Raos

Zappone

. Tuning adhesion and energy dissipation in polymer films between solid surfaces via grafting and cross-linking. Macromolecules 2024; 57(14): 6502–6513.

26.

Solar

Qin

Buehler

. Molecular mechanics and performance of crosslinked amorphous polymer adhesives. J Mater Res 2014; 29(9): 1077–1085.

27.

Langeloth

Sugii

Böhm

, et al. The glass transition in cured epoxy thermosets: a comparative molecular dynamics study in coarse-grained and atomistic resolution. J Chem Phys 2015; 143(24): 243158.

28.

Langeloth

Sugii

Böhm

, et al. Formation of the interphase of a cured epoxy resin near a metal surface: reactive coarse-grained molecular dynamics simulations. Soft Mater 2014; 12(Suppl. 1): S71–S79.

29.

Dötschel

Richter

Possart

, et al. Reactive coarse-grained MD models to capture interphase formation in epoxy-based structural adhesive joints. Eur J Mech A Solid 2026; 116: 105801.

30.

Laio

Parrinello

. Escaping free-energy minima. Proc Natl Acad Sci 2002; 99(20): 12562–12566.

31.

Lau

Büyüköztürk

Buehler

. Characterization of the intrinsic strength between epoxy and silica using a multiscale approach. J Mater Res 2012; 27(14): 1787–1796.

32.

Keten

Buehler

. Asymptotic strength limit of hydrogen-bond assemblies in proteins at vanishing pulling rates. Phys Rev Lett 2008; 100(19): 198301.

33.

Yang

. An investigation of the tensile deformation and failure of an epoxy/cu interface using coarse-grained molecular dynamics simulations. Model Simul Mater Sc 2014; 22(6): 065011.

34.

Giuntoli

Hansoge

van Beek

, et al. Systematic coarse-graining of epoxy resins with machine learning-informed energy renormalization. NPJ Comput Mater 2021; 7(1): 168.

35.

Xia

Song

Jeong

, et al. Energy-renormalization for achieving temperature transferable coarse-graining of polymer dynamics. Macromolecules 2017; 50(21): 8787–8796.

36.

Schulz

Speekenbrink

Krause

. A tutorial on Gaussian process regression: modelling, exploring, and exploiting functions. J Math Psychol 2018; 85: 1–16.

37.

Pal

Dansuk

Giuntoli

, et al. Predicting the effect of hardener composition on the mechanical and fracture properties of epoxy resins using molecular modeling. Macromolecules 2023; 56(12): 4447–4456.

38.

Ries

Tiefenthäler

Losert

, et al. Deciphering elastoplastic properties from atomistic structure: reactive coarse-grained MD for epoxies. Int J Adhes Adhes 2026; 145: 104196.

39.

Noid

. Perspective: coarse-grained models for biomolecular systems. J Chem Phys 2013; 139(9): 090901.

40.

Yang

Gao

. A molecular dynamics study of tensile strength between a highly-crosslinked epoxy molding compound and a copper substrate. Polymer 2013; 54(18): 5064–5074.

41.

Chu

Voth

. Coarse-grained modeling of the actin filament derived from atomistic-scale simulations. Biophys J 2006; 90(5): 1572–1582.

42.

Dans

Zeida

Machado

, et al. A coarse grained model for atomic-detailed DNA simulations with explicit electrostatics. J Chem Theory Comput 2010; 6(5): 1711–1725.

43.

Ding

Sharma

Chalasani

, et al. Ab initio RNA folding by discrete molecular dynamics: from structure prediction to folding mechanisms. RNA 2008; 14(6): 1164–1173.

44.

Noid

Szukalo

Kidder

, et al. Rigorous progress in coarse-graining. Annu Rev Phys Chem 2024; 75(1): 21–45.

45.

Jin

Pak

Durumeric

, et al. Bottom-up coarse-graining: principles and perspectives. J Chem Theory Comput 2022; 18(10): 5759–5791.

46.

Laio

Rodriguez-Fortea

Gervasio

, et al. Assessing the accuracy of metadynamics. J Phys Chem B 2005; 109(14): 6714–6721.

47.

Laio

Gervasio

. Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep Prog Phys 2008; 71(12): 126601.

48.

consortium

. Promoting transparency and reproducibility in enhanced molecular simulations. Nat Methods 2019; 16(8): 670–673.

49.

Barducci

Bussi

Parrinello

. Well-tempered metadynamics: a smoothly converging and tunable free-energy method. Phys Rev Lett 2008; 100(2): 020603.

50.

Micheletti

Laio

Parrinello

. Reconstructing the density of states by history-dependent metadynamics. Phys Rev Lett 2004; 92(17): 170601.

51.

Plimpton

. Fast parallel algorithms for short-range molecular dynamics. J Comput Phys 1995; 117(1): 1–19.

52.

Thompson

Plimpton

Mattson

. General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions. J Chem Phys 2009; 131(15): 154107.

53.

Tribello

Bonomi

Branduardi

, et al. PLUMED 2: new feathers for an old bird. Comput Phys Commun 2014; 185(2): 604–613.

54.

Evans

Holian

. The Nose–Hoover thermostat. J Chem Phys 1985; 83(8): 4069–4074.

55.

Sun

Han

Swanson

JMJ

, et al. Molecular transport through membranes: accurate permeability coefficients from multidimensional potentials of mean force and local diffusion constants. J Chem Phys 2018; 149(7): 072310.

56.

Gravelle

Alvares

CMS

Gissinger

, et al. A set of tutorials for the LAMMPS simulation package [article v1.0]. Living J Computat Mol Sci 2025; 6(1): 3027.

57.

Kumar

Rosenberg

Bouzida

, et al. The weighted histogram analysis method for free-energy calculations on biomolecules. i. the method. J Comput Chem 1992; 13(8): 1011–1021.

58.

Kumar

Rosenberg

Bouzida

, et al. Multidimensional free-energy calculations using the weighted histogram analysis method. J Comput Chem 1995; 16(11): 1339–1350.

59.

Grossfield

. An implementation of wham: the weighted histogram analysis method, version 2.0. 9. 2014. http://membrane.urmc.rochester.edu/sites/default/files/wham/doc.pdf.

60.

Savitzky

Golay

MJE

. Smoothing and differentiation of data by simplified least squares procedures. Anal Chem 1964; 36(8): 1627–1639.

61.

Ries

. How to evaluate the interaction between epoxy and a free surface via metadynamics: a coarse-grained example - dataset. 2025. https://doi.org/10.5281/zenodo.16282488

62.

Stukowski

. Visualization and analysis of atomistic simulation data with OVITO–the open visualization tool. Model Simul Mater Sci Eng 2009; 18(1): 015012.

63.

Marrink

De Vries

Mark

. Coarse grained model for semiquantitative lipid simulations. J Phys Chem B 2004; 108(2): 750–760.

How to evaluate the interaction between epoxy and a free surface via metadynamics: A coarse-grained example

Abstract

Keywords

1. Introduction

2. Methods

2.1. Coarse-grained model

2.1.1. Epoxy

2.1.2. Substrate

2.1.3. Transferability

2.2. Metadynamics

2.3. Sample configurations

2.4. Simulation parameters

3. Results and discussions

3.1. Gaussian width δs

3.2. Convergence time t C

3.3. Gaussian height ω

3.4. Impact of non-bonded interactions

3.4.1. DGEBA

3.4.2. PACM

4. Conclusion and outlook

Simulation software and data

Footnotes

Appendix 1

Acknowledgements

ORCID iD

Funding

Declaration of conflicting interests

References

3.2. Convergence time $t_{C}$