Size-dependent polarizabilities and van der Waals dispersion coefficients of fullerenes from large-scale complex polarization propagator calculations

Ever since its most prominent member C₆₀ was first characterized in 1985,¹ the remarkable family of carbon fullerenes has gained increasing interest among both the experimental and theoretical communities.^2–4 Located in the realm between molecular and nanostructural materials, the properties of fullerenes vary with size, shape, and structure, making them appealing for a wide range of nanotechnological applications such as sensing, photovoltaics, and electronics. Dedicated to the functionalization of fullerenes, the emerged sub-field of fullerene chemistry enhances their potential and increases the scope of applications even further.^5–8 

The cage-like structure with all sp²-hybridized carbon atoms located at the surface leads to delocalized π-electron densities with significant curvature effects, in particular, for fullerenes with smaller diameter. Due to their closed structure, the attractive interaction between neutral fullerene pairs is dominated by dispersive long-range contributions.^9,10 At large intermolecular separation distances R where charge density overlaps are negligible but retardation effects are not yet coming into play, i.e., the van der Waals (vdW) region, the interaction energy exhibits an R⁻⁶ dependence, as expressed by the C₆ dispersion coefficient. As introduced by Casimir and Polder,¹¹ the interaction potential for two polarizable systems connects the C₆ coefficients to the associated electric-dipole polarizabilities, which are the subject of previous experimental and theoretical studies.^12–17 

An efficient first-principles computational scheme to obtain dispersion coefficients has been established using the linear complex polarization propagator (CPP) method,^18–20 where the molecular polarizabilities entering into the Casimir–Polder integral are directly evaluated on the imaginary frequency axis. Application of this procedure to a set of carbon fullerenes ranging from C₆₀ to C₈₄ revealed anomalous size-dependencies and scalings with respect to the number of carbon atoms N. For the static polarizabilities and C₆ coefficients, scalings were found to be equal to N^1.2 and N^2.2, respectively, at the level of density functional theory (DFT) in conjunction with the B3LYP functional.²¹ These results deviate from the prediction made from a classical model using a spherical-shell approximation, obtaining a scaling of N^2.75 for the van der Waals coefficients.²² Later, an extension of the scope of the investigated carbon fullerenes was made as to include up to a size of 720 atoms. For computational reasons, this study was forced to use an atomistic electrodynamics model based on the capacitance–polarizability interaction (CPI), and it resulted in scaling exponents of 1.46 and 2.8.¹⁷ Agreement with the classical model was in this case found to be close, but the accuracy of the CPI model was difficult to assess. The fact that quantum size effects play an important role in nanoscale materials is widely accepted and previously reported,^23–26 so it is arguably worthwhile to revisit the study of the property-to-size scalings of fullerenes but at a first-principles level of theory.

Theoretical and methodological developments and hardware advancements at the high-performance computing (HPC) centers enable new applications in theoretical chemistry. Recently, we have launched the VeloxChem project aimed at the development of a modern, object-oriented, Python/C++ software for DFT modeling of spectroscopic properties on contemporary and future HPC platforms.²⁷ Key components for this study are the highly efficient hybrid OpenMP/MPI parallel Fock-matrix construction and an efficient multi-frequency linear CPP equation solver.^28,29 In the present work, we have further improved the solver algorithm/implementation so that the storage of reduced-space vectors and the associated vector–vector operations are carried out across the available cluster nodes, enabling linear-response calculations to be performed on large-scale systems as several terabytes of distributed memory are easily accessible.

As a demonstration of the capability of this software, we extend the scope of systems to be treated with the CPP method by performing calculations of static polarizabilities and C₆ dispersion coefficients for carbon fullerenes up to C₅₄₀. We first validate our computational procedure by comparing our results for smaller systems to available experimental and theoretical benchmark values. After affirmation of the integrity of the applied approach, we establish reference data for the scalings of these properties with respect to system size.

The long-range isotropic average interaction energy between two separated systems A and B is given by the Casimir–Polder potential,¹¹ which, neglecting retardation effects, simplifies to

where R is the distance between the systems and C₆ is the dipole–dipole dispersion coefficient. The latter can be expressed by means of dipole polarizabilities $α¯$⁠,

More specifically, $α¯A(iωI)$ refers to the orientation-averaged dipole polarizability of system A at an imaginary frequency iω^I.

With the polarizability being the first-order response of the molecular dipole moment to an external electric field, its response function for an incident frequency ω may be written in terms of a sum-over-state (SOS) expression, which can be evaluated at any general complex-frequency argument ω = ω^R + iω^I by means of the polarization propagator approach³⁰ and its extension to the complex-frequency domain (CPP).¹⁹ 

While the spectral representation is valid for exact states, the application to approximate state methods leads to matrix equations instead. In the framework of single determinant Hartree–Fock (HF) or Kohn–Sham DFT approximations, the linear-response functions for purely imaginary frequency arguments take the following form:³¹ 

where E^[2] and S^[2] are the so-called Hessian and metric matrices, respectively, and A^[1] and B^[1] are the property gradients corresponding to the components of the linear polarizability. The computational effort to perform explicit inversion of the involved matrix expressions exceeds feasibility/practicality even for rather small systems due to their large dimension. The response functions are instead evaluated via an iterative subspace algorithm using symmetrized trial vectors,²⁸ in which solving the matrix equation

for the complex response vector X = X^R + iX^I is central. The given equation results as a reduction in a more general form involving a 4 × 4 matrix, making use of the facts that the dipole gradient is purely real and of ungerade symmetry by its nature and only imaginary frequencies are considered in the case of α(iω^I).

An efficient hybrid OpenMP/MPI implementation of this algorithm enabling parallel handling of multiple frequencies and gradients in a common reduced space including the evaluation of the entire set of two-electron integrals in every iteration has been reported recently.²⁷ Additional improvements of the implementation address the storage and distribution of the reduced-space vectors. Arrangement of those column vectors in matrices, which are then distributed row-wise, enables the necessary algebraic operations of the vectors to be carried out in a parallel manner over the available computational nodes while making use of the entirety of cluster-aggregated memory. Inter-nodal data communication was thereby observed to account for less than 0.1% of the total computation time spent in the solver routine.

The C₆ dispersion coefficients were obtained by approximating the integral in Eq. (2) using a Gauss–Legendre quadrature after substituting the integration variables according to

where a transformation factor of ω₀ = 0.3 a.u. was used.³² 

The CPP approach has been used to determine dispersion coefficients for noble gases and n-alkanes,¹⁸ polyacenes,³³ sodium clusters,³⁴ and fullerenes up to a system size of 84 carbon atoms.²¹ 

The geometries were taken from the fullerene structure library created by Yoshida³⁵ and optimized with the approximate normal coordinate rational function optimization (ANCOPT) scheme based on the semi-empirical density functional tight-binding (DFTB)³⁶ method. The geometry optimizations were carried out using the xtb program (version 6.2.3) employing the GFN2-xTB parameterization.³⁷ 

The molecular properties were calculated at the Hartree–Fock and DFT levels of theory. For the latter, the B3LYP³⁸ functional was used as well as the local density approximation (LDA).

The augmented def2-SVPD basis set³⁹ optimized with respect to atomic static polarizabilities was applied for all property calculations carried out in this work. The contraction scheme [8s4p2d|4s2p2d] is applied to carbon atoms in the def2-SVPD basis set, resulting in a total of 1200 to 10 800 contracted basis functions for the investigated fullerenes.

The threshold for the relative residual norm used to solve the response equations was 10⁻³. Compared to a threshold of 10⁻⁴, the static polarizability of C₆₀ changed insignificantly (4.7 × 10⁻⁵ a.u.). A seven-point Gauss–Legendre integration scheme was employed to evaluate the Casimir–Polder integral for the C₆ dispersion coefficients since it was found to yield the C₆ coefficient for C₆₀ with four-digit accuracy compared to a 12-point integration scheme.

All property calculations have been carried out with the VeloxChem program (modified version 1.0).²⁷ 

The foundation of the property calculations carried out in this work is formed by the molecular geometries. Therefore, their accuracy had to be assessed by means of a comparison with both experimental data and theoretical approaches used in other studies. Table I contains average bond lengths obtained from calculations and experiments. Although conventional DFT calculations are evidently in better agreement with the experimentally obtained values for C₆₀, the small deviations of the DFTB optimized geometries of around 0.01 Å are still acceptable. Notably, both DFT and DFTB approaches appear to describe the structure of C₇₀ equally well, differing for every type of bond (Fig. 1). Our need for a method applicable to large systems—supported by the reported deployment of the tight-binding DFT scheme to the calculation of geometries and ¹³C-NMR shielding constants of carbon fullerenes⁴⁰—prompted us to employ the DFTB method for all geometry optimizations in this study.

Although the Sadlej basis set⁴³ has been shown to accurately describe polarizabilities of carbon fullerenes,²¹ the smaller def2-SVPD basis set was used after the evaluation for the property calculations in order to reduce the computational cost. The static polarizabilities of C₆₀ for both basis sets deviate by 1.1% and 0.9% for the HF and B3LYP level of theory, respectively. For C₇₀, the values differ by 1.0% and 0.8% for the same methods. A comparison with the experimentally obtained dipole polarizabilities and theoretical benchmark values (Table II) confirms that the basis sets perform similarly. Existing experimental values obtained by deflectometry methods are accompanied by rather large statistical and systematic uncertainties,^12,44 which was recently reduced significantly by Fein et al.¹³ for both C₆₀ and C₇₀. Additionally, these measurements include vibrational contributions not accounted for in our calculations, hampering a direct comparison of the values. Accurate theoretically calculated values were obtained using dipole oscillator strength distribution (DOSD)⁴⁵ and the approximated coupled cluster singles and doubles model RICC2.⁴⁶ Another ab initio value for the polarizability was calculated with the linear-response formalism of the coupled cluster singles doubles (LR-CCSD) method,⁴⁷ all providing reliable benchmark values in the range of 555.3–559.6 a.u. These values are underestimated by the CPP approach applied herein by 6%–7% and 3%–4% for the def2-SVPD basis set at the HF and B3LYP levels of theory, which has to be compared with the slightly smaller deviation of 5%–6% and 2%–3% found for the Sadlej basis set. We consider the results for both basis sets to be in very good agreement with the theoretical benchmark.

As the next larger stable carbon fullerene, C₇₀ and its properties were experimentally addressed extensively. In line with the values for C₆₀, the results found for the static polarizability in various deflectometry studies are afflicted with large uncertainties.^44,48 Although a higher accuracy was achieved in more recent experiments,¹³ a direct comparison with our calculated results is left out for the reasons mentioned above. A calculated benchmark value using methods at a sufficiently high level of theory is not known until this point. While the lack of a highly accurate value does not allow us to compare our results for C₇₀ directly, an indirect assessment can be done considering the ratio of the polarizabilities of C₇₀ and C₆₀, which has been measured experimentally with the benefit of having much smaller uncertainties than the individual polarizabilities themselves. The measured ratios of 1.22 ± 0.03⁴⁴ and 1.33 ± 0.03⁴⁸ can be found in the literature, while a value of 1.22 was extracted from Ref. 13 with the uncertainty expected to be in the same magnitude as the ones given explicitly. The ratios found for the CPP approach are 1.23 and 1.24 for HF and DFT/B3LYP. In the given precision, the ratios are identical for the def2-SVPD and Sadlej basis sets. With the ratios of $α¯(0)$ being in excellent agreement with the experimental ones, we conclude that the quality of the values for the polarizabilities of C₇₀ is equally good as for C₆₀. We further assume that the capability of our computational protocol to describe $α¯(0)$ of those two systems is comparably high for larger carbon fullerenes. Whereas the result obtained with the Sadlej basis set is in marginally better agreement with the C₆₀ benchmark value, the slightly larger underestimation of the values for $α¯(0)$ is a sacrifice we are willing to make in order to reduce the computational cost to make much larger systems up to 540 carbon atoms accessible with the CPP approach.

The ratio of the static polarizabilities of C₇₀ and C₆₀ already suggests that the scaling of this property is not behaving in a linear fashion with respect to the number of carbon atoms. With the assumption of an N^η scaling, different values for η were found previously for different computational treatments of the systems. Limited by the computational costs, a completely ab initio approach has not exceeded the size of C₈₄ until this point, leaving a question mark behind the scaling of not only the polarizability but other size-dependent properties as well.

Therefore, the established computational protocol was used to calculate the static polarizabilities and C₆ dispersion coefficients of the carbon fullerenes C₁₀₀, C₁₈₀, C₂₄₀, and C₅₄₀. The selection of the systems was mainly motivated by an interest in the group of fullerenes with N = 60n atoms, with n being a positive integer. This group of fullerenes is expected to have a higher energy gap between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO)⁴⁹ associated with a higher kinetic stability.⁵⁰ However, C₁₀₀ was chosen as an intermediate step between C₇₀ and C₁₈₀, as no geometry for C₁₂₀ was contained in the fullerene structure library.³⁵ 

Table III shows the obtained static polarizabilities for C₆₀–C₅₄₀ for the different levels of theory along with their values per carbon atom. The latter directly reveal the non-additive scaling, indicating an increase in the per-atom-polarizability with growing system size.

The homomolecular C₆ coefficients for the investigated systems are summarized in Table IV. Since the polarizability is considered quadratically in the calculation of the coefficients, the values are also divided by the square of the number of atoms N to yield the C₆ coefficients per atom. Those can also be seen as the contribution per single carbon–carbon interaction. Evidently, the dispersion also behaves in a non-additive manner, as reported before for fullerenes^17,21,22 and other carbon-based nanomaterials.⁵¹ Consequently, all calculated values except for C₆₀ are higher than the recommended DFT-D3 values 27.4 a.u. and 22.4 a.u. for sp²-and sp³-hybridized carbon–carbon dispersions, respectively.⁵² 

Presuming an N^η dependence, the values for $α¯(0)$ were plotted over the number of atoms N applying an logarithmic scale on both axes (Fig. 2). The scaling exponents η can be read directly from this plot as the slopes of the linear regressions. They were obtained as η = 1.184 and 1.222 at the HF and B3LYP levels of theory, respectively. Compared to the values reported earlier by Kauczor et al.²¹ of η = 1.230 and 1.246 for the same methods, but smaller scope of fullerenes, the scalings found herein are slightly lower, whereas the trend in between methods is preserved. The decrease in the scaling with an expansion of the system size might seem surprising in the light of the significantly higher value of η = 1.46 found in Ref. 17. Although fullerenes up to a system size of 720 atoms were investigated therein and the semi-empirical parameters were fit to ab initio results, the accuracy of the used CPI model can by no means compete with a fully ab initio approach, as its application has led to—what now has to be considered as—a less precise value for the static polarizability of C₆₀ (502.7 a.u.).⁵³ We, therefore, believe that the found scaling of N^1.2 provides the most accurate representation of today of the size-dependency for $α¯(0)$ of carbon fullerenes. Calculations for the intermediate-sized fullerenes C₇₂, C₇₆, C₈₀, C₈₄, and C₉₀ confirm this trend (see S1 of the supplementary material), providing a direct link between the different size domains. In an attempt to provide a closer description of the electronic structure to that of a metallic model system with a homogeneous electron gas, the LDA exchange functional was used for the calculation of $α¯(0)$ for C₆₀ to C₁₈₀. The obtained scaling exponent of η = 1.190 showed that our results are stable with respect to functional variations. Calculations for C₂₄₀ and C₅₄₀ were omitted as no further physical insights can be expected.

Following the same C₆ ∝ N^η ansatz, the logarithmic plot of the dispersion coefficients over the system size (Fig. 3) gives rise to values of η = 2.190 and 2.218 for HF and B3LYP. In contrast to the scaling for the polarizability, the corresponding exponents found in Ref. 21 of 2.176 and 2.188 are marginally lower than our results, although not significantly. As for the polarizability, calculations for the smaller systems employing the LDA exchange functional result in a scaling exponent of η = 2.185.

While the exponents obtained with different methods vary only within a range of 0.03 and can be condensed to η = 2.2, the divergence from studies investigating fullerenes up to 720 and 3840 carbon atoms predicting values of 2.8¹⁷ and 2.75²² is more striking. However, as discussed above, our approach used herein provides a higher expected accuracy in comparison with Ref. 22, treating the fullerenes simply as classical spherical shells. Using a quantum harmonic oscillator (QHO) approach,⁵⁴ Gobre and Tkatchenko⁵¹ reported a N^2.35 dependence by looking at C₂₀, C₆₀, C₈₀, and C₅₄₀ and reducing to N^2.25 for only considering the first three systems. Although the latter dependence matches quite well with our findings, the increase in the exponent deviates significantly. This might be caused by the sparse selection of systems or simply by the poor performance of the chosen method involving classical electrodynamics not benchmarked for fullerenes specifically. However, all predictions agree in the non-additive behavior for the van der Waals dispersion coefficients, while their size-scaling was significantly overestimated using (semi-)classical methods. Our ab initio treatment confirms the trend of N^2.2 reported previously for a smaller scope of carbon fullerenes. This confirmation is backed up by the calculation of fullerenes in the range of C₇₂ to C₉₀, supporting the found scaling (see S2 of the supplementary material). Agreement of these values with results reported previously reveals once more that our computational protocol for obtaining the structures yields reasonable results.

Consideration of isomeric effects on the properties was taken into account by optimizing all 450 isomers included in the library for C₁₀₀ and calculating $α¯(0)$ and C₆ coefficients for the five isomers with the lowest energy. The reported values for this system (Table III or Table IV) are averaged over these five isomers. The difference between the average value and the lowest/highest values of the isomers was introduced as an “isomeric spread” to determine a systematic uncertainty of the property scaling. However, the uncertainty obtained with this approach was found to be insignificant within the given precision.

With an improved and efficient distributed memory handling in the complex linear-response equation solver, stable calculations of dynamic electric-dipole polarizabilities for systems involving up to and beyond 10 000 contracted and partly diffuse basis functions are made available.

The complex linear polarization propagator method was used to compute static and dynamic polarizabilities for fullerenes C_N with N ≤ 540. Employing the Hartree–Fock and DFT/B3LYP levels of theory, the obtained static polarizabilities of C₆₀ underestimate accurate benchmark values by 6%–7% and 3%–4%, respectively, whereas the obtained ratio of the static polarizabilities C₇₀/C₆₀ is in excellent agreement with the experimentally measured one. Given the available experimental data together with their large uncertainties, there is no evidence other than that the introduced computational protocol is reliable for the set of carbon fullerenes.

The scalings with the number of carbon atoms of the static polarizabilities and the C₆ van der Waals dispersion coefficients, calculated from the dynamic polarizabilities, were found to be N^1.2 and N^2.2, respectively. These results confirm the findings of a previous study using the same approach but limited to smaller system sizes.²¹ Predictions of the scaling exponent for C₆ coefficients made by other studies are 2.75,²² 2.8,¹⁷ and 2.35.⁵¹ Although the scope of those investigations was sufficiently large for a comparison, all the methods applied therein were of classical, semi-classical, or semi-empirical nature.

We are, therefore, confident that our fully analytical approach—made applicable to those system sizes by an efficient implementation of the underlying algorithm combined with the access to high-performance computing resources—provides a more accurate description of the size-dependency of polarizabilities and C₆ dispersion coefficients.

See the supplementary material for static polarizabilities and C₆ dispersion coefficients of system sizes C₇₂–C₁₀₀.

Financial support from the European Commission in the form of the ITN titled “Computational Spectroscopy in Natural Sciences and Engineering” (COSINE) (Grant No. 765739), the Swedish e-Science Research Centre (SeRC), and the Swedish Research Council (Grant No. 2018-4343) is acknowledged. Computational resources are provided by the Swedish National Infrastructure for Computing (SNIC).

The data that support the findings of this study are available from the corresponding author upon reasonable request.