# Impact of the Force Field on the Calculation of Density and Surface Tension of Epoxy-Resins

Mathilde Orselly, Cécile Richard, Julien Devémy, Agathe Bouvet-Marchand, Alain Dequidt, Cédric Loubat, Patrice Malfreyt

► To cite this version:

Mathilde Orselly, Cécile Richard, Julien Devémy, Agathe Bouvet-Marchand, Alain Dequidt, et al. Impact of the Force Field on the Calculation of Density and Surface Tension of Epoxy-Resins. Journal of Physical Chemistry B, 2023, 127 (11), pp.2617-2628. 10.1021/acs.jpcb.2c09087. [hal-04049687](https://uca.hal.science/hal-04049687)

HAL Id: hal-04049687

https://uca.hal.science/hal-04049687

Submitted on 28 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Impact of the Force Field on the Calculation of Density and Surface Tension of Epoxy-Resins

Mathilde Orselly,† Cécile Richard,† Julien Devémy,‡ Agathe Bouvet-Marchand,† Alain Dequidt,‡ Cédric Loubat,† and Patrice Malfreyt*,‡

† *Specific Polymers, 150 Avenue des Cocardières, 34160, Castries, France*

‡*Université Clermont Auvergne, CNRS, Clermont Auvergne INP, Institut de Chimie de Clermont-Ferrand, F-63000 Clermont-Ferrand, France*

E-mail: Patrice.Malfreyt@uca.fr

## Abstract

The molecular simulation of interfacial systems is a matter of debate because of the choice of many input parameters that can affect significantly the performance of the force field of reproducing the surface tension and the coexisting densities. After developing a robust methodology for the calculation of the surface tension on a Lennard-Jones fluid, we apply it with different force fields to calculate the density and surface tension of pure constituents of epoxy resins. By using the model that best reproduces the experimental density and surface tension, we investigate the impact of composition in mass fraction on uncured epoxy resins and the effects of degree of cross-linking on cured resins.

## 1 Introduction

Thermosetting polymers are matrices of choice in many industrial fields such as aerospace, automotive, construction sectors. For example, they are increasingly used in aircraft manufacturing for weight reduction and energy saving purposes. Epoxy polymers are some of the most prominent thermosetting polymers valued nowadays due to their good adhesion to many substrates, high stiffness, strength, creep resistance and thermal resistance when compared

with thermoplastic polymers.¹,²

Epoxy resins are formed from a liquid solution that evolves irreversibly into a solid material during the curing process through polymerization reactions between epoxy monomers and curing agents (often called hardeners). The resulting atomic scale structure can be assimilated to a giant hyper-branched molecule³ through a three-dimensional cross-linked network. Numerous combinations between resins and curing agents is possible, but the formulation of these resins is often guided by the target property through a mostly empirical approach.

The knowledge of the interfacial properties of epoxy resins is critical to understand and predict the adhesion properties of these resins in composite systems. A key-property is the surface tension, but very few data of surface tension are available in the literature.⁴,⁵ This lack of surface tension data can be explained in large part by the complexity of the nature of these multi-component materials including binder, hardener and additives and an increasing viscosity due to the cross-linking reaction. The experimental measurements⁶ of an epoxy mixture formed by diglycidyl ether of bisphenol-A (DGEBA) diluted with diglycidyl ether of 1,4 butanediol (DGE) at a weight ratio of 80/20 (DGEBA/DGE) show a weak decrease of the surface tension (γ) from about 41.5 to 38.5 mN m⁻¹ over a temperature range of 40 K. Other experiments⁷ report the surface tensions
of DGEBA prepolymers with different degrees of polymerization and observe that the surface tension increases by 1.3 mN m⁻¹ as the degree of polymerization increases from 0.03 to 0.47. Another DGEBA-based epoxy system constituted by a DGEBA prepolymer modified with polypropylene glycol diglycidyl ether and mixed with an isophorone diamine hardener shows a decreasing function⁵ of $γ = 37.5 - 0.05 \times T$ with temperature. An experimental work⁸ measured the surface tensions of the epoxy resin and the hardener components to conclude that the surface tension of the curing epoxy resin is not necessarily linked to those of the uncured epoxy resin and hardener. This study also reported very little changes in the surface tension as the curing reaction proceeds.

Since the interfacial tension is required for the end-use properties of these multi-component materials, it is possible to use molecular simulations that can compensate for this lack of surface tension. One of the major benefits of molecular simulations is the ability to rationalize and interpret macroscopic properties from molecular interactions.⁹ Indeed, the order of magnitude of surface tension depend on the strength of interactions. For example at 298 K, the surface tension of the *n*-pentane is $15 \text{ mN m}^{-1}$¹⁰ whereas it is equal to $72 \text{ mN m}^{-1}$ for water.¹¹⁻¹³ The strong electrostatic interactions in water characterized by the formation of hydrogen bonds are in the origin of this higher surface tension, but the surface tension of water looses about $20 \text{ mN m}^{-1}$ over the 300..400 K temperature range. The increasingly realistic hope is then to design the resin with the targeted property based on the nature of the interactions. The reader is redirected to ref. 14 for a comprehensive review about the prediction of thermo-mechanical properties of curing thermoset polymers. A number of molecular models¹⁴ has been applied to the simulation of polymers, knowing that the performance and quality of the model may vary depending on the property to be calculated and the systems to be modelled.

Concerning epoxy systems, extensive research in molecular simulations was found in resins resulted from the cross-linking be-

tween diglycidyl ether of bisphenol F (DGEBF) and diethyl toluenediamine (DETDA),¹⁵⁻¹⁹ triethylenetetramine (TETA)²⁰⁻²³ or diethylenetriamine (DETA),²⁴ especially in the fields of modern aeronautics. The crosslinking of the diglycidyl ether of bisphenol A (DGEBA) resin was also widely studied with different hardeners, such as isophorone diamine (IPDA),²⁵⁻²⁷ trimethylene glycol dip-aminobenzoate (TMAB),²⁸ diethyltoluenediamine (DETDA),²⁹⁻³² triethylenetetramine (TETA),²⁹ ethylenediamine (EDA),³³ diaminodiphenyl sulfone (DDS),³⁴ methylenedianiline (MDA),³⁵ poly(oxopropylene) diamines (POP),³⁶⁻³⁸ 4,4'-methylenebis (cyclohexylamine) (MCA),³⁷ diethylenetriamine (DETA),³⁹⁻⁴² and polyetheramine JEFFAMINE D-230.⁴³⁻⁴⁶ The quality of the prediction of properties such as density, glass transition temperature, elastic constants, and heat capacity depends mainly on the accuracy of the molecular models but not too much on the methodology used for the calculation of the property.

The same does not necessarily apply to the calculation of the surface tension of liquids.⁴⁷ Of course, the quality of the model contributes to the reproduction of the surface tension, but a set of parameters related to the methodology can affect significantly the performance of the prediction⁹,⁴⁷ of this interfacial property. Indeed, calculating the surface tension of a liquid requires to model a two-phase system with an explicit interface between the vapor and liquid bulk phases, leading then to a heterogeneous system in density along the direction normal to the interface. This density gradient along this direction means that certain assumptions that are valid in a homogeneous system are no longer verified in a heterogeneous environment. More specifically, the truncation of the Lennard-Jones potential is not too problematic in bulk phases, since it is corrected by adding long range corrections to the energy and forces equations. In contrast, a significant number of studies have focused on the long range corrections to be applied to macroscopic properties⁴⁷⁻⁵³ and some of them proposed to apply local specific long range corrections on energy, forces and related properties.⁴⁹,⁵⁰,⁵⁴,⁵⁵ The im-
pact of the truncation of the force and energy$^{56,57}$ was investigated through the deviations observed on surface tension values calculated by Monte Carlo and Molecular Dynamics simulations. Surface tension dependencies were established on system-sizes$^{58–61}$ and some recommendations were proposed to avoid strong system-size effects.$^{61}$ For a comprehensive review about the calculation of the surface tension, the reader is directed to Ref. 47.

Here, we propose to investigate how general force fields such as GAFF2, OPLS, PCFF and CGenFF perform in the prediction of the densities and surface tensions of DGEBA, DGEBU, DGEBF and DGEVA prepolymers and IPDA, DETA, MDA and TEPA hardener liquids (see Figure 1). We extend the investigation to mixtures of polymers and hardeners and cured epoxy resins. Since these force fields differ in the way of truncating the non-bonded interactions, we start the study by establishing a methodology that should allow us to make a quantitative comparison between the different molecular models on both surface tension and density of epoxy resins.

# 2 Computational methods

## 2.1 All-Atom Force fields

The general expression of a force field sums intramolecular and intermolecular energy contributions of a molecular system to yield the energy of the configuration. The intramolecular interactions (see Eq. 1) normally consists of bonds stretching, angle bending, dihedral and improper torsions and non-bonded energy contributions. We used the All-Atom (AA) version of the classical force fields indicating that all the atoms of the molecules are treated explicitly. In the class I force fields such as General Amber Force Field (GAFF2),$^{62–64}$ Chemistry at Harvard Macromolecular Mechanics General Force Field (CGenFF),$^{65–68}$ Optimized Potentials for Liquid Simulations (OPLS)$^{69,70}$ molecular models, the bond stretching and angle bending energy are described by harmonic potentials. In class II force fields such as Polymer Consis-

tent Force Field (PCFF),$^{71}$ cross terms represented by $U_{coupling}$ in Eq. 1 can be added to the intramolecular interactions to model couplings between stretching, bending and torsion.

$$ U_{\text{intra}} = U_{\text{bonds}} + U_{\text{angles}} + U_{\text{torsions}} + U_{\text{impropers}} + U_{\text{coupling}} + U_{\text{nb}} \quad (1) $$

The non-bonded interactions can occur between atoms in the same molecule thus contributing to intramolecular interactions and between atoms of different molecules giving rise in that case to intermolecular interactions (see Eq. 2). In most cases, these intramolecular non-bonded interactions take place between atoms separated by more than three bonds, with a possible scaling of the 1-4 interactions whose value depends on the force field.$^{65,69}$

$$ U_{\text{inter}} = U_{\text{nb}} = U_{\text{elect}} + U_{\text{LJ}} \quad (2) $$

As shown by Eq. 2, the intermolecular interactions result from repulsion-dispersion and electrostatic interactions represented by Lennard-Jones and Coulomb potentials, respectively.

The electrostatic interactions between two charges $q_i$ and $q_j$ at a distance $r_{ij}$ is represented by the Coulomb potential as

$$ U_{\text{elect}} = \sum_{i=1}^{N-1} \sum_{j=i+1}^{N} S(r_{ij}) \frac{q_i q_j}{4\pi \epsilon_0 r_{ij}} \quad (3) $$

where $\epsilon_0$ represents the dielectric constant of the vacuum. These electrostatic interactions can be calculated by using the Ewald$^{72}$ and PPPM$^{73}$ methods. The latter consists of making an approximate calculation of the reciprocal space contribution by mapping the system on a mesh and using fast Fourier transform.

The van der Waals interactions are mostly described by the 12-6 Lennard-Jones (LJ) potential

$$ U_{\text{LJ}}^{12,6} = \sum_{i=1}^{N-1} \sum_{j=i+1}^{N} S(r_{ij}) 4\epsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right] \quad (4) $$

where $\epsilon_{ij}$ and $\sigma_{ij}$ correspond to the energy
parameter of the interaction and the Lennard-Jones core diameter. The $S(r_{ij})$ function in Eq. 3 and Eq. 4 indicates how are truncated the interactions. For the GAFF2, OPLS and PCFF force fields, the interactions are truncated by using a cutoff radius $r_c$ such as

$$S(r_{ij}) = \begin{cases} 1, & \text{if } r_{ij} \le r_c \\ 0, & \text{if } r_{ij} > r_c \end{cases} \qquad (5)$$

For the CGenFF model, the electrostatic and LJ interactions are modified by a switching function in order to make energy and forces equations continuous at the cutoff radii $r_{c1}$ and $r_{c2}$. The switching function used in the CGenFF force field is defined by

$$S(r_{ij}) = \begin{cases} 1, & \text{if } r_{ij} \le r_{c1} \\ \frac{(r_{c2}^2 - r_{ij})^2 (r_{c2}^2 + 2r_{ij}^2 - 3r_{c1}^2)}{(r_{c2}^2 - r_{c1}^2)^3}, & \text{if } r_{c1} < r_{ij} \le r_{c2} \end{cases} \qquad (6)$$

In the PCFF model, the repulsion term $r_{ij}^{12}$ of the LJ potential is replaced by a softer repulsion term $r_{ij}^9$ as in the following equation

$$U_{\text{LJ}}^{9,6} = \sum_{i=1}^{N-1} \sum_{j=i+1}^{N} S(r_{ij}) \varepsilon_{ij} \left[ 2 \left( \frac{\sigma_{ij}}{r_{ij}} \right)^9 - 3 \left( \frac{\sigma_{ij}}{r_{ij}} \right)^6 \right] \qquad (7)$$

The LJ parameters for the interactions between unlike sites in the CGenFF and GAFF2 force fields were calculated by using the Lorentz-Berthelot combing rules$^{74,75}$ which are geometric for $\varepsilon_{ij} = (\varepsilon_{ij}\varepsilon_{jj})^{1/2}$ and arithmetic for $\sigma_{ij} = 1/2(\sigma_{ii} + \sigma_{jj})$. In contrast, OPLS model uses a geometric mean mixing rule$^{76}$ for both $\sigma$ and $\epsilon$ whereas PCFF applies specific Waldman-Hagler rules$^{77}$ defined as follows for $\sigma_{ij}$ and $\varepsilon_{ij}$, respectively.

$$\sigma_{ij} = \left( \frac{\sigma_{ii}^6 + \sigma_{jj}^6}{2} \right)^{1/6} \qquad (8)$$

$$\varepsilon_{ij} = 2 \left( \frac{\sigma_{ii}^3 \sigma_{jj}^3}{\sigma_{ii}^6 + \sigma_{jj}^6} \right) \sqrt{\varepsilon_{ij} \varepsilon_{jj}} \qquad (9)$$

## 2.2 Dispersion term

Since the second term in the LJ potential is short-ranged, the dispersion interactions are

normally truncated beyond a cutoff radius $r_c$. The truncation of the LJ potential has no impact on energy and structural properties of homogeneous systems due to the compensation of effects in an isotropic environment. This is no longer valid for heterogeneous systems, especially at the interfacial region characterized by a strong gradient of concentration.$^{56,57}$ One solution widely used is to add long-range corrections to the calculated thermodynamic properties to compensate for the missing long-range part of the potential. An alternative to the calculation of the long-range corrections is to include them explicitly by applying the Ewald summation method to the following dispersion term.$^{48,51,52,78,79}$

$$U_{\text{dispersion}} = - \sum_{i=1}^{N-1} \sum_{j>i}^{N} \frac{B_{ij}}{r_{ij}^6} \qquad (10)$$

For example, Eq.(11) expresses the force on the atom $i$ due to the dispersion term of Eq.(10)

$$\mathbf{f}_i = \sum_{j=1}^{N} B_{ij} \left( 6 + 6\beta^2 r_{ij}^2 + 3\beta^4 r_{ij}^4 + \beta^6 r_{ij}^6 \right) \times \exp(-\beta^2 r_{ij}^2) \frac{\mathbf{r}_{ij}}{r_{ij}^8} + \frac{\pi^{1/2}}{12V} \operatorname{Im} \left( \sum_{\mathbf{h} \neq \mathbf{0}} ib_i \exp(-i\mathbf{h} \cdot \mathbf{r}_i) h^3 \right) \times \left[ \pi^{1/2} \operatorname{erfc}(b) + \left( \frac{1}{2b^3} - \frac{1}{b} \right) \exp(-b^2) \right] S_6(\mathbf{h})\mathbf{h} \qquad (11)$$

where $\beta$ is the Ewald parameter for dispersion interactions. $b$ is defined as $|\mathbf{k}|/2\beta$ where $\mathbf{k}$ represents the vectors from the discrete $2\pi\mathbf{n}/L$ with $L$ the length of the box vectors. $h = |\mathbf{h}|$ is the reciprocal lattice vector. $S_6(\mathbf{h})$ is a complex number describing the structure factor and defined as

$$S_6(\mathbf{h}) = \sum_j b_j \exp(-i\mathbf{h} \cdot \mathbf{r}_j) \qquad (12)$$

When PPPM is used, the expression of the force $\mathbf{f}_i$ can be found elsewhere.$^{51}$
## 2.3 Surface tension

The surface tension $\gamma$, originally given by Kirkwood and Buff,⁸⁰ is defined by

$$ \gamma = \frac{1}{2} \langle p_N - p_T \rangle L_z \quad (13) $$

where $p_N$ and $p_T$ are the normal and tangential components of the pressure and $L_z$ is the length of the simulation cell in the $z$ direction. Since a two phase simulation (see Figure 1.i) with periodic boundary conditions consists of two interfaces, the surface tensions calculated from Eq.(13) is divided by 2 to calculate $\gamma$ for a single interface. For a planar interface, $p_N$ is given by $p_{zz}$, whereas the tangential component $p_T$ is given by $\frac{1}{2}(p_{xx} + p_{yy})$.

Irving and Kirkwood⁸⁰⁻⁸⁵ have shown that $\gamma$ can be calculated from $p_N(z)$ and $p_T(z)$, the components of the pressure tensor as a function of $z$;

$$ \gamma = \frac{1}{2} \int_{-L_z/2}^{L_z/2} (p_N(z) - p_T(z)) dz \quad (14) $$

This local definition of $\gamma$, based on the mechanical road, uses the force acting across a unit area in the $z$-plane for one interface. There is no unique way of calculating the forces across a particular area, since it is unclear which atoms contribute to this force. This has no effect on $p_N(z)$ but different choices of the contour can affect the definition of $p_T(z)$. However, these choices have no effect on the integral in Eq.(14). We use here the Harasima definition⁸⁰,⁸²,⁸⁶ for the normal and tangential pressure components. The normal component of the pressure tensor can be written as

$$ \begin{aligned} p_N(z) = & \langle \rho(z) \rangle k_B T \\ & + \frac{1}{2 A \Delta z} \left\langle \sum_{i=1}^{N-1} \sum_{j=i+1}^{N} \left( z_{ij} (\mathbf{f}_{ij})_z \right) \right. \\ & \left. \times (\delta(z - z_i) + \delta(z - z_j)) \right\rangle \end{aligned} \quad (15) $$

where $A$ is the surface area, $\Delta z$ is the thickness of the slab and $\rho(z)$ is the local number density. The simulation box is divided into $N_z$

slabs of thickness $\delta z$. $\mathbf{f}_{ij}$ is the force between atoms $i$ and $j$ defined as:

$$ \mathbf{f}_{ij} = - \frac{\mathbf{r}_{ij}}{r_{ij}} \frac{d\mathcal{U}(r_{ij})}{dr_{ij}} \quad (16) $$

where $\mathcal{U}$ represents all the intramolecular and intermolecular energy contributions described in Eqs.(1) and (2). The tangential component of the pressure tensor is then expressed as follows

$$ \begin{aligned} p_T(z) = & \langle \rho(z) \rangle k_B T \\ & + \frac{1}{4 A \Delta z} \left\langle \sum_{i=1}^{N-1} \sum_{j=i+1}^{N} \left( x_{ij} (\mathbf{f}_{ij})_x + y_{ij} (\mathbf{f}_{ij})_y \right) \right. \\ & \left. \times (\delta(z - z_i) + \delta(z - z_j)) \right\rangle \end{aligned} \quad (17) $$

These equations indicate that half of the virial contributions arising between atom $i$ and atom $j$ is assigned to the slab where $i$ is located and the other half to the slab where $j$ is located. For CGenFF, since intermolecular energies and forces are modified by a switching function that makes energy and forces equations decrease smoothly to zero at the cutoff, no long-range corrections due to the truncated potentials need to be applied to the pressure components and surface tension.

## 2.4 Simulations details

The liquid-vapor interface of the single LJ center was simulated at T = 100 K using an united model with $\sigma = 3.7327$ Å and $\epsilon = 1246$ kJ mol$^{-1}$ and 5000 molecules. The dimensions of the cell are $L_x = L_y = 50$ Å = 13.4 $\sigma$ while $L_z = 300$ Å = 80.4 $\sigma$. These dimensions are significantly larger than those recommended⁵⁶,⁶¹ to avoid any dependence of the surface tension on the system size. The liquid-vapor interface of water, modelled with TIP4P/2005,⁸⁷ was simulated at T = 300 K using an all-atom model. The dimensions of the cell are similar to the ones used with the single LJ center. Likewise with pre-polymers and hardeners, the dimensions of the cell are only changed for $L_z = 180$ Å
and the simulation takes place at 300 K. Each cell contains about 20.000 atoms. The dimension of the liquid phase is thus close to 85 Å. When treating the mixture of pre-polymers and hardeners, a global of 400 molecules is used inside the simulation box. The following mole fractions are used 0.10, 0.20, 0.30, 0.40, 0.50, 0.60 and 0.80. These mole fractions are then computed as a weight fractions to compare with the experimental data.

More generally, the periodic boundary conditions were applied in all three directions. The MD simulations were performed in the constant-NVT statistical ensemble. The Velocity-Verlet integrator was used to integrate the equations of motion using a time step of 1 fs. In order to control the temperature, a Nosé-Hoover thermostat was applied.⁸⁸ The thermodynamic interfacial properties were averaged over 10 independent simulations and the statistical fluctuations were calculated using these 10 superblock averages. The equilibration period of each independent simulation was performed over 1 nanosecond, and data were collected over additional 10 nanoseconds. When the dispersion term of the Lennard-Jones potential is calculated with the PPPM method, we used meshes of size 30 × 30 × 180Å. The Ewald parameter was set to $\beta = 0.2894 \AA^{-1}$.

## 2.5 Experimental details

A common setup, based on pendant drop tensiometry, was used to determine the liquid-vapor surface tension experimentally.⁸⁹ The OCA 50 DATAPHYSICS INSTRUMENT apparatus consists in a light source, a needle and a camera. A known volume of liquid is pushed through the needle in order to obtain a pendant drop, then a digital image is captured. Numerically, the edge of the drop is computed and surface tension is acquired via the Young-Laplace equation :

$$ \gamma(\frac{1}{R_1} + \frac{1}{R_2}) = \Delta P = \Delta P_0 - \Delta \rho g z \quad (18) $$

<figure><img src="image_2.png" /><figcaption>Figure 1: Typical configurations of epoxy polymers a) DGEBA, b) DGEBF, c) DGEBU, d) DGEVA and curing agents e) DETA, f) TEPA, g) IPDA h) MDA and i) a liquid-vapor interface of a DGEBA+IPDA mixture.</figcaption></figure>

# 3 Results and discussions

## 3.1 Methodological issues

In order to avoid any troubles with the methodology and to verify that the thermodynamic equilibrium is respected whatever the form of the model used, we plot in Figure 2 the differences between the normal and tangential components of the pressure tensor calculated as a function of $z$ by using Eqs.(15) and (17), respectively. As expected for a planar interface in the NVT ensemble, the normal and tangential pressures must be equal in vapor and liquid phases, i.e the difference should be zero. Figure 2 confirms that $p_N(z) - p_T(z) = 0$ in the bulk phases. Since $p_T(z)$ should be negative and $p_N(z)$ constant at the interface, the positive peaks of $p_N(z) - p_T(z)$ refers to the tension at the interface. The two peaks are approximately symmetric indicating that our two-phase system is at equilibrium. The surface tension is then measured as a function of $z$ by integrating $p_N(z) - p_T(z)$. As shown in Figure 2, this property is constant in the bulk phases and the contribution for both interphases is the same. The mechanical equilibrium is valid whatever the method used for the calculation of the dispersion term.
<figure><img src="image_3.png" /><figcaption>Figure 2: Profiles of $p_N(z) - p_T(z)$ (left axis) and of its integral defined by $\frac{1}{2}\gamma(z)$ (right axis) for the liquid-vapor of methane at T= 100 K for both spherically truncated and spherically truncated and shifted LJ potential models and for the dispersion term modeled by the PPPM method.</figcaption></figure>

We report in Figure 3 the values of surface tension of methane calculated with different models for the dispersion term for a cutoff radius changing from 8 to 24 Å. As the dispersion term of the LJ potential is truncated at the cutoff radius (see Eqs.(5) and (6)), a long range correction is needed to account for the missing part. This long-range contribution was calculated by using the following expression$^{90}$

$$ \gamma_{\text{LRC}} = \frac{\pi}{2}(\rho_l - \rho_v)^2 \int_0^1 ds \int_{r_c}^{+\infty} dr \coth\left(\frac{2rs}{d}\right) \\ \times \frac{dU_{\text{LJ}}^{12,6}}{dr} r^4 (3s^3 - s) \quad (19) $$

where $\rho_l, \rho_v$ represent the densities of the liquid and vapor phases, respectively. d is an estimation of the thickness of the interface and s is a parameter defined as $s = (z_i - z_j)/r_{ij}$. The values of $\rho_l, \rho_v$ and d are obtained by assuming that the density profile $\rho(z)$ can be fitted to a hyperbolic tangent function of width d.

$$ \rho(z) = \frac{1}{2}(\rho_l + \rho_g) - \frac{1}{2}(\rho_l - \rho_g) \tanh(2(z - z_g)/d) \quad (20) $$

where $z_g$ is the position of the Gibbs dividing surface. Figure S1 of the Supporting Information

shows that the atomistic density profile can be fitted accurately with Eq.(20). With the dispersion term calculated with the switching function of Eq.(6), no long range correction is needed since the potential and force equations decrease smoothly to zero. For the PPPM and Ewald methods, the long range interactions are already explicitly treated in Fourier space and no correction needs to be applied to the surface tension. Figure S2 of the Supporting Information shows the long-range correction to the surface tension calculated with Eq.(19) at different cutoffs along with the intrinsic or short range part of the surface tension. We observe that this tail contribution decreases from about 8 to 1 mN$^{-1}$ as the cutoff increases from 8 to 24 Å. For a standard value of cutoff of 12 Å, the long range correction to the surface tension contributes by 25% to the total value of $\gamma$.

Figure 3 reports the cutoff-dependence of both surface tensions and liquid densities of the LJ fluid when the dispersion term is calculated with a truncated and a truncated and shifted LJ potential but also with the PPPM and Ewald methods. First, the truncated and shifted LJ potential of Eqs. (4) and (6) shows in Figure 3a a strong dependence of $\gamma$ on the cutoff with a value that more than doubles over the 8..24 Å cutoff range. Adding a long-range correction to the surface tension (see Eq.(19)) strongly mitigates the cutoff-dependence of $\gamma$. From $r_c = 16$ Å, we may consider that the surface tension becomes independent of the cutoff. Interestingly, when the dispersion term is calculated by Ewald and PPPM methods, we no longer observe a dependence of the surface tension on the cutoff radius with values that varies by no more than 2% with respect to the limit value obtained with the largest cutoff of 24 Å.

Figure 3b shows how the liquid densities of the liquid-vapor equilibrium changes in the 8..24 Å cutoff range. Since the liquid density cannot be modified by tail contributions, we observe the same cutoff-dependence of liquid densities for truncated LJ potentials. The Ewald and PPPM methods used for treating the dispersion part avoid any dependence of the liquid density on cutoff radius. We next sought to investigate in Figure 4 whether the calculation
<figure><img src="image_4.png" /><figcaption>Figure 3: a) Surface tensions and b) liquid densities of the liquid-vapor equilibrium of methane at T= 100 K calculated at different cutoffs ($r_c$) with several approaches for the dispersion term as indicated in the legend</figcaption></figure>

of electrostatic interactions through the PPPM method in the liquid-vapor interface of water could modify the cutoff-dependence of the surface tension. This is not the case and applying electrostatic interaction via the PPPM method does not change the conclusion obtained with LJ fluids. The tail correction to the surface tension of water was also calculated by Eq.(19).

The main conclusion we can draw from Figures 3 and 4 is that it is possible to avoid an impact of the cutoff radius by applying the PPPM method for the calculation of the dispersion term. With truncated potentials, we show here

<figure><img src="image_5.png" /><figcaption>Figure 4: Surface tensions of water at T= 300 K calculated with different cutoffs ($r_c$). The dispersion term is calculated using the PPPM method and a standard truncated LJ potential. Water was modeled using the TIP4P/2005 model.⁸⁷</figcaption></figure>

that the values of surface tension modified by the addition of a correction term match very well the values obtained by using the PPPM method provided that you apply a cutoff of 16 Å. Taking into account the prohibitive calculation time required with the PPPM method and the need of modelling molecular systems much more complex than a single LJ center, we therefore retain the option of adding a long range contribution calculated with Eq.19 by using a cutoff of 16 Å to correct the surface tension.

## 3.2 Quality of the force field

### 3.2.1 Prepolymers and hardeners

Before evaluating the performance of the models on surface tension, it may be useful to test it on the density which requires a more straightforward calculation. In addition, for most liquids, we observe that the higher the density the greater the surface tension. This observation suggests a relationship between the density and the surface tension. Macleod⁹¹ established that the surface tension is proportional to the power
four of the difference in the densities between the liquid and vapor phases. However, like any theoretical model, it is difficult to transfer from one family of molecules to another and requires numerous reparametrizations.

Figure 5 shows the correlation between experimental and simulated densities of prepolymers (DGEBA, DGEBF, DGEBU, DGEVA) and curing agents (IPDA, DETA, MDA, TEPA) for the GAFF2, OPLS, PCFF and CGenFF force fields. We can reasonably conclude that the prediction is good when the deviation from the experimental data remains below 5%.

<figure><img src="image_6.png" /></figure>

Figure 5: Correlation between experimental and calculated densities at T=300 K. The densities are calculated with different force fields and for various prepolymers and hardeners as indicated in the legend. The dotted lines represent simulated densities that deviate by ± 5% from experiments. The error bars are less than symbols' size.

The analysis of Figure 5 informs about the ability of the different force fields to reproduce the experimental densities at 300 K. First, the best prediction is obtained with the PCFF force field with an overall average absolute deviation of 2.7%. The following models OPLS, CGenFF, GAFF2 reproduce the experimental densities within overall absolute deviations of 4.5, 4.7 and 7.2 %. The worst prediction is obtained with MDA for which no model is able to predict the density within 9%. For the DETA prepoly-

mer, only the PCFF model successfully reproduces the density at less than 5%. Finally, we observe that the maximum deviation from experiments is then less than 10% for all considered resins. From this comparison of densities with experiments, we can deduce that PCFF is the most successful model for all the polymers and hardeners used here. As a result, we propose to investigate its transferability on the surface tension.

Figure 6 shows the liquid-vapor surface tension of pure polymers and hardeners at 300 K and 400 K. The experimental surface tensions are given when available at 300 K. For completeness, Table 1 summarizes the available experimental surface tensions at 300 K along with the simulated surface tensions at both temperatures. The atomistic density profiles of the liquid-vapor interfaces of DGEBA and IPDA are given in Figure S3 in the Supporting Information for completeness.

<figure><img src="image_7.png" /></figure>

Figure 6: Surface tensions of the liquid-vapor equilibrium of pure polymers and hardeners calculated with the PCFF model. The experimental surface tensions are only available at 300 K. These values are either taken from the literature, from our own experiments or deduced from $-\frac{d\gamma}{dT}$. The origin of these data is specified in Table 1. The simulated surface tensions at 400 K are given for completeness.

Figure 6 illustrates that the comparison with experiments shows significant deviations for
DGEBF (34%) and DGEBU (54%) prepoly-
mers and reasonable agreements within a max-
imum deviation of 10% for IPDA, DETA,
DGEBA and DGEVA systems. In some cases,
such as for DGEBA, we observe deviations of
the order of 25% between the experimental val-
ues alone as shown in Table 1. The worst
prediction of the surface tension observed for
DGEBF and DGEBU cannot be related to a
poor prediction of the density. Indeed, the den-
sities of DGEBF and DGEBU polymers devi-
ate by less of 1% from experiments. This also
shows the difficulty of measuring this property
for epoxy resin. The same can be applied to
the experimental aspect of this work. The var-
ious methodological bias can lead to uncertain-
ties and endangers the quantitative results, but
the simulation could guide new experiments in
this field. We also shown in Figure 6 the val-
ues of γ<sub>LV</sub> calculated at 400 K. These values
inform about the temperature dependence of
the surface tension represented in Table 1 by
the slope dγ/dT. Concerning the temperature
dependence of γ, the molecular simulation pre-
dicts slopes ranging from 0.15 to 0.32 whereas
experiments give coefficients in the range of 0.07
to 0.19 with substantial differences between the
experimental values. Finally, the average devi-
ation on the surface tension is less than 10%
with the PCFF model. The prediction of the
surface tension is much more sensitive in the
same way as the experimental determination.

One of the advantages of molecular model-
ing lies in the ability to highlight a correlation
between the molecular structure and a macro-
scopic property. Figure 7 shows the surface ten-
sion with respect to density for epoxy polymers
and curing agents. The compounds are sepa-
rated into three families. A family composed of
DGEBU, DGEBA, DGEBF and DGEVA poly-
mers, a second one of the three curing agents
IPDA, DETA and TEPA and the third one
composed only of the MDA curing agent. MDA
has been singled out in a family because this
compound contains both an amine function and
two hydroxyl groups. We observe in Figure
7 that the density and surface tension of cur-
ing agents increase with the number of amine
groups in line with the increase in the number

of hydrogen bonds and a more associated liq-
uid. For polymers, the surface tension and den-
sity are highest for the compound that has an
aromatic group, two epoxy groups, three ether
groups and a hydroxyl group. The smallest den-
sity and surface tension was found for DGEBU
that only has two ethers and two epoxy groups
with no aromatic group.

<figure><img src="image_8.png" /><figcaption>Figure 7: Calculated surface tension of poly-
mers and hardeners as a function of the simu-
lated density.</figcaption></figure>

3.2.2 Uncured epoxy resins

After studying the pure compounds of the
epoxy resins, we propose now to investigate the
density and surface tension of mixtures of these
compounds with respect to the composition in
mass fraction. In this section, we do not con-
sider any polymerization. We investigate three
mixtures DGEBA + IPDA, DGEBU + IPDA
and DGEVA+DETA mixtures. The simulated
densities with the PCFF model are represented
in Figure 8a along with available experimental
densities of the pure compounds. We also add
for comparison the density of the mixture calcu-
lated by assuming that the mixture has an ideal
behavior with a zero enthalpy change of mixing
and ΔV = 0. In the case of DGEBA + IPDA
<table><caption>Table 1: Compilation of experimental and simulated data for the liquid-vapor surface tensions of epoxy resins and associated hardeners at 300K. The values marked with * were made experimentally in this work. The subscripts indicate the precision of the last decimal, e.g. 40.5$_{30}$ means 40.5 ± 3.0.</caption><thead><tr><th rowspan="2">Systems</th><th colspan="2">γ<sub>LV,300K</sub> (mN m<sup>-1</sup>)</th><th colspan="2">γ<sub>LV,400K</sub> (mN m<sup>-1</sup>)</th><th colspan="2">-dγ/dT (mN m<sup>-1</sup> K<sup>-1</sup>)</th></tr><tr><th>Exp</th><th>Simu</th><th>Exp</th><th>Simu</th><th>Exp</th><th>Simu</th></tr></thead><tbody><tr><td>DGEBA</td><td>38.1*<br>47.2<sup>92</sup><br>47.2<sup>93</sup><br>46.4<sup>94</sup><br>46<sup>95</sup><br>43.9<sup>a 7</sup><br>36.3<sup>b 5</sup><br>44.1<sub>2</sub><sup>8</sup></td><td>50.7<sub>30</sub></td><td></td><td>34.3<sub>14</sub></td><td></td><td>0.164</td></tr><tr><td>DGEBF</td><td>41<sub>2</sub>*</td><td>55.0<sub>10</sub></td><td></td><td>27.4<sup>a</sup></td><td>0.186<sup>95</sup></td><td></td></tr><tr><td>DGEVA</td><td>50<sub>3</sub>*</td><td>56.0<sub>27</sub></td><td></td><td>31.8<sup>a</sup></td><td>0.098<sup>7</sup></td><td></td></tr><tr><td>DGEBU</td><td>45<sub>2</sub>*<br>28.7<sup>a 7</sup><br>41.9<sup>93</sup></td><td>44.4<sub>4</sub></td><td>22.4<sup>a</sup></td><td>31.3<sup>a,b</sup></td><td>0.05<sup>5</sup></td><td></td></tr><tr><td>IPDA</td><td>34<sub>2</sub>*<br>34.0<sup>a 7</sup></td><td>31.0<sub>4</sub></td><td></td><td>16.8<sub>3</sub></td><td></td><td>0.142</td></tr><tr><td>DETA</td><td>47<sub>4</sub>*<br>42.0<sup>93</sup></td><td>45.8<sub>13</sub></td><td>21<sup>a</sup></td><td>24.6<sub>5</sub></td><td>0.106<sup>7</sup></td><td>0.212</td></tr><tr><td>TEPA</td><td>—</td><td>61.6<sub>42</sub></td><td></td><td>29.3<sub>4</sub></td><td></td><td>0.323</td></tr><tr><td>MDA</td><td>—</td><td>48.4<sub>38</sub></td><td></td><td>25.3<sub>8</sub></td><td></td><td>0.231</td></tr></tbody></table>

mixture, we have the following relationship

$$ \rho_{\text{mixture}} = \frac{m_{\text{DGEBA}} + m_{\text{IPDA}}}{V_{\text{DGEBA}} + V_{\text{DETA}}} = \frac{1}{\frac{\chi_{\text{DGEBA}}}{\rho_{\text{DGEBA}}} + \frac{\chi_{\text{IPDA}}}{\rho_{\text{IPDA}}}} \quad (21) $$

where $\chi_{\text{DGEBA}}$ is the mass fraction of DGEBA component in the mixture and $\rho_{\text{DGEBA}}$ is the density of the pure DGEBA. These calculated densities are reported in Figure 8a with dotted lines. First, Figure 8 shows that the simulated densities of both mixtures match very well with the densities predicted by Eq. 21. Second, the simulated surface tensions of the DGEBA+IPDA mixture in Figure 8b show an increase from 30 to 50 mN m$^{-1}$ with a reasonable agreement with the experimental surface tensions measured in this study. The increase of the surface tension with the composition is well-

reproduced by our model in line with the increase of the density with the mass fraction. In the case of the DGEBA+DETA mixture, if we take the evolution of the density with respect to the composition as a basis, the surface tension is expected to increase. Indeed, we observe an increase of $\gamma_{LV}$ by 10 mN m$^{-1}$ over the composition range. Similarly, when studying the effect of the resin's nature on the mixture, we note an lower increase of $\gamma_{LV}$ for DGEBU+IPDA than for DGEBA+IPDA in tune with the surface tension of the pure resins. The surface tensions of the mixture were proven to be quite difficult to obtain experimentally due to various parameters such as viscosity for instance, simulations could thus aid in determining these surface tensions.

<footnote>a estimated values from -dγ/dT</footnote>

<footnote>b mixed with polypropylene glycol diglycidyl ether</footnote>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: a) Simulated densities and b) surface tensions obtained with the PCFF force field for DGEBA+IPDA; DGEBU+IPDA and DGEVA+DETA mixtures as a function of the composition in mass fraction in DGEBA, DGEBU and DGEVA prepolymers. The experimental densities of pure components are given for comparison.</figcaption></figure>

### 3.2.3 Cured epoxy resins

We now investigate the properties of a cured epoxy resin composed of DGEBA and curing agent IPDA as a function of the degree of cross-linking. We calculate the conversion rate or degree of cross-linking as the percent conversion which gives the ratio between the number of cure reactions carried out and the maximum number of theoretically possible reactions. We

restrict the degree of cross-linking to 50% because beyond this value, we observe a solidification (gelation/vitrification) at 300 K during the curing process through the analysis of the local surface tension (see Figure S4 of the Supporting Information). This investigation of the impact of degree of cross-linking is possible by molecular simulation but is not experimentally controlled. The methodology used for building a three dimensional cross-linked polymer network has been published in ref. 96. Figure 9a shows the density of the cured DGEBA+IPDA epoxy resin as a function of the conversion rate at 300 K and 400 K. For both temperatures, molecular simulations reproduce the increase of density with increasing degree of crosslinking as expected from the volume shrinkage due to the increasing formation of covalent bonds. Over the range 0..50 % of conversion rate, the density increases by about 1% and 5% at 300 K and 400 K, respectively. This is the order of magnitude reported in the literature<sup>17,97,98</sup> for other epoxy resins. Additionally, a plateau is observed for a conversion rate of 20% at 300 K whereas a monotonous increase of density upon increasing degree of cross-linking is found at 400 K, confirming a solidification process at the lowest temperature.

We now turn to the study of the liquid-vapor surface tension of the cross-linked epoxy polymers formed by DGEBA and IPDA. We do not retain the values of $\gamma_{LV}$ at 300 K since the resins exhibit the first steps of a solidification. Indeed, in the case of liquid-vapor equilibrium, the local surface tension defined by $\gamma(z) = \frac{1}{2} \int_{-L_z/2}^{z} (p_N(z) - p_T(z))dz$ should be constant in the bulk liquid phase. Figure S4a of the Supporting information shows that $\gamma(z)$ is no longer constant from a degree of cross-linking of 20%. It means that the simulation is not long enough for the system to relax to its equilibrium state. The solution is either to perform much longer simulations or to perform the simulation at a higher temperature. We choose the second option by running simulations at 400 K. Indeed, Figure S4b of the Supporting Information shows the profiles of the surface tensions at 400 K. These profiles confirm that the DGEBA+IPDA remains liquid
<figure><img src="image_11.png" /><figcaption>Figure 9: a) Simulated densities of the DGBA+IPDA cured epoxy resin at both temperatures as a function of the conversion rate; b) liquid-vapor surface tension of the epoxy resin calculated only at 400 K.</figcaption></figure>

within the range 0..50% of conversion rate. As a result, Figure 9b reports the surface tensions with respect to the cross-linking density. We note an increase of the surface tension of about $7\ \text{mN}\ \text{m}^{-1}$ upon increasing conversion rate corresponding to about an increase of $\gamma_{LV}$ of 25% with respect to the uncured resin. From a theoretical viewpoint, the increase of the surface tension upon increasing cross-linking density is expected from the increase of number of covalent bonds and as a consequence of stronger intramolecular and intermolecular energy contributions.

## 4 Conclusions

Used as part of a robust simulation methodology, the force field is the key element for reproducing or predicting material properties. While the calculation of the density does not require a specific methodology, this is not the case for the calculation of surface tension. Indeed, this interfacial property is sensitive to a number of factors such as the truncation of the energy equation, size-effects and the definition used for the calculation.

For this reason, the first part of the work was devoted to methodological work. We have indeed demonstrated that the PPPM method should be used to model the dispersion term of the Lennard-Jones potential. By using this option on a single LJ center, the simulated surface tensions were shown to be independent of the cutoff radius. However, this alternative is very time-consuming and could be not routinely applied to molecular systems involving van der Waals and electrostatic interactions. The use of a truncated potential with the addition of long range correction to the surface tension approaches the value calculated with the PPPM method from a cutoff radius of 16 Å with a much more reasonable computing time.

We then showed that PCFF was the best force field for reproducing the density of pure liquids formed by prepolymers and hardeners. We then tested the transferability of the PCFF model on the surface tension of different polymers and hardeners by using the methodology developed in this work. In this case, the deviation from experiments can only be assigned to the quality of the model. We obtained an average maximum deviation of 10%, with two polymers (DGEBF and DGEBU) showing dramatic deviations from experiments. These results must also consider that there are many discrepancies between the experimental values due to the complexity of obtaining such properties.

We completed this study by investigating uncured and cured epoxy resins. We studied the effects of increasing the mass fraction of prepolymers in epoxy resins on the density and surface tension properties. The values of surface
tensions were found in a reasonable agreement
with the available experimental data. Finally,
we investigate the impact of the cross-linking
density of a cured DGEBA+IPDA epoxy resin
at a higher temperature to preserve the liq-
uid state of the resin. The surface tension was
found to increase upon increasing the reaction
conversion.

This study confirms molecular simulation as
an essential tools for calculating the surface ten-
sions of epoxy resins where experimental mea-
surements are rare and difficult to perform. In
addition, the PCFF force field is a good candi-
date for the study of bulk and interfacial prop-
erties of these three dimensional polymer net-
works.

## Supporting Information Available

Fit to the LJ density profile by using a tangent
hyperbolic function, importance of the long-
range corrections to the surface tension as a
function of the cutoff radius, density profiles of
the liquid-vapor interface of DGEBA and IPDA
molecules, profiles of the surface tension along
the direction perpendicular to the interface at
300 K and 400 K.

**Acknowledgments** MO, CR, ABM and CL
would like to thank all the members of Simat-
Lab for stimulating discussions about this work.
SimatLab is a joint public-private laboratory
dedicated to the multi-scale modelling of poly-
mer materials. This laboratory is supported
by Michelin, Clermont Auvergne University
(UCA), CHU of Clermont-Ferrand and CNRS.
We are also grateful to the Mesocenter Cler-
mont Auvergne University for providing com-
puting and storage resources.
References

(1) Pascault, J.-P.; Williams, R. J. J. *Epoxy Polymers*; Wiley, 2010; Chapter 1, pp 1–12.

(2) Auvergne, R.; Caillol, S.; David, G.; Boutevin, B.; Pascault, J. P. Biobased Thermosetting Epoxy: Present and Future. *Chem. Rev.* **2014**, *114*, 1080–1115.

(3) Gavrilov, A. A.; Komarov, P. V.; Khafatur, P. G. Thermal Properties and Topology of Epoxy Networks: A Multiscale Simulation Methodology. *Macromolecules* **2015**, *48*, 206–212.

(4) Tyson, W. R.; Miller, W. A. Surface Free Energies of Solid Metals: Estimation From Liquid Surface Tension Measurements. *Surf. Sci.* **1977**, *62*, 267–276.

(5) Page, S. A.; Berg, J. C.; Manson, J. A. A. Characterization of epoxy resin surface energetics Characterization of epoxy resin surface energetics. *J. Adhesion Sci. Technol.* **2001**, *15*, 153–170.

(6) Abbott, J. R.; Higgins, B. G. Surface Tension of a Curing Epoxy. *J. Polym. Sci. Polym. Chem. Ed.* **1988**, *26*, 1985–1988.

(7) Chen, D.; Pascault, J. P. Influence of the Chemical Structure on Surface Energy of Monomers, Comonomers and Additives. *Makromol. Chem.* **1991**, *192*, 867–882.

(8) Synytska, A.; Michel, S.; Pleul, D.; Bellmann, C.; Schinner, R.; Eichhorn, K.-J.; Grundke, K.; Neumann, A. W.; Stamm, M. Monitoring the Surface Tension of Reactive Epoxy-Amine Systems Under Different Environmental Conditions. *J. Adhes.* **2004**, *80*, 667–683.

(9) Allen, M. P.; Tildesley, D. J. *Computer Simulation of Liquids : Second Edition*; Oxford: Clarendon Press, 2017.

(10) Somayajulu, G. R. A Generalized Equation for Surface Tension from the Triple Point to the Critical Point. *Int. J. Thermophys.* **1988**, *9*, 559–566.

(11) Gittens, G. J. Variation of Surface Tension of Water with Temperature. *J. Colloid Interface Sci.* **1968**, *30*, 406–412.

(12) Jasper, J. J. The Surface Tension of Pure Liquid Compounds. *J. Phys. Chem. Ref. Data* **1972**, *1*, 841–1009.

(13) Vasquez, G.; Alvarez, E.; Navaza, J. M. Surface Tension of Alcohol + Water from 20 to 50 °C. *J. Chem. Eng. Data* **1995**, *40*, 611–614.

(14) Li, C.; Strachan, A. Molecular Scale Simulations on Thermoset Polymers: A Review. *J. Polym. Sci., Part B : Polym. Phys.* **2015**, *53*, 103–122.

(15) Tack, J. L.; Ford, D. M. Thermodynamic and Mechanical Properties of Epoxy Resin DGEBF Crosslinked with DETDA by Molecular Dynamics. *J. Mol. Graph. Model.* **2008**, *26*, 1269–1275.

(16) Varshney, V.; Patnaik, S. S.; Roy, A. K.; Farmer, B. L. A Molecular Dynamics Study of Epoxy-Based Networks: Cross-Linking Procedure and Prediction of Molecular and Material Properties. *Macromolecules* **2008**, *41*, 6837–6842.

(17) Li, C.; Strachan, A. Molecular Simulations of Crosslinking Process of Thermosetting Polymers. *Polymer* **2010**, *51*, 6058–6070.

(18) Nouri, N.; Ziaei-Rad, S. A Molecular Dynamics Investigation on Mechanical Properties of Cross-Linked Polymer Networks. *Macromolecules* **2011**, *44*, 5481–5489.

(19) Bandyopadhyay, A.; Valavala, P. K.; Clancy, T. C.; Wise, K. E.; Odegard, G. M. Molecular Modeling of Crosslinked Epoxy Polymers: The Effect of Crosslink Density on Thermomechanical Properties. *Polymer* **2011**, *52*, 2445–2452.

(20) Fan, H. B.; Yuen, M. M. Material Properties of the Cross-Linked Epoxy Resin Compound Predicted by Molecular Dynamics Simulation. *Polymer* **2007**, *48*, 2174–2178.
(21) Yu, S.; Yang, S.; Cho, M. Multiscale Modeling of Cross-Linked Epoxy Nanocomposites to Characterize the Effect of Particle Size on Thermal Conductivity. *J. Appl. Phys.* **2011**, *110*, 124302.

(22) Choi, J.; Yu, S.; Yang, S.; Cho, M. The Glass Transition and Thermoelastic Behavior of Epoxy-Based Nanocomposites: A molecular Dynamics Study. *Polymer* **2011**, *52*, 5197–5203.

(23) Yu, S.; Yang, S.; Cho, M. Multi-Scale Modeling of Cross-Linked Epoxy Nanocomposites. *Polymer* **2009**, *50*, 945–952.

(24) Schichtel, J. J.; Chattopadhyay, A. Modeling Thermoset Polymers Using an Improved Molecular Dynamics Crosslinking Methodology. *Comput. Mater. Sci.* **2020**, *174*, 109469.

(25) Wu, C.; Xu, W. Atomistic molecular modeling of crosslinked epoxy resin. *Polymer* **2006**, *47*, 6004–6009.

(26) Wu, C.; Xu, W. Atomistic Molecular Simulations of Structure and Dynamics of Crosslinked Epoxy Resin. *Polymer* **2007**, *48*, 5802–5812.

(27) Wang, Z.; Lv, Q.; Chen, S.; Li, C.; Sun, S.; Hu, S. Glass Transition Investigations on Highly Crosslinked Epoxy Resins by Molecular Dynamics Simulations. *Mol. Simul.* **2015**, *41*, 1515–1527.

(28) Lin, P.-H.; Khare, R. Molecular Simulation of Cross-Linked Epoxy and Epoxy-POSS Nanocomposite. *Macromolecules* **2009**, *42*, 4319–4327.

(29) Jeyranpour, F.; Alahyarizadeh, G.; Arab, B. Comparative Investigation of Thermal and Mechanical Properties of Cross-Linked Epoxy Polymers With Different Curing Agents by Molecular Dynamics Simulation. *J. Mol. Graph. Model.* **2015**, *62*, 157–164.

(30) Clancy, T.; Frankland, S.; Hinkley, J.; Gates, T. Molecular Modeling for Calculation of Mechanical Properties of Epoxies with Moisture Ingress. *Polymer* **2009**, *50*, 2736–2742.

(31) Gavrilov, A. A.; Komarov, P. V.; Khatur, P. G. Thermal Properties and Topology of Epoxy Networks: A Multiscale Simulation Methodology. *Macromolecules* **2015**, *48*, 206–212.

(32) Demir, B.; Walsh, T. R. A Robust and Reproducible Procedure for Cross-Linking Thermoset Polymers using Molecular Simulation. *Soft Matter* **2016**, *12*, 2453–2464.

(33) Gavrielides, A.; Duguet, T.; Aufray-Moi, M.; Lacaze-Dufaure, C. Model of the DGEBA-EDA Epoxy Polymer: Experiments and Simulation Using Classical Molecular Dynamics. *Int. J. Polym. Sci.* **2019**, *2019*, 9604714.

(34) Liu, H.; Li, M.; Lu, Z.-Y.; Zhang, Z.-G.; Sun, C.-C.; Cui, T. Multiscale Simulation Study on the Curing Reaction and the Network Structure in a Typical Epoxy System. *Macromolecules* **2011**, *44*, 8650–8660.

(35) Chang, S.-H.; Kim, H.-S. Investigation of Hydroscopic Properties in Electronic Packages Using Molecular Dynamics Simulation. *Polymer* **2011**, *52*, 3437–3442.

(36) Soni, N. J.; Lin, P.-H.; Khare, R. Effect of Cross-Linker Length on the Thermal and Volumetric Properties of Cross-Linked Epoxy Networks: A Molecular Simulation Study. *Polymer* **2012**, *53*, 1015–1019.

(37) Sirk, T. W.; Karim, M.; Khare, K. S.; Lenhart, J. L.; Andzelm, J. W.; Khare, R. Bi-Modal Polymer Networks: Composition-Dependent Trends in Thermal, Volumetric and Structural Properties from Molecular Dynamics Simulation. *Polymer* **2015**, *58*, 199–208.
(38) Jang, C.; Sirk, T. W.; Andzelm, J. W.; Abrams, C. F. Comparison of Crosslinking Algorithms in Molecular Dynamics Simulation of Thermosetting Polymers. *Macromol Theory Simul* **2015**, *24*, 260–270.

(39) Shokuhfar, A.; Arab, B. The Effect of Cross Linking Density on the Mechanical Properties and Structure of the Epoxy Polymers: Molecular Dynamics Simulation. *J. Mol. Model.* **2013**, *19*, 3719–3731.

(40) Arab, B.; Shokuhfar, A.; Ebrahimi-Nejad, S. Glass Transition Temperature of Cross-Linked Epoxy Polymers: A Molecular Dynamics Study. *Proceedings of the International Conference Nanomaterials: Applications and Properties*. 2012; pp 1–4.

(41) Arab, B.; Shokuhfar, A. Molecular Dynamics Simulation of Cross-Linked Epoxy Polymers: the Effect of Force Field on the Estimation of Properties. *J. Nano-Electron. Phys.* **2013**, *5*, 01013–1–01013–5.

(42) Langeloth, M.; Sugii, T.; Böhm, M. C.; Müller-Plathe, F. The glass Transition in Cured Epoxy Thermosets: A Comparative Molecular Dynamics Study In Coarse-Grained And Atomistic Resolution. *J. Chem. Phys.* **2015**, *143*, 243158.

(43) Masoumi, S.; Arab, B.; Valipour, H. A Study of Thermo-Mechanical Properties of the Cross-Linked Epoxy: An Atomistic Simulation. *Polymer* **2015**, *70*, 351–360.

(44) Fan, J.; Anastassiou, A.; Macosko, C. W.; Tadmor, E. B. Molecular Dynamics Predictions of Thermomechanical Properties of an Epoxy Thermosetting Polymer. *Polymer* **2020**, *196*, 122477.

(45) Chowdhury, S. C.; Elder, R. M.; Sirk, T. W.; Gillespie, J. W. Epoxy Resin Thermo-Mechanics and Failure Modes: Effects of Cure and Cross-Linker Length. *Compos. B. Eng.* **2020**, *186*, 107814.

(46) Patil, P. N.; Rath, S. K.; Sharma, S. K.; Sudarshan, K.; Maheshwari, P.; Patri, M.; Praveen, S.; Khandelwal, P.; Pujari, P. K. Free Volumes and Structural Relaxations in Diglycidyl Ether of Bisphenol-A Based Epoxy-Polyether Amine Networks. *Soft Matter* **2013**, *9*, 3589.

(47) Ghoufi, A.; Malfreyt, P.; Tildesley, D. J. Computer Modelling of the Surface Tension of the Gas-Liquid and Liquid-Liquid Interface. *Chem. Soc. Rev.* **2016**, *45*, 1387–1409.

(48) Lopez-Lemus, J.; Alejandre, J. Thermodynamic and Transport Properties of Simple Fluids using Lattice Sums: Bulk Phases and Liquid-Vapour Interface. *Mol. Phys.* **2002**, *100*, 2983–2992.

(49) Janeček, J. Long Range Corrections In Inhomogeneous Simulations. *J. Chem. Phys.* **2006**, *131*, 6264–6269.

(50) Janeček, J.; Krienke, H.; Schmeer, G. Inhomogeneous Monte Carlo simulation of the vapor-liquid equilibrium of benzene between 300 K and 530 K. *Condens. Matter Phys.* **2007**, *10*, 415–423.

(51) Isele-Holder, R. E.; Mitchell, W.; Ismail, A. E. Development and Application of a Particle-Particle Particle-Mesh Ewald Method for Dispersion Interactions. *J. Chem. Phys.* **2012**, *157*, 174107.

(52) Isele-Holder, R. E.; Mitchell, W.; Hammond, J. R.; Kohlmeyer, A.; Ismail, A. E. Reconsidering Dispersion Potentials: Reduced Cutoffs in Mesh-Based Ewald Solvers Can Be Faster Than Truncation. *J. Chem. Theory. Comput.* **2013**, *9*, 5412–5420.

(53) Goujon, F.; Ghoufi, A.; Malfreyt, P.; Tildesley, D. J. Controlling the long-range corrections in atomistic Monte Carlo simulations of two-phase systems. *J. Chem. Theory Comput.* **2015**, *11*, 4575–4585.
(54) Guo, M.; Lu, B. Long Range Corrections to Thermodynamic Properties of Inhomogeneous Systems with Planar Interfaces. *J. Chem. Phys.* **1997**, *106*, 3688–3695.

(55) Guo, M.; Lu, B. C. Y. Long Range Corrections to Mixture Properties of Inhomogeneous Systems. *J. Chem. Phys.* **1998**, *109*, 1134.

(56) Trokhymchuk, A.; Alejandre, J. Computer Simulations of Liquid/Vapor Interface in Lennard-Jones Fluids: Some Questions and Answers. *J. Chem. Phys.* **1999**, *111*, 8510–8523.

(57) Goujon, F.; Malfreyt, P.; Simon, J. M.; Boutin, A.; Rousseau, B.; Fuchs, A. H. Monte Carlo Versus Molecular Dynamics Simulations in Heterogeneous Systems: An Application to The n-Pentane Liquid-Vapor Interface. *J. Chem. Phys.* **2004**, *121*, 12559–12571.

(58) Errington, J. R. Evaluating Surface Tension using Grand-Canonical Transition-Matrix Monte Carlo Simulation and Finite-Size Scaling. *J. Chem Phys.* **2003**, *67*, 012102.

(59) Orea, P.; Lopez-Lemus, J.; Alejandre, J. Oscillatory Surface Tension due to Finite-Size Effects. *J. Chem. Phys.* **2005**, *123*, 114702.

(60) Gonzalez-Melchor, M.; Orea, P.; Lopez-Lemus, J.; Bresme, F.; Alejandre, J. Stress Anisotropy Induced by Periodic Boundary Conditions. *J. Chem. Phys.* **2005**, *122*, 094503.

(61) Biscay, F.; Ghoufi, A.; Goujon, F.; Lachet, V.; Malfreyt, P. Calculation of the Surface Tension from Monte Carlo Simulations: Does the Model Impact on the Finite-Size Effects? *J. Chem. Phys.* **2009**, *130*, 184710.

(62) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Merz, I. R. G. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second

Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. *J. Am. Chem. Soc.* **1995**, *117*, 5179–5197.

(63) Wang, J.; Wolf, R.; Caldwell, J.; Kollman, P.; Case, D. Development and Testing of a General Amber Force Field. *J. Comput. Chem.* **2004**, *25*, 1157–1174.

(64) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. *J. Mol. Graph. Model.* **2006**, *25*, 247–260.

(65) A. D. MacKerell, J.; Bashford, D.; Bellott, M.; R. L. Dunbrack, J.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; W. E. Reiher, I.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. *J. Phys. Chem. B* **1998**, *102*, 3586–3616.

(66) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; Mackerell, A. D. CHARMM General Force Field: A Force Field For Drug-Like Molecules Compatible with The CHARMM All-Atom Additive Biological Force Fields. *J. Comput. Chem.* **2010**, *31*, 671–690.

(67) Vanommeslaeghe, K.; MacKerell, A. D. Automation of the CHARMM General Force Field (CGenFF) I: Bond Perception and Atom Typing. *J. Chem. Inf. Model.* **2012**, *52*, 3144–3154.

(68) Vanommeslaeghe, K.; Raman, E. P.; MacKerell, A. D. Automation of the CHARMM General Force Field
(CGenFF) II: Assignment of Bonded Parameters and Partial Atomic Charges. *J. Chem. Inf. Model.* **2012**, *52*, 3155–3168.

(69) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. *J. Am. Chem. Soc.* **1996**, *118*, 11225–11236.

(70) Jorgensen, W. L.; Tirado-Rives, J. Potential Energy Functions for Atomic-Level Simulations of Water and Organic and Biomolecular Systems. *Proc. Natl. Acad. Sci. U. S. A.* **2005**, *102*, 6665–6670.

(71) Sun, H. Ab Initio Calculations And Force Field Development For Computer Simulation of Polysilanes. *Macromolecules* **1995**, *28*, 701–712.

(72) Ewald, P. P. Die Berechnung Optischer und Elektrostatischer Gitterpotentiale. *Ann. Phys.* **1921**, *369*, 253–287.

(73) Eastwood, J. W.; Hockney, R. W.; Lawrence, D. N. P3M3DP - the 3-Dimensional Periodic Particle-Particle-Particle-Mesh Program. *Comput. Phys. Commun.* **1980**, *19*, 215–261.

(74) Lorentz, H. A. Ueber die Anwendung des Satzes vom Virial in der Kinetischen Theorie der Gase. *Ann. Phys.* **1881**, *12*, 127–136.

(75) Berthelot, D. C. R. Sur le Mélange des Gaz. *Hebd. Acad. Sci.* **1898**, *126*, 1703–1855.

(76) Good, R. J.; Hope, C. J. New Combining Rule for Intermolecular Distances in. Intermolecular Potential Functions. *J. Chem. Phys.* **1970**, *53*, 540–543.

(77) Waldman, M.; Hagler, A. T. New Combining Rules for Rare Gas van der Waals Parameters. *J. Comput. Chem.* **1993**, *14*, 1077–1084.

(78) in 't Veld, P. J.; Ismail, A. E.; Grest, G. S. Application of Ewald Summations to Long-Range Dispersion Forces. *J. Chem. Phys.* **2007**, *127*, 144711.

(79) Tameling, D.; Springer, P.; Bientinesi, P.; Ismail, A. E. Multilevel Summation for Dispersion: A Linear-Time Algorithm for $r^{-6}$ Potentials. *J. Chem. Phys.* **2014**, *140*, 024105.

(80) Kirkwood, J. G.; Buff, F. P. The Statistical Mechanical Theory of Surface Tension. *J. Chem. Phys.* **1949**, *17*, 338–343.

(81) Irving, J. H.; Kirkwood, J. The Statistical Mechanical Theory of Transport Processes IV. The Equations of Hydrodynamics. *J. Chem. Phys.* **1950**, *18*, 817–829.

(82) Rowlinson, J. S.; Widom, B. *Molecular Theory of Capillarity*; Clarendon Press: Oxford, 1982.

(83) Walton, J. P. R. B.; Tildesley, D. J.; Rowlinson, J. S.; Henderson, J. R. The Pressure Tensor at the Planar Surface of a Liquid. *Mol. Phys.* **1983**, *48*, 1357–1368.

(84) Walton, J. P. R. B.; Tildesley, D. J.; Rowlinson, J. S.; Henderson, J. R. The Pressure Tensor at the Planar Surface of a Liquid. *Mol. Phys.* **1983**, *50*, 1381, erratum.

(85) Shi, K.; Santiso, E. E.; Gubbins, K. E. Can We Define a Unique Microscopic Pressure In Inhomogeneous Fluids ? *J. Chem. Phys.* **2021**, *154*, 084502.

(86) Harasima, A. *Advances in Chemical Physcis*; Wiley & Sons, 1958; pp 203–237.

(87) Abascal, J. L. F.; Vega, C. A General Purpose Model For The Condensed Phases of Water: TIP4P/2005. *J. Chem Phys.* **2005**, *123*, 234505.

(88) Nosé, S. A Molecular Dynamics Method for Simulations in The Canonical Ensemble. *Mol. Phys.* **1984**, *52*, 255–268.
(89) Berry, J. D.; Neeson, M. J.; Dagastine, R. R.; Chan, D. Y.; Tabor, R. F. Measurement of surface and interfacial tension using pendant drop tensiometry. *J. Colloid Interface Sci.* **2015**, *454*, 226–237.

(90) Blokhuis, E. M.; Bedaux, D.; Holcomb, C. D.; Zollweg, J. A. Tail Corrections to the Surface Tension of a Lennard-Jones Liquid-Vapour Interface. *Molec. Phys.* **1995**, *85*, 665–669.

(91) Macleod, D. B. Relation Between Surface Tension and Density. *Trans. Faraday Soc.* **1923**, *19*, 38.

(92) Sharpe, L. H.; Schonhorn, H. Surface Ennergetics, Adhesion, and Adhesive Joints. *Advances in Chem. Series 43* **1964**, 189–201.

(93) Dearlove, T. J. Diluted Epoxy Adhesives. I. Physical Properties and Lap Shear Strengths. *J. Appl. Polym. Sci.* **1978**, *22*, 2509–2521.

(94) Foister, R. T. Dynamic Surface Properties Due to Amine Migration And Chemical Reaction in Primary Amine/Epoxide Systems. *J. Colloid. Interface Sci.* **1984**, *99*, 568–585.

(95) Sohn, J. E.; Emerson, J. A.; Thompson, P. A.; Koberstein, J. T. Rubber-modified epoxies: Interfacial tension and morphology. *J. Appl. Polym. Sci.* **1989**, *37*, 2627–2636.

(96) Orselly, M.; Devemy, J.; Bouvet-Marchand, A.; Dequidt, A.; Loubat, C.; Malfreyt, P. Molecular Simulations of Thermomechanical Properties of Epoxy-Amine Resins. *ACS Omega* **2022**, *7*, 30040–30050.

(97) Yarovsky, I.; Evans, E. Computer Simulation of Structure and Properties of Crosslinked Polymers: Application to Epoxy Resins. *Polymer* **2022**, *43*, 963–969.

(98) Guseva, D. V.; Rudyak, V. Y.; Kamarov, P. V.; Sulimov, A. V.; Bulgakov, B. A.; Chertovich, A. V. Crosslinking Mechanisms, Structure and Glass Transition in Phthalonitrile Resins: Insight from Computer Multiscale Simulations and Experiments. *J Polym Sci B Polym Phys* **2018**, *56*, 362–374.
TOC Graphic

<figure><img src="image_12.png" /></figure>