Adhesion of Epoxy Resin with Hexagonal Boron Nitride and Graphite.

Adhesion interaction of epoxy resin with the basal surfaces of h-BN and graphite is investigated with the first-principles density functional theory calculations in conjunction with the dispersion correction. The h-BN/epoxy and graphite/epoxy interfaces play an important role in producing nanocomposite materials with excellent thermal dissipation properties. The epoxy resin structure is simulated by using four kinds of fragmentary models. Their structures are optimized on the h-BN and graphite surfaces after an annealing simulation. The distance between the epoxy fragment and the surface is about 3 Å. At the interface between h-BN and epoxy resin, no H-bonding formation is observed, though one could expect that the active functional groups of epoxy resin, such as hydroxyl (-OH) group, would be involved in a hydrogen-bonding interaction with nitrogen atoms of the h-BN surface. The adhesion energies for the two interfaces are calculated, showing that these two interfaces are characterized by almost the same strength of adhesion interaction. To obtain the adhesion force-separation curve for the two interfaces, the potential energy surface associated with the detachment of the epoxy fragment from the surface is calculated with the help of the nudged elastic band method and then the adhesion force is obtained by using either the Morse-potential approximation or the Hellmann-Feynman force calculation. The results from both methods agree with each other. The maximum adhesion force for the h-BN/epoxy interface is as high as that for the graphite/epoxy interface. To better understand this result, a force-decomposition analysis is carried out, and it has been disclosed that the adhesion forces working at both interfaces mainly come from the dispersion force. The trend of increase in the C 6 parameters used for the dispersion correction for the atoms included in the h-BN or graphite surface is in the order: N < C < B, which reasonably explains why the strengths of the dispersion forces operating at the two interfaces are similar. Also, the electron localization function analysis can explain why the h-BN surface cannot form an H bond with the hydroxyl group in epoxy resin.

Introduction

Polymer nanocomposites technology has found a broad range of applications in industry, including structural materials,[1] coating,[2] medical products,[3] and electronic devices.[4] They are generated by mixing two components, namely, polymer matrix and nanofiller. Epoxy resin is one of the most commonly used polymer matrices of nanocomposites because of its high modulus, chemical resistance, and good processability.[5] When nanofillers are homogeneously dispersed in the epoxy matrix, the mechanical, optical, thermal, electronic, or magnetic properties of the resin will be improved dramatically.[4,6] Though researchers speculate that the size, shape, aspect ratio, chemical composition, and surface modification of the fillers should have a significant impact on the properties of nanocomposite materials, the detailed mechanism of how the filler properties influence the properties of the epoxy matrix through the interaction between them remains unclear.[6] Generally, the compatibility between the host matrix and the filler can be regarded as a key factor that should be taken into account when one designs a new composite material. Thus, the community has been active in scrutinizing the interfacial interaction between the constitutive phases of nanocomposite materials.[7,8]

In this paper, we present a comparative study on the interfacial interaction in two types of epoxy nanocomposites: one uses hexagonal boron nitride (h-BN) as a nanofiller and the other adopts graphite or graphene layers.

The h-BN/epoxy composites can pave the way for thermally conductive composites with excellent electronically insulating behavior.[9] Owing to the recent advancement in the field of the electronic information technology, there has been a constant demand for heat-dissipating polymer materials that can be applied to the substrate or package for integrated electronic devices that emit large amount of heat.[10] So high thermal conductivity and low electrical conductivity are the necessary conditions. And, the h-BN/epoxy composites can meet the requirements.[11,12] In addition, by the addition of h-BN, one can also improve the mechanical properties of epoxy resin, such as elastic modulus, ultimate tensile strength, and toughness.[6]

For the case of graphite/epoxy composites (or graphene/epoxy composites), similarly to the h-BN/epoxy composites, one can observe the improvement of thermal conductivity[13,14] and mechanical strength.[15,16] An important difference between the two composites is that the graphite/epoxy composites show an increasing electrical conductivity as the filler load increases.[17] Besides, the electromagnetic interference shielding behavior[18] and electric heating behavior[19] emerge in the composites.

To enhance the thermal conductivity of h-BN/epoxy and graphite/epoxy nanocomposites, one needs to improve the dispersion and affinity between the epoxy matrix and the nanofillers. For this purpose, it is indispensable to understand the dominant interacting force between them on a molecular scale. Herein, we intend to provide insight into the adhesive interaction at the graphite/epoxy and h-BN/epoxy interfaces in nanocomposites by means of first-principles quantum chemical calculations.

Recently, the theoretical modeling of the interface between adhesive and adherend has become an important research area because it can provide us with a molecular-level insight into the interfacial adhesion interaction, which is of significant importance in industrial applications of adhesives, such as an anticorrosive coating for metallic surfaces. Lee et al. studied the interactions of epoxy resin with pure iron and Fe2O3 surfaces using the density functional theory (DFT) computation,[20,21] and Ramezanzadeh et al. studied the interaction of epoxy resin with various iron oxides, cerium oxide, lanthanum oxide, and praseodymium oxide using a molecular dynamics simulation.[22−25]

Both h-BN and graphite are two-dimensional materials consisting of hexagonal ring layers (see Figure ). Although they are isoelectronic and share common structural features, their optical and electronic properties are very different. Owing to the difference in electronegativity between B and N, electrons are expected to be localized on the N atoms, leading to polarized B–N bonds having partial ionic character.[26] That is why h-BN is an insulator. Since the B atoms are positively charged, whereas the N atoms are negatively charged, the B atoms interact with the N atoms directly above and below them in adjacent layers (so-called AA′ stacking mode). The interlayer interactions are likely to come from electrostatic forces as well as van der Waals forces.[27]

Figure 1

Crystal structures of (a) graphite and (b) h-BN.

In graphite, on the other hand, the honeycomb layers are stacked in the AB stacking mode, in which half of the C atoms interact with the C atoms directly above and below them in adjacent layers. The C–C bond in the graphite layer is nonpolar and the overlap between adjacent C 2p orbitals makes π-electrons delocalized over the entire layer. The interlayer bonding can work mainly through van der Waals forces. Also, orbital interactions between the adjacent layers should contribute to the interlayer bonding to some extent.[28] However, the electrostatic forces are not important.

When one considers the interaction of a polymer matrix with h-BN or graphite sheets at the molecular level, it is reasonable to assume the interaction can essentially be explained by functional groups in the side chains or in the polymer–matrix main chain and the basal or (001) surface of h-BN or graphite. Chen and co-workers studied the interactions between graphene sheets and a polymer matrix by measuring the Raman spectra, finding a signature that is consistent with the formation of noncovalent π–π stacking between phenyl rings included in the polymer matrix and the basal planes of graphene.[29] In a molecular dynamics (MD) study on interfacial thermal conductance of h-BN–polymer composites, Luo and co-workers adopted a model of the basal plane of h-BN because they assumed that the majority of the interfaces are between the basal plane of h-BN and the polymer matrices.[30]

Yoshizawa and co-workers have carried out an energy decomposition analysis to an epoxy/graphene interface modeled by a fragmentary structure of epoxy resin and the basal plane of graphene.[31] It was found that the electrostatic and charge-transfer interactions are canceled by the exchange–repulsion interactions, and hence only dispersion interactions work at the interface. Also, the study demonstrated that the dispersion interaction does not vary with the surface type. So, it is reasonable to speculate that the epoxy resin would be more strongly bound by the h-BN surface than to the graphite surface because dispersion interaction is expected to exist at both epoxy/graphene and epoxy/h-BN interfaces, but the electrostatic one would only work at the epoxy/h-BN interface. In this paper, we intend to clarify this point theoretically.

Results and Discussion

Bulk Structure Optimization

The crystallographic data for h-BN and graphite bulk structures were taken from the literature.[32,33] We need to optimize the bulk structures using DFT. However, a general drawback of widely used DFT functionals is that they cannot properly describe dispersion interactions,[34] which should play an important role in determining the interlayer distance of h-BN and graphite. To resolve this problem, one can add a correction term, Edisp, to the conventional Kohn–Sham DFT energy, EDFT, as follows:[35]We checked the reproducibility of the interlayer distances of h-BN and graphite using several dispersion correction methods available in VASP, finding that the Tkatchenko–Scheffler (TS) dispersion correction[36] is a suitable for reproducing the experimentally observed interlayer distances. The detailed results of the benchmark test are summarized in the Supporting Information (SI) of this paper.

Slab Model Construction

To construct slab models for the basal planes of h-BN and graphite, one needs to decide how many layers should be included in the model. For this purpose, we calculated the surface energy, γ, which is the energy needed to cleave the bulk crystal and is obtained from a slab calculation using the following equation:[37]where Esurf is the total energy of the surface slab, Ebulk is the total energy of bulk, and n is the ratio of the number of atoms in bulk, Nbulk, to that in the surface slab, Nslab, so n = Nslab/Nbulk. The surface energy is related to the affinity of surfaces to adsorbates, describing the adhesion properties of a surface to a range of adhesives.[38,39] Also, since the higher the surface energy of a surface, the more energy is gained upon bringing the surface into contact with adhesive materials,[39] the surface energy should be a good measure of the adhesion energy. We then calculated the surface energies of two-monolayer (2-ML), 4-ML, and 6-ML slab models for h-BN and graphite. Their plots are shown in the SI. Owing to the layered nature of the structures and a large interlayer separation, the surface energies obtained are not so large when compared with metal surfaces.[40] The surface energy seems not to be affected by the number of layers when one moves from the 2-ML slab to larger-numbered layers (see the SI). Thus, we adopted the 2-ML slab models of h-BN and graphite.

In our previous studies, we have investigated the effects of surface water on the adhesion properties because we have dealt with the hydrophilic surfaces, such as alumina[41] and silica.[42] Since the h-BN and graphite surfaces are hydrophobic and water slippage has been theoretically predicted on such surfaces,[43] as a first approximation, we neglect surface water molecules in our model.

Modeling of Epoxy Resin

As shown in Figure a, epoxy resin is usually formed from the polymerization of the diglycidylether of bisphenol A (DGEBA), which is the product of the reaction between epichlorohydrin and bisphenol-A.[44] To model epoxy resin, we cut out four different units of the polymer chain from the structure shown in Figure a. The resultant fragment models are shown in Figure b. When the fragment is cut out, we need to break a covalent bond, so newly formed terminal C and O atoms are passivated by H and CH3, respectively.

Figure 2

Chemical structure of (a) epoxy resin and (b) its fragment models.

Determining Adhesion Structures of Fragment Models on h-BN and Graphite Basal Surfaces

Given the size of epoxy fragments, we estimated that a 5 × 5 supercell of the slab model constructed in Section is large enough to accommodate the fragment structure on the surface in the unit cell. A top view of the surface model used is presented in the SI. The polymer fragments were randomly placed above the h-BN and graphite surfaces and then optimized. But there is no guarantee that the structures obtained from the initial optimization correspond to the global minimum. One should notice that the initial configuration of the optimization could affect the extent of attachment to the surface as well as adhesion properties. We are blessed with a wide variety of optimization algorithms, and among them the simulated annealing method is thought to be a powerful tool to find the global minimum structure of complicated systems with high probability.[45,46] Here, to find the energetically most stable adhesion structure, we employed the annealing method combined with ab initio Born–Oppenheimer molecular dynamics, which is implemented in VASP. Starting at 300 K, the temperature was gradually decreased down to 0 K. By scaling the velocities of the atoms, a continuous decrease of the kinetic energy is achieved. The interval between velocity scaling events is simulated in a microcanonical ensemble. For each annealing simulation, a total of 5000 molecular dynamics steps with a time step of 1 fs were conducted. We carried out further static geometry optimizations at 0 K for the final annealed structures.

The adhesion structures of the fragments (1–4) on the h-BN surface obtained by using the annealing simulation were found to be over 10 kcal/mol more stable than before the annealing. On the other hand, though the annealing simulation was carried out two times for the fragments 1 and 2 on the graphite surface, we could not obtain a more stable structure. Thus, we assume that the adhesion structures of the fragments 1 and 2 on the graphite surface have already been located at the global minima. The energies of the structures of the fragments 3 and 4 adhered on the graphite surface became lower than before the annealing by 1.9 and 7.9 kcal/mol, respectively. In Figure , we show the determined adhesion structures of fragments 1, 2, 3, and 4 on the basal planes of h-BN and graphite.

Figure 3

Top and side views of the adhesion structures of the fragments 1, 2, 3, and 4 on the basal plane of (a) h-BN or (b) graphite. Ivory, blue, gray, red, and white balls represent B, N, C, O, and H atoms, respectively. The height of an atom in each epoxy fragment that is the closest to the surface is shown in units of Å. The fragment models are shown by a thick ball-and-stick model, whereas the surfaces are represented as a thin ball-and-stick model. Only the topmost layer is shown.

Our way of regarding each fragment as relaxed freely from other fragments might be an ideal limit. The fragments can be assumed to mutually interact with each other through the covalent bond as well as the intra/intermolecular steric, vdW, and long-range Coulomb interactions. The adhesion strength could be affected substantially by such mutual interactions in the real system, as demonstrated recently by a large-scale DFT-based calculation.[47] However, the primary purpose of this study is to gain a qualitative insight into the governing interaction working at the interface between epoxy resin and the h-BN or graphite surface rather than the interaction working inside the epoxy resin or the one between the fragmentary structures. Thus, we believe our small, freely relaxed fragments will do.

As shown in Figure , an H atom bonded to a C atom is located in the closest proximity to the h-BN or graphite surface in all the adhesion structures. This observation hints at the presence of the so-called CH/π interaction.[48] The distances between the H atom and the surface are close to those found in molecular calculations for CH/benzene (C6H6) or CH/borazine (B3N3H6) complexes.[49,50] Usually, experimentally observed CH/π distances have been found to be shorter than 3.05 Å.[51] Given that the distance between the H atom of the CH bond and the ring of benzene or borazine falls in a range between 2.25 and 2.44 Å,[49] we can expect that there should be the CH/π interaction in the optimized adhesion structures.

Except the fragment 3, the benzene rings in the epoxy fragment are not parallel to the surface so that we cannot expect the existence of the π–π stacking interaction. This is because the benzene rings, in the fragments 1, 2, and 4, are adjacent to the isopropyl group, which imposes steric hindrance and prevents the formation of the π–π stacking interaction with the surface. One benzene ring in the fragment 3 is not adjacent to the isopropyl group, so it can form the π–π stacking interaction[52] and one can see the parallel alignment of the benzene ring to the surface. However, it comes from the artificial fragmentation of the polymer into a monomeric unit. We may not need to take this observation seriously for now.

There may be some ways to analyze the adhesion structures, though they are inherently complicated. In this paper, we adopt two of them: one is the mass density profile and the other is the radial distribution function.

Figure shows the mass density profile normal to the h-BN or graphite (001) surface along the z axis for the four fragment models. The horizontal axis of the graph corresponds to the height measured from the topmost layer of h-BN or graphite. On the left side of the graphs, one can see the tails of the peaks, which correspond to the mass density due to the atoms included in the topmost layer of h-BN or graphite. Of course, the shape of profiles differs from model to model, but a few common features can be found. For example, in both h-BN and graphite cases, the peaks corresponding to the epoxy fragment start to grow at around 2.5 Å and their first high peaks are located between 3 and 4 Å.

Figure 4

Mass density profiles of the fragments 1, 2, 3, and 4 along the z-axis on (a) the h-BN or (b) graphite (001) surface. The bin size of the horizontal axis, which determines the resolution of the profiles, is set to 0.25 Å.

Let us move on to a scrutiny of atom–atom separations at the interface between the adhesive and adherend. For this purpose, we employ the radial distribution function (RDF), g(r), as shown in Figure . Since Figure implies that the variation from fragment to fragment is not so significant, the RDF plots for the adhesion interface shown in Figure are represented by those calculated for the fragment 1/h-BN or the fragment 1/graphite interface. The RDF plots for the other fragments are shown in the SI.

Figure 5

Radial distribution function (RDF), g(r), calculated for a single optimized structure of the fragment 1/h-BN interface (a) and that of the fragment 1/graphite interface (b). In (a), the separations of the C, H, and O atoms in the epoxy fragment from the B or N atoms of the h-BN surface are sampled, whereas in (b), the separations of the C, H, and O atoms in the epoxy fragment from the C atoms of the graphite surface are sampled, where C(g) and C(e) mean the C atoms of graphite and those of the epoxy fragment, respectively. The upper limit of sampled distances is set to 6.5 Å (but here they are shown up to 4 Å). The resolution of the histogram is set to 0.05 Å.

The RDF plots for the B–H, N–H, and C(g)–H atomic pairs are characterized by relatively broaden peaks, which start to grow at around r = 2.5 Å. Here, C(g) means the C atoms of graphite. This observation agrees well with the distances measured in Figure . By contrast, the RDF plots for the B–O, N–O, and C(g)–O atomic pairs are characterized by sharp peaks. This feature may have roots in the fact that there are limited number of O atoms in the structures. The first peaks for the B–O, N–O, and C(g)–O atomic pairs are located between r = 3 and 3.5 Å. Of course, the peaks associated with the separations of the C atoms of the epoxy fragment from the surface atoms are found at further distances.

In the literature,[53−55] one can find some intermolecular N–H and B–O interactions as exemplified by the ones shown in Scheme . An acid–base interaction can be found between the B atom and the O atom of alcohol or ether. Such an interaction is characterized by a relatively short intermolecular distance of around 1.5 Å. A hydrogen bond can be found between the N atom and the hydroxyl group. Such an interaction is characterized by a much shorter intermolecular distance of 1.25 Å.

Scheme 1

Structures of 1:1 Molecular Complexes Bonded by a B–O or N–H Intermolecular Interaction

The B–O and N–H distances are indicated in Angstrom. The structures shown in (a), (b), and (c) are taken from the crystal structures presented in refs[60, 61], and (62), respectively.

Given the intermolecular B–O and N–H distances found in Scheme , one cannot expect any specific intermolecular interactions at the interface between the epoxy fragment and the h-BN(001) surface on the basis of the RDF plots shown in Figure . Also, the RDF plots for the epoxy/h-BN structure look very much akin to those for the epoxy/graphite structure. This observation makes us think that one might not be able to find any significant difference in the adhesion energy or force between the epoxy/h-BN and epoxy/graphite interfaces. This motivated us to theoretically probe the adhesion force acting on the interfaces as well as adhesion energies. This is what we will present in the subsequent sections.

Adhesion Energy

One can gauge the adhesion energy between the epoxy fragments and the h-BN or graphite surface using the following equation:[56,57]where E(adhesive+adherend) is the total energy of the surface–resin complex shown in Figure and Eadhesive and Eadherend are, respectively, the energies of the epoxy fragments and the h-BN or graphite surface optimized separately. Calculated adhesion energies are shown in Table . Mostly, the adhesion energies of the fragments on the h-BN surface are slightly larger than those on the graphite surface, but the difference is very tiny (about 0.1 eV). This difference may be traced back to the difference in charge separation between the h-BN and graphite surfaces. Our Hirshfeld-charge analysis[58] for the surface models shows that the B and N atoms of the h-BN surface bear charges of +0.21|e| and −0.21|e|, respectively, whereas the C atoms of the graphite surface have no charge. Thus, the electrostatic interaction between the epoxy fragments and the h-BN surface is expected to exist, resulting in a little bit larger adhesion energy.

Table 1

Adhesion Energies of the Fragments 1, 2, 3, and 4 on the h-BN or Graphite Surface in Units of eV

	fragment 1	fragment 2	fragment 3	fragment 4
h-BN	1.53	1.57	1.73	1.51
graphite	1.42	1.45	1.76	1.43

As for fragment 3, one can see the opposite trend. The contribution from the π–π stacking to the adhesion energy in the case of fragment 3 might make a difference. We are planning to theoretically design an epoxy adhesive whose affinity to the graphite or h-BN surface is reinforced through π–π stacking.

The adhesion energy can be decomposed into the interaction energy and the distortion energies of the adherend and the surface.[59] Such a scheme of energy decomposition can be traced back to the distortion/interaction model of Houk and co-workers[60,61] or the activation strain model of Bickelhaupt and co-workers.[62,63] The distortion energies of the epoxy fragments and the h-BN or graphite surface are presented in the SI. The distortion energies of the h-BN and graphite surfaces are at most 0.02 eV, whereas those of the fragment models are at most 0.05 eV. This difference may reflect the flexibility of the epoxy fragments and the robust nature of h-BN and graphite. However, these energies can generally be regarded as a tiny number. So, we can expect that the adhesion energies tabulated in Table mainly come from the interaction between the surface and the adherend, and the effects of distortion can be negligible.

Adhesion Force

Here, we turn to how large is the adhesion force during the detachment of the epoxy fragments from the h-BN or graphite surface. In our previous studies, we proposed a way of estimating the adhesion force.[41,56,57] The main points of the methodology are summarized as follows: First, one needs to obtain the interaction potential energy curve with respect to separation between the adhesive molecule and the adherend surface, and then one can convert the potential curve to the adhesion–force curve by using the derivative of the potential curve with respect to the displacement of the adhesive molecule from the position of stable equilibrium.

There may be some ways to probe the potential energy curve associated with the detachment of the adhesive molecule; for example, one could increase the displacement of the adhesive molecule in the direction perpendicular to the surface step by step and perform a single-point calculation or a partial optimization at each step, resulting in a potential-energy plot as a function of the displacement, i.e., energy–distance curve.[41,56,57]

The last 2 decades have witnessed a significant development of the computational methods for the potential-surface search or the search for the minimum-energy path (MEP), such as the nudged-elastic-band (NEB) method,[64] which is what we adopted for the purpose of the calculation of the potential-energy curve, which, thus obtained, is sometimes termed the NEB chain.

The NEB method is implemented in VASP and available, though the computational cost for the NEB calculation is usually large. But our computational resource is enough for conducting it. An advantage of the NEB calculations is that we do not need to impose artificial constraints on the adhesive structure. We need such constraints in the MEP search with the partial optimization. Therefore, the MEP or the distance–energy curve obtained from the NEB computation is likely to be more plausible.

To carry out an NEB calculation, one needs to locate two local minima: the initial state and the final state. In our simulation, the initial state corresponds to the structure in which the adhesive molecule is attached on the adherend surface. The structures shown in Figure were used as the initial state. As for the final-state structure, we used a structure in which the adhesive molecule is well separated from the surface. To obtain such a structure, the epoxy fragments were displaced in the direction perpendicular to the surface by 3 Å, and then the whole structure was optimized without any constraints except for the lowermost layer of the h-BN or graphite slab.[65] Using these initial-state and final-state structures, we performed an NEB calculation with 16 images for the fragments 1, 2, 3, and 4 on the h-BN and graphite surfaces.

Figure shows the MEP or energy–distance curve associated with the detachment of the epoxy fragments 1, 2, 3, and 4 from the h-BN or graphite surface. Notice that there are 18 points on each potential curve: two of them correspond to the initial and the final states and the other 16 points correspond to the intermediate images of the NEB chain.

Figure 6

Energy–distance curves or the NEB chains associated with the detachment of the epoxy fragments from (a) the h-BN or (b) graphite surface. The horizontal axis indicates the displacement of the center of gravity of the epoxy fragment from the position of stable equilibrium. The lines are least-squares fitting to the Morse potential (see the text for more detail).

The calculated plots of the energy versus distance were approximated by a Morse potential curve by using the least-squares method.[41,56,57] The Morse potential used is written as followswhere De is the binding energy or the adhesive energy, a is a constant specific to a system, and Δr is the displacement of the epoxy fragment. The approximated Morse potential curves are also shown in Figure . It is clear from this figure that the potential curves are perfectly fitted by the Morse function, where the R2 value of the fitting is no smaller than 0.98.

The fitting parameters De and a are tabulated in the SI. One could get something out of them, as can be seen in the literature,[42] but here the variation from fragment to fragment or that from surface to surface appears not so large. Nevertheless, we note that there is a good correspondence between the De values and the adhesion energies.

The derivative of eq with respect to the displacement produces a force–displacement curve as followsFigure shows calculated adhesion forces exerted on the fragments 1, 2, 3, and 4 upon the detachment from the h-BN or graphite surface. On both surfaces, the adhesion forces take their maximum values, Fmax, when the displacement approaches 0.6 Å. The Fmax values are summarized in Table . The Fmax value of the fragment 3 is about 1.6 nN in both cases, and this is the largest of the four fragments. The other three fragments have almost the same Fmax value, which falls in a range from 1.3 to 1.5 nN. This trend is almost consistent with that found in the adhesion energy (see Table ).

Figure 7

Adhesion force–displacement curves for the four fragments on (a) the h-BN or (b) graphite surface calculated from eq .

Table 2

Maximum Adhesion Forces of the Fragments 1, 2, 3, and 4 onto the h-BN or Graphite Surface in Units of nN

	fragment 1	fragment 2	fragment 3	fragment 4
h-BN	1.44	1.46	1.58	1.39
graphite	1.36	1.42	1.62	1.37

Hellmann–Feynman Force as Adhesion Force

When it comes to interatomic forces in general, putting aside our convention of the derivative of the Morse potential, one may apply a remarkable theorem named after Hellmann and Feynman.[66,67] From this theorem, the z-component of the force on the Ith nucleus can be calculated as[68,69]where Z is the z-coordinate of the Ith nucleus, Ψ is the normalized wave function, and V is the potential energy that consists of electronic repulsion, nuclear repulsion, and electron–nuclear interaction. The magnitude of the sum of the z-component of the Hellmann–Feynman force on a nucleus over all atoms included in the epoxy fragment leads to another measure of the adhesion force as followsVASP provides us with the value of F on each atom so that we can evaluate the adhesion force working on the fragments using this equation. Here, one should notice that since we turn on the dispersion correction in our computation, the van der Waals force acting on each atom is also included in F.

Figure shows the adhesion force–displacement plots calculated from eq (Hellmann–Feynman force). For comparison, the adhesion force–displacement curves obtained from the Morse-potential approximation are also shown. We do not see any significant inconsistency between them. The difference in the maximum adhesion force between these two methods is no larger than 0.2 nN. Since the two different methods predict almost the same adhesion force, the numbers listed in Table can be viewed as reliable. However, it is a pity that we do not find any experimental values to be compared in the literature.

Figure 8

Hellmann–Feynman adhesion forces acting on the fragments 1 (blue diamonds), 2 (orange square), 3 (gray triangle), and 4 (yellow square) evaluated at each image of the NEB chain using eq on the h-BN (a) or graphite (b) surface. For comparison, they are plotted on the adhesion force–displacement curves obtained from the Morse-potential approximation, which are the same as those shown in Figure .

Force-Decomposition Analysis

In previous sections, we have calculated the adhesion energies and the adhesion forces for the epoxy/h-BN and epoxy/graphite interfaces; we have not perceived any notable difference between them. This observation is actually a bit of a surprise to us. In this section, we will probe the origin of the adhesion forces in both interfaces and try to understand why both surfaces have almost the same affinity to the epoxy fragments. To this end, we conduct a force–decomposition analysis.

In this study, we adopt the dispersion correction scheme of Tkatchenko–Scheffler. As shown in eq , the whole energy can be separated into the contributions from the DFT term and the dispersion–correction one. The derivative of each term with respect to the displacement results in the decomposition of the whole adhesion force–displacement curve into two force curves: one is due to the PBE-functional-based force, denoted as FDFT, and the other is due to the dispersion force, denoted as Fdisp, as expressed byUp to this point, we have not seen any significant variation from fragment to fragment, and hence we only perform the force–decomposition analysis for the fragment 1 on the h-BN or graphite surface.

As shown in Figure , we decompose the energies of the images in the NEB chain corresponding to the detachment process of fragment 1 from the h-BN or graphite surface into the contributions from the DFT energy and the dispersion–correction one. In contrast with Figure , where the energy of the initial image of the NEB chain is set to the origin of energy, the energy of the final image of the NEB chain is set to the origin of energy in this figure. As such, here the positive energy corresponds to repulsion, whereas the negative one corresponds to attraction.

Figure 9

For all images in the NEB chain used for the simulation of the detachment process of fragment 1 from (a) the h-BN or (b) graphite surface, the potential energy at each point is decomposed into the contributions from the DFT energy (denoted by red diamonds) and dispersion–correction one (denoted by blue squares). The points of the DFT and dispersion energies are well fitted to 6th order polynomials, which are drawn by the red and blue curves, respectively.

Comparing Figures a,b, we notice that both graphs show a similar trend that the contribution from the PBE functional of DFT is very small but the effect of dispersion force is significant. As the epoxy fragment gets closer to the surface, the potential energy coming from the dispersion force goes down, whereas the DFT energy goes up. Thus, we come to a conclusion that epoxy resin can be adhered to both surfaces of h-BN and graphite mainly through dispersion interaction, so the inclusion of the dispersion correction in the DFT calculation of adhesion is essential.

Let us convert Figure to the adhesion force–displacement curve. To this end, we need to calculate the derivatives of the curves in this figure. When we made Figure , the energy–displacement plots were nicely approximated by the Morese potential, but now we cannot use it for the dispersion energy–displacement curves because the shapes of the curves appear different from the Morse potential. Therefore, for now, we use a polynomial regression: two plots in each panel in Figure are fitted to 6th order polynomials with an R2 value of 1.00. The exact forms of the approximated polynomial equations and their derivatives are shown in the SI.

The derivatives of the approximated polynomial functions are plotted in Figure . It is clear from this figure that the main attractive force acting at the two interfaces is the dispersion force, whereas the DFT contribution acts as a repulsive force in a region where the separation of the resin from the surface is small. The attractive contribution of the DFT force can be found in a range of displacement from 0.8 to 2.4 Å, but its magnitude is very small.

Figure 10

Adhesion force–displacement curve for fragment 1 on (a) the h-BN or (b) graphite surface is decomposed into contributions from the DFT (red) and dispersion (blue) parts.

Origin of Dispersion Force

As shown in Figure , the dispersion force working at the epoxy/h-BN interface is as large as that at the epoxy/graphite interface, and the dispersion force is the dominant intermolecular force resulting in the adhesion force. In this section, we intend to pursue the reason why the magnitudes of the dispersion forces at the two interfaces are almost the same. To this end, we need to trace back to the original formulation of the dispersion correction by Tkatchenko and Scheffler (TS),[36] which was used in our study. The general formulation used in this method is formally identical to that of Grimme’s D2 method,[70] where a pairwise interatomic C6R–6 term is added to the DFT energy. The added term reads[35]where C6 is the C6 coefficient for the atomic pair (i and j), r is the distance between the ith and jth atoms, and f(r) is the short-ranged damping function that scales the force field so as to avoid the singularity. The summations run over all N atoms. Also, all translations of the unit cell have to be included in the summation, but this is not indicated explicitly in eq . From the form of this equation, the dispersion energy is significantly affected by the C6 parameter.

In Grimme’s formulation, C6 is calculated from the geometric mean of the element-specific C6 parameters for the ith and jth atoms. Grimme proposed a simple computational scheme for the C6 parameters based on the London formula for dispersion, where C6 can be regarded to be proportional to the atomic ionization potential and the static dipole polarizability. In this method, however, the parameter is only dependent on the element, fixed during the calculation, and not sensitive to the particular chemical environment.

On the other hand, in the TS method, the C6 parameter is dependent on the atomic charge density so that variation in the local chemical environment can be taken into account. The C6 parameter for the ith atom can be calculated fromwhere C6free is the C6 parameter for the free atom and ν is the effective atomic volume estimated from the partitioning of the total electron density proposed by Hirshfeld.[58] ν can be written aswhere r is the distance from the nucleus of the atom i, ρ is the total electron density, ρfree is the electron density of the neutral free atom i, and w is the Hirshfeld atomic partitioning weight for the atom i. The partitioning weight w can be estimated fromOnce the parameters C6 and C6 are determined using eqs –12, one can evaluate the pairwise parameter C6 using the following equationwhere α is the free-atomic reference value for the static polarizability of the ith atom.

As far as the interfacial dispersion interaction goes (cf. Figure ), the formulation of the dispersion correction (eq ) may be rewritten aswhere i runs over all the atoms included in the surface, whereas j runs over all the atoms in the epoxy fragment. The contributions of the dispersion energy coming from the atom pairs inside the surface or the epoxy fragment also exist, but the amplitudes of such contributions are not expected to change in the process of detachment more or less because the structural deformation of the adherend surface and the epoxy resin was found to be small; their contributions are likely to be canceled when one calculates the energy difference ΔEdisp.

When compared between the epoxy/h-BN and epoxy/graphite interfaces, the epoxy fragment is common. Thus, the C6 parameters for the atoms included in the h-BN or graphite surface should be investigated closely. Table summarizes the average values of the C6 parameters in the TS scheme for the B and N atoms in the h-BN surface and the one for the C atoms in the graphite surface. Note that the C6 parameter in the TS scheme varies from atom to atom even if they are the same element because of different chemical environment. Thus, we calculated the average value as listed in Table .

Table 3

Average Values of the C6 Parameters for the B and N Atoms in the h-BN Surface and the One for the C Atoms in the Graphite Surface in Units of Jnm6 mol–1

B	N	C
3.68	1.07	1.97

The C6 parameter increases in the order: N < C < B. Hence, we come to a conclusion that at the epoxy/h-BN interface, the B atoms result in a strong dispersion force, but the N atoms result in a weak dispersion force; therefore, the total dispersion force is moderate, whereas the C6 parameter for the C atoms in the epoxy/graphite interface has a middle value between the C6 parameters for the B and N atoms, so the total dispersion force for this interface is also moderate. We believe this should be a possible reason why the two interfacial interactions end up with almost the same dispersion force.

Missing Hydrogen Bond

We had expected a hydrogen bond formation between the OH group in the epoxy fragment and the N atoms on the h-BN surface, but such a hydrogen bond is not observed. When it comes to the N···H–O type H-bonding, the two-orbital–two-electron interaction is expected to exist between the lone-pair orbital of the N atom and the unoccupied antibonding orbital of the O–H bond (σO-H*).[31] We hypothesized that a lone-pair orbital could be found on the N atoms of the h-BN surface, but that has proven not to be true, as delineated below.

We carried out the calculation of the electron localization function (ELF)[71,72] for the h-BN surface to identify whether or not the lone pair exists on the N atoms in the h-BN surface. For comparison, the ELF for the graphite surface was also calculated. They are plotted in Figure . The ELF takes a value ranging from 0 to 1. High ELF values, approaching 1, correspond to paired electrons, such as lone pairs, bonds, and core electrons, whereas the ELF values around 0.5 correspond to a uniform electron gas.[73]

Figure 11

Isosurface (yellow) of the electron localization function (ELF) with the value of 0.85 for (a) the h-BN and (b) graphite (001) surfaces. Gray balls indicate N, green balls indicate B, and brown balls indicate C.

In Figure a, we cannot find any ELF isosurface corresponding to N’s lone pair, which could be a potential H acceptor of H bond. This finding should be important for understanding the absence of H-bonding in epoxy/h-BN interface. If the ELF is calculated for an NH3 molecular crystal, a clear lobe corresponding to the lone pair can be observed on the N atom in the ELF isosurface plot (see the SI). Thus, the absence of N’s lone pair is likely to be an inherent feature of h-BN. This may also be a key to the inertness of the h-BN surface.

Instead, high ELF values are found in the middle of the B–N bonds and the C–C bonds in Figures a,b, respectively. These two plots look alike in that the high-ELF regions are centered on the bonds, implying that electrons are mainly localized in bonds. Such ELF plots reflect the covalent nature of the B–N and C–C bonds. If h-BN was ionic, high ELF values would be observed only near N atoms, but that is not true.

Summary and Conclusions

In this paper, we have presented a theoretical study on the interfacial interaction between epoxy resin and h-BN or graphite with applications to nanocomposite materials in mind. h-BN and graphite share a structural feature that they consist of the stacking of hexagonal honeycomb layers; however, the C–C bonds in graphite and the B–N bonds in h-BN have quite a different electronic character. For example, the former has delocalized π-electrons and is characterized as a conductor, whereas the latter is an insulator due to the polarity of the B–N bonds. We expected that such a difference in the electronic structure would result in different adhesion interaction with epoxy resin. But that has not proven to be true. Both h-BN and graphite surfaces have almost the same affinity to epoxy resin, which is found to be governed by dispersion interaction.

To investigate the interfacial interaction of the basal surfaces of h-BN and graphite with epoxy resin, we used fragment models of epoxy resin, which can reproduce a local interaction between the surface and resin. The adhesion structures formed between the fragment structure and the surface slab were obtained from the first-principles simulated annealing method followed by geometry optimization. The adhesion structures were analyzed by using the mass density profile and the radial distribution function, but we have not seen a significant variation between the epoxy/h-BN and epoxy/graphite models. No H-bonding was observed in the two kinds of interfaces. However, we have found a structure that hints at the importance of the π–π stacking between the benzene ring of the epoxy fragment and the h-BN or graphite surface for the reinforcement of adhesion. The above-mentioned observation is almost consistent with the calculated adhesion energies of the investigated interfaces.

We have simulated the adhesion force–separation distance curves in two different methods: one is based on the Morse-potential approximation and the other is evaluated from the Hellmann–Feynman force. In both of them, the potential-energy surface associated with the detachment of the epoxy fragment from the surface is sampled by using the NEB method. The results of the calculations for the two interfaces coincide with each other. Very similar adhesion–force curves were obtained for the two interfaces. To understand this feature, we conducted a force-decomposition analysis, in which the total adhesion force is separated into the contributions from the DFT force and the dispersion force due to the Tkatchenko–Scheffler dispersion correction. As a result, we found that most of the attractive force working at both epoxy/h-BN and epoxy/graphite interfaces is dominated by the dispersion force and its magnitude is comparable between them. The contribution from the DFT force to attraction was found to be negligible, and hence it is inevitable to include the dispersion correction when one wants to calculate adhesion interfaces with DFT.

The fact that the adhesion force between the epoxy/h-BN and epoxy/graphite interfaces is the same originates from the C6 parameters of the B, N, and C atoms used for the dispersion correction. At the epoxy/h-BN interface, the C6 parameter of the B atoms is larger but that of the N atoms is smaller, resulting in a moderate dispersion force. In the case of the epoxy/graphite interface, the C6 parameter of the C atoms of graphite is larger than that of the N atoms but smaller than that of the B atoms, resulting a moderate dispersion force. This is why the two interfaces end up with almost the same adhesion interaction. As for the absence of the H-bonding in epoxy/h-BN interfaces, which we had anticipated, our ELF calculation has revealed that there is no lone-pair orbital on the N atoms, so the N atoms in the h-BN surface cannot act as an H acceptor of H bond.

Computational Methods

In this study, all of the periodic density functional theory (DFT) calculations were carried out by using the Vienna Ab initio Simulation Package (VASP 5.4.1).[74−76] The generalized gradient approximation of Perdew–Burke–Ernzerhof (GGA-PBE)[77] was adopted for the exchange–correlation functional. The electron–ion interactions were treated within the projector-augmented wave scheme.[78,79] A plane-wave basis set cutoff of 600 eV, self-consistent field tolerance of 1.0 × 10–5 eV, Brillouin zone sampling on a grid of spacing 2π × 0.05 Å–1, and 0.05 eV/Å threshold of forces on atoms guaranteed good convergence. The structure graphics were produced by using VESTA.[80]