Methodology

In this work, we address this issue by going beyond the dipole approximation in $$ {\widehat{\mathcal{H}}}_{\mathrm{DC}}$$ while retaining an efficient atomistic formalism. Given that the equivalent of Equation (1) with the full Coulomb interaction does not allow for a straightforward closed-form solution, we derive and evaluate the correction toward full Coulomb coupling, $$ {\widehat{V}}^{\prime }=f_{\mathrm{damp}}\left({\widehat{V}}_{\mathrm{Coul}}-{\widehat{T}}_{\mathbf{pp}}\right)$$, to first order in perturbation theory,

Typically, the electronic Coulomb energy is divided into its classical and correlation parts in electronic-structure theory. Likewise, we also divide the oscillator Coulomb energy into its classical component, J[ρ], and the correlation energy, E_corr[Ψ], and denote $$ E_{\mathrm{dip}}\left[\mathrm{\Psi }\right]=⟨\mathrm{\Psi }\mid {\widehat{T}}_{\mathbf{pp}}\mid \mathrm{\Psi }⟩$$. As E_corr[Ψ₀] = E_dip[Ψ₀] = 0, the first-order full Coulomb contribution can be written as

Fig. 1

Schematic representation of dipole-correlated Coulomb singles (DCS).

a Green arrows represent dipole coupling between electronic fragments. First-order perturbation theory (PT) on top of the many-body dispersion formalism (MBD) captures the interaction energy, E_DCS, between δρ_A and δρ_B, depicted by field lines. b E_DCS vanishes in 3D isotropic vacuum because of symmetry. c Under rotational symmetry-breaking confinement, electric field lines between electronic fragments deform, which leads to E_DCS ≠ 0. d Further inclusion of higher-order terms leads to full Coulomb-coupled vdW interaction.

We note that in the presence of a polarizable environment, the DCS interaction between two bodies can decay asymptotically slower than the zeroth-order MBD energy. For example, the dipole correlation in the MBD ground-state wavefunction between a finite body and its environment induces permanent quadrupole moments on the oscillators, causing the resulting interaction to decay as R⁻⁵. In contrast, the MBD interaction energy between two bodies in such a system decays asymptotically as R⁻⁶. The physical origin of R⁻⁵-dependent repulsive interactions between oscillators with reduced dimensionality as reported in ref. ³² has been controversially debated^46,47. The generalized explanations and results for realistic systems reported in this publication finally resolve this discussion and clearly show that this leading-order behavior is not owing to purely electrostatic interactions, but originates from dispersion-induced electron density polarization effects. From this discussion, it is also clear that the DCS contribution is a leading-order term, which in general could be comparable to or even more important than the renowned higher-order multipolar terms (dipole–quadrupole and higher terms⁴⁴) in the vdW dispersion energy. Higher-order terms in the perturbation theory used here are naturally lower in magnitude and show a comparably quick decay with interatomic separations³². This justifies the limitation to the more-slowly decaying first-order DCS term.

The accurate prediction of interaction energies for non-covalently bound materials is an ongoing challenge for researchers and workhorse density-functional theory (DFT) methods are at the forefront of development efforts. The complex balance of intermolecular interactions is particularly challenging to predict and what is more, the target accuracy is generally in the range of a few kJ/mol in interaction energies. Our approach to better accuracy is to improve the physical basis of theoretical methods—here, by incorporating the DCS contributions. First, we combine DCS with MBD in MBD+DCS to compute interaction energies of small molecular dimers. Following this, MBD+DCS is used to compute the binding energies of supramolecular host–guest complexes and confined Xe dimers in capped CNTs. These larger confined systems reveal the impact of E_DCS on binding energies, as well as the length scale and character of the emergent changes to long-range interactions.

DCS in small molecular dimers

We begin by applying the first-order perturbation term (3) to the S66 data set⁴⁸ of small, unconfined molecular dimers. For such systems, semi-local or hybrid density-functional approximations in conjunction with the MBD formalism are designed to provide excellent agreement with accurate reference results from coupled-cluster calculations with single, double, and perturbative triple excitations (CCSD(T))³³. The S66 set contains non-covalently interacting dimers in 3D isotropic vacuum and in accordance with Fig. 1, we expect minuscule first-order Coulomb corrections in this case. Indeed, we find that the DCS contributions for all systems in S66 are very small and they can have both positive as well as negative values (see Supplementary Fig. 2). As a result, the accuracy of vdW-inclusive DFT remains equally good upon account for DCS contributions to the interaction energies of small molecular dimers as contained in S66.

Coulomb corrections for host–guest complexes

Host–guest molecular systems are significantly more complex than the S66 dimers, but are still tractable with accurate benchmark methods such as diffusion quantum Monte-Carlo (DQMC). Here, we first ascertain the impact of DCS on the binding energies of host–guest complexes and demonstrate the accuracy of PBE0+MBD+DCS with respect to DQMC reference interaction energies from previous works. Fig. 2 shows two examples of host–guest complexes. For both systems, the guest molecule is a C₇₀-fullerene (buckyball) while the host is either [6]-cycloparaphenyleneacetylene (6-CPPA, Fig. 2 left) or the “buckyball-catcher” molecule (Fig. 2 right)⁴⁹. In the framework of Fig. 1, the host molecule serves as both confinement and the interaction partner. The 6-CPPA and catcher molecules provide a different confining environment for the buckyball by virtue of their geometry. We therefore focus on these systems to showcase the contribution of E_DCS in confined systems.

Fig. 2

Binding energies and dispersion–polarization of a C₇₀ fullerene different host molecules.

Host molecules: 6-CPPA (left), “buckyball-catcher” (right). PBE0+MBD results (orange), diffusion quantum Monte-Carlo reference (DQMC, blue line; error bars shown as boxes), PBE0+MBD including dipole-correlated Coulomb singles (DCS, black line). DQMC reference data were taken from ref. ¹⁶. The depiction of the complexes includes iso-surfaces at ±0.003 (a.u.) of the change in the density of electronic fluctuations with respect to the isolated monomers (red: decrease, blue: increase).

Within the dipole approximation, it has already been shown that many-body effects have an important role in the description of binding energies in such host–guest complexes¹⁶. Here, we show that also interactions beyond the dipolar MBD, in the form of DCS, have a significant effect and that inclusion of DCS to PBE0+MBD yields excellent agreement with the reference results from DQMC. The DCS contribution for C₇₀ in 6-CPPA and in the buckyball-catcher is 6.4 kJ/mol and 4.9 kJ/mol, respectively. Considering our findings for the S66 data set above, this clearly highlights the importance of DCS corrections once the 3D isotropy of vacuum is substantially perturbed.

The relative contribution of DCS to the total binding energy for C₇₀ in 6-CPPA and in the buckyball-catcher is 6.2% and 3.9%, respectively. As can be seen from Fig. 2 (and more clearly from Fig. 3 below), the DCS contribution does not correlate with the system size nor the vdW interaction within the dipolar approximation. However, the different E_DCS contributions can be interpreted in terms of the physics described in the dipole-coupled state, which represents the starting point (unperturbed state) for the calculation of the DCS contribution. One important factor is how the dipolar coupling changes the density of electronic fluctuations (i.e., δρ shown in Fig. 1), which can be obtained as the expectation value of the charge density operator via the wavefunction of the dipole-coupled, Ψ_DC, and non-interacting set, Ψ₀, of quantum harmonic oscillators. Analysis of the difference of δρ in the host–guest complex with respect to the isolated monomers gives us a measure of how much the density deforms upon the host–guest interaction and, therefore, forms a connection between confinement and electronic fluctuations and polarizability. The density difference shown in Fig. 2 shows that, upon dipole coupling, the density of electronic fluctuations for C₇₀ in 6-CPPA is more strongly deformed into the plane of 6-CPPA. Furthermore, we find that the overall displaced charge, i.e., the integral over the absolute density difference, can serve as a descriptor for the DCS contribution to the interaction energy: with increasing displaced charge, we observe an increased DCS contribution to the interaction energy. We point out that this also applies to the systems considered below. Hence, the electronic properties obtained for the dipole-coupled state can serve as a qualitative rule-of-thumb to estimate the magnitude of E_DCS.

Fig. 3

Dipole-correlated Coulomb singles (DCS) contributions to binding energies of ring–C₇₀ complexes and correlation to structural features.

A Binding energies for four ring–C₇₀ host–guest complexes (R1–R4): PBE0+MBD results (orange), diffusion quantum Monte-Carlo reference (DQMC, blue line; error bars as boxes), PBE0+MBD including DCS (black line). DQMC reference data taken from ref. ¹⁶. The hosts for R1–R4 are 8-CPPA rings. B Measure of axial–radial asymmetry (f_a), proximity measure (f_d), and DCS contribution (f_e) to binding energy (all values normalized to results for R4). Definition of axial and radial phenyl units of 8-CPPA via P_v plane shown in C.

Asymmetry and steric effects

To explore the connection between E_DCS and confinement, we analyze a set of geometrically similar ring–C₇₀ complexes as depicted in Fig. 3: the four complexes are C₇₀ hosted by four different conformations of 8-CPPA. In previous work, PBE0+MBD has been shown to provide reasonably accurate binding energies with respect to DQMC¹⁶. As can be seen from Fig. 3A, the addition of the DCS contribution to PBE0+MBD further improves the binding energies of all four complexes. However, the individual DCS contributions vary significantly across these conformations (see relative E_DCS as shown in Fig. 3B).

In order to study the potential role of asymmetry and steric effects for E_DCS, we define two geometrical measures: one for proximity (f_d), which is given by the sum of inverse distances between the atoms of the fullerene guest molecule and the CPPA-host, and one for the asymmetry of the system (f_a). For the latter, we define a plane along the elongated axis of C₇₀ that is perpendicular to the CPPA ring (labeled P_v plane in Fig. 3C). The four phenyl units closest to this plane are considered axial and the remaining four radial. Based on this classification, we define axial vicinity (A_∥) and radial vicinity (A_⊥) by summing the inverse distances between all fullerene atoms and atoms of the axial and radial phenyl rings, respectively. Our measure of (axial–radial) asymmetry is then given by f_a = (A_∥ − A_⊥)/(A_∥ + A_⊥). Fig. 3B summarizes the results for the proximity and asymmetry measures and the ratio of E_DCS of each system and that of R4 (f_e = E_DCS(R_i)/E_DCS(R₄), i = {1, 2, 3, 4}). It is clear that in principle, proximity has a role in electronic confinement (cf. proximity and DCS contributions for C₇₀ in 6-CPPA and R1). As can be seen from the detailed analysis in Fig. 3B, however, purely geometric considerations do not correspond directly to the trends in E_DCS. First, the proximity measure, f_d, is insensitive to the different confining environments and remains almost constant among all four conformations, whereas the DCS contribution varies significantly. Furthermore, the asymmetry measure, f_a, has no correlation with E_DCS (see R2 versus R3 and R4 in Fig. 3B). Thus, also a pairwise description of asymmetry between atomic positions is insufficient to predict the qualitative trend of the contribution of DCS.

The failure to capture the behavior of E_DCS in terms of simple geometric characteristics stems from the fact that the DCS contribution is a quantum-mechanical effect arising from long-range electron correlation, which shows a non-trivial dependence on the geometrical features of a system. A considerable part of E_DCS represents charge polarization effects due to long-range electron correlation, cf. Equation (4).

As discussed for the previous complexes, the displaced charge within MBD (as depicted in Fig. 2) can provide a measure for the dispersion–polarization-like term and indeed tracks the qualitative trend in the relative DCS interaction energies for all supramolecular complexes treated here. The best geometry-based metric found in this work is a sum of inverse distances to the power of five, which resembles an interaction of quadrupoles induced by long-range correlation (vide supra). Further information on qualitative descriptors can be found in Supplementary Fig. 3.

DCS in nanotubes

Having established the importance of E_DCS for confined host–guest systems, we now employ our methodology to investigate the effects of confinement for a Xe dimer inside CNTs. In particular, we answer the question how the presence and strength of confinement changes the importance of the DCS correction relative to dipolar vdW interactions. Xe does not possess permanent multipoles and has substantial polarizability. As a result, the long-range Xe–Xe interaction has pure vdW character.

CNTs can be generally classified according to their chiral indices (m, n), where armchair CNTs with m = n are metallic in nature. As the DCS contribution becomes especially important in the presence of metallic screening³², we analyze the Xe–Xe interaction inside two armchair, hydrogen-capped CNTs with m = n = 5 and m = n = 6. The diameters of the (5, 5)- and (6, 6)-CNTs are 6.78 Å and 8.13 Å, respectively. The length of each nanotube was chosen to be 30 Å, which is sufficient to avoid any significant edge effects. The binding energies of the Xe dimer inside the nanotubes are calculated as $$ E_{\mathrm{int}}=E_{{\mathrm{Xe}}_2\left(\mathrm{NT}\right)}+E_{\mathrm{CNT}}-E_{{\mathrm{Xe}}_{\mathrm{A}}\left(\mathrm{NT}\right)}-E_{{\mathrm{Xe}}_{\mathrm{B}}\left(\mathrm{NT}\right)}$$, where $$ E_{{\mathrm{Xe}}_2\left(\mathrm{NT}\right)}$$, $$ E_{{\mathrm{Xe}}_{\mathrm{A}/\mathrm{B}}\left(\mathrm{NT}\right)}$$, and E_CNT are the energies of Xe₂ inside the CNT, the two single Xe atoms hosted by the nanotube, and the bare CNT, respectively.

Fig. 4 summarizes the effect of the confining potential of capped CNTs on the Xe dimer binding energy. We focus on the variation of the dipole-coupled MBD and the corresponding DCS contributions as a function of the inter-Xe distance, R. In Fig. 4A, we show the effect of confinement on the individual contributions by comparing a Xe dimer inside the (6, 6)-CNT and in gas-phase. One clearly sees that both E_MBD as well as E_DCS become less attractive owing to confinement. In the case of the MBD interaction energy this can be attributed to (i) decreased Xe polarizabilities due to the screening by the CNT and (ii) the restriction of electronic fluctuations on the Xe atoms due to correlation with fluctuations in the CNT. This reduction of the vdW interaction as a result of many-body correlation has been observed and detailed in a number of previous works^12,13,16,33. It is notable that the bare presence of the confinement affects E_MBD more strongly than the DCS interaction energy.

Fig. 4

MBD and dipole-correlated Coulomb singles (DCS) contributions to the Xe–Xe interaction inside carbon nanotubes (CNTs).

A Comparing the MBD and DCS contributions inside a (6, 6)-CNT and in gas-phase as a function of the Xe–Xe separation, R. B Effect of the different confinements of a (5, 5)- and (6, 6)-CNT on E_MBD and E_DCS. C Two Xe atoms (violet) encapsulated in a CNT. Total van der Waals interaction energy given as sum of E_MBD and E_DCS including the results when increasing the Xe polarizability by 50% (black). The inset shows the variation of the absolute value of the ratio of E_DCS and E_MBD.

Fig. 4B then shows the MBD and DCS components for the two CNTs with different radii. In contrast to the bare presence of confinement, the type and strength of the confinement, as represented by the different nanotubes, has a larger effect on the DCS interaction than on the MBD contribution. We also note that the effect of the different environment on the MBD interaction is negligible after 6 Å and the binding curves follow the same behavior, whereas the effect is more long-ranged for E_DCS. DCS are more sensitive to the characteristics of the confinement compared to E_MBD. As expected, the destabilization of the Xe dimer owing to screening and many-body correlation effects is less pronounced inside the (6, 6)-CNT.

Fig. 4C shows the cumulative vdW binding energy E_vdW = E_DCS + E_MBD. In total, confinement in a CNT leads to a substantial decrease of the Xe–Xe vdW interaction. Mostly owing to the DCS contribution, this destabilization is strongly dependent on the confining environment. This shows that the total long-range vdW interaction can be in fact substantially altered by (nano-)confinement, whereas the bare MBD treatment would predict the environment to have no effect beyond interatomic distances of 6 Å. For Xe₂ inside a (5, 5)-CNT, the interplay of the repulsive DCS contribution and the attractive MBD interaction interestingly leads to a near-linear behavior for separations of 6 Å to ~8 Å. To explore the balance between the repulsive E_DCS and the attractive E_MBD more clearly, we show the absolute value of their ratio as a function of R in the inset of Fig. 4C. In all cases, the ratio, i.e., the relative importance of the DCS contribution, increases with larger inter-Xe distances and reaches a maximum at a distance of ~5.8 Å before converging to zero. The interatomic distance at which we observe the maximum is surprisingly independent of the presence and strength of the confinement. The ratio of E_DCS and E_MBD and its maximum value, on the other side, strongly depends on the environment of the interacting particles.

In order to highlight the critical role played by the response properties of the objects interacting under confinement, we have performed the analysis for Xe₂ inside a (5, 5)-CNT while increasing the Xe polarizability by 50% (black curve). The results indicate a pivotal role of the polarizability in the total vdW interaction in confined systems: at shorter interatomic distances, the overall interaction is increased as the attractive MBD contribution is affected more strongly. In the very long-range limit, the interaction converges to the same behavior as for normal Xe₂. In the intermediate region, however, the repulsive contribution from DCS increases more strongly than its MBD counterpart, which leads to a substantial destabilization and eventually repulsive interaction energy. The interplay between E_DCS and E_MBD in this intermediate region gives rise to a maximum followed by a very shallow minimum creating a small barrier of ~0.1 kJ/mol in the binding curve. All this can be explained by a much higher sensitivity of E_DCS to the Xe polarizability compared to the MBD interaction energy. Accordingly, changing the polarizability has a strong effect on the ratio of E_DCS and E_MBD and its maximum (see Fig. 4C). The ratio surpasses 1, meaning that E_DCS supersedes the MBD contribution in magnitude and introduces a region of repulsive interaction at ~7 Å. The position of the maximum is thus increased by almost 1 Å compared to normal Xe₂. Altogether, we can conclude that with increasing polarizability, the relative importance of the DCS contribution increases, becomes more long-ranged, and can lead to non-trivial qualitative changes in the overall vdW interaction. For the considered system of Xe₂ inside CNTs, the vdW interaction thereby fully governs the total long-range interaction. The corresponding PBE-DFT interaction energies are of negligible magnitude at inter-Xe distances beyond ~6 Å and only introduce the well-known repulsive contributions at shorter separations. As a result, the qualitative changes owing to DCS reported above remain unaltered by inclusion of contributions captured in semi-local DFT. One particularly interesting aspect of the total PBE+MBD+DCS interaction is that PBE contributions together with DCS cancel out the meta-stable state of Xe₂ in the (5,5)-CNT. On the level of PBE+MBD, one observes a local minimum at a inter-Xe distance of ~5.5 Å, whereas PBE+MBD+DCS predicts only repulsive or negligible interaction at all distances. The corresponding PBE-DFT and total interaction energies are summarized in Supplementary Fig. 4.