To study the contribution of three-body dispersion to crystal lattice energies, we compute the three-body contributions to the lattice energies for crystalline benzene, carbon dioxide, and triazine using various computational methods. We show that these contributions converge quickly as the intermolecular distances between the monomers grow. In particular, the smallest value among the three pairwise intermonomer closest-contact distances, Rmin, shows a strong correlation with the three-body contribution to the lattice energy, and, here, the largest of the closest-contact distances, Rmax, serves as a cutoff criterion to limit the number of trimers to be considered. We considered all trimers up to Rmax=15Å. The trimers with Rmin<4Å contribute 90.4%, 90.6%, and 93.9% of the total three-body contributions for crystalline benzene, carbon dioxide, and triazine, respectively, for the coupled-cluster singles, doubles, and perturbative triples [CCSD(T)] method. For trimers with Rmin>4Å, the second-order Møller–Plesset perturbation theory (MP2) supplemented with the Axilrod–Teller–Muto (ATM) three-body dispersion correction reproduces the CCSD(T) values for the cumulative three-body contributions with errors of less than 0.1 kJ mol−1. Moreover, three-body contributions are converged within 0.15 kJ mol−1 by Rmax=10Å. From these results, it appears that in molecular crystals where dispersion dominates the three-body contribution to the lattice energy, the trimers with Rmin>4Å can be computed with the MP2+ATM method to reduce the computational cost, and those with Rmax>10Å appear to be basically negligible.

The many-body expansion (MBE) is a widely used computational method for computing the interaction energies of non-covalent molecular aggregates, such as molecular clusters and molecular crystals.1,2 The MBE expresses the interaction energy as the sum of n-body (n = 2, 3, 4, …) terms, each of which involves interactions between n molecules simultaneously. Because the MBE series converges quickly as the number of molecules involved n increases, in approximate computations, it may be truncated at a low order. The simplest MBE approach would be to truncate the series at n = 2, i.e., approximate the total interaction energy of the many-body system with the sum of all pairwise interactions. Such an approximation benefits from its low computational cost, which stems from the relatively modest number of possible dimer combinations, N(N1)2, as opposed to 2N − 1 possible monomers, dimers, trimers, etc., for the full many-body expansion, where N is the total number of molecules in the system. However, the many-body (n ≥ 3) effect can contribute ∼10%–20% to the total interaction energy, e.g., in crystalline benzene3,4 or water clusters,5 and is thus not negligible. Such examples point to the importance of an effective method to compute many-body interaction energies, especially the three-body ones, which usually make up the major contribution of many-body interaction energies.

In 2015, Beran et al. constructed the 3B-69 benchmark set, which consists of the three-body interaction energies of 69 trimers with their geometries extracted from 23 different molecular crystals,6 filling the vacancy of a benchmark dataset for three-body interactions. These workers demonstrated the crucial fact that DFT-D3 does not perform as well for the nonadditive three-body contribution to interaction energies as it does for two-body contributions, mainly because of the delocalization error that prevents density functional theory (DFT) from predicting the three-body polarization energies correctly. This work has also shown that the popular exchange functionals seem to have difficulty describing three-body exchange effects. For the wavefunction methods, they recommended MP2.5 (average of MP2 and MP3 energies)7 and spin-component scaled coupled-cluster singles doubles (CCSD)8 for noncovalent interactions [SCS(MI)-CCSD]9 for computing three-body interaction energies, but both of them scale as O(N6), restricting their usefulness for important practical applications like computing the lattice energies of molecular crystals. They have also used an ab initio force field4,10 to estimate the polarization and dispersion contribution to the three-body interaction energy, with the latter based on the Axilrod–Teller–Muto (ATM) formula.11,12 From these results, they reconfirmed that MP2 is a good approximation to three-body polarization effects, but it lacks three-body dispersion.

In addition to the wavefunction and DFT supermolecular approaches, an alternative approach to model three-body interactions was introduced by Podeszwa and Szalewicz in 2007.13 In this approach, the symmetry-adapted perturbation theory based on the density functional theory [SAPT(DFT)]14,15 was generalized to three-body interactions. In addition to demonstrating how to express several of the components of the three-body interaction in the SAPT(DFT) framework, they also proposed a hybrid SAPT/supermolecular model that they called MP2+SDFT, which adds the third-order nonadditive three-body dispersion energy from the coupled Kohn–Sham (CKS) frequency-dependent density susceptibility (FDDS), denoted as Edisp(3)(CKS), to the supermolecular MP2 interaction energy. This model is inspired by the fact that the supermolecular MP2 intrinsically contains two-body and three-body interactions with a partial accounting of electron correlation, but it misses three-body dispersion. Although the MP2+SDFT model has a low scaling of O(N5), their results seems to show that the Edisp(3)(CKS) correction might not be an ideal description of the interactions missing in the supermolecular MP2 approach. This model obtains accurate interaction energies for helium and water trimers, but the results are less accurate for argon and benzene trimers, especially at closer intermolecular distances. Such contrasting results between different systems implies the possibility that the MP2+SDFT model could be less promising for trimer systems with larger dispersion effects.

In 2015, Huang and Beran attempted to replace the CKS dispersion in the MP2+SDFT model of Podeszwa and Szalewicz with a damped ATM dispersion term, and added it to the supermolecular MP2 three-body interaction energy.16 This is called MP2+ATM in their paper. They have shown that MP2+SDFT model (or, as Huang and Beran refer to it, MP2+CKS) overestimates the repulsive nonadditive three-body interaction energy, and they have attributed this discrepancy to some missed higher-order attractive exchange terms in this model. On the other hand, the ATM dispersion is effectively the triple-dipole term of the multipole expansion of the three-body CKS dispersion.17 The three-body dispersion energy becomes less repulsive by dropping out terms with quadrupoles or higher, and its magnitude is further diminished by the empirical damping. As a result of fortuitous error cancellation, the MP2+ATM model provides much better agreement with respect to CCSD(T) than MP2+CKS for various systems including 3B-69.16 The potential energy surface of pyrazole trimer also shows that the error cancellation behaves consistently well along different intermolecular distances. In the end of their work, they suggested that MP2+ATM would be an effective model for treating three-body intermolecular interactions. This study also examined the MP2+ATM model for 65 symmetry-unique trimers in crystalline benzene, and it indicated that MP2+ATM is promising, with a small error of 0.18 kcal mol−1, for the three-body contribution to the crystal lattice energy. Nevertheless, an in-depth analysis of the effectiveness of MP2+ATM, MP2+CKS, and other computational methods over trimers of different intermolecular distances and geometric configurations still remains lacking.

In this work, we provide a more thorough assessment of several approximate methods for computing the three-body contributions to the crystal lattice energy of small organic molecules, i.e., benzene, carbon dioxide, and 1,3,5-triazine, expected to feature significant three-body dispersion effects. We note that all of these molecules lack permanent dipole moments, so that three-body induction/polarization is small in these systems (although, in the case of triazine, it is almost as large as the three-body dispersion; see below). Future work examining additional crystals having significant three-body dispersion and large three-body polarization energies would also be beneficial. The approximate methods considered here include CCSD(T), MP2+ATM, MP2+CKS, MP2.5, and SCS(MI)-CCSD. To assess these methods, we obtain benchmark-quality CCSD(T) three-body interaction energies for a large number of trimers in the crystals, up to 15 Å cutoff for the closest contact distance between monomers. This includes 1977, 1948, and 750 trimers for the benzene, carbon dioxide, and triazine crystals, respectively. Since interaction energies of all significant trimers within these crystals are computed, we are able to analyze the contribution to the three-body lattice energy for various geometrical subsets of the trimers within each crystal based on intermolecular distances between monomers. This analysis reveals that the trimers with at least two monomers with a close intermolecular distance contribute the vast majority of the three-body lattice energies. In turn, this demonstrates that one may reduce computational costs dramatically by ignoring the trimers with monomers far apart from each other and/or by using a cheaper computational method to compute their contributions. Such strategies are quantified for the approximate theoretical methods considered here.

The three-body dispersion energy of a trimer can be approximated as a summation over triples of atoms in the system, using the Axilrod–Teller–Muto formula for three-body triple-dipole dispersion,11,12

EATMabc=C9abc1+3cosθacosθbcosθc(RabRbcRca)3,
(1)

where Rab is the distance between atoms a and b (and similarly for Rbc and Rca); θa, θb, and θc are the angles of the triangle formed by the three atoms, each symbol representing the angle with the corresponding atom as the vertex. C9abc is the dispersion coefficient for the atom triplet abc and can be obtained from the frequency-dependent dipole polarizabilities of the atoms,

C9abc=3π0dωαa(iω)αb(iω)αc(iω).
(2)

In practice, the C9abc coefficient can be approximated using the two-body C6 coefficients as18 

C9abcC6abC6bcC6ca.
(3)

In this work, we utilize C6 coefficients computed for each pair of atoms within each trimer. These C6 coefficients are computed with the DFT-D4 program of Grimme et al.,19–21 and they are based on interpolations of tabulated atomic imaginary frequency-dependent polarizabilities, using information for each atom’s chemical environment.

The total ATM three-body dispersion energy of trimer ABC is obtained by summing EATMabc over all atom triplets from the trimer, in which one atom is on each monomer (terms with two or more atoms on each monomer would contribute to dimer or monomer energies and thus, cancel out when computing the nonadditive three-body energy),

EATM=aAbBcCf(9)EATMabc,
(4)

where f(9) is an empirical damping function that damps the ATM dispersion energy at short interatomic distances. In this work, we implement two different damping functions that are previously mentioned in other studies. In the work of Beran et al. on three-body interactions for 3B-69 and other systems,6,16 they wrote the damping function as a product of three two-body Tang–Toennies (TT) damping functions,22,23

f9abc(β)=f6ab(Rab,β)f6bc(Rbc,β)f6ca(Rca,β),
(5)

where the two-body damping function f6(R, β) is given as

f6(R,β)=1k=06βRkk!eβR.
(6)

β is an empirical parameter that depends on the van der Waals radii rvdW of the atoms involved in f6(R, β),16 

β=0.31(ravdW+rbvdW)+3.43(Bohrs1),
(7)

and the values of rvdW are taken from the representative CCSD van der Waals radii, which are 2.63, 3.34, 3.18, and 3.07 Bohrs for hydrogen, carbon, nitrogen, and oxygen, respectively.22 

Another form of damping function, employed in the DFT-D4 program package, uses a zero-damping scheme proposed by Chai and Head-Gordon (CHG),24 

fdamp(9)R̄abc=11+6R̄abc16,
(8)

where R̄abc is a measure of the average interatomic distance,

R̄abc=RabRbcRca/R0,BJabR0,BJbcR0,BJca1/3.
(9)

In this damping scheme, the quantity R0,BJab is defined as it is in the Becke–Johnson (BJ) damping function, i.e.,25 

R0,BJab=a1R0ab+a2,
(10)

where the R0ab is the cutoff distance of each atom pair,

R0ab=C8abC6ab,
(11)

and the parameters a1 and a2 are dependent on the density functional chosen for a DFT-D calculation. In our work, we use the values for Hartree–Fock, which are a1 = 0.4496 and a2 = 3.3574.

An alternative way to compute the three-body interaction energy arises from the third-order dispersion energy in SAPT(DFT). In two-body SAPT(DFT), the leading contribution of the dispersion energy is of the second order, which corresponds to the product of the FDDS’s of two monomers. Using the Casimir–Polder identity, the two-body dispersion energy can be written as an integration over frequency, ω, of the FDDS product,26,27

Edisp(2)=12π0dωdrAdrAdrBdrB1rArB1rArB×χArA,rA|iωχBrB,rB|iω,
(12)

and transforming the FDDS in the position space into the density-fitting auxiliary basis representation yields28 

Edisp(2)=12π0dωTrS1χAS1χB,
(13)

where S is the Coulomb metric in the auxiliary basis set,

SPQ=P|1r12|Q.
(14)

As discussed in Ref. 14, the coupled FDDS needs to be computed rather than the uncoupled one, in order to obtain accurate dispersion energies. The coupled FDDS can be solved with the time-dependent coupled Kohn–Sham (TD-CKS) theory.14 We recently reported an improved algorithm to compute the dispersion contribution in SAPT(DFT) using hybrid kernels, based on a QR factorization of certain intermediates and using a density-fitting representation of the FDDS and other intermediate quantities. Further details are provided in Ref. 29.

For the three-body case, the leading contribution of the nonadditive three-body dispersion energy is of the third order. The formula of three-body Edisp(3) is similar to that of two-body Edisp(2),13 

Edisp(3)=1π0dωdrAdrAdrBdrBdrCdrC1rArB×1rArC1rBrCχArA,rA|iωχBrB,rB|iω×χBrC,rC|iω,
(15)

or in the auxiliary basis representation,16 

Edisp(3)=1π0dωTrS1χAS1χBS1χC.
(16)

In this work, we investigate the three-body contribution to the lattice energy of crystalline benzene, carbon dioxide, and triazine. The crystal structures are obtained from the Cambridge Structural Database (CSD),30 and we have used an updated development version of the CrystaLattE program package of Borca et al.31 to generate sets of trimer geometries from the CIF files for each crystal structure. In the CrystaLattE program, the crystal lattice energies are computed per monomer, i.e., by selecting a “central” or “reference” monomer and only counting trimers that include this monomer. This set of trimers may include symmetry-equivalent replicas, and CrystaLattE is able to determine the symmetry-unique trimers and how many replicas RM there are of each. The contribution (per monomer) of the symmetry-unique trimer M to the total three-body lattice energy CM is thus computed as

CM=RM×ΔEMint,33,
(17)

where ΔEMint,3 is the nonadditive three-body interaction energy of trimer M.

We compute the ΔEMint,3 for each trimer with various methods, including CCSD(T), MP2, MP2.5, and SCS(MI)-CCSD. We also employ the MP2+ATM and MP2+CKS dispersion correction schemes in addition to these wavefunction methods, where the ATM dispersion term is computed both undamped or damped (with either TT or CHG as the damping function). For the CKS dispersion correction, we use the PBE0 density functional32,33 combined with the gradient-regulated asymptotic correction (GRAC),34 as discussed in Ref. 29. The MP2 interaction energy is computed by performing a two-point CBS extrapolation35 for the MP2 correlation energy using aug-cc-pVTZ and aug-cc-pVQZ basis sets. The Hartree–Fock reference part of the MP2 interaction energy is computed with the aug-cc-pVQZ basis set. For the purpose of readability, in the rest of the paper, we will refer to this scheme for computing MP2 interaction energy as MP2/CBS(a[TQ]Z).

For other wavefunction methods, a focal point approach is used, in which the interaction energy is approximated by adding to the MP2/CBS(a[TQ]Z) interaction energy a correction for higher-order electron correlation, computed as the difference between a higher-level method and MP2, using a relatively small basis set, aug-cc-pVDZ in this work.36,37 By employing the focal point approach, one avoids expensive computations with high-level theory like CCSD(T) with a large basis set. For example, the CCSD(T)/CBS interaction energy is approximated by

E(CCSD(T)/CBS)E(MP2/CBS)+δMP2CCSD(T)/aDZ,
(18)

where the difference between CCSD(T) and MP2 energies is computed under the aug-cc-pVDZ basis set,

δMP2CCSD(T)/aDZ=E(CCSD(T)/aDZ)E(MP2/aDZ).
(19)

The focal point approach exploits the fact that higher-order correlation effects beyond MP2 are usually captured adequately in smaller basis sets,38,39 as demonstrated specifically for intermolecular interactions.36 The interaction energies for MP2.5 and SCS(MI)-CCSD can be computed in a similar way. In the rest of the paper, we will use a shorthand notation for such focal point schemes for convenience; the CCSD(T) interaction energy computed from the focal point approach as shown above will be represented as CCSD(T)/CBS(a[TQ]Z; δ:aDZ), and other methods can be abbreviated similarly. For benzene, we are able to re-use the CCSD(T) results at this level of theory from our recent benchmark study, where we showed excellent agreement with even higher-level CCSD(T)/CBS(a[Q5]Z; δ:aTZ) results for the four closest trimers, where basis set effects should be the most significant (total three-body contribution of 2.66 vs 2.61 kJ mol−1 for these four trimers using the smaller and larger basis sets, respectively).40 The CCSD(T) reference values for carbon dioxide and triazine trimers are newly reported here. The results were obtained using a development version of Psi4 1.641 and utilized density fitting, the frozen core approximation, and the counterpoise correction of Boys and Bernardi.42 There is more than one way to define counterpoise correction for three-body interactions, and we are using the one suggested by Wells and Wilson, which is to use the trimer-centered basis set for computing the energies of all fragments.43 Energies were converged tightly to 10−10 hartree.

When computing the crystal lattice energy using the many-body expansion, an infinite number of n-mers could be generated due to the periodic nature of the system. Therefore, one needs to apply a cutoff criterion based on intermolecular separations, assuming that the interaction energies for n-mers will decay reasonably quickly as the molecules grow further apart. In this work, we follow some definitions and procedures from Ref. 31. First, the closest-contact separation for a dimer is defined as the minimum interatomic distance between two atoms, one of which is on each monomer. Second, the trimer cutoff criterion computes the longest of the three pairwise closest-contact separations within a trimer and ensures that it is shorter than a given cutoff distance. In this work, we will apply a trimer cutoff of 15 Å for the trimers considered in the computation of crystal lattice energies. For convenience, we will use Rmax to represent the largest of the closest separations between the three pairs of monomers in a trimer, and we also define Rmin as the smallest of the three closest-contact separations. To describe the spatial arrangement of the trimer, we also define θmax of trimer as the largest angle of the triangle formed by centers of mass (COMs) of its three monomers.

For each crystal (benzene, carbon dioxide, and triazine), we generated a single set of trimers consistent with a “trimer cutoff” of 15 Å. We will analyze the crystal lattice energy from various subsets based on geometrical parameters defined in Subsection II C. By evaluating the contribution to the three-body lattice energies for different subsets, we can study how the energies depend on the geometrical parameters of the trimers involved.

We first look at how the three-body lattice energy grows as the Rmin cutoff threshold increases. In Figs. 13, the accumulative three-body lattice energy is plotted against the Rmin cutoff value, for which any trimers with Rmin greater than the cutoff are ignored. Each curve corresponds to a level of theory introduced in Subsection II C, where shorthand notations are used to represent different damping schemes for the ATM dispersion correction in the MP2+ATM method: (u) for undamped, (TT) for the Tang–Toennies damping function,22,23 and (CHG) for the Chai–Head-Gordon damping function.24 From these graphs, we observe a few trends among the three molecular crystals.

FIG. 1.

Accumulated three-body lattice energy for crystalline benzene as a function of Rmin cutoff value in kJ mol−1. MP2+CKS represents MP2 corrected by the three-body CKS FDDS dispersion term, and MP2+ATM represents MP2 corrected by ATM dispersion correction. The suffixes (u), (TT), and (CHG) correspond to undamped, damped with Tang–Toennies damping function,22,23 and damped with Chai–Head-Gordon damping function,24 respectively.

FIG. 1.

Accumulated three-body lattice energy for crystalline benzene as a function of Rmin cutoff value in kJ mol−1. MP2+CKS represents MP2 corrected by the three-body CKS FDDS dispersion term, and MP2+ATM represents MP2 corrected by ATM dispersion correction. The suffixes (u), (TT), and (CHG) correspond to undamped, damped with Tang–Toennies damping function,22,23 and damped with Chai–Head-Gordon damping function,24 respectively.

Close modal
FIG. 2.

Accumulated three-body lattice energy for crystalline carbon dioxide as a function of Rmin cutoff value, in kJ mol−1. MP2+CKS represents MP2 corrected by the three-body CKS FDDS dispersion term, and MP2+ATM represents MP2 corrected by ATM dispersion correction. The suffixes (u), (TT), and (CHG) correspond to undamped, damped with Tang–Toennies damping function,22,23 and damped with Chai–Head-Gordon damping function,24 respectively.

FIG. 2.

Accumulated three-body lattice energy for crystalline carbon dioxide as a function of Rmin cutoff value, in kJ mol−1. MP2+CKS represents MP2 corrected by the three-body CKS FDDS dispersion term, and MP2+ATM represents MP2 corrected by ATM dispersion correction. The suffixes (u), (TT), and (CHG) correspond to undamped, damped with Tang–Toennies damping function,22,23 and damped with Chai–Head-Gordon damping function,24 respectively.

Close modal
FIG. 3.

Accumulated three-body lattice energy for crystalline triazine as a function of Rmin cutoff value, in kJ mol−1. MP2+CKS represents MP2 corrected by the three-body CKS FDDS dispersion term, and MP2+ATM represents MP2 corrected by ATM dispersion correction. The suffixes (u), (TT), and (CHG) correspond to undamped, damped with Tang–Toennies damping function,22,23 and damped with Chai–Head-Gordon damping function,24 respectively.

FIG. 3.

Accumulated three-body lattice energy for crystalline triazine as a function of Rmin cutoff value, in kJ mol−1. MP2+CKS represents MP2 corrected by the three-body CKS FDDS dispersion term, and MP2+ATM represents MP2 corrected by ATM dispersion correction. The suffixes (u), (TT), and (CHG) correspond to undamped, damped with Tang–Toennies damping function,22,23 and damped with Chai–Head-Gordon damping function,24 respectively.

Close modal

For all of these crystals, the trimers with Rmin less than 4 Å contribute most of the three-body lattice energy, for all methods involved, and the cumulative lattice energies converge quickly as the Rmin cutoff increases. MP2, which lacks three-body dispersion effects, only accounts for a very small portion of the true three-body contribution to the lattice energy, especially for crystalline benzene, as it only captures 14% of the reference value from CCSD(T)/CBS(a[TQ]Z; δ:aDZ). This indicates the importance of three-body dispersion for the lattice energies of the crystals studied in this work. The CKS or undamped ATM dispersion corrections are shown to significantly overestimate the three-body dispersion for all three systems, and the damping of ATM dispersion results in fortuitous but consistent error cancellation, which agrees with the conclusion of Huang and Beran.16 The SCS(MI)-CCSD lattice energies match the CCSD(T) references very well for the three crystals studied, proving that it is effective as an alternative to CCSD(T) for studying three-body interactions, not only for gas-phase trimers at equilibrium distance, as shown in the 3B-69 benchmark study,6 but also for molecular crystals. On the other hand, the MP2.5 method, which performs well for the 3B-69 systems, lacks consistency in predicting the three-body lattice energies accurately; the MP2.5 three-body lattice energy for crystalline triazine reproduces the CCSD(T) result well, but it overestimates the three-body contribution to the lattice energy for crystalline benzene and underestimates it for crystalline carbon dioxide. Possibly, the MP3 scaling factor of 0.5 in the MP2.5 approach might only be optimal for trimers at near-equilibrium intermolecular distances.

It is also worth mentioning that although Huang and Beran have stated that the ATM formula underestimates the CKS dispersion in Ref. 16, the results for carbon dioxide crystal show the opposite, with the value of ATM dispersion being slightly greater than that of CKS dispersion. As pointed out by Huang and Beran, the damping function for the ATM model partially cancels the error of CKS dispersion due to ignoring higher-order exchange and correlation terms. From our results, the Tang–Toennies damping function damps the ATM dispersion more heavily than the Chai–Head-Gordon damping function and produces more accurate three-body lattice energy as a result of better error cancellation. MP2+ATM(TT) and MP2.5 exhibit similar errors (with different signs) vs CCSD(T) for crystalline carbon dioxide, and MP2+ATM(TT) is even slightly more accurate than MP2.5 for crystalline benzene (although it is not as good for crystalline triazine). Given the cheap computational cost of the ATM dispersion correction, the MP2+ATM approach with the Tang–Toennies damping function seems to be a good approximation to CCSD(T) for computing three-body contributions to crystal lattice energies.

In Tables IIII, the three-body lattice energy contributions for each Rmin domain are presented for each computational method. The trimers within each crystal structure are split into three domains, or subsets, based on their Rmin values: short- (Rmin4Å), medium- (4Å<Rmin7Å), and long-range (7Å<Rmin15Å). In these tables, we also list the percentage contributions of each subset to the total three-body lattice energy for each method. For most computational methods in all three crystals, the short-range subset accounts for more than 90% of the three-body lattice energies, with less than 30% of the number of symmetry-unique trimers considered. In contrast, the long-range subset only contributes 1%–2% of the total three-body lattice energy. We have also listed the errors of the subset contributions with respect to the CCSD(T)/CBS reference values for each method in the parentheses. From these errors, one can conclude that any method that involves proper treatment of three-body dispersion effects, i.e., any method except for MP2, has its error almost completely in the short-range domain. For each method except for MP2 and for all three crystal systems, the total error in the medium- and long-range is less than 0.1 kJ mol−1. From the results above, we can conclude that the three-body dispersion effect at medium- and long-ranges, although being less significant in its magnitude, can be accurately modeled by any of the approximate methods studied in this work (except for MP2), including direct dispersion corrections like ATM and CKS and wavefunction methods like MP2.5 and SCS(MI)-CCSD. This fact agrees with the comments on CKS dispersion correction by Huang and Beran16 that the discrepancy between MP2+CKS and CCSD(T) for three-body interaction energies is because the MP2+CKS model lacks higher-order exchange terms. For trimers with Rmin>4Å, the molecules are apparently far apart enough from each other for any exchange-repulsion effects to vanish. Therefore, the CKS dispersion term, as well as ATM, which is an approximation to it, describes the three-body dispersion effect very accurately at long-range where exchange-repulsion effects do not exist. The percentages and errors of three-body lattice energy contributions in different ranges of Rmin also suggest that it is sufficient to perform CCSD(T) calculations only for the trimers in the short-range domain and compute the contribution in the medium- and long-range domains with a cheaper approximate method such as MP2+CKS or MP2+ATM, while still obtaining very high-accuracy three-body lattice energies with less than 0.1 kJ mol−1 error compared with computing all trimers with CCSD(T), allowing one to avoid more than 70% of the expensive CCSD(T) computations that are not necessary in order to obtain high-quality lattice energies. Even larger savings may be obtained by adding an additional geometric closeness constraint on the trimers to be computed by CCSD(T): we may also require that the largest of the closest contact distances between the monomers, Rmax, is less than some threshold value. Dependencies of the energies based on this cutoff are considered next.

TABLE I.

Benchmark CCSD(T)/CBS(a[TQ]Z; δaDZ) values and errors for approximate methods for three-body contributions to the lattice energy of crystalline benzene (kJ mol−1), split by the separation domain. Cutoff distances (in Å) are determined by the smallest of three intermolecular separations, Rmin.

Range
TotalShortMediumLong
MethodRmin ≤ 15Rmin ≤ 44 < Rmin ≤ 77 < Rmin ≤ 15
Reference values 
CCSD(T)/CBS(a[TQ]Z; δ:aDZ) 3.63 3.28 0.30 0.04 
Relative contribution (%) 100.0 90.4 8.4 1.2 
Errors 
MP2/CBS(a[TQ]Z) −3.12 −2.78 −0.28 −0.06 
MP2.5/CBS(a[TQ]Z; δ:aDZ) 0.78 0.85 −0.05 −0.01 
SCS(MI)-CCSD/CBS(a[TQ]Z; δ:aDZ) 0.34 0.35 −0.00 0.00 
MP2 + CKS disp 2.46 2.44 0.01 0.00 
MP2 + undamped ATM disp 1.32 1.33 −0.01 −0.00 
MP2 + damped(TT) ATM disp 0.60 0.60 −0.01 −0.00 
MP2 + damped(CHG) ATM disp 1.06 1.07 −0.01 −0.00 
Number of unique trimers 1977 494 710 773 
 100.0% 25.0% 35.9% 39.1% 
Range
TotalShortMediumLong
MethodRmin ≤ 15Rmin ≤ 44 < Rmin ≤ 77 < Rmin ≤ 15
Reference values 
CCSD(T)/CBS(a[TQ]Z; δ:aDZ) 3.63 3.28 0.30 0.04 
Relative contribution (%) 100.0 90.4 8.4 1.2 
Errors 
MP2/CBS(a[TQ]Z) −3.12 −2.78 −0.28 −0.06 
MP2.5/CBS(a[TQ]Z; δ:aDZ) 0.78 0.85 −0.05 −0.01 
SCS(MI)-CCSD/CBS(a[TQ]Z; δ:aDZ) 0.34 0.35 −0.00 0.00 
MP2 + CKS disp 2.46 2.44 0.01 0.00 
MP2 + undamped ATM disp 1.32 1.33 −0.01 −0.00 
MP2 + damped(TT) ATM disp 0.60 0.60 −0.01 −0.00 
MP2 + damped(CHG) ATM disp 1.06 1.07 −0.01 −0.00 
Number of unique trimers 1977 494 710 773 
 100.0% 25.0% 35.9% 39.1% 
TABLE II.

Benchmark CCSD(T)/CBS(a[TQ]Z; δaDZ) values and errors for approximate methods for three-body contributions to the lattice energy of crystalline carbon dioxide (kJ mol−1), split by the separation domain. Cutoff distances (in Å) are determined by the smallest of three intermolecular separations, Rmin.

Range
TotalShortMediumLong
MethodRmin ≤ 15Rmin ≤ 44 < Rmin ≤ 77 < Rmin ≤ 15
Reference values 
CCSD(T)/CBS(a[TQ]Z; δ:aDZ) 1.43 1.30 0.13 0.01 
Relative contribution (%) 100.0 90.6 9.0 0.4 
Errors 
MP2/CBS(a[TQ]Z) −1.01 −0.89 −0.11 −0.01 
MP2.5/CBS(a[TQ]Z; δ:aDZ) −0.31 −0.27 −0.04 0.01 
SCS(MI)-CCSD/CBS(a[TQ]Z; δ:aDZ) 0.04 0.03 −0.00 0.02 
MP2 + CKS disp 0.79 0.76 0.01 0.02 
MP2 + undamped ATM disp 0.99 0.91 0.05 0.03 
MP2 + damped(TT) ATM disp 0.39 0.32 0.05 0.03 
MP2 + damped(CHG) ATM disp 0.60 0.52 0.05 0.03 
Number of unique trimers 1948 309 637 1002 
 100.0% 15.9% 32.7% 51.4% 
Range
TotalShortMediumLong
MethodRmin ≤ 15Rmin ≤ 44 < Rmin ≤ 77 < Rmin ≤ 15
Reference values 
CCSD(T)/CBS(a[TQ]Z; δ:aDZ) 1.43 1.30 0.13 0.01 
Relative contribution (%) 100.0 90.6 9.0 0.4 
Errors 
MP2/CBS(a[TQ]Z) −1.01 −0.89 −0.11 −0.01 
MP2.5/CBS(a[TQ]Z; δ:aDZ) −0.31 −0.27 −0.04 0.01 
SCS(MI)-CCSD/CBS(a[TQ]Z; δ:aDZ) 0.04 0.03 −0.00 0.02 
MP2 + CKS disp 0.79 0.76 0.01 0.02 
MP2 + undamped ATM disp 0.99 0.91 0.05 0.03 
MP2 + damped(TT) ATM disp 0.39 0.32 0.05 0.03 
MP2 + damped(CHG) ATM disp 0.60 0.52 0.05 0.03 
Number of unique trimers 1948 309 637 1002 
 100.0% 15.9% 32.7% 51.4% 
TABLE III.

Benchmark CCSD(T)/CBS(a[TQ]Z; δaDZ) values and errors for approximate methods for three-body contributions to the lattice energy of crystalline triazine (kJ mol−1), split by the separation domain. Cutoff distances (in Å) are determined by the smallest of three intermolecular separations, Rmin.

Range
TotalShortMediumLong
MethodRmin ≤ 15Rmin ≤ 44 < Rmin ≤ 77 < Rmin ≤ 15
Reference values 
CCSD(T)/CBS(a[TQ]Z; δ:aDZ) 4.75 4.46 0.21 0.08 
Relative contribution (%) 100.0 93.9 4.3 1.7 
Errors 
MP2/CBS(a[TQ]Z) −2.59 −2.37 −0.16 −0.06 
MP2.5/CBS(a[TQ]Z; δ:aDZ) −0.16 −0.08 −0.05 −0.02 
SCS(MI)-CCSD/CBS(a[TQ]Z; δ:aDZ) −0.01 0.00 −0.01 0.00 
MP2 + CKS disp 2.41 2.40 0.01 0.00 
MP2 + undamped ATM disp 1.42 1.38 0.04 0.01 
MP2 + damped(TT) ATM disp 0.77 0.73 0.03 0.01 
MP2 + damped(CHG) ATM disp 1.25 1.21 0.04 0.01 
Number of unique trimers 750 211 210 329 
 100.0% 28.1% 28.0% 43.9% 
Range
TotalShortMediumLong
MethodRmin ≤ 15Rmin ≤ 44 < Rmin ≤ 77 < Rmin ≤ 15
Reference values 
CCSD(T)/CBS(a[TQ]Z; δ:aDZ) 4.75 4.46 0.21 0.08 
Relative contribution (%) 100.0 93.9 4.3 1.7 
Errors 
MP2/CBS(a[TQ]Z) −2.59 −2.37 −0.16 −0.06 
MP2.5/CBS(a[TQ]Z; δ:aDZ) −0.16 −0.08 −0.05 −0.02 
SCS(MI)-CCSD/CBS(a[TQ]Z; δ:aDZ) −0.01 0.00 −0.01 0.00 
MP2 + CKS disp 2.41 2.40 0.01 0.00 
MP2 + undamped ATM disp 1.42 1.38 0.04 0.01 
MP2 + damped(TT) ATM disp 0.77 0.73 0.03 0.01 
MP2 + damped(CHG) ATM disp 1.25 1.21 0.04 0.01 
Number of unique trimers 750 211 210 329 
 100.0% 28.1% 28.0% 43.9% 

In Subsection III A, we have shown that Rmin is a very good descriptor for the contribution to the three-body crystal lattice energies. However, as stated in Subsection II C, the number of trimers in a periodic crystal system is infinite, and a cutoff solely based on Rmin would not be able to generate a subset with finite number of trimers. Therefore, one also needs to set a cutoff value for Rmax to ensure that one only has to compute the three-body interaction energy for a finite number of trimers to obtain the three-body crystal lattice energy. Similar to the analysis for Rmin, in Figs. 46, we plot the accumulative three-body lattice energy against the Rmax cutoff value.

FIG. 4.

Accumulated three-body lattice energy for crystalline benzene as a function of Rmax cutoff value, in kJ mol−1. MP2+CKS represents MP2 corrected by the three-body CKS FDDS dispersion term, and MP2+ATM represents MP2 corrected by ATM dispersion correction. The suffixes (u), (TT), and (CHG) correspond to undamped, damped with Tang–Toennies damping function,22,23 and damped with Chai–Head-Gordon damping function,24 respectively.

FIG. 4.

Accumulated three-body lattice energy for crystalline benzene as a function of Rmax cutoff value, in kJ mol−1. MP2+CKS represents MP2 corrected by the three-body CKS FDDS dispersion term, and MP2+ATM represents MP2 corrected by ATM dispersion correction. The suffixes (u), (TT), and (CHG) correspond to undamped, damped with Tang–Toennies damping function,22,23 and damped with Chai–Head-Gordon damping function,24 respectively.

Close modal
FIG. 5.

Accumulated three-body lattice energy for crystalline carbon dioxide as a function of Rmax cutoff value, in kJ mol−1. MP2+CKS represents MP2 corrected by the three-body CKS FDDS dispersion term, and MP2+ATM represents MP2 corrected by ATM dispersion correction. The suffixes (u), (TT), and (CHG) correspond to undamped, damped with Tang–Toennies damping function,22,23 and damped with Chai–Head-Gordon damping function,24 respectively.

FIG. 5.

Accumulated three-body lattice energy for crystalline carbon dioxide as a function of Rmax cutoff value, in kJ mol−1. MP2+CKS represents MP2 corrected by the three-body CKS FDDS dispersion term, and MP2+ATM represents MP2 corrected by ATM dispersion correction. The suffixes (u), (TT), and (CHG) correspond to undamped, damped with Tang–Toennies damping function,22,23 and damped with Chai–Head-Gordon damping function,24 respectively.

Close modal
FIG. 6.

Accumulated three-body lattice energy for crystalline triazine as a function of Rmax cutoff value, in kJ mol−1. MP2+CKS represents MP2 corrected by the three-body CKS FDDS dispersion term, and MP2+ATM represents MP2 corrected by ATM dispersion correction. The suffixes (u), (TT), and (CHG) correspond to undamped, damped with Tang–Toennies damping function,22,23 and damped with Chai–Head-Gordon damping function,24 respectively.

FIG. 6.

Accumulated three-body lattice energy for crystalline triazine as a function of Rmax cutoff value, in kJ mol−1. MP2+CKS represents MP2 corrected by the three-body CKS FDDS dispersion term, and MP2+ATM represents MP2 corrected by ATM dispersion correction. The suffixes (u), (TT), and (CHG) correspond to undamped, damped with Tang–Toennies damping function,22,23 and damped with Chai–Head-Gordon damping function,24 respectively.

Close modal

As opposed to the quick convergence for the crystal lattice energies with respect to Rmin, the crystal lattice energies do not show similar quick convergence along different cutoff values of Rmax. This suggests that the long-range interactions are crucial in computing three-body lattice energies accurately, and one needs to consider trimers with relatively large Rmax. In other words, although one of the intermonomer distances (Rmin) needs to be relatively small for the three-body interaction energy to be significant, the longest intermonomer distance (Rmax) can be relatively large. As discussed above, CKS and ATM dispersion corrections overestimate the dispersion contribution to the three-body lattice energies, and the Tang–Toennies damping function provides consistent error cancellations and allows the MP2+ATM model to predict the three-body lattice energies with good accuracy.

The accumulative three-body lattice energies for all three crystals appear to be essentially converged by 10 Å. In order to observe if 10 Å is a good cutoff value of Rmax, we tabulate the contributions to three-body crystal lattice energies, for various computational methods, over different domains of Rmax values in Tables IVVI. For the domains of Rmax, we split the “long-range” domain in two: medium-long-range (7Å<Rmax10Å) and long-range (10Å<Rmax15Å), resulting in four Rmax domains. The larger errors from medium-, medium-long-, and long-range bins in Tables IVVI, compared with those seen in Tables IIII, indicate that errors correlate more strongly with the smallest of three pairwise closest contact distances, Rmin, than they do with the largest one Rmax. From these tables, ignoring trimers with Rmax greater than 10 Å introduces errors around 0.05–0.10 kJ mol−1 (maximum 0.15 kJ mol−1) in the three-body lattice energy at all levels of theory, including the reference method CCSD(T)/CBS. On the other hand, the long-range domain includes more than 80% of the unique-symmetry trimers considered; hence, ignoring them would greatly speed up the computation. Therefore, at least for present systems under consideration, where three-body dispersion is expected to dominate, it appears that three-body dispersion contributions are essentially converged for a cutoff Rmax of 10 Å. Using our trimer energy data, we examined errors caused if we (a) neglected all trimers with Rmax > 10 Å, (b) used CCSD(T)/CBS for the remaining close trimers with Rmin < 4 Å, and (c) used MP2+ATM for the remaining trimers with Rmin > 4 Å. Such a series of approximations yields errors vs the full CCSD(T)/CBS benchmark three-body contribution to the lattice energy of 0.13, 0.14, and −0.03 kJ mol−1 for crystalline benzene, carbon dioxide, and triazine, respectively (the approximations employed to result in some error cancellation). Note that the damping function in the ATM dispersion correction is irrelevant for Rmin>4Å; hence, undamped and damped ATM would produce essentially the same results.

TABLE IV.

Benchmark CCSD(T)/CBS(a[TQ]Z; δaDZ) values and errors for approximate methods for three-body contributions to the lattice energy of crystalline benzene (kJ mol−1), split by the separation domain. Cutoff distances (in Å) are determined by the largest of the three intermolecular separations, Rmax.

Range
TotalShortMediumMedium-longLong
MethodRmax ≤ 15Rmax ≤ 44 < Rmax ≤ 77 < Rmax ≤ 1010 < Rmax ≤ 15
Reference values 
CCSD(T)/CBS(a[TQ]Z; δ:aDZ) 3.63 2.66 0.75 0.35 −0.14 
Relative contribution (%) 100.0 73.3 20.8 9.6 −3.8 
Errors 
MP2/CBS(a[TQ]Z) −3.12 −2.34 −0.77 −0.02 0.02 
MP2.5/CBS(a[TQ]Z; δ:aDZ) 0.78 −0.10 0.50 0.27 0.12 
SCS(MI)-CCSD/CBS(a[TQ]Z; δ:aDZ) 0.34 0.08 0.17 0.07 0.03 
MP2 + CKS disp 2.46 1.16 0.87 0.29 0.14 
MP2 + undamped ATM disp 1.32 0.55 0.49 0.19 0.09 
MP2 + damped(TT) ATM disp 0.60 −0.11 0.40 0.20 0.10 
MP2 + damped(CHG) ATM disp 1.06 0.29 0.49 0.19 0.09 
Number of unique trimers 1977 52 173 1748 
 100.0% 0.2% 2.6% 8.8% 88.4% 
Range
TotalShortMediumMedium-longLong
MethodRmax ≤ 15Rmax ≤ 44 < Rmax ≤ 77 < Rmax ≤ 1010 < Rmax ≤ 15
Reference values 
CCSD(T)/CBS(a[TQ]Z; δ:aDZ) 3.63 2.66 0.75 0.35 −0.14 
Relative contribution (%) 100.0 73.3 20.8 9.6 −3.8 
Errors 
MP2/CBS(a[TQ]Z) −3.12 −2.34 −0.77 −0.02 0.02 
MP2.5/CBS(a[TQ]Z; δ:aDZ) 0.78 −0.10 0.50 0.27 0.12 
SCS(MI)-CCSD/CBS(a[TQ]Z; δ:aDZ) 0.34 0.08 0.17 0.07 0.03 
MP2 + CKS disp 2.46 1.16 0.87 0.29 0.14 
MP2 + undamped ATM disp 1.32 0.55 0.49 0.19 0.09 
MP2 + damped(TT) ATM disp 0.60 −0.11 0.40 0.20 0.10 
MP2 + damped(CHG) ATM disp 1.06 0.29 0.49 0.19 0.09 
Number of unique trimers 1977 52 173 1748 
 100.0% 0.2% 2.6% 8.8% 88.4% 
TABLE V.

Benchmark CCSD(T)/CBS(a[TQ]Z; δaDZ) values and errors for approximate methods for three-body contributions to the lattice energy of crystalline carbon dioxide (kJ mol−1), split by the separation domain. Cutoff distances (in Å) are determined by the largest of the three intermolecular separations, Rmax.

Range
TotalShortMediumMedium-longLong
MethodRmax ≤ 15Rmax ≤ 44 < Rmax ≤ 77 < Rmax ≤ 1010 < Rmax ≤ 15
Reference values 
CCSD(T)/CBS(a[TQ]Z; δ:aDZ) 1.43 0.86 −0.05 0.72 −0.09 
Relative contribution (%) 100.0 60.1 −3.5 49.9 −6.6 
Errors 
MP2/CBS(a[TQ]Z) −1.01 −0.79 −0.17 −0.05 0.00 
MP2.5/CBS(a[TQ]Z; δ:aDZ) −0.31 −0.34 −0.02 0.03 0.01 
SCS(MI)-CCSD/CBS(a[TQ]Z; δ:aDZ) 0.04 −0.03 0.04 0.02 0.01 
MP2 + CKS disp 0.79 0.27 0.44 0.05 0.03 
MP2 + undamped ATM disp 0.99 0.46 0.46 0.03 0.04 
MP2 + damped(TT) ATM disp 0.39 −0.04 0.36 0.04 0.04 
MP2 + damped(CHG) ATM disp 0.60 0.09 0.45 0.03 0.04 
Number of unique trimers 1948 27 159 1760 
 100.0% 0.1% 1.4% 8.2% 90.3% 
Range
TotalShortMediumMedium-longLong
MethodRmax ≤ 15Rmax ≤ 44 < Rmax ≤ 77 < Rmax ≤ 1010 < Rmax ≤ 15
Reference values 
CCSD(T)/CBS(a[TQ]Z; δ:aDZ) 1.43 0.86 −0.05 0.72 −0.09 
Relative contribution (%) 100.0 60.1 −3.5 49.9 −6.6 
Errors 
MP2/CBS(a[TQ]Z) −1.01 −0.79 −0.17 −0.05 0.00 
MP2.5/CBS(a[TQ]Z; δ:aDZ) −0.31 −0.34 −0.02 0.03 0.01 
SCS(MI)-CCSD/CBS(a[TQ]Z; δ:aDZ) 0.04 −0.03 0.04 0.02 0.01 
MP2 + CKS disp 0.79 0.27 0.44 0.05 0.03 
MP2 + undamped ATM disp 0.99 0.46 0.46 0.03 0.04 
MP2 + damped(TT) ATM disp 0.39 −0.04 0.36 0.04 0.04 
MP2 + damped(CHG) ATM disp 0.60 0.09 0.45 0.03 0.04 
Number of unique trimers 1948 27 159 1760 
 100.0% 0.1% 1.4% 8.2% 90.3% 
TABLE VI.

Benchmark CCSD(T)/CBS(a[TQ]Z; δaDZ) values and errors for approximate methods for three-body contributions to the lattice energy of crystalline triazine (kJ mol−1), split by the separation domain. Cutoff distances (in Å) are determined by the largest of the three intermolecular separations, Rmax.

Range
TotalShortMediumMedium-longLong
MethodRmax ≤ 15Rmax ≤ 44 < Rmax ≤ 77 < Rmax ≤ 1010 < Rmax ≤ 15
Reference values 
CCSD(T)/CBS(a[TQ]Z; δ:aDZ) 4.75 3.62 1.17 −0.11 0.07 
Relative contribution (%) 100.0 76.3 24.6 −2.3 1.4 
Errors 
MP2/CBS(a[TQ]Z) −2.59 −2.27 −0.29 −0.03 −0.01 
MP2.5/CBS(a[TQ]Z; δ:aDZ) −0.16 −0.57 0.21 0.16 0.04 
SCS(MI)-CCSD/CBS(a[TQ]Z; δ:aDZ) −0.01 −0.09 0.03 0.03 0.01 
MP2 + CKS disp 2.41 1.21 0.75 0.38 0.08 
MP2 + undamped ATM disp 1.42 0.77 0.44 0.16 0.05 
MP2 + damped(TT) ATM disp 0.77 0.12 0.40 0.20 0.05 
MP2 + damped(CHG) ATM disp 1.25 0.60 0.44 0.16 0.05 
Number of unique trimers 750 15 98 634 
 100.0% 0.4% 2.0% 13.1% 84.5% 
Range
TotalShortMediumMedium-longLong
MethodRmax ≤ 15Rmax ≤ 44 < Rmax ≤ 77 < Rmax ≤ 1010 < Rmax ≤ 15
Reference values 
CCSD(T)/CBS(a[TQ]Z; δ:aDZ) 4.75 3.62 1.17 −0.11 0.07 
Relative contribution (%) 100.0 76.3 24.6 −2.3 1.4 
Errors 
MP2/CBS(a[TQ]Z) −2.59 −2.27 −0.29 −0.03 −0.01 
MP2.5/CBS(a[TQ]Z; δ:aDZ) −0.16 −0.57 0.21 0.16 0.04 
SCS(MI)-CCSD/CBS(a[TQ]Z; δ:aDZ) −0.01 −0.09 0.03 0.03 0.01 
MP2 + CKS disp 2.41 1.21 0.75 0.38 0.08 
MP2 + undamped ATM disp 1.42 0.77 0.44 0.16 0.05 
MP2 + damped(TT) ATM disp 0.77 0.12 0.40 0.20 0.05 
MP2 + damped(CHG) ATM disp 1.25 0.60 0.44 0.16 0.05 
Number of unique trimers 750 15 98 634 
 100.0% 0.4% 2.0% 13.1% 84.5% 

We have computed the three-body contributions to the lattice energies of crystalline benzene, carbon dioxide, and triazine, using various computational methods, including CCSD(T), MP2, MP2.5, and SCS(MI)-CCSD. In these crystals, the three-body dispersion is crucial in estimating the crystal lattice energies, and MP2 is known to be lacking three-body dispersion. Therefore, we have also employed the MP2+CKS and MP2+ATM methods, which add the coupled Kohn–Sham and Axilrod–Teller–Muto dispersion corrections (with the Tang–Toennies22,23 or Chai–Head-Gordon24 damping function) to the MP2 interaction energies, to provide the three-body dispersion effects missing in MP2. Our results show that MP2+CKS significantly overestimates the three-body dispersion contribution to the lattice energies for these crystals, and MP2+ATM with the Tang–Toennies damping scheme provides a fortuitous but consistent error cancellation and matches well with the reference values from CCSD(T). Interestingly enough, the outstanding performance of MP2+ATM(TT) was also observed in predicting liquid state properties,44,45 indicating that the Tang–Toennies damping scheme might be a good universal choice to combine with the ATM dispersion formula when estimating three-body dispersion energies. SCS(MI)-CCSD retains its outstanding performance in estimating the three-body interaction energies for the gas-phase trimer systems6 in computing the three-body lattice energies in molecular crystals, and MP2.5 becomes less robust in crystal systems despite its good performance in equilibrium-geometry gas-phase trimers.

We also showed that the three-body contributions to lattice energies converge rapidly as the smallest of three pairwise intermonomer closest contact distances Rmin grows. The trimers with Rmin<4Å contribute to more than 90% of the total three-body lattice energies for all three crystals and with almost all computational methods studied. For trimers with Rmin > 4 Å, the approximate methods considered, except for MP2, become much more accurate (in part due to the reduced magnitude of the interactions).

The three-body contributions to lattice energies, especially the dispersion contribution, also essentially converge as the largest of the three intermonomer closest contact distances Rmax exceeds 10 Å in the crystals studied. Such an approximation drastically reduces the computational resources needed to compute the three-body lattice energies, which involves three-body interaction energy calculations for a massive set of trimers. The errors introduced by this approximation are less than 0.15 kJ mol−1 for the crystals studied here at all levels of theory considered. In systems with larger three-body induction effects, which we might expect in molecular crystals of molecules with non-zero dipole moments, larger cutoff values for Rmin and Rmax could be needed. Nevertheless, we would expect MP2 to be satisfactory in describing the long-range induction interactions.

The supplementary material contains csv files with the raw data analyzed in this study (three-body lattice energy contributions computed with various theoretical methods) and geometrical information about the trimers. It also contains Cartesian coordinates (as xyz files) for all trimers considered in this study.

The authors gratefully acknowledge financial support from the U.S. National Science Foundation through Grant No. CHE-1955940.

The authors have no conflicts to disclose.

Yi Xie: Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (supporting); Software (supporting); Visualization (lead); Writing – original draft (lead); Writing – review & editing (supporting). Zachary L. Glick: Formal analysis (supporting); Methodology (equal); Software (lead); Validation (lead). C. David Sherrill: Conceptualization (lead); Formal analysis (supporting); Funding acquisition (lead); Investigation (supporting); Methodology (equal); Project administration (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available with in the article and its supplementary material.

1.
D.
Hankins
,
J. W.
Moskowitz
, and
F. H.
Stillinger
,
J. Chem. Phys.
53
,
4544
(
1970
).
2.
H.
Stoll
and
H.
Preuß
,
Theor. Chim. Acta
46
,
11
(
1977
).
3.
R.
Podeszwa
,
B. M.
Rice
, and
K.
Szalewicz
,
Phys. Rev. Lett.
101
,
115503
(
2008
).
4.
S.
Wen
and
G. J. O.
Beran
,
J. Chem. Theory Comput.
7
,
3733
(
2011
).
5.
M.
Kamiya
,
S.
Hirata
, and
M.
Valiev
,
J. Chem. Phys.
128
,
074103
(
2008
).
6.
J.
Řezáč
,
Y.
Huang
,
P.
Hobza
, and
G. J. O.
Beran
,
J. Chem. Theory Comput.
11
,
3065
(
2015
).
7.
M.
Pitoňák
,
P.
Neogrády
,
J.
Černý
,
S.
Grimme
, and
P.
Hobza
,
Chem. Phys. Chem.
10
,
282
(
2009
).
8.
T.
Takatani
,
E. G.
Hohenstein
, and
C. D.
Sherrill
,
J. Chem. Phys.
128
,
124111
(
2008
).
9.
M.
Pitoňák
,
J.
Řezác
, and
P.
Hobza
,
Phys. Chem. Chem. Phys.
12
,
9611
(
2010
).
10.
A.
Sebetci
and
G. J. O.
Beran
,
J. Chem. Theory Comput.
6
,
155
(
2010
).
11.
B. M.
Axilrod
and
E.
Teller
,
J. Chem. Phys.
11
,
299
(
1943
).
12.
Y.
Muto
,
Proc. Phys. Math. Soc. Jpn.
17
,
629
(
1943
).
13.
R.
Podeszwa
and
K.
Szalewicz
,
J. Chem. Phys.
126
,
194101
(
2007
).
14.
A. J.
Misquitta
,
R.
Podeszwa
,
B.
Jeziorski
, and
K.
Szalewicz
,
J. Chem. Phys.
123
,
214103
(
2005
).
15.
A.
Heßelmann
,
G.
Jansen
, and
M.
Schütz
,
J. Chem. Phys.
122
,
014103
(
2005
).
16.
Y.
Huang
and
G. J. O.
Beran
,
J. Chem. Phys.
143
,
044113
(
2015
).
17.
R. J.
Bell
,
J. Phys. B: At. Mol. Phys.
3
,
751
(
1970
).
18.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
,
J. Chem. Phys.
132
,
154104
(
2010
).
19.
E.
Caldeweyher
,
C.
Bannwarth
, and
S.
Grimme
,
J. Chem. Phys.
147
,
034112
(
2017
).
20.
E.
Caldeweyher
 et al,
J. Chem. Phys.
150
,
154122
(
2019
).
21.
E.
Caldeweyher
,
J.-M.
Mewes
,
S.
Ehlert
, and
S.
Grimme
,
Phys. Chem. Chem. Phys.
22
,
8499
(
2020
).
22.
A.
Tkatchenko
,
R. A.
DiStasio
,
R.
Car
, and
M.
Scheffler
,
Phys. Rev. Lett.
108
,
236402
(
2012
).
23.
K. T.
Tang
and
J. P.
Toennies
,
J. Chem. Phys.
80
,
3726
(
1984
).
24.
J.-D.
Chai
and
M.
Head-Gordon
,
Phys. Chem. Chem. Phys.
10
,
6615
(
2008
).
25.
S.
Grimme
,
S.
Ehrlich
, and
L.
Goerigk
,
J. Comput. Chem.
32
,
1456
(
2011
).
26.
M.
Jaszunski
and
R.
McWeeny
,
Mol. Phys.
55
,
1275
(
1985
).
27.
P. J.
Knowles
and
W. J.
Meath
,
Mol. Phys.
59
,
965
(
1986
).
28.
M.
Pitoňák
and
A.
Heßelmann
,
J. Chem. Theory Comput.
6
,
168
(
2010
).
29.
Y.
Xie
,
D. G. A.
Smith
, and
C. D.
Sherrill
,
J. Chem. Phys.
157
,
024801
(
2022
).
30.
F. H.
Allen
,
Acta Crystallogr., Sect. B: Struct. Sci.
58
,
380
(
2002
).
31.
C. H.
Borca
,
B. W.
Bakr
,
L. A.
Burns
, and
C. D.
Sherrill
,
J. Chem. Phys.
151
,
144103
(
2019
).
32.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
33.
C.
Adamo
and
V.
Barone
,
J. Chem. Phys.
110
,
6158
(
1999
).
34.
M.
Grüning
,
O. V.
Gritsenko
,
S. J. A.
van Gisbergen
, and
E. J.
Baerends
,
J. Chem. Phys.
114
,
652
(
2001
).
35.
A.
Halkier
 et al,
Chem. Phys. Lett.
286
,
243
(
1998
).
36.
L. A.
Burns
,
M. S.
Marshall
, and
C. D.
Sherrill
,
J. Chem. Theory Comput.
10
,
49
(
2014
).
37.
M. O.
Sinnokrot
,
E. F.
Valeev
, and
C. D.
Sherrill
,
J. Am. Chem. Soc.
124
,
10887
(
2002
).
38.
A. L. L.
East
and
W. D.
Allen
,
J. Chem. Phys.
99
,
4638
(
1993
).
39.
A. G.
Császár
,
W. D.
Allen
, and
H. F.
Schaefer
,
J. Chem. Phys.
108
,
9751
(
1998
).
40.
C. H.
Borca
,
Z. L.
Glick
,
D. P.
Metcalf
,
L. A.
Burns
, and
C. D.
Sherrill
, “
Benchmark coupled-cluster lattice energy of crystalline benzene, and assessment of multi-level approximations in the many-body expansion
,”
J. Chem. Phys.
(submitted).
41.
D. G. A.
Smith
 et al,
J. Chem. Phys.
152
,
184108
(
2020
).
42.
S. F.
Boys
and
F.
Bernardi
,
Mol. Phys.
19
,
553
(
1970
).
43.
B. H.
Wells
and
S.
Wilson
,
Chem. Phys. Lett.
101
,
429
(
1983
).
44.
R.
Bukowski
and
K.
Szalewicz
,
J. Chem. Phys.
114
,
9518
(
2001
).
45.
J. G.
McDaniel
and
J. R.
Schmidt
,
J. Phys. Chem. B
118
,
8042
(
2014
).

Supplementary Material