Research Article  Open Access
Kewei Tang, Weihong Qi, Yaru Wei, Guoliang Ru, Weimin Liu, "HighThroughput Calculation of Interlayer van der Waals Forces Validated with Experimental Measurements", Research, vol. 2022, Article ID 9765121, 11 pages, 2022. https://doi.org/10.34133/2022/9765121
HighThroughput Calculation of Interlayer van der Waals Forces Validated with Experimental Measurements
Abstract
Interlayer van der Waals interactions play an important role in twodimensional (2D) materials on various occasions. The interlayer binding force is often directly measured and is considered more closely related to the exfoliation condition. However, a binding force database from accurate theoretical calculations does not yet exist. In this work, the critical interlayer binding force and energy are directly calculated for 230 2D materials, which exhibit divergent trends. A linear relationship that links the two quantities with the equilibrium interlayer distance is found and checked. Experiments are carried out for three different materials using atomic force microscopy. The measured forces show a consistent trend with the calculated results, and the estimated binding strengths are of the same order of magnitude as the predicted values. Our work can provide a reliable reference for interlayer adhesion studies and help establish accurate models of exfoliation processes.
1. Introduction
Interest in twodimensional (2D) materials continues to grow due to their promising properties and the ability to arbitrarily integrate themselves [1, 2]. The latter takes advantage of comparatively weak interlayer van der Waals (vdW) interactions. The interlayer vdW interaction is far weaker than ionic bonding [3] but strong enough to hold layers together and enable a certain degree of electronic coupling at the interface. This characteristic of interlayer vdW interactions has enabled many fascinating properties and novel physical phenomena, such as flatband and vanHove singularities [4–6]. From perturbation theory, the vdW potential is found to be a series expansion: (where is the distance) that should be terminated at different terms for different materials to achieve a certain accuracy [7, 8]. While taking only the first term of the series expansion as the vdW potential and modeling the exchange potential with an inverse power term, one obtains the famous Lennard–Jones (LJ) potential. Similarly, if the exchange potential is modeled more accurately with the Born–Mayer term [9], then the wellknown Buckingham potential is obtained. Moreover, there are still many other potential models extended from these two models [10–13]. These models are very simple in form. However, in complex systems, the vdW interaction can be affected by the local bonding environment and screening effect [14, 15], thus making it harder to correctly characterize it with semiempirical methods. In general, the most acknowledged method to characterize accurate vdW interactions is still perturbation theory (PT) calculations [16, 17].
One of the basic concerns related to the interlayer vdW interaction is exfoliation. Examining all the available exfoliation technologies, one would find that exfoliation is always a complicated process in the sense of mechanics [18]. The detailed mechanical process is very important but rarely visited. To obtain a general sense of the difficulty level of exfoliation, exfoliation energy is widely used, which is defined as the energy to extract one layer from a bulk material [19]. In calculations or simulations, interlayer binding energies are used instead since they are more computationally convenient (the two concepts are closely related and differ slightly in value [17]). While the binding energy carries the information for how much energy is needed to complete exfoliation when exfoliation is known to succeed, it cannot directly give the critical mechanical condition. However, while the interlayer binding energy (or cohesive energy in some cases) was used, calculated, and measured in many cases, the interlayer binding force usually only receives its attention in fundamental studies focusing on interlayer interaction [14, 20].
For 2D materials, although adhesion energy measurement is viable with fracture mechanics approaches [21], adhesion force measurement is often more direct and intuitive. As a result, many works study interlayer binding by measuring the force and converting it to energy with approximate models. For example, a recent study measured the critical force to split a square mesa made of 2D materials using atomic force microscopy (AFM), and the interfacial adhesion energy was obtained by integrating the retraction curves [22]. This can be seen as an approximation to the integration of the real interlayer force curve, but it can be significantly overestimated (especially when the spring constant of the probe is low). More reliable measurements often rely on classical continuum mechanics models such as JKR (JohnsonKendallRoberts) [23], DMT (DerjaguinMullerToporov) [24], MaugisDugdale [25], and MYD (MullerYushchenkoDerjaguin) [26] to extract adhesion energy from pulloff force [27–29]. The JKR model includes only the “inside adhesion” (within a defined area) and is only applicable for soft substrates. The DMT model uses a Hertzian profile for repulsion, which is far from realistic and is only applicable for very hard substrates. The MaugisDugdale model uses a simple Dugdale model and has a variable coefficient that can be tuned for different circumstances. The MaugisDugdale model can be reduced to either the JKR or the DMT models by tuning this coefficient and can be viewed as a transitional model between these two models. Obviously, the determination of the coefficient is vital to the accuracy of the model. Since it cannot be derived solely from forcedistance curves, many studies obtain this coefficient by fitting an empirical equation proposed by Carpick et al. [30]. The MYD model is another step taken for filling the gaps between the JKR and DMT models by using the Lennard–Jones potential, which is considered more realistic. Similarly, there are parameters that must be predetermined for MYDlike models, such as the equilibrium distance. Different models can be chosen for different cases to reach optimal accuracy. However, these models are approximations and can sometimes induce unexpected errors. Adhesion energies measured using these models are often not comparable to other similar studies due to the differences in experimental details (such as tip curvature), and the directly measured forces are not compared to reliable theoretical references.
Reliable and abundant interlayer force data from accurate theoretical predictions are urgently needed. In this study, we aim to provide such a reliable data set with reliable firstprinciple calculations. Targeting 230 different common 2D materials provided by Mounet’s 2D material database [19], a simple effective descent algorithm is designed and implemented to work with the Vienna Ab initio Simulation Package (VASP) code to perform firstprinciple calculations for the critical interlayer binding force and energy . Both the optB88vdW and vdWDF2 functionals are used, and the results are compared with accurate random phase approximation (RPA) calculations. Experimental measurements are carried out by coating 2D materials onto homemade silica probes and measuring the force between homogeneous 2D material interfaces in AFM. The predicted results show significant divergence between the trend of binding force and binding energy, which emphasizes the irreducibility between the two quantities and hence the meaning for providing the binding force data. The measured data show excellent agreement with the predictions.
2. Results
2.1. Prediction of Interlayer Binding Energy and Force
To locate the minimum interlayer binding force and energy, we present a simple descent algorithm, as shown briefly in Figure 1(a). The binding energy and force are calculated in bilayers, which is different from many other works since the binding force in the bulk structure cannot be calculated with periodic boundary conditions. To capture the dispersive interaction [17], two wellknown nonlocal vdW functionals are adopted: optB88 (an optimized vdWDF functional) [31] and the vdWDF2 [32] functional. An example output with the vdWDF2 functional of the algorithm implementation is shown in Figure 1(b).
Figure 1(c) shows the and calculated with the vdWDF2 functional for 230 different 2D materials provided by Mounet’s 2D database [19], which is marked as “easily exfoliable.” The data points are located within a band around a straight line, which indicates loose linearlike dependencies between and with notable divergence. The calculated and for many materials share a similar trend because they are gathered in the center of the band where the data points severely overlap with each other. However, significant divergence can be found for many systems located away from the center. For example, two different 2D materials with the same or similar binding energies can possess a difference in binding forces as large as 5 eV Å^{3}. For quite a few of the systems, the trend of can even go against the trend of (C and As e.g.). The results calculated with the optB88vdW functional show a similar divergence (which can be found in Supplementary Figure 1). The subplots at the top and right sides of Figure 1(c) show the distribution of the two quantities. More than 80% of the calculated 2D materials possess interlayer binding energy within the range of 10 to 20 eV Å^{2}; approximately 90% of the calculated 2D materials possess interlayer binding force within the range of 5 to 10 eV Å^{3}. with the equilibrium distance and with the critical distance calculated from the vdWDF2 functional for a few common 2D materials are listed in Table 1, and the full results calculated with the vdWDF2 and optB88vdW functionals can be found in Supplementary Table 1 and Supplementary Table 2, respectively. The distribution of and is shown in Supplementary Figure 4.

It is worth noting from the calculated results that the equilibrium interlayer distance calculated with vdWDF2 might be slightly overestimated, as it has been reported as one of the drawbacks of the vdWDF methods [17, 33, 34]. Taking graphene as an example, the experimentally observed interlayer distance for graphite is approximately 3.35 Å, which is significantly smaller than the predicted value of 3.66 Å. However, one should be aware that the experimental value is given for the AB stacking mode, while the predicted value is for the AA stacking mode, and the AA stacking mode usually results in a discernably larger distance due to stronger repulsion. This difference caused by the stacking mode has been demonstrated by other studies [33, 34]. A reliable reference for the AA stacking interlayer distance is 3.42 Å from the RPA calculation for bulk graphite [33]. Considering that the interlayer distance increases with the decrease in the number of layers [35], the actual overestimation of the predicted value in this study should be notably less than 7.0%. These assessments suggest that the distances calculated can be slightly overestimated, but not significantly.
The calculated from the two different vdW functionals are compared with each other and with the results calculated from an RPA calculation [17] for some systems. The comparison is shown in Figure 1(d). Obviously, both results share a similar trend with that of the RPA calculations. However, the calculated with the vdWDF2 functional is closer to the RPA results, while the calculated with the optB88vdW functional is notably larger than that of the vdWDF2 functional. This drawback of the optB88vdW functional has also been reported by other studies [36, 37]. The vdWDF2 results are slightly lower than the RPA references since the RPA results are for bulk systems. Figure 1(e) shows the comparison between the forces calculated with the two functionals. It is obvious that optB88vdW tends to overestimate as well, but the overestimation is slighter.
2.2. The Loose Relationship between , , and
Assuming the interlayer binding energy is modeled by the traditional n6 potential where the repulsion term takes the form of , the interlayer binding force is then represented as its derivative. By solving , a simple relationship between and is then achieved (details can be found in the Supplementary Information), where is a constant determined by . Although the n6 potential might not fit equally well for all kinds of 2D materials, this expression might still be good for rough estimations. To verify this idea, the calculated and are fitted with a line through the origin, as shown in Figure 2(a) and Figure 2(b) for the vdWDF2 and optB88vdW results, respectively. During the fitting of the data, it was found that the data points with small tend to stray away from the main trend, while other data points fit comparatively well with the line. These abnormal data points with exceptionally small values, as labeled in the graph, are mostly 2D hydroxides, which in their natural state should have hydrates in between layers. It should be noted that the results presented here do not involve hydrates between layers. To obtain a clearer trend of the data, only the data points with larger than 2 Å are included in the fitting. For the vdWDF2 results, most of the data points are gathered closely (with the standard deviation being just 1.12 eV Å^{3}) around the fitting line with a slope of 1.56. Similar trends can be seen for the optB88vdW results, but the slope is much smaller (1.29), which can be attributed to the significantly overestimated compared to the slight overestimation of , as mentioned above. The good agreement with the calculated results suggests that equation (1) can be used to roughly predict the interlayer binding force from the equilibrium distance and binding energy.
As the expression is derived solely from the n6 potential, the deviation might infer the inherent inaccuracy of the potential model. To determine the details, the interlayer binding energies produced by the descent algorithm are fitted with various potential models. Here, we choose Ca(OH)_{2} (, which cannot fit in Figure 1(a)) and graphene (which fits well in Figure 1(a)) from the vdWDF2 calculations to fit three potential models: the n6 model, the Buckingham model, and the modified Buckingham model [38]. The expressions of the three potential models are shown in the following equations: where , , , , _{2,} and are constant parameters to be determined. Global optimization is used to search for the best parameters. The fitting results of C and Ca(OH)_{2} are shown in Figure 2(c) and Figure 2(d), respectively. The n6 potential fits well for graphene but is comparatively poor for Ca(OH)_{2}, which may be the reason why the Ca(OH)_{2} strays from the main trend in Figure 1(a). The failure of the n6 potential in Ca(OH)_{2} is common in systems with small . In such cases, the variation in the exchange energy can be very drastic, whereas the Born–Mayer term has an exponential form, can model the shortrange exchange effect more reasonably, thus fitting better for many systems [9, 39–41]. The Buckingham potential is just a direct modification from the n6 potential by replacing the term with the Born–Mayer term. The fitting with the Buckingham potential in Figure 2(d) shows a notable improvement. Moreover, having more dispersive terms, such as the modified Buckingham potential in equation (4), can be advantageous, as shown in Figure 2(d), but the benefits are less significant. These analyses demonstrate that the deviation of the calculated results from the proposed relationship (equation (1)) is mainly due to the inherent inaccuracy brought by the inversepower repelling term. For this reason, models based on n6like potentials (such as the LJ potential) for converting to or vice versa (the Maugis model [25], e.g.) should not be considered accurate methods. This applies to the loose relationship presented in equation (1), which should only be used for rough estimations.
2.3. Experimental Measurements of Interlayer Cohesive Force
Experimental evidence is needed to confirm the validity of the predictions. Since corresponding measurements are few and some might be inaccurate, an experiment is designed by using atomic force microscopy (AFM) with 2Dmaterialcoated probes. Three different 2D materials with chemical formulas of C, BN, and In_{2}Se_{3} are chosen for measurements. The key point of this method is to fabricate AFM probes with appropriate tip curvature (silica balls are used to achieve this purpose) and wrap the 2D materials onto them. The probe fabrication process using glue for welding resembles the method reported by Li et al. [42], but a highstrength type of glue and a homemade micromanipulator are used, which enable us to fabricate the probes more efficiently and precisely. The results are shown in Figure 3. The fabrication process is shown in Figure 3(a), which is further illustrated in the methods section. The advantage of this fabrication process is that it will not contaminate the cantilevers with extra glue since all glue is applied precisely to where it needs to be using the micromanipulator, and the upward coating of 2D materials with soft substrates can make them wrap more seamlessly onto the ball.
The design of the measurement and the analysis of feasibility can be found in Section 2 of the Supplementary Information. Figures 3(e)–3(g) show the distribution of the measured data. The data fit very well against the Gaussian distribution, with the expected values of , , and . These values are then corrected for radius differences and converted to strength (perunitarea forces) with an effective contact area of ~89.00 nm^{2} (details in Section 2.4 of the Supplementary Information), which is shown as a comparison to the calculated results in Table 2. The statistical results show a trend: , which is consistent with the theoretical prediction, and the values are of the same magnitude as those of other similar studies [20]. The corrected forces are very close to the uncorrected values since the radii of the silica balls are rather uniform. The estimated strengths are of the same magnitude as the calculation results, proving that the measured results are quite reliable and can correctly reflect the trend of the actual interlayer binding force.
 
This is only the result of a rough correction based on the LJ model. 
Since the relationship predicted between the binding energies of the three different materials, , is consistent with other theoretical studies [17, 19], together with the force relationship measured as , the inconsistency between and is also proven. The of BN and In_{2}Se_{3} are similar, but the differs quite notably. This inconsistency between and emphasizes that attention should be given to binding forces.
3. Discussion
In this study, highthroughput calculations are carried out for 230 2D materials to determine their maximum binding energies and forces. Three materials are chosen to perform experimental interlayer cohesion force measurements. The calculation results show excellent agreement with both accurate RPA calculations and the measured results. The predicted critical binding forces for 230 different 2D materials validated with experimental measurements can serve as a reliable source of references for future studies related to interlayer adhesion. It can be especially helpful for indirect measurements of interlayer adhesion energy from adhesion force since the data provided by this work can be used to check with the directly measured forces before conversion to ensure the correctness of the direct data in the first place. The proposed method for highthroughput calculation and experimental measurement of interlayer binding forces can be easily extended to other 2D materials. In a more general sense, we hope the results of this work can help future studies to better understand the characteristics of vdW interactions in 2D materials and to build accurate mechanical models for exfoliation.
4. Materials and Methods
4.1. Calculation Parameters for VASP Calls
The firstprinciple calculations were performed with the VASP [43] code using the projectoraugmented planewave (PAW) method. The PBE [44] functional was adopted to address the exchange interaction, while the optB88vdW [31] and vdWDF2 [32] functionals were adopted to address the correlation interaction. For all the calculations, the energy cutoff of the planewave basis was set to 550 eV, the first Brillouin zone was sampled with a Monkhorst Pack grid, and the convergence criteria for selfconsistent field loops were set to 10^{5} eV. For monolayer relaxation, only the inplane adjustments were allowed with an atomic force tolerance of 0.01 eV Å^{1}. For static calculations of monolayers and bilayers, the highest accuracy has been enabled to achieve more accurate atomic forces.
4.2. Descent Algorithm for Locating and
The algorithm assumes that the interlayer energy or force profile is a singlevalley curve and searches for the minimum by ensuring that every step reduces the energy or force. Upon convergence, the maximum binding energy is yielded at the equilibrium distance , while the critical binding force is yielded at the critical distance . The algorithm is implemented in Python to work with VASP. The results yielded for demonstration are shown in Figure 1(b). For most cases, the algorithm is capable of locating and with approximately 15 VASP calls, which is efficient enough for this work. The details of how the routine was called are explained in the next subsection.
4.3. HighThroughput Calculation
The whole calculation was automated with Python (3.6.3) scripting, and the descent algorithm was implemented with Python. Bilayer systems are constructed with a custom Python library PyCrystal, which provides various lattice operations. For each 2D system, monolayer relaxation and monolayer static calculations were carried out first. Then, the bilayer systems were built with relaxed structures, and the descent algorithm was started to locate the minimum binding energy and atomic force. During the minimum search, both the energy and force data were produced. is searched first, the algorithm starts from and searches toward , with a convergence criterion of 10^{4} eV per atom. Then, the produced data are fed to the search for , which starts from the deepest valley of the provided data. The criterion for convergence of was set to 10^{3} eV Å^{1} per atom. To avoid local convergence while there are multiple valleys in the or profiles (there are a few such cases in the optB88vdW calculations), a check for global convergence is performed before real convergence. To ensure sufficient numerical accuracy, the convergence criterion was set to 10^{4} eV per atom for the binding energy and 10^{3} eV Å^{1} per atom for the binding force.
4.4. Global Search for Fitting Parameters
The simulated annealing algorithm for potential fitting is implemented in C. Exponential temperature management is used with a base value of 0.999, the initial temperature is 300 K, and the critical temperature is 10^{5} K. The algorithm searches 10000 times in every temperature step and breaks at the critical temperature for divergence. Convergence is reached when the optimized values remain the same for 100 runs. During the fitting of the potentials, both the global fitting technology and the traditional trustregionbased local search are used together to ensure the best fitting results.
4.5. Fabrication of 2DMaterialCoated Probes
Before starting, a homemade micromanipulator with a tip size of ~10 μm was equipped onto the 2D material transfer system. Using the micromanipulator, a small amount of epoxy glue was precisely applied on the pointed end of a tipless cantilever. Then, the micromanipulator was replaced with a clean one to take a silica ball with a radius of ~4.5 μm from the glass slide (this is possible because the adhesion between the silica balls and the glass slide is much weaker than the adhesion between the micromanipulator and the silica balls). Aligning the silica ball on the micromanipulator with the previously glued cantilever, the silica ball was carefully attached to the cantilever. After the glue was dried (>24 h), a small amount of glue was precisely applied onto the top of the previously glued silica ball. Then, the 2D material flake prepared on polydimethylsiloxane (PDMS) was loaded in whole in the transfer system and slowly pressed against the ball probe while keeping the alignment at the same time. It is pressed down to a certain depth after initial contact to ensure full contact between the 2D material flake and silica ball. In this way, extra glue will be squeezed out, and the residual glue between the interface will be very thin so that its influence on the curvature of the 2D material flake can be neglected. Finally, the glass slide carrying the PDMS was immediately lifted to detach the 2D flake and leave it on the silica ball. The 2D flakes can detach from the PDMS since the glue has a very high tensile strength (sticky), even when it has not yet dried.
4.6. Interlayer Cohesion Force Measurements
The 2D material transfer system is E1T, manufactured by MetaTest Co. Highly oriented pyrolytic graphite (HOPG) is produced by the Institute of Metals, Chinese Academy of Sciences; highly oriented hexagonal BN is purchased from Onway Technology Co.; αIn_{2}Se_{3} (2H phase) is produced by HQ Graphene Co. The tipless cantilever is TL_CONT, Nanosensors. The spring constant given by the manufacturer is 0.020.77 N/m. Highquality silica balls were acquired from NanoMicro, Co. All three different materials were characterized by Xray diffraction (XRD). The results are shown in Supplementary Figure 8. The 2D material flakes were prepared with mechanical exfoliation using Scotch tape. The size (slightly smaller than the size of the silica balls) and thickness (>20 nm as demanded) of the flakes were filtered by AFM morphology characterization to find the flakes appropriate for coating. Smooth large flakes are selected and transferred onto silicon substrates using the 2D material transfer system. The epoxy glue used was T200 manufactured by Ergo Co., Ltd. The fabricated probes were characterized with scanning electron microscopy (SEM), and more images can be found in Supplementary Figure 9 and Supplementary Figure 10. The measurements were carried out using Asylum Research MFP3D Origin+ AFM in tapping mode. For each kind of 2D material, 2 to 3 large flakes that are the smoothest were chosen as the substrate. On these selected flakes, several different areas (each with a size of ) were sampled in the measurements. The sample points without a successful jumpoff were not accounted for as valid data and were not included in the statistics. All measurements were carried out under ambient conditions. The temperature was 2127°C, and the relative humidity was 20%25%. The humidity is very low, and the capillary effects are assumed to be negligible, as found in many studies [45–47]. In addition, no longtail characteristic caused by water bridges [48] is found in the forcedistance curves, which supports this assumption. Force measurements and SEM characterizations were performed alternately. After the first measurement and SEM characterization, a second measurement (for verification only) and SEM characterization were carried out. The measurement results before and after the first SEM characterization are highly identical, which supports that charge pollution induced by the electron beam is negligible in our case [47]. The SEM characterizations before and after the measurements are shown in Supplementary Figure 11.
Data Availability
The main data not shown in the main text can be found in the Supplementary Information. Other data, including the initial structures of the 2D materials, the data produced by each step of the algorithm, and the measurement results, were published in the Materials Cloud Archive (doi:10.24435/materialscloud:t0hn). The implementation of the descent algorithm is a single file named ‘descender.py’ in the data published in the Materials Cloud Archive. The code that automated the calculation process is not publicly available. The custom code PyCrystal is licensed under GPL v3.0, which can be found at https://github.com/CementMixer/pycrystal.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this article.
Authors’ Contributions
K. T. and W. Q. conceptualized the study. K. T. designed and carried out the calculations, designed and implemented the algorithm, and performed the scripting work; K. T. and Y. W. designed and performed the experiments; K. T. and G. R. carried out the characterizations; and W. L. and W. Q. supervised the work and revised the manuscript.
Acknowledgments
We acknowledge the support from the National Natural Science Foundation of China (Grant No. 52072308), the Fundamental Research Funds for the Central Universities (Grant Nos. 3102021MS0404 and 3102019JC001) and Open Testing Foundation of the Analytical & Testing Center of Northwestern Polytechnical University.
Supplementary Materials
Supplementary Figure 1:  calculated with the optB88vdW functional. Supplementary Figure 2: distribution of energy error . Supplementary Figure 3: distribution of . Supplementary Figure 4: distribution of the calculated and . Supplementary Figure 5: diagram of the geometry for the integration model. Supplementary Figure 6: visualization of Supplementary equation (13). Supplementary Figure 7: visualization of with . Supplementary Figure 8: XRD characterization of the three different 2D materials. Supplementary Figure 9: SEM characterization (with BSE) of the three fabricated probes. Supplementary Figure 10: highresolution SEM characterization (with SE) of the fabricated probes. Supplementary Figure 11: comparison of the coated 2D material flakes before and after measurements. Supplementary Table 1: the binding energy and binding force calculated with the vdWDF2 functional. Supplementary Table 2: the binding energy and binding force calculated with the optB88vdW functional. (Supplementary Materials)
References
 R. Frisenda, E. NavarroMoratalla, P. Gant et al., “Recent progress in the assembly of nanodevices and van der Waals heterostructures by deterministic placement of 2D materials,” Chemical Society Reviews, vol. 47, no. 1, pp. 53–68, 2018. View at: Publisher Site  Google Scholar
 K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. Castro Neto, “2D materials and van der Waals heterostructures,” Science, vol. 353, no. 6298, article aac9439, 2016. View at: Publisher Site  Google Scholar
 Y. Liu, Y. Huang, and X. Duan, “Van der Waals integration before and beyond twodimensional materials,” Nature, vol. 567, no. 7748, pp. 323–333, 2019. View at: Publisher Site  Google Scholar
 K. Tang and W. Qi, “Moirépatterntuned electronic structures of van der Waals heterostructures,” Advanced Functional Materials, vol. 30, no. 32, article 2002672, 2020. View at: Publisher Site  Google Scholar
 R. Bistritzer and A. H. MacDonald, “Moire bands in twisted doublelayer graphene,” Proceedings of the National Academy of Sciences, vol. 108, no. 30, pp. 12233–12237, 2011. View at: Publisher Site  Google Scholar
 I. Brihuega, P. Mallet, H. GonzálezHerrero et al., “Unraveling the intrinsic and robust nature of van Hove singularities in twisted bilayer graphene by scanning tunneling microscopy and theoretical analysis,” Physical Review Letters, vol. 109, no. 19, article 196802, 2012. View at: Publisher Site  Google Scholar
 A. Dalgarno and J. Lewis, “The representation of long range forces by series expansions I: the divergence of the series II: the complete perturbation calculation of long range forces,” Proceedings of the Physical Society, vol. 69, no. 1, pp. 57–64, 1956. View at: Publisher Site  Google Scholar
 S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, “A consistent and accurate ab initio parametrization of density functional dispersion correction (DFTD) for the 94 elements HPu,” The Journal of Chemical Physics, vol. 132, no. 15, article 154104, 2010. View at: Publisher Site  Google Scholar
 A. A. Abrahamson, “BornMayertype interatomic potential for neutral groundstate atoms with Z=2 to Z=105,” Physics Review, vol. 178, no. 1, pp. 76–79, 1969. View at: Publisher Site  Google Scholar
 K. T. Tang and J. P. Toennies, “A simple theoretical model for the van der Waals potential at intermediate distances. I. Spherically symmetric potentials,” The Journal of Chemical Physics, vol. 66, no. 4, pp. 1496–1506, 1977. View at: Publisher Site  Google Scholar
 U. Buck, “Inversion of molecular scattering data,” Reviews of Modern Physics, vol. 46, no. 2, pp. 369–389, 1974. View at: Publisher Site  Google Scholar
 A. F. Wagner, G. Das, and A. C. Wahl, “Calculated longrange interactions and low energy scattering of Ar–H,” The Journal of Chemical Physics, vol. 60, no. 5, pp. 1885–1891, 1974. View at: Publisher Site  Google Scholar
 Z. Kozioł, G. Gawlik, and J. Jagielski, “Van der Waals interlayer potential of graphitic structures: From Lennard–Jones to Kolmogorov–Crespy and Lebedeva models,” Chinese Physics B, vol. 28, no. 9, article 096101, 2019. View at: Publisher Site  Google Scholar
 X. Liu, J. Yang, and W. Guo, “Semiempirical van der Waals method for twodimensional materials with incorporated dielectric functions,” Physical Review B, vol. 101, no. 4, article 045428, 2020. View at: Publisher Site  Google Scholar
 X. Liu, Z. Zhang, and W. Guo, “van der Waals screening by graphenelike monolayers,” Physical Review B, vol. 97, no. 24, article 241411, 2018. View at: Publisher Site  Google Scholar
 B. Jeziorski, R. Moszynski, and K. Szalewicz, “Perturbation theory approach to intermolecular potential energy surfaces of van der Waals complexes,” Chemical Reviews, vol. 94, no. 7, pp. 1887–1930, 1994. View at: Publisher Site  Google Scholar
 T. Björkman, A. Gulans, A. V. Krasheninnikov, and R. M. Nieminen, “van der Waals bonding in layered compounds from advanced densityfunctional firstprinciples calculations,” Physical Review Letters, vol. 108, no. 23, article 235502, 2012. View at: Publisher Site  Google Scholar
 E. Gao, S. Z. Lin, Z. Qin, M. J. Buehler, X. Q. Feng, and Z. Xu, “Mechanical exfoliation of twodimensional materials,” Journal of the Mechanics and Physics of Solids, vol. 115, pp. 248–262, 2018. View at: Publisher Site  Google Scholar
 N. Mounet, M. Gibertini, P. Schwaller et al., “Twodimensional materials from highthroughput computational exfoliation of experimentally known compounds,” Nature Nanotechnology, vol. 13, no. 3, pp. 246–252, 2018. View at: Publisher Site  Google Scholar
 B. Li, J. Yin, X. Liu et al., “Probing van der Waals interactions at twodimensional heterointerfaces,” Nature Nanotechnology, vol. 14, no. 6, pp. 567–572, 2019. View at: Publisher Site  Google Scholar
 K. M. Liechti, “Characterizing the interfacial behavior of 2D materials: a review,” Experimental Mechanics, vol. 59, no. 3, pp. 395–412, 2019. View at: Publisher Site  Google Scholar
 H. Rokni and W. Lu, “Direct measurements of interfacial adhesion in 2D materials and van der Waals heterostructures in ambient air,” Nature Communications, vol. 11, no. 1, article 5607, 2020. View at: Publisher Site  Google Scholar
 K. L. Johnson, K. Kendall, A. D. Roberts, and D. Tabor, “Surface energy and the contact of elastic solids,” Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, vol. 324, no. 1558, pp. 301–313, 1971. View at: Publisher Site  Google Scholar
 B. V. Derjaguin, V. M. Muller, and Y. P. Toporov, “Effect of contact deformations on the adhesion of particles,” Progress in Surface Science, vol. 45, no. 14, pp. 131–143, 1994. View at: Publisher Site  Google Scholar
 D. Maugis, “Adhesion of spheres: the JKRDMT transition using a dugdale model,” Journal of colloid and interface science, vol. 150, no. 1, pp. 243–269, 1992. View at: Publisher Site  Google Scholar
 V. M. Muller, V. S. Yushchenko, and B. V. Derjaguin, “On the influence of molecular forces on the deformation of an elastic sphere and its sticking to a rigid plane,” Journal of Colloid and Interface Science, vol. 77, no. 1, pp. 91–101, 1980. View at: Publisher Site  Google Scholar
 Y. Jiang and K. T. Turner, “Measurement of the strength and range of adhesion using atomic force microscopy,” Extreme Mechanics Letters, vol. 9, pp. 119–126, 2016. View at: Publisher Site  Google Scholar
 Y. Li, S. Huang, C. Wei, C. Wu, and V. N. Mochalin, “Adhesion of twodimensional titanium carbides (MXenes) and graphene to silicon,” Nature Communications, vol. 10, no. 1, p. 3014, 2019. View at: Publisher Site  Google Scholar
 T. Jiang and Y. Zhu, “Measuring graphene adhesion using atomic force microscopy with a microsphere tip,” Nanoscale, vol. 7, no. 24, pp. 10760–10766, 2015. View at: Publisher Site  Google Scholar
 R. W. Carpick, D. F. Ogletree, and M. Salmeron, “A general equation for fitting contact area and friction vs load measurements,” Journal of colloid and interface science, vol. 211, no. 2, pp. 395–400, 1999. View at: Publisher Site  Google Scholar
 T. Thonhauser, V. R. Cooper, S. Li, A. Puzder, P. Hyldgaard, and D. C. Langreth, “Van der Waals density functional: selfconsistent potential and the nature of the van der Waals bond,” Physical Review B, vol. 76, no. 12, article 125112, 2007. View at: Publisher Site  Google Scholar
 K. Lee, É. D. Murray, L. Kong, B. I. Lundqvist, and D. C. Langreth, “Higheraccuracy van der Waals density functional,” Physical Review B, vol. 82, no. 8, article 081101, 2010. View at: Publisher Site  Google Scholar
 S. Lebègue, J. Harl, T. Gould, J. G. Ángyán, G. Kresse, and J. F. Dobson, “Cohesive properties and asymptotics of the dispersion interaction in graphite by the random phase Approximation,” Physical Review Letters, vol. 105, no. 19, article 196401, 2010. View at: Publisher Site  Google Scholar
 M. S. Alam, J. Lin, and M. Saito, “Firstprinciples calculation of the interlayer distance of the twolayer graphene,” Japanese Journal of Applied Physics, vol. 50, no. 8, article 080213, 2011. View at: Publisher Site  Google Scholar
 Y. Saito, T. Yoshikawa, S. Bandow, M. Tomita, and T. Hayashi, “Interlayer spacings in carbon nanotubes,” Physical Review B, vol. 48, no. 3, pp. 1907–1909, 1993. View at: Publisher Site  Google Scholar
 Z. Wang, S. M. Selbach, and T. Grande, “van der Waals density functional study of the energetics of alkali metal intercalation in graphite,” RSC Advances, vol. 4, no. 8, pp. 3973–3983, 2014. View at: Publisher Site  Google Scholar
 I. Hamada, “van der Waals density functional made accurate,” Physical Review B, vol. 89, no. 12, article 121103, 2014. View at: Publisher Site  Google Scholar
 R. Heller, “Theory of some van der Waals molecules,” The Journal of Chemical Physics, vol. 9, no. 2, pp. 154–163, 1941. View at: Publisher Site  Google Scholar
 V. I. Gaydaenko and V. K. Nikulin, “BornMayer interatomic potential for atoms with to ,” Chemical Physics Letters, vol. 7, no. 3, pp. 360–362, 1970. View at: Publisher Site  Google Scholar
 K. Tang, W. Qi, Y. Li, and T. Wang, “Electronic properties of van der Waals heterostructure of black phosphorus and MoS_{2},” Journal of Physical Chemistry C, vol. 122, no. 12, pp. 7027–7032, 2018. View at: Publisher Site  Google Scholar
 K. Tang, W. Qi, Y. Li, and T. Wang, “Tuning the electronic properties of van der Waals heterostructures composed of black phosphorus and graphitic SiC,” Physical Chemistry Chemical Physics, vol. 20, no. 46, pp. 29333–29340, 2018. View at: Publisher Site  Google Scholar
 J. Li, J. Li, L. Jiang, and J. Luo, “Fabrication of a graphene layer probe to measure force interactions in layered heterojunctions,” Nanoscale, vol. 12, no. 9, pp. 5435–5443, 2020. View at: Publisher Site  Google Scholar
 G. Kresse and J. Furthmuller, “Efficient iterative schemes for ab initio totalenergy calculations using a planewave basis set,” Physical Review B, vol. 54, no. 16, pp. 11169–11186, 1996. View at: Publisher Site  Google Scholar
 J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple,” Physical Review Letters, vol. 77, no. 18, pp. 3865–3868, 1997. View at: Publisher Site  Google Scholar
 F. Podczeck, J. M. Newton, and M. B. James, “Influence of relative humidity of storage air on the adhesion and autoadhesion of micronized particles to particulate and compacted powder surfaces,” Journal of colloid and interface science, vol. 187, no. 2, pp. 484–491, 1997. View at: Publisher Site  Google Scholar
 L.O. Heim, J. Blum, M. Preuss, and H.J. Butt, “Adhesion and friction forces between spherical micrometersized particles,” Physical Review Letters, vol. 83, no. 16, pp. 3328–3331, 1999. View at: Publisher Site  Google Scholar
 B. Li, X. Liu, and W. Guo, “Probing interactions at twodimensional heterointerfaces by boron nitridewrapped tip,” Nano Research, vol. 14, no. 3, pp. 692–698, 2021. View at: Publisher Site  Google Scholar
 J. W. Suk, S. R. Na, R. J. Stromberg et al., “Probing the adhesion interactions of graphene on silicon oxide by nanoindentation,” Carbon, vol. 103, pp. 63–72, 2016. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2022 Kewei Tang et al. Exclusive Licensee Science and Technology Review Publishing House. Distributed under a Creative Commons Attribution License (CC BY 4.0).