Adhesion of Epoxy Resin with Hexagonal Boron Nitride and Graphite

5 days ago - Adhesion interaction of epoxy resin with the basal surfaces of h-BN and graphite is investigated with the first-principles density functi...
1 downloads 0 Views 4MB Size
This is an open access article published under an ACS AuthorChoice License, which permits copying and redistribution of the article or any adaptations for non-commercial purposes.

Article Cite This: ACS Omega 2019, 4, 4491−4504

http://pubs.acs.org/journal/acsodf

Adhesion of Epoxy Resin with Hexagonal Boron Nitride and Graphite Yuta Tsuji,† Yasuhiro Kitamura,† Masao Someya,‡ Toshihiko Takano,‡ Michio Yaginuma,‡ Kohei Nakanishi,‡ and Kazunari Yoshizawa*,† †

Institute for Materials Chemistry and Engineering and IRCCS, Kyushu University, Nishi-ku, Fukuoka 819-0395, Japan Mitsubishi Gas Chemical Company Inc., Chiyoda-ku, Tokyo 100-8324, Japan



Downloaded via 5.101.222.87 on March 5, 2019 at 15:04:04 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: 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 forcedecomposition 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 C6 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.

1. 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 © 2019 American Chemical Society

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 Received: January 14, 2019 Accepted: February 11, 2019 Published: March 1, 2019 4491

DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504

ACS Omega

Article

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 conductivity13,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 behavior18 and electric heating behavior19 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 1). Although

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 2pz 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.

2. RESULTS AND DISCUSSION 2.1. 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 E DFT ‐ disp = E DFT + Edisp (1)

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

We checked the reproducibility of the interlayer distances of hBN and graphite using several dispersion correction methods available in VASP, finding that the Tkatchenko−Scheffler (TS) dispersion correction36 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. 2.2. 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

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 4492

DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504

ACS Omega

Article

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

2.4. 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 2.2 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 hBN 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 hBN 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 3, we show the determined adhesion

energy needed to cleave the bulk crystal and is obtained from a slab calculation using the following equation:37 γ = (Esurf − nE bulk )/2A

(2)

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 2ML 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 alumina41 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. 2.3. Modeling of Epoxy Resin. As shown in Figure 2a, 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 2a. The resultant fragment models are shown in Figure 2b. 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. 4493

DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504

ACS Omega

Article

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 balland-stick model. Only the topmost layer is shown.

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 Å.

structures of fragments 1, 2, 3, and 4 on the basal planes of hBN and graphite. 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 DFTbased 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 3, 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 4494

DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504

ACS Omega

Article

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 Å.

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 1. An acid−base interaction can be found between the

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 interaction52 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 4 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 Å. 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 5. Since Figure 4 implies that the variation from fragment to fragment is not so significant, the RDF plots for the adhesion interface shown in Figure 5 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. 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 3. 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

Scheme 1. Structures of 1:1 Molecular Complexes Bonded by a B−O or N−H Intermolecular Interactiona

a

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 refs60, 61, and 62, respectively.

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 Å. Given the intermolecular B−O and N−H distances found in Scheme 1, one cannot expect any specific intermolecular interactions at the interface between the epoxy fragment and 4495

DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504

ACS Omega

Article

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-workers60,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 1 mainly come from the interaction between the surface and the adherend, and the effects of distortion can be negligible. 2.6. 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

the h-BN(001) surface on the basis of the RDF plots shown in Figure 5. 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. 2.5. 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 Eadhesion = Eadhesive + Eadherend − E(adhesive + adherend)

(3)

where E(adhesive+adherend) is the total energy of the surface−resin complex shown in Figure 3 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 1. Mostly, the adhesion energies of Table 1. Adhesion Energies of the Fragments 1, 2, 3, and 4 on the h-BN or Graphite Surface in Units of eV h-BN graphite

fragment 1

fragment 2

fragment 3

fragment 4

1.53 1.42

1.57 1.45

1.73 1.76

1.51 1.43

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 analysis58 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. 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

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). 4496

DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504

ACS Omega

Article

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

also shown in Figure 6. 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 4 with respect to the displacement produces a force−displacement curve as follows

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 potentialenergy 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 3 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 6 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. The calculated plots of the energy versus distance were approximated by a Morse potential curve by using the leastsquares method.41,56,57 The Morse potential used is written as follows E = De(1 − e−aΔr )2

F(Δr ) =

dE = 2aDe e−aΔr (1 − e−aΔr ) dΔr

(5)

Figure 7 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, F max , when the displacement approaches 0.6 Å. The Fmax values are summarized in Table 2. The Fmax value of the fragment 3 is about 1.6 nN in both Table 2. Maximum Adhesion Forces of the Fragments 1, 2, 3, and 4 onto the h-BN or Graphite Surface in Units of nN h-BN graphite

fragment 1

fragment 2

fragment 3

fragment 4

1.44 1.36

1.46 1.42

1.58 1.62

1.39 1.37

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 1). 2.7. 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 as68,69 FIz = −

(4)

∂E ∂V =− Ψ Ψ ∂ZI ∂ZI

(6)

where ZI 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

where 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 4497

DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504

ACS Omega

Article

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 7 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 7.

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.

that we do not find any experimental values to be compared in the literature. 2.8. 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 1, 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 by

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 follows F=

∑ I ∈ epoxy

FIz (7)

VASP provides us with the value of FIZ 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 FIZ. Figure 8 shows the adhesion force−displacement plots calculated from eq 7 (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 2 can be viewed as reliable. However, it is a pity

F = FDFT + Fdisp 4498

(8) DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504

ACS Omega

Article

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.

2.9. Origin of Dispersion Force. As shown in Figure 10, 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 reads35

Up 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 9, 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 6, 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. Comparing Figures 9a,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 9 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 7, 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 9 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 10. 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.

Edisp = −

N N C 1 6ij ∑ ∑ 6 fd ,6 (rij) 2 i = 1 j = 1 rij

(9)

where C6ij is the C6 coefficient for the atomic pair (i and j), rij is the distance between the ith and jth atoms, and fd,6(rij) 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 9. From the form of this equation, the dispersion energy is significantly affected by the C6 parameter. In Grimme’s formulation, C6ij 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 from C6i = νi2C6free i

(10)

where Cfree 6i is the C6 parameter for the free atom and νi is the effective atomic volume estimated from the partitioning of the 4499

DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504

ACS Omega

Article

total electron density proposed by Hirshfeld.58 νi can be written as νi =

∫ r 3wi(r)ρ(r) d3r ∫ r 3ρi free (r) d3r

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. 2.10. 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 11. The ELF takes a

(11)

where r is the distance from the nucleus of the atom i, ρ is the total electron density, ρfree is the electron density of the neutral i free atom i, and wi is the Hirshfeld atomic partitioning weight for the atom i. The partitioning weight wi can be estimated from wi =

ρi free (r) N

∑ j = 1 ρjfree (r)

(12)

Once the parameters C6i and C6j are determined using eqs 10−12, one can evaluate the pairwise parameter C6ij using the following equation C6ij =

2C6iC6j αj αi

C6i +

αi C αj 6j

(13)

where αi is the free-atomic reference value for the static polarizability of the ith atom. As far as the interfacial dispersion interaction goes (cf. Figure 9), the formulation of the dispersion correction (eq 9) may be rewritten as N

ΔEdisp ≈ −

1 ∑ 2 i ∈ surf.

N

∑ j ∈ epoxy

C6ij rij6

fd ,6 (rij)

(14)

where 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 3 summarizes the

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.

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

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 In Figure 11a, 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 hBN. 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 11a,b, respectively.

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 3. 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 4500

DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504

ACS Omega

Article

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/hBN 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.

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.

3. 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 Morsepotential approximation and the other is evaluated from the Hellmann−Feynman force. In both of them, the potentialenergy 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

4. 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 (GGAPBE)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



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acsomega.9b00129. Benchmark test of the dispersion−correction methods, effects of layers included in the slab model on the surface energy, top view of the surface slab model, radial distribution function plots, the distortion energies of the epoxy fragments and the h-BN or graphite surface, fitting parameters for the Morse potential, exact forms of the approximated polynomial equations, ELF isosurface for an NH3 crystal structure, coordinates of the optimized structures, and supplementary references (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Yuta Tsuji: 0000-0003-4224-4532 Kazunari Yoshizawa: 0000-0002-6279-9722 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by KAKENHI grant numbers JP15K13710, JP17K14440, and JP17H03117 from Japan Society for the Promotion of Science (JSPS) and the Ministry of Education, Culture, Sports, Science and Technology of Japan (MEXT), the MEXT Projects of “Integrated Research Consortium on Chemical Sciences”, “Cooperative Research Program of Network Joint Research Center for Materials and Devices”, “Elements Strategy Initiative to Form Core Research Center”, and JST-CREST “Innovative Catalysts”. The computation was mainly carried out using the computer 4501

DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504

ACS Omega

Article

(18) Liang, J.; Wang, Y.; Huang, Y.; Ma, Y.; Liu, Z.; Cai, J.; Zhang, C.; Gao, H.; Chen, Y. Electromagnetic Interference Shielding of Graphene/Epoxy Composites. Carbon 2009, 47, 922−925. (19) An, J. E.; Jeong, Y. G. Structure and Electric Heating Performance of Graphene/Epoxy Composite Films. Eur. Polym. J. 2013, 49, 1322−1330. (20) Lee, J. H.; Kang, S. G.; Choe, Y.; Lee, S. G. Mechanism of Adhesion of the Diglycidyl Ether of Bisphenol A (DGEBA) to the Fe(100) Surface. Compos. Sci. Technol. 2016, 126, 9−16. (21) Lee, J. H.; Choe, Y.; Lee, S. G. Adhesion Mechanism of Bisphenol A Diglycidyl Ether (BADGE) on an α-Fe2O3 (0001) Surface. J. Ind. Eng. Chem. 2017, 53, 62−67. (22) Bahlakeh, G.; Ghaffari, M.; Saeb, M. R.; Ramezanzadeh, B.; De Proft, F.; Terryn, H. A Close-up of the Effect of Iron Oxide Type on the Interfacial Interaction between Epoxy and Carbon Steel: Combined Molecular Dynamics Simulations and Quantum Mechanics. J. Phys. Chem. C 2016, 120, 11014−11026. (23) Bahlakeh, G.; Ramezanzadeh, B. A Detailed Molecular Dynamics Simulation and Experimental Investigation on the Interfacial Bonding Mechanism of an Epoxy Adhesive on Carbon Steel Sheets Decorated with a Novel Cerium−Lanthanum Nanofilm. ACS Appl. Mater. Interfaces 2017, 9, 17536−17551. (24) Bahlakeh, G.; Ramezanzadeh, B.; Saeb, M. R.; Terryn, H.; Ghaffari, M. Corrosion Protection Properties and Interfacial Adhesion Mechanism of an Epoxy/Polyamide Coating Applied on the Steel Surface Decorated with Cerium Oxide Nanofilm: Complementary Experimental, Molecular Dynamics (MD) and First Principle Quantum Mechanics (QM) Simulation Methods. Appl. Surf. Sci. 2017, 419, 650−669. (25) Bahlakeh, G.; Ramezanzadeh, B.; Ramezanzadeh, M. New Detailed Insights on the Role of a Novel Praseodymium Nanofilm on the Polymer/Steel Interfacial Adhesion Bonds in Dry and Wet Conditions: An Integrated Molecular Dynamics Simulation and Experimental Study. J. Taiwan Inst. Chem. Eng. 2018, 85, 221−236. (26) Chen, Y. I. Ed. Nanotubes and Nanosheets: Functionalization and Applications of Boron Nitride and Other Nanomaterials; CRC Press: Boca Raton, 2015. (27) Marom, N.; Bernstein, J.; Garel, J.; Tkatchenko, A.; Joselevich, E.; Kronik, L.; Hod, O. Stacking and Registry Effects in Layered Materials: The Case of Hexagonal Boron Nitride. Phys. Rev. Lett. 2010, 105, No. 046801. (28) Yoshizawa, K.; Yumura, T.; Yamabe, T.; Bandow, S. The Role of Orbital Interactions in Determining the Interlayer Spacing in Graphite Slabs. J. Am. Chem. Soc. 2000, 122, 11871−11875. (29) Wu, H.; Zhao, W.; Hu, H.; Chen, G. One-Step In Situ Ball Milling Synthesis of Polymer-Functionalized Graphene Nanocomposites. J. Mater. Chem. 2011, 21, 8626−8632. (30) Ma, R.; Wan, X.; Zhang, T.; Yang, N.; Luo, T. Interfacial Thermal Transport in Boron Nitride-Polymer Nanocomposite. arXiv preprint arXiv:1711.11201, 2017. (31) Yoshizawa, K.; Semoto, T.; Hitaoka, S.; Higuchi, C.; Shiota, Y.; Tanaka, H. Synergy of Electrostatic and van der Waals Interactions in the Adhesion of Epoxy Resin with Carbon-Fiber and Glass Surfaces. Bull. Chem. Soc. Jpn. 2017, 90, 500−505. (32) Klesnar, H. P.; Rogl, P. Phase Relations in the Ternary Systems Rare-Earth Metal (RE)-Boron-Nitrogen, Where RE = Tb, Dy, Ho, Er, Tm, Lu, Sc, and Y. High Temp. High Press. 1990, 22, 453−457. (33) Zhao, Y. X.; Spain, I. L. X-Ray Diffraction Data for Graphite to 20 GPa. Phys. Rev. B 1989, 40, 993. (34) Tuma, C.; Sauer, J. Treating Dispersion Effects in Extended Systems by Hybrid MP2: DFT CalculationsProtonation of Isobutene in Zeolite Ferrierite. Phys. Chem. Chem. Phys. 2006, 8, 3955−3965. (35) Kresse, G.; Furthmüller, J. VASP the Guide. http://cms.mpi. univie.ac.at/vasp/vasp/vasp.html. (36) Tkatchenko, A.; Scheffler, M. Accurate Molecular van der Waals Interactions from Ground-State Electron Density and FreeAtom Reference Data. Phys. Rev. Lett. 2009, 102, No. 073005.

facilities at Research Institute for Information Technology, Kyushu University. Y.T. thanks JSPS Grant-in-Aid for Scientific Research on Innovative Areas “Discrete Geometric Analysis for Materials Design”: Grant Number JP18H04488.



REFERENCES

(1) Mai, Y.-W.; Yu, Z.-Z. Eds. Polymer Nanocomposites; CRC: Boca Raton, 2006. (2) Thakur, V. K.; Thakur, M. K. Eco-Friendly Polymer Nanocomposites: Chemistry and Applications. Advanced Structured Materials; Springer: New Delhi, 2015. (3) Pandey, J. K.; Takagi, H.; Nakagaito, A. N.; Kim, H.-J. Processing, Performance and Application: Polymer Nanocomposites of Cellulose Nanoparticles. Handbook of Polymer Nanocomposites; Springer: Berlin, 2014; Vol. C. (4) Han, Z.; Wood, J. W.; Herman, H.; Zhang, C.; Stevens, G. C. Thermal Properties of Composites Filled with Different Fillers; IEEE International Symposium on Electrical Insulation: Vancouver, 2008. (5) Pissis, P.; Kotsilkova, R. Thermoset Nanocomposites for Engineering Applications; Smithers Rapra Technology Limited: Shawbury, 2007. (6) Lee, D.; Song, S. H.; Hwang, J.; Jin, S. H.; Park, K. H.; Kim, B. H.; Hong, S. H.; Jeon, S. Enhanced Mechanical Properties of Epoxy Nanocomposites by Mixing Noncovalently Functionalized Boron Nitride Nanoflakes. Small 2013, 9, 2602−2610. (7) Rittigstein, P.; Torkelson, J. M. Polymer−Nanoparticle Interfacial Interactions in Polymer Nanocomposites: Confinement Effects on Glass Transition Temperature and Suppression of Physical Aging. J. Polym. Sci., Part B: Polym. Phys. 2006, 44, 2935−2943. (8) Zare, Y.; Garmabi, H. Modeling of Interfacial Bonding between Two Nanofillers (Montmorillonite and CaCO3) and a Polymer Matrix (PP) in a Ternary Polymer Nanocomposite. Appl. Surf. Sci. 2014, 321, 219−225. (9) Yung, K. C.; Liem, H. Enhanced Thermal Conductivity of Boron Nitride Epoxy-Matrix Composite through Multi-Modal Particle Size Mixing. J. Appl. Polym. Sci. 2007, 106, 3587−3591. (10) Barzic, A. I., Ioan, S., Eds. Multiphase Polymer Systems: Micro-to Nanostructural Evolution in Advanced Technologies; CRC Press: Boca Raton, 2016. (11) Zhu, B. L.; Ma, J.; Wu, J.; Yung, K. C.; Xie, C. S. Study on the Properties of the Epoxy-Matrix Composites Filled with Thermally Conductive AlN and BN Ceramic Particles. J. Appl. Polym. Sci. 2010, 118, 2754−2764. (12) Lin, Z. Y.; McNamara, A.; Liu, Y.; Moon, K. S.; Wong, C. P. Exfoliated Hexagonal Boron Nitride-Based Polymer Nanocomposite with Enhanced Thermal Conductivity for Electronic Encapsulation. Compos. Sci. Technol. 2014, 90, 123−128. (13) Teng, C.-C.; Ma, C.-C. M.; Lu, C.-H.; Yang, S.-Y.; Lee, S.-H.; Hsiao, M.-C.; Yen, M.-Y.; Chiou, K.-C.; Lee, T.-M. Thermal Conductivity and Structure of Non-Covalent Functionalized Graphene/Epoxy Composites. Carbon 2011, 49, 5107−5116. (14) Li, Q.; Guo, Y.; Li, W.; Qiu, S.; Zhu, C.; Wei, X.; Chen, M.; Liu, C.; Liao, S.; Gong, Y.; Mishra, A. K.; Liu, L. Ultrahigh Thermal Conductivity of Assembled Aligned Multilayer Graphene/Epoxy Composite. Chem. Mater. 2014, 26, 4459−4465. (15) Tang, L. C.; Wan, Y. J.; Yan, D.; Pei, Y. B.; Zhao, L.; Li, Y. B.; Wu, L. B.; Jiang, J. X.; Lai, G. Q. The Effect of Graphene Dispersion on the Mechanical Properties of Graphene/Epoxy Composites. Carbon 2013, 60, 16−27. (16) Naebe, M.; Wang, J.; Amini, A.; Khayyam, H.; Hameed, N.; Li, L. H.; Chen, Y.; Fox, B. Mechanical Property and Structure of Covalent Functionalised Graphene/Epoxy Nanocomposites. Sci. Rep. 2014, 4, No. 4375. (17) Jović, N.; Dudic, D.; Montone, A.; Antisari, M. V.; Mitric, M.; Djokovic, V. Temperature Dependence of the Electrical Conductivity of Epoxy/Expanded Graphite Nanosheet Composites. Scr. Mater. 2008, 58, 846−849. 4502

DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504

ACS Omega

Article

(37) Wu, Z.; Overbury, S. H.; Eds. Catalysis by Materials with WellDefined Structures, 1st ed.; Elsevier: Amsterdam, 2015. (38) Cherkaoui, M.; Capolungo, L. Atomistic and Continuum Modeling of Nanocrystalline Materials; Springer-Verlag: New York, 2009. (39) Chung, J.; Ryu, B. Crystallization Behavior of Phosphate Glasses with Hydrophobic Coating Materials. J. Nanomater. 2015, 2015, No. 981720. (40) Sholl, D. S.; Steckel, J. A. Density Functional Theory: A Practical Introduction; Wiley: Hoboken, 2009. (41) Semoto, T.; Tsuji, Y.; Yoshizawa, K. Molecular Understanding of the Adhesive Force between a Metal Oxide Surface and an Epoxy Resin: Effects of Surface Water. Bull. Chem. Soc. Jpn. 2012, 85, 672− 678. (42) Higuchi, C.; Tanaka, H.; Yoshizawa, K. Molecular Understanding of the Adhesive Interactions between Silica Surface and Epoxy Resin: Effects of Interfacial Water. J. Comput. Chem. 2019, 40, 164−171. (43) Tocci, G.; Joly, L.; Michaelides, A. Friction of Water on Graphene and Hexagonal Boron Nitride from Ab Initio Methods: Very Different Slippage despite Very Similar Interface Structures. Nano Lett. 2014, 14, 6872−6877. (44) Bai, J.; Ed. Advanced Fibre-Reinforced Polymer (FRP) Composites for Structural Applications; Woodhead Publ. Oxford, 2013. (45) Prodhomme, P. Y.; Raybaud, P.; Toulhoat, H. Free-Energy Profiles Along Reduction Pathways of MoS2 M-Edge and S-Edge by Dihydrogen: A First-Principles Study. J. Catal. 2011, 280, 178−195. (46) Bedolla, P. O.; Feldbauer, G.; Wolloch, M.; Gruber, C.; Eder, S. J.; Dörr, N.; Mohn, P.; Redinger, J.; Vernes, A. Density Functional Investigation of the Adsorption of Isooctane, Ethanol, and Acetic Acid on a Water-Covered Fe(100) Surface. J. Phys. Chem. C 2014, 118, 21428−21437. (47) Ogata, S.; Uranagase, M. Unveiling the Chemical Reactions Involved in Moisture-Induced Weakening of Adhesion between Aluminum and Epoxy Resin. J. Phys. Chem. C 2018, 122, 17748− 17755. (48) Nishio, M.; Hirota, M.; Umezawa, Y. The CH/π Interaction. Evidence, Nature and Consequences; Wiley: New York, 1998. (49) Khanh, P. N.; Ngan, V. T.; Man, N. T. H.; Nhung, N. T. A.; Chandra, A. K.; Trung, N. T. An Insight into Csp−H··· π Hydrogen Bonds and Stability of Complexes Formed by Acetylene and Its Substituted Derivatives with Benzene and Borazine. RSC Adv. 2016, 6, No. 106662. (50) Tsuzuki, S.; Fujii, A. Nature and Physical Origin of CH/π Interaction: Significant Difference from Conventional Hydrogen Bonds. Phys. Chem. Chem. Phys. 2008, 10, 2584−2594. (51) Novoa, J., Ed. Intermolecular Interactions in Crystals: Fundamentals of Crystal Engineering; The Royal Society of Chemistry: Cambridge, 2017. (52) The distances measured from the benzene ring of the fragment 3 to the h-BN and graphite surfaces are 3.3 and 3.4 Å, respectively. These values are very close to the typical value of the π-π stacking separation and the interlayer distances of h-BN and graphite shown in Figure 1. (53) Mootz, D.; Steffen, M. Fluoride und Fluorosäuren. IV. Kristallstrukturen von Bortrifluorid und seinen 1: 1-Verbindungen mit Wasser und Methanol, Hydroxo-und Methoxotrifluoroborsäure. Z. Anorg. Allg. Chem. 1981, 483, 171−180. (54) Wagner, T.; Eigendorf, U.; Herberich, G. E.; Englert, U. From Crystal to Crystal: A Reconstructive Solid-State Reaction. Struct. Chem. 1994, 5, 233−237. (55) Steiner, T.; Majerz, I.; Wilson, C. C. First O−H···N Hydrogen Bond with a Centered Proton Position Obtained by Thermally Induced Proton Migration. Angew. Chem., Int. Ed. Engl. 2001, 40, 2651−2654. (56) Semoto, T.; Tsuji, Y.; Yoshizawa, K. Molecular Understanding of the Adhesive Force between a Metal Oxide Surface and an Epoxy Resin. J. Phys. Chem. C 2011, 115, 11701−11708.

(57) Semoto, T.; Tsuji, Y.; Tanaka, H.; Yoshizawa, K. Role of Edge Oxygen Atoms on the Adhesive Interaction between Carbon Fiber and Epoxy Resin. J. Phys. Chem. C 2013, 117, 24830−24835. (58) Hirshfeld, F. L. Bonded-Atom Fragments for Describing Molecular Charge Densities. Theor. Chim. Acta 1977, 44, 129−138. (59) Tsuji, Y.; Yoshizawa, K. Adsorption and Activation of Methane on the (110) Surface of Rutile-Type Metal Dioxides. J. Phys. Chem. C 2018, 122, 15359−15381. (60) Ess, D. H.; Houk, K. N. Distortion/Interaction Energy Control of 1,3-Dipolar Cycloaddition Reactivity. J. Am. Chem. Soc. 2007, 129, 10646−10647. (61) Liu, F.; Paton, R. S.; Kim, S.; Liang, Y.; Houk, K. N. Diels− Alder Reactivities of Strained and Unstrained Cycloalkenes with Normal and Inverse-Electron-Demand Dienes: Activation Barriers and Distortion/Interaction Analysis. J. Am. Chem. Soc. 2013, 135, 15642−15649. (62) Bickelhaupt, F. M. Understanding Reactivity with Kohn−Sham Molecular Orbital Theory: E2−SN2 Mechanistic Spectrum and Other Concepts. J. Comput. Chem. 1999, 20, 114−128. (63) Diefenbach, A.; Bickelhaupt, F. M. Oxidative Addition of Pd to C−H, C−C and C−Cl Bonds: Importance of Relativistic Effects in DFT Calculations. J. Chem. Phys. 2001, 115, 4030−4040. (64) Mills, G.; Jonsson, H. Quantum and Thermal Effects in H2 Dissociative Adsorption: Evaluation of Free Energy Barriers in Multidimensional Quantum Systems. Phys. Rev. Lett. 1994, 72, 1124. (65) If the displacement is 1 or 2 Å, the subsequent optimization leads to the structure which is the same as the initial-state structure. However, such a situation does not happen if the displacement is 3 Å. So, we adopted a displacement of 3 Å. (66) Hellmann, H. Einführung in die Quantenchemie; F. Deuticke: Leipzig, 1937. (67) Feynman, R. P. Forces in Molecules. Phys. Rev. 1939, 56, 340− 343. (68) Chandra, A. K. Introductory Quantum Chemistry; Tata McGraw Hill New Delhi, 1994. (69) Finnis, M. Interatomic Forces in Condensed Matter; Oxford University Press: New York, 2003. (70) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (71) Becke, A. D.; Edgecombe, K. E. A Simple Measure of Electron Localization in Atomic and Molecular Systems. J. Chem. Phys. 1990, 92, 5397−5403. (72) Silvi, B.; Savin, A. Classification of Chemical Bonds Based on Topological Analysis of Electron Localization Functions. Nature 1994, 371, 683−686. (73) Kurzydłowski, D.; Ejgierd-Zaleski, P.; Grochala, W.; Hoffmann, R. Freezing in Resonance Structures for Better Packing: XeF2 Becomes (XeF+)(F−) at Large Compression. Inorg. Chem. 2011, 50, 3832−3840. (74) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47, 558−561. (75) Kresse, G.; Hafner, J. Ab Initio Molecular-Dynamics Simulation of the Liquid-Metal-Amorphous-Semiconductor-Transition in Germanium. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 49, 14251− 14269. (76) Kresse, G.; Furthmüller, J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15−50. (77) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (78) Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. Phys. Rev. B 1999, 59, 1758− 1775. (79) Adolph, B.; Furthmüller, J.; Beckstedt, F. Optical Properties of Semiconductors Using Projector-Augmented Waves. Phys. Rev. B 2001, 63, No. 125108. 4503

DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504

ACS Omega

Article

(80) Momma, K.; Izumi, F. VESTA 3 for Three-Dimensional Visualization of Crystal, Volumetric and Morphology Data. J. Appl. Crystallogr. 2011, 44, 1272−1276.

4504

DOI: 10.1021/acsomega.9b00129 ACS Omega 2019, 4, 4491−4504