Generalized Energy-Based Fragmentation DLPNO-CCSD(T) Approach at Complete Basis Set Limit and Its Application to Benzene Clusters

I.   INTRODUCTION

The π-π interactions play a crucial role in biomolecules and molecular crystals [1, 2]. Understanding these interactions is essential for elucidating the properties and behaviors of complex systems [3−6], including benzene cluster systems [7] as well. Previous studies have calculated the energies of benzene dimers [8−11]. However, the computational cost increases significantly as the number of benzene molecules increases. To address this issue, extensive research on benzene clusters has been conducted [12−19]. Chakrabarti and coworkers [15] used symmetry-adapted perturbation theory (SAPT) based on density functional theory (DFT) to study the energy landscape of a benzene cluster, (C₆H₆)₁₃. Takeuchi [16] investigated benzene clusters (C₆H₆)_n (n ≤ 30) using OPLS-AA potential, revealing that the local structures deviate from regular icosahedrons and exhibit anisotropic interactions, and the clusters with n = 13, 19, 23, 26, and 29 (magic numbers) are relatively stable. Gonzalez et al. [18] applied the second-order Møller-Plesset perturbation theory (MP2) method to calculate the energies of benzene clusters (n=3, 4, 5). Furthermore, Ahirwar and coworkers [19] employed the molecular tailoring approach (MTA), a linear scaling method, to calculate the energies of benzene clusters (n=3–7) at the MP2.5/aug-cc-pVDZ level. However, high-accuracy quantum chemistry calculations for complex benzene clusters remain relatively rare‎ in current research‎ [20, 21]‎.

Various methods have been developed to study the average binding energies (ABEs) of molecule clusters, with the density functional theory (DFT) method being particularly notable [22, 23]. However, DFT struggles to accurately describe dispersion interactions. To address this, the empirical D3 dispersion correction [24, 25] proposed by Grimme’s group enhances the accuracy of DFT calculations by mitigating these limitations. Despite this improvement, DFT still faces challenges when applied to clusters containing a large number of atoms. Moreover, ab initio electron correlation methods, such as MP2, can provide more accurate calculations of the ABEs of benzene clusters, but they exhibit computational complexity that scales as N⁵ (N is the number of basis functions). Higher-level electron correlation methods, such as coupled cluster singles and doubles with noniterative triple excitations ‎‎(CCSD(T)) [26−28], particularly at the complete basis set (CBS) limit (known as CCSD(T)/CBS) [29] scales as N⁷. Therefore, the high computational cost makes it difficult to directly use the CCSD(T)/CBS method to compute large molecular clusters.

To uphold computational accuracy while managing memory and disk resources, researchers have explored methods that ensure resource requirements scale linearly with system size. Among these linear scaling methods, the energy-based fragmentation (EBF) methods [30−41] have gained broader acceptance. EBF methods represent the energy (or gradient) of a large system as a combination of energies (or gradients) of small subsystems, offering a concise and widely applicable approach. Notably, the generalized energy-based fragmentation (GEBF) approach [33, 34, 37, 42] has been successfully applied to perform high-level ab initio calculations on complex systems. For instance, Li and coworkers used the GEBF method to examine the ABEs of large water clusters, specifically (H₂O)_n (n=32 and 64), at the CCSD(T)/CBS levels [43]. The domain-based local pair natural orbital‎ (DLPNO) method is a “direct” local correlation approach developed by Neese and coworkers [44, 45]. This approach simultaneously solves a set of equations for the amplitudes of strongly interacting pairs. The DLPNO approach has been combined with the highly accurate CCSD(T) method to efficiently calculate the energies of large molecules [46−48]. Li and coworkers [49−51] proposed a cluster-based local correlation method, the cluster-in-molecule (CIM) approach, which was combined with the DLPNO-CCSD(T) method to obtain CBS energies for benzene clusters of various sizes [49]. However, there are no reports of using the GEBF method to calculate the ABEs of benzene clusters at the CCSD(T)/CBS level of accuracy.

In this work, we have combined the GEBF approach with the DLPNO-CCSD(T) method to accurately compute the ABEs and relative stability of benzene clusters using the aug-cc-pVXZ (X=D,T) basis set and extrapolated to the CBS limit. These energies served as reference values and were compared with those obtained using GEBF-DFT methods with various functionals. The work also identifies the most accurate DFT functionals for calculating ABEs and discusses the impact of empirical D3 correction on these functionals, as well as changes in binding energies and stabilities across different sizes of benzene clusters.

II.   METHODOLOGY AND COMPUTATIONAL DETAILS

A   GEBF approach

In the GEBF approach, we calculate the total ground-state energy of a large system by combining the ground-state energies of a series of subsystems, each of which is electrostatically embedded (EE) in background point charges outside the atoms of subsystem. Then, the total energy using the GEBF method is expressed as follows [37]:

Here, M represents the total number of subsystems. ${\stackrel{~}{E}}_{i}$ is the ground-state energy of the ith EE subsystem, with its coefficient denoted by ${C}_{i}$. ${Q}_{\rm{A}}$ represents the natural charge of atom A, while ${R}_{\rm{AB}}$ denotes the interatomic distance between atoms A and B.

The basic procedure of GEBF approach is summarized as follows using FIG. 1 as an illustration: (1) Divide a target system into various fragments (see FIG. 1(a)). (2) For each fragment (denoted as central fragment), construct a primitive subsystem by adding its spatially close neighboring fragments within a distance threshold (ξ), adding hydrogen atoms for valence saturation when needed (see FIG. 1(b)), although this is not necessary for benzene clusters. The size of the subsystem is limited by a maximum number of fragments (η). (3) Additionally, generate extra primitive subsystems for pairs of fragments separated by less than 2ξ that are not included together in any other primitive subsystem (see FIG. 1(c)). (4) Obtain derivative subsystems (see FIG. 1(d)) using inclusion-exclusion principle to address overlaps among the primitive subsystems. (5) Embed each subsystem in background point charges generated by all atoms outside the subsystem. These charges are iteratively derived from natural population analysis (NPA) [52, 53] at HF level (for HF or post-HF calculations) or DFT level (for DFT calculations) on all primitive subsystems and electrostatically embedded (EE) primitive subsystems. (6) Finally, perform conventional ab initio calculations on all EE subsystems, and obtain the total energy of the target system according to Eq.(1).

Figure  1.  The illustration of constructing GEBF subsystems: (a) divide target system into seven fragments; (b) add environmental fragments to central fragment to form primitive subsystems; (c) construct extra subsystems for pairs of fragments with distance being less than 2ξ but not simultaneously appearing in other primitive subsystems; (d) obtain derivative subsystems by inclusion-exclusion principle.

DownLoad: Full-Size Img PowerPoint

B   Basis set extrapolation

The CCSD(T) energies calculated with the aug-cc-pVXZ (X=D, T) basis sets are extrapolated to the CBS limit by using the following equations [29]:

for Hartree-Fock (HF) energy

and for CCSD(T) correlation energy.

When using the aug-cc-pVDZ (X = 2) and aug-cc-pVTZ (Y = 3) basis sets [54], the recommended values for α and β are 4.42 and 2.46, respectively. Then, the total CCSD(T) energy at the CBS limit is obtained by adding the HF energy in Eq.(2) and CCSD(T) correlation energy in Eq.(3). In our approach, we compute the GEBF-DLPNO-CCSD(T) correlation energies with different basis set and then extrapolate the GEBF-DLPNO-CCSD(T)/CBS correlation energy. An alternative way is to combine the DLPNO-CCSD(T)/CBS correlation energies of all subsystems to obtain the total correlation energies at the CBS limit. Since the extrapolation is linear with respective to ${E}_{\text{corr}}^{\left(X\right)}$ and ${E}_{\text{corr}}^{\left(Y\right)}$, the two extrapolation ways are equivalent.

C   Computational details

In this work, the ABEs of ten benzene clusters (C₆H₆)₁₃ (displayed in FIG. 2) are computed at the GEBF-DLPNO-CCSD(T) and GEBF-DFT levels. The GEBF-DLPNO-CCSD(T) energies are obtained with the aug-cc-pVXZ (X=D, T) basis sets and extrapolated to CBS limit using the basis set extrapolation. For DFT calculations, the selected functionals, including the PBE [55], PBE-D3(BJ) [25], B3LYP [56], B3LYP-D3(BJ), M06-2X [57], M06-2X-D3 [24], MN15 [58], and ωB97M-V [59], all utilize the aug-cc-pVTZ basis set. Then the ABEs (per benzene molecule) and relative stabilities are obtained for (C₆H₆)_n (n=11, 12, 13, 14, 15) cluster with the aug-cc-pVTZ basis set. All structures of the benzene clusters are taken out from the work of Takeuchi [16]. The coordinates of these structures are available in Supplementary materials (SM). The GEBF calculations are performed using the LSQC 2.5 program [60, 61], and each benzene molecule is considered as a fragment. For GEBF subsystems, the DLPNO-CCSD(T) and DFT calculations are carried using ORCA 5.0.4 program [62−67] and Gaussian 16 package [68], respectively. In the DLPNO calculations, the “TightPNO” settings are utilized [69]. For CCSD, the PNO is from non-iterative semi-local MP2. For (T), iterative DLPNO-(T) is employed [70]. The DFTD3 (version 3.1) program [24, 25] is used to compute the dispersion correction for some functionals.

Figure  2.  The structures of the 10 isomers of benzene clusters (C₆H₆)₁₃ under study.

DownLoad: Full-Size Img PowerPoint

III.   RESULTS AND DISCUSSION

A   Benchmark calculations for the accuracy of GEBF and basis set effects

In GEBF calculations, selecting the values of two parameters, ξ and η, can balance the computational costs and accuracy. Based on the experience of previous GEBF calculations, we set the value of ξ as 4 Å and tested different η values of 4, 5, 6, and 7 for benzene clusters. The GEBF(4, η)-PBE (η=4, 5, 6, 7) ABEs (per benzene molecule) $\Delta {E_n} = ({E_n} - n{E_1})/n$ for 10 isomers of benzene cluster (C₆H₆)₁₃ with the aug-cc-pVTZ basis set are compared in Table I. It shows that with increasing η, the GEBF ABEs converge with the largest deviation being only 0.02 kcal/mol between η=6 and 7. It indicates that our GEBF(4, 6) computations could achieve computational accuracy of ABEs below kilojoule per mole. Therefore, we choose η = 6 in all the subsequent GEBF-DLPNO-CCSD(T) and GEBF-DFT calculations to balance the computational expenses and precision of the results.

Table  I.  The ABEs per benzene molecule for ten (C₆H₆)₁₃ isomers calculated at the GEBF(4, η)-PBE (η=4, 5, 6, 7) level using the aug-cc-pVTZ basis set.

Isomer	ABE/(kcal/mol)
	η=4	η=5	η=6	η=7
1	−0.31	−0.12	−0.07	−0.06
2	−0.35	−0.09	−0.03	−0.02
3	−0.25	0.01	0.09	0.09
4	−0.29	−0.09	0.05	0.05
5	−0.17	0.00	0.09	0.09
6	−0.24	0.03	0.12	0.12
7	−0.27	−0.01	0.07	0.09
8	−0.25	−0.06	0.03	0.03
9	−0.30	0.00	0.12	0.12
10	−0.13	0.08	0.23	0.24

 | Show Table

DownLoad: CSV

In this study, GEBF-HF energies and GEBF-DLPNO-CCSD(T) correlation energies were obtained using the aug-cc-pVXZ (X=D, T) basis sets and extrapolated to the CBS limit. Subsequently, the GEBF-DLPNO-CCSD(T) ABEs of each (C₆H₆)₁₃ cluster were compared across the aug-cc-pVDZ, aug-cc-pVTZ, and complete basis sets, as shown in Table II. This comparison demonstrates that both the aug-cc-pVDZ and aug-cc-pVTZ basis sets correctly order the ten isomers average binding to the CBS results. The maximum unsigned error (MUE) and the average unsigned error (AUE) for the aug-cc-pVDZ results are 5.72 and 5.58 kcal/mol, respectively, while those for the aug-cc-pVTZ results are 1.80 and 1.76 kcal/mol, respectively. It indicates that the DLPNO-CCSD(T) ABEs using the aug-cc-pVTZ basis set are still not accurate enough compared to CBS results. In addition, the ABEs (kcal/mol) and the extrapolated ABEs calculated using the M06-2X-D3 method among 10 isomers of (C₆H₆)₁₃ are compared in Table S4 (SM). It shows that the MUE and AUE for the M06-2X-D3/aug-cc-pVTZ results are only 0.69 and 0.67 kcal/mol, respectively, when compared with the CBS results. Therefore, the DFT results with the aug-cc-pVTZ basis set could be employed to compute the ABEs for the benzene clusters.

Table  II.  The ABEs and the extrapolated ABEs calculated using the GEBF-DLPNO-CCSD(T) method based on the aug-cc-pVXZ (X=D, T) basis set among 10 isomers of (C₆H₆)₁₃.

Isomer	ABE/(kcal/mol)
	aug-cc-pVDZ	aug-cc-pVTZ	CBS
1	−12.48	−8.69	−6.96
2	−12.36	−8.54	−6.77
3	−12.26	−8.45	−6.70
4	−12.33	−8.42	−6.61
5	−12.19	−8.38	−6.63
6	−12.18	−8.34	−6.57
7	−12.16	−8.40	−6.67
8	−12.17	−8.37	−6.63
9	−12.08	−8.23	−6.46
10	−12.17	−8.34	−6.57
MUEa	5.72	1.80	0.00
AUEb	5.58	1.76	0.00
a The maximum unsigned error relative to the GEBF-DLPNO-CCSD(T)/CBS results. b The average unsigned error relative to the GEBF-DLPNO-CCSD(T)/CBS results.

 | Show Table

DownLoad: CSV

B   The average binding energies of benzene clusters (C₆H₆)₁₃

Using the GEBF-DLPNO-CCSD(T)/CBS ABEs as reference values, the ABEs calculated from various DFT functionals with the aug-cc-pVTZ basis set are compared in FIG. 3 (for selected functionals only) and Table III. These functionals include PBE, PBE-D3(BJ), B3LYP, B3LYP-D3(BJ), M06-2X, M06-2X-D3, MN15, and ωB97M-V. According to the GEBF-DLPNO-CCSD(T)/CBS results, the ABEs gradually increase from isomer 1 to 10 within the benzene clusters, indicating that isomer 1 is the most stable. Both FIG. 3 and Table III demonstrate that incorporating the D3 dispersion correction generally improves the accuracy of the DFT calculations. For example, the MUEs for the PBE, B3LYP and M06-2X functionals without D3 corrections are 6.89, 9.08 and 1.62 kcal/mol, respectively. In contrast, the results for these functionals with D3 corrections are 0.51, 0.24 and 0.21 kcal/mol, respectively. The AUE values of these DFT calculations support the same conclusion. The PBE-D3(BJ), B3LYP-D3(BJ), and M06-2X-D3 functionals perform very well, with MUEs of ABEs no greater than 0.51 kcal/mol and AUEs of only 0.40 kcal/mol. Notably, the MN15 functional without empirical D3 corrections already yields satisfactory results. Therefore, the PBE-D3(BJ), B3LYP-D3(BJ), and M06-2X-D3 functionals are suitable for treating the noncovalent interactions in benzene clusters.

Figure  3.  Comparison of ABEs calculated by GEBF-DLPNO-CCSD(T)/CBS and several DFT functionals for 10 isomers of (C₆H₆)₁₃ cluster. The aug-cc-pVTZ basis set is employed for DFT calculations.

DownLoad: Full-Size Img PowerPoint

Table  III.  A comparison of ABEs, calculated using GEBF-DLPNO-CCSD(T)/CBS and various DFT functionals with the aug-cc-pVTZ basis set for 10 isomers of the (C₆H₆)₁₃ cluster.

Isomer	ABE/(kcal/mol)
	GEBF-DLPNO-CCSD(T)	PBE (-D3(BJ))	B3LYP (-D3(BJ))	M06-2X (-D3)	MN15	ωB97M-V
1	−6.96	−0.07 (−6.44)	2.13 (−7.04)	−5.35 (−6.79)	−7.65	−5.67
2	−6.77	−0.03 (−6.35)	2.17 (−6.92)	−5.23 (−6.65)	−7.53	−5.65
3	−6.70	0.09 (−6.25)	2.29 (−6.83)	−5.08 (−6.50)	−7.38	−5.55
4	−6.67	0.05 (−6.25)	2.23 (−6.83)	−5.13 (−6.56)	−7.46	−5.55
5	−6.63	0.09 (−6.25)	2.30 (−6.82)	−5.10 (−6.51)	−7.41	−5.45
6	−6.63	0.12 (−6.22)	2.31 (−6.80)	−5.02 (−6.46)	−7.36	−5.48
7	−6.61	0.07 (−6.28)	2.30 (−6.85)	−5.11 (−6.52)	−7.44	−5.52
8	−6.57	0.03 (−6.16)	2.19 (−6.72)	−5.08 (−6.46)	−7.27	−5.53
9	−6.57	0.12 (−6.22)	2.32 (−6.80)	−5.04 (−6.47)	−7.33	−5.43
10	−6.46	0.23 (−6.09)	2.44 (−6.65)	−4.89 (−6.31)	−7.21	−5.26
MUE	0.00	6.89 (0.51)	9.08 (0.24)	1.62 (0.21)	0.82	1.28
AUE	0.00	6.73 (0.40)	8.93 (0.17)	1.55 (0.13)	0.75	1.15
Note: The values in parentheses represent the ABEs obtained by adding the D3 or D3(BJ) dispersion correction to the functional.

 | Show Table

DownLoad: CSV

Although the DFT results exhibited some deviations from the reference values, a notable observation was the varying effects of the D3 correction on each functional’s performance. To explore the impact of the D3 correction on the accuracy of various DFT functionals, the AUEs of three functionals, both with and without the empirical D3 dispersion correction, for 10 isomers of (C₆H₆)₁₃ clusters were compared, as shown in FIG. 4. The data revealed that the addition of the D3 correction generally reduced the error in ABEs calculated by each functional. Specifically, the AUEs decreased by 6.33 and 8.76 kcal/mol for the PBE and B3LYP functionals, respectively. However, the AUE for the M06-2X functional decreased only by 1.42 kcal/mol, indicating that this functional already accounts for most dispersion interactions. This comparative study elucidates the effects of D3 correction on the performance of common DFT functionals, providing valuable insights for future computational research.

Figure  4.  The average unsigned errors (AUEs) of three DFT functionals with and without empirical D3 correction for 10 isomers of (C₆H₆)₁₃ clusters with the aug-cc-pVTZ basis set.

DownLoad: Full-Size Img PowerPoint

C   Binding energies and relative stability of benzene cluster (C₆H₆)_n

Based on prior calculations, we found that the M06-2X functional achieves high accuracy in computing benzene clusters without significant improvements from the D3 correction. Consequently, we employed the GEBF-M06-2X method with the aug-cc-pVTZ basis set to calculate the binding energies and relative stabilities of five benzene clusters, (C₆H₆)_n for n=11, 12, 13, 14, and 15, derived from stable structures reported by Takeuchi [16]. The specific structures of each cluster are depicted in FIG. 5, with detailed coordinates provided in SM. Here, the relative stability index (${S_n}$) is defined as ${E_{n - 1}} + {E_{n + 1}} - 2{E_n}$ [16], where ${E_n}$ is the energy of (C₆H₆)_n. Higher relative stability indicates a more stable benzene cluster (C₆H₆)_n. The ABEs of (C₆H₆)_n with n=10, 11, 12, 13, 14, 15 and 16 are compared in FIG. 6. It shows that the ABEs for (C₆H₆)₁₁ and (C₆H₆)₁₂ are −6.07 and −6.20 kcal/mol, respectively, with the changes being only 0.13 kcal/mol. However, the ABEs for (C₆H₆)₁₃ rapidly decreases to −6.79 kcal/mol. The ABEs of (C₆H₆)_n (n=14, 15, 16) gradually increase by only 0.01, −0.07 and −0.03 kcal/mol by comparing with those (C₆H₆)_n–1 cluster, respectively. Meanwhile, the relative stabilities of (C₆H₆)_n with n=11, 12, 13, 14 and 15 are compared in FIG. 7. A positive value of relative stability indicates that the cluster is more stable than its neighbors, while a negative value indicates that it is less stable. Therefore, (C₆H₆)₁₃ is the most stable among these five benzene clusters with its relative stability being 7.05 kcal/mol. These findings align with prior calculations and experimental observations, where each molecule in liquid benzene typically interacts with approximately 12 other benzene molecules, suggesting that n=13 is a “magic number” for benzene clusters (C₆H₆)_n.

Figure  5.  The stable structures of the (C₆H₆)_n (n=10, 11, 12, 13, 14, 15, 16).

DownLoad: Full-Size Img PowerPoint

Figure  6.  ABEs of benzene clusters (C₆H₆)_n with different n.

DownLoad: Full-Size Img PowerPoint

Figure  7.  Relative stability index of benzene clusters (C₆H₆)_n with different n.

DownLoad: Full-Size Img PowerPoint

IV.   CONCLUSION

In this study, we employed the GEBF method to investigate the average binding energies of (C₆H₆)₁₃ clusters at the DLPNO-CCSD(T)/CBS level. The results demonstrate that the GEBF method can accurately calculate the ABEs of large benzene clusters, showcasing its potential as a powerful tool for studying complex molecular systems. Furthermore, we extrapolated the CBS limit by calculating the energies of GEBF-DLPNO-CCSD(T)/aug-cc-pVXZ (X=D, T), providing a robust benchmark for evaluating the performance of various DFT methods. By using the calculated GEBF-DLPNO-CCSD(T)/CBS ABEs as reference data, we assessed the accuracy of a series of DFT methods in computing the ABEs of benzene clusters. Our results show that the PBE-D3(BJ), B3LYP-D3(BJ), and M06-2X-D3 functionals provide accurate ABEs of the benzene clusters. Moreover, adding empirical D3 corrections to some DFT functionals (PBE and B3LYP) can improve the calculation accuracy, providing simple suggestions for future computational work. Finally, we verify that (C₆H₆)₁₃ is the most stable cluster through calculations of ABEs and relative stabilities of benzene clusters (C₆H₆)_n (n=10–16), which is consistent with experimental results.

Supplementary materials: Details of total energies using the GEBF method with selected parameters, HF energies and correlation energies (including extrapolated results), total energies calculated using different DFT functionals, total energies, and the Cartesian coordinates of (C₆H₆)_n (n=1, 10–16) (DOCX) are shown.