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, *R*_{min}, shows a strong correlation with the three-body contribution to the lattice energy, and, here, the largest of the closest-contact distances, *R*_{max}, serves as a cutoff criterion to limit the number of trimers to be considered. We considered all trimers up to $Rmax=15A\u030a$. The trimers with $Rmin<4A\u030a$ 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>4A\u030a$, 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=10A\u030a$. From these results, it appears that in molecular crystals where dispersion dominates the three-body contribution to the lattice energy, the trimers with $Rmin>4A\u030a$ can be computed with the MP2+ATM method to reduce the computational cost, and those with $Rmax>10A\u030a$ appear to be basically negligible.

## I. INTRODUCTION

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(N\u22121)2$, as opposed to 2^{N} − 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 benzene^{3,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*(*N*^{6}), restricting their usefulness for important practical applications like computing the lattice energies of molecular crystals. They have also used an *ab initio* force field^{4,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*(*N*^{5}), 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.

## II. THEORY AND METHOD

### A. Axilrod–Teller–Muto three-body dispersion

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}

where *R*_{ab} is the distance between atoms a and b (and similarly for *R*_{bc} and *R*_{ca}); *θ*_{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,

In practice, the $C9abc$ coefficient can be approximated using the two-body *C*_{6} coefficients as^{18}

In this work, we utilize *C*_{6} coefficients computed for each pair of atoms within each trimer. These *C*_{6} 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),

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}

where the two-body damping function *f*_{6}(*R*, *β*) is given as

*β* is an empirical parameter that depends on the van der Waals radii *r*^{vdW} of the atoms involved in *f*_{6}(*R*, *β*),^{16}

and the values of *r*^{vdW} 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}

where $R\u0304abc$ is a measure of the average interatomic distance,

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

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

and the parameters *a*_{1} and *a*_{2} are dependent on the density functional chosen for a DFT-D calculation. In our work, we use the values for Hartree–Fock, which are *a*_{1} = 0.4496 and *a*_{2} = 3.3574.

### B. Three-body CKS FDDS dispersion

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}

and transforming the FDDS in the position space into the density-fitting auxiliary basis representation yields^{28}

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

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}

or in the auxiliary basis representation,^{16}

### C. Three-body crystal lattice energies

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

where $\Delta EMint,3$ is the nonadditive three-body interaction energy of trimer *M*.

We compute the $\Delta 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 functional^{32,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 extrapolation^{35} 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

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

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.6^{41} 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 *R*_{max} to represent the largest of the closest separations between the three pairs of monomers in a trimer, and we also define *R*_{min} 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.

## III. RESULTS

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.

### A. Dependence of energies on the minimum intermonomer distance, *R*_{min}

We first look at how the three-body lattice energy grows as the *R*_{min} cutoff threshold increases. In Figs. 1–3, the accumulative three-body lattice energy is plotted against the *R*_{min} cutoff value, for which any trimers with *R*_{min} 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.

For all of these crystals, the trimers with *R*_{min} less than 4 Å contribute most of the three-body lattice energy, for all methods involved, and the cumulative lattice energies converge quickly as the *R*_{min} 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 I–III, the three-body lattice energy contributions for each *R*_{min} domain are presented for each computational method. The trimers within each crystal structure are split into three domains, or subsets, based on their *R*_{min} values: short- $(Rmin\u22644A\u030a)$, medium- $(4A\u030a<Rmin\u22647A\u030a)$, and long-range $(7A\u030a<Rmin\u226415A\u030a)$. 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 Beran^{16} 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>4A\u030a$, 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 *R*_{min} 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, *R*_{max}, is less than some threshold value. Dependencies of the energies based on this cutoff are considered next.

. | Range . | |||
---|---|---|---|---|

. | Total . | Short . | Medium . | Long . |

Method . | R_{min} ≤ 15
. | R_{min} ≤ 4
. | 4 < R_{min} ≤ 7
. | 7 < R_{min} ≤ 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 . | |||
---|---|---|---|---|

. | Total . | Short . | Medium . | Long . |

Method . | R_{min} ≤ 15
. | R_{min} ≤ 4
. | 4 < R_{min} ≤ 7
. | 7 < R_{min} ≤ 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 . | |||
---|---|---|---|---|

. | Total . | Short . | Medium . | Long . |

Method . | R_{min} ≤ 15
. | R_{min} ≤ 4
. | 4 < R_{min} ≤ 7
. | 7 < R_{min} ≤ 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 . | |||
---|---|---|---|---|

. | Total . | Short . | Medium . | Long . |

Method . | R_{min} ≤ 15
. | R_{min} ≤ 4
. | 4 < R_{min} ≤ 7
. | 7 < R_{min} ≤ 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 . | |||
---|---|---|---|---|

. | Total . | Short . | Medium . | Long . |

Method . | R_{min} ≤ 15
. | R_{min} ≤ 4
. | 4 < R_{min} ≤ 7
. | 7 < R_{min} ≤ 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 . | |||
---|---|---|---|---|

. | Total . | Short . | Medium . | Long . |

Method . | R_{min} ≤ 15
. | R_{min} ≤ 4
. | 4 < R_{min} ≤ 7
. | 7 < R_{min} ≤ 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% |

### B. Dependence of energies on the maximum intermonomer distance, *R*_{max}

In Subsection III A, we have shown that *R*_{min} 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 *R*_{min} would not be able to generate a subset with finite number of trimers. Therefore, one also needs to set a cutoff value for *R*_{max} 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 *R*_{min}, in Figs. 4–6, we plot the accumulative three-body lattice energy against the *R*_{max} cutoff value.

As opposed to the quick convergence for the crystal lattice energies with respect to *R*_{min}, the crystal lattice energies do not show similar quick convergence along different cutoff values of *R*_{max}. 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 *R*_{max}. In other words, although one of the intermonomer distances (*R*_{min}) needs to be relatively small for the three-body interaction energy to be significant, the longest intermonomer distance (*R*_{max}) 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 *R*_{max}, we tabulate the contributions to three-body crystal lattice energies, for various computational methods, over different domains of *R*_{max} values in Tables IV–VI. For the domains of *R*_{max}, we split the “long-range” domain in two: medium-long-range ($7A\u030a<Rmax\u226410A\u030a$) and long-range ($10A\u030a<Rmax\u226415A\u030a$), resulting in four *R*_{max} domains. The larger errors from medium-, medium-long-, and long-range bins in Tables IV–VI, compared with those seen in Tables I–III, indicate that errors correlate more strongly with the smallest of three pairwise closest contact distances, *R*_{min}, than they do with the largest one *R*_{max}. From these tables, ignoring trimers with *R*_{max} 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 *R*_{max} of 10 Å. Using our trimer energy data, we examined errors caused if we (a) neglected all trimers with *R*_{max} > 10 Å, (b) used CCSD(T)/CBS for the remaining close trimers with *R*_{min} < 4 Å, and (c) used MP2+ATM for the remaining trimers with *R*_{min} > 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>4A\u030a$; hence, undamped and damped ATM would produce essentially the same results.

. | Range . | ||||
---|---|---|---|---|---|

. | Total . | Short . | Medium . | Medium-long . | Long . |

Method . | R_{max} ≤ 15
. | R_{max} ≤ 4
. | 4 < R_{max} ≤ 7
. | 7 < R_{max} ≤ 10
. | 10 < R_{max} ≤ 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 | 4 | 52 | 173 | 1748 |

100.0% | 0.2% | 2.6% | 8.8% | 88.4% |

. | Range . | ||||
---|---|---|---|---|---|

. | Total . | Short . | Medium . | Medium-long . | Long . |

Method . | R_{max} ≤ 15
. | R_{max} ≤ 4
. | 4 < R_{max} ≤ 7
. | 7 < R_{max} ≤ 10
. | 10 < R_{max} ≤ 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 | 4 | 52 | 173 | 1748 |

100.0% | 0.2% | 2.6% | 8.8% | 88.4% |

. | Range . | ||||
---|---|---|---|---|---|

. | Total . | Short . | Medium . | Medium-long . | Long . |

Method . | R_{max} ≤ 15
. | R_{max} ≤ 4
. | 4 < R_{max} ≤ 7
. | 7 < R_{max} ≤ 10
. | 10 < R_{max} ≤ 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 | 2 | 27 | 159 | 1760 |

100.0% | 0.1% | 1.4% | 8.2% | 90.3% |

. | Range . | ||||
---|---|---|---|---|---|

. | Total . | Short . | Medium . | Medium-long . | Long . |

Method . | R_{max} ≤ 15
. | R_{max} ≤ 4
. | 4 < R_{max} ≤ 7
. | 7 < R_{max} ≤ 10
. | 10 < R_{max} ≤ 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 | 2 | 27 | 159 | 1760 |

100.0% | 0.1% | 1.4% | 8.2% | 90.3% |

. | Range . | ||||
---|---|---|---|---|---|

. | Total . | Short . | Medium . | Medium-long . | Long . |

Method . | R_{max} ≤ 15
. | R_{max} ≤ 4
. | 4 < R_{max} ≤ 7
. | 7 < R_{max} ≤ 10
. | 10 < R_{max} ≤ 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 | 3 | 15 | 98 | 634 |

100.0% | 0.4% | 2.0% | 13.1% | 84.5% |

. | Range . | ||||
---|---|---|---|---|---|

. | Total . | Short . | Medium . | Medium-long . | Long . |

Method . | R_{max} ≤ 15
. | R_{max} ≤ 4
. | 4 < R_{max} ≤ 7
. | 7 < R_{max} ≤ 10
. | 10 < R_{max} ≤ 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 | 3 | 15 | 98 | 634 |

100.0% | 0.4% | 2.0% | 13.1% | 84.5% |

## IV. CONCLUSIONS

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–Toennies^{22,23} or Chai–Head-Gordon^{24} 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 systems^{6} 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 *R*_{min} grows. The trimers with $Rmin<4A\u030a$ 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 *R*_{min} > 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 *R*_{max} 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 *R*_{min} and *R*_{max} could be needed. Nevertheless, we would expect MP2 to be satisfactory in describing the long-range induction interactions.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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