research papers
A transferable quantum mechanical energy model for intermolecular interactions using a single empirical parameter
^{a}School of Molecular and Life Sciences, Curtin University, Perth, Western Australia 6845, Australia, and ^{b}School of Molecular Sciences, University of Western Australia, Perth, Western Australia 6009, Australia
^{*}Correspondence email: peter.spackman@curtin.edu.au
The calculation of intermolecular interactions in molecular crystals using model energies provides a unified route to understanding the complex interplay of driving forces in crystallization, elastic properties and more. Presented here is a new singleparameter interaction energy model (CE1p), extending the previous CrystalExplorer energy model and calibrated using density functional theory (DFT) calculations at the ωB97MV/def2QZVP level over 1157 intermolecular interactions from 147 crystal structures. The new model incorporates an improved treatment of dispersion interactions and polarizabilities using the exchangehole dipole model (XDM), along with the use of effective core potentials (ECPs), facilitating application to molecules containing elements across the periodic table (from H to Rn). This new model is validated against highlevel reference data with outstanding performance, comparable to stateoftheart DFT methods for molecular energies over the X23 set (mean absolute deviation 3.6 kJ mol^{−1}) and for intermolecular interactions in the S66x8 benchmark set (root meansquare deviation 3.3 kJ mol^{−1}). The performance of this model is further examined compared to the GFN2xTB tightbinding model, providing recommendations for the evaluation of intermolecular interactions in molecular crystal systems.
1. Introduction
A detailed quantitative description of intermolecular interactions – at low computational cost – is essential to the modelling of molecular crystals and many other chemical problems. This includes rationalization of crystal packing, the factors driving crystal growth and dissolution, mechanical properties of crystals and more. The study of these interactions is well established, but the prediction and understanding of the relative energetic contributions and terms driving these interactions is as important as ever. Moreover, there is an everexpanding taxonomy of intermolecular interactions within the study of molecular crystals (and more generally), which now include not only terms such as hydrogen bonding and π–π stacking but also halogen, chalcogen, pnictogen and even aerogen bonding. Thus, accessible and accurate methods to compute both intermolecular interaction energies and their substituent factors are extremely useful, not only for their predictive power but also for their capacity to provide a unified view of intermolecular interactions as not simply a collection of categorized interaction types or intermolecular `bonds', but in terms of the fundamental underlying physical forces and concepts.
Energetic approaches to understanding and rationalizing intermolecular interactions can be broken down broadly into two categories, composition and decomposition. The former encompasses models that build up the total interaction energy from prediction of separate components, and the latter encompasses those that start from a total energy and break it down into components after the fact. There are a myriad approaches to this end, with some components (e.g. Coulomb or electrostatic interactions) being more universally present than others. Further, the boundary between these two approaches is often unclear – some terms may be naturally more separable than others due to their mathematical/physical derivation and, likewise, it is always possible to produce decompositions that may be more artificial than meaningful. In either case, though, the motivation is not only to predict the total energy of a given system, but also to understand the components which lead to this energy and therefore, in the context of a molecular crystal, rationalize why a particular interaction may be present (or not) in the experimentally observed structure.
For molecular crystalline solids, a natural composition (or decomposition) approach is to consider the total energy as a sum over all pairwise dimer energies within a given range. Calculations of conventional quantum mechanical (QM) dimer interaction energies involve three full selfconsistent field (SCF) calculations, where the interaction energy of the system is computed as the difference between the total energy of the dimer, AB, and the sum of the total energies of the monomers, A and B, i.e. E_{int} = E_{AB} − (E_{A} + E_{B}). In the absence of a complete (or nearcomplete) basis set this can be complicated by the basisset superposition error (BSSE), necessitating a counterpoise correction using either the full dimer basis set (Boys & Bernardi, 1970) or approximations to this correction such as the geometric counterpoise correction (GCP) (Kruse & Grimme, 2012).
We have previously demonstrated the value of an alternative approach to computing dimer energies that involves fitting scale factors to components from a QM method based on the Hartree–Fock (HF) interaction energy (Su & Li, 2009) to accelerate the determination of these interaction energies (Turner et al., 2014; Mackenzie et al., 2017) and made this available through the CrystalExplorer software (Spackman et al., 2021), where it has found widespread use, particularly in the field of crystal engineering. This model avoids the full dimer (AB) SCF through approximating the change in molecular orbitals (MOs) from A interacting with B via superposition of the two wavefunctions, along with a single orthogonalization step for the orbitals to respond to one another. This model, when coupled with scaling coefficients and force fieldlike terms for polarization and dispersion, has proved to be an effective and useful model for intermolecular interactions to aid in rationalizing molecular packing (Tan et al., 2019), mechanical slip planes (Wang & Sun, 2018), halogen bonding (Brammer, 2017) and even prediction of crystal growth (Spackman et al., 2023).
Despite their overall success, the previous models (CEHF and CEB3LYP) suffered from several primary (technical) shortcomings:
(i) The choice of the 631G(d,p) basis for the more accurate CEB3LYP model prevented its application to many heavier elements.
(ii) The global dispersion model D2 (Grimme, 2006) does not depend on the monomer wavefunctions, which limits its accuracy.
(iii) The use of free atomic polarizabilities was a more extreme approximation, as such quantities change significantly when the atom is part of a molecule.
It is not only timely to address some shortcomings of the previous method, but also to highlight the value that the CrystalExplorer model energies, and other models such as the density functional tightbinding (DFTB) method GFN2xTB (Bannwarth et al., 2019), can provide in our classification and understanding of intermolecular interactions. Characterization of the strengths, weaknesses and range of applicability of these models, and making available a cohesive implementation to use both in the context of CrystalExplorer and on highperformance computing systems, is our fundamental purpose here. Thus in this work we aim to address these issues, while producing a model that is more accurate and transferable wherever possible.
To the above end, the interaction model, along with its software implementation, has been updated to incorporate the use of effective core potentials (ECPs, see Section S2.2 in the supporting information) which model core electrons in heavier elements. This not only increases the speed of the calculations for heavier elements, but enables the use of basis sets parameterized for use with ECPs, for example the Ahlrichs def2 basis sets (Weigend & Ahlrichs, 2005), which we make use of in this work. Furthermore, the previously used D2 dispersion model (Grimme, 2006) has been replaced with the theoretically motivated exchangehole dipole model (XDM) (Johnson & Becke, 2006; Becke & Johnson, 2007; OterodelaRoza & Johnson, 2012b). This yields moleculedependent dispersion parameters while also facilitating an elegant solution for atominmolecule polarizabilities, and is itself a drastic improvement in the treatment of dispersion.
2. Methods
2.1. Reference values used in the fitting procedure
In order to fit and evaluate new models, a modified version of the training set of molecular crystals used in the previous CE model fitting (Mackenzie et al., 2017) was chosen. The revised set consists of 1157 molecule/ion pairs extracted from 147 organic, inorganic and organometallic/metal–organic molecular crystal structures, incorporating atoms up to and including bromine, and also iodine and xenon. For each molecular molecule/ion pairs were obtained by generating a cluster of nearestneighbour interactions around each symmetryunique molecule/ion in the crystal. A complete list of the crystal structures, including Cambridge Structural Database (CSD) refcodes (Groom et al., 2016) and Inorganic Database (ICSD) identifiers (Levin, 2020), is available in the supporting information.
Reference ωB97MV/def2QZVP interaction energies were then calculated via evaluation of the full dimer SCF energy, subtracting the monomer energies (in the monomer basis). The ωB97MV method (Mardirossian & HeadGordon, 2016) was selected because of its general accuracy and transferability across many problems, including noncovalent interactions (Najibi & Goerigk, 2018). The large def2QZVP (Weigend et al., 2003; Weigend & Ahlrichs, 2005; Peterson et al., 2003) basis set was chosen to minimize BSSE where possible, while still being computationally feasible for dimers comprising hundreds of atoms. All calculations were performed using the ORCA software suite (Version 5; Neese, 2022). After evaluating reference energies, pairs with a significant amount of charge transfer were removed from the training set. This was deemed necessary since charge transfer is not an aspect we are attempting to model in this work, and so inclusion of values where this is a significant contribution to the interaction energy would only serve to confound the results of the fitting procedure(s) and error estimates. Furthermore, particularly in the context of molecular crystals, the charge transfer occurring in a gasphase dimer does not necessarily correspond well to that present when the dimer is embedded in a crystalline environment with competing chargetransfer possibilities. Significant charge transfer, where a large portion of an electron has moved from one monomer to another, was, for our purposes, defined as
where q_{A} is the resulting Löwdin net charge on monomer A in the dimer calculation and is the formal charge on monomer A in the monomer calculation (along with the corresponding definitions for monomer B). This resulted in the removal of 24 interaction pairs, all from the groups of organic salts and organic salt solvates.
The final set consists of 221 neutral organic pairs, 276 neutral closedshell organometallic/metal–organic pairs, 341 pairs from organic salts and 319 organic salt solvates; a full listing of the crystal structures used is available in Table S1 in the supporting information. As in our previous work, it is our belief that the training set is well balanced, incorporating a wide range of typical atom environments and interaction types, and is robust enough that removal or addition of a few structures will have a minimal effect on the outcomes.
2.2. Fitting procedure
The energy model utilized throughout this work can be expressed as follows:
where the subscripts tot, coul, rep, exch, pol and disp indicate total, Coulomb, repulsion, exchange, polarization and dispersion components of the energy, respectively. Additional details of the evaluation of these terms are provided in the supporting information in Section S2. In previous work, the term E_{exchrep} = E_{rep} + E_{exch} was used in place of expressing the two components separately. In that work, the E_{exchrep} term was denoted simply as E_{rep}, but the notation has been changed here to better express that it is a combination of the HF exchange term and repulsion due to MO changes, rather than a pure repulsive term. Likewise, the notation E_{coul} has been utilized here rather than E_{ele}, although they are the same term.
The previous CE energy models utilized four empirical parameters, one scale factor for each term E_{coul}, E_{exchrep}, E_{pol} and E_{disp}. Although we can simply fit the same parameters for the model(s) pursued in this work, we have instead examined the possibility of separately scaling the exchange and repulsion terms E_{exch} and E_{rep}.
Likewise, there are free parameters in the XDM model that may be fitted (see Section 2.4). As preliminary results and examination showed that fitting these XDM parameters, in addition to other scale factors, did not significantly improve the quality of the fit, we elected not to fit them, instead choosing simply to use values of a_{1} = 0.65 and a_{2} = 1.70 Å, which are relatively close to those utilized for the B86bPBE25 hybrid by Price et al. (2023). Similarly, preliminary exploration of other terms outlined in Section 2.5, with (or without) their own scale factors, demonstrated little improvement to the fit.
Since the introduction of more parameters (and terms), as described above, was not found to give a systematic improvement to the quality of the fit to our training data – a strong indication that overfitting is occurring – three possible models were evaluated:
(i) CE1p, a global singleparameter fit for the repulsion term and polarization term (i.e. k_{rep} = k_{pol}, using the same parameter regardless of wavefunction source).
(ii) CE2p, a twoparameter fit for the exchangerepulsion term and polarization term, fitting k_{exchrep} and k_{pol} separately, and using the same parameters regardless of wavefunction source.
(iii) CE5p, a fiveparameter fit with separate scale factors k_{coul}, k_{rep}, k_{exch}, k_{pol} and k_{disp}.
The main motivation for examining these particular fits were the following theoretical arguments. The longrange behaviour for E_{coul} should be exact, along with it typically being the dominant term (especially for charged systems). Likewise, the XDM dispersion term is well defined and theoretically sound, while the weakest terms theoretically are the polarization treatment (which is only exact for spherical atoms, which is not the case here) and the E_{rep} term which would be exact if the full dimer SCF were performed. However, since we are only orthonormalizing the MOs for the monomers in the dimer systems, this ought to overestimate the repulsion as we are not allowing the MOs to relax after the initial change in the dimer environment.
2.3. Monomer wavefunctions
There are two primary considerations when selecting a source of monomer wavefunctions for the proposed model energies: choice of method, be it HF, density functional theory (DFT) or otherwise, and choice of basis set. All of the options considered in this work utilize the def2SVP (Weigend & Ahlrichs, 2005) basis set. This was chosen because of its moderate cost (no diffuse functions, but retaining polarization functions) and wide support for elements across the periodic table (H to Rn), with corresponding ECPs (Peterson et al., 2003) for heavier elements.
When considering candidate computational methods, we would intuitively expect to see systematic improvement from better molecular wavefunction sources, corresponding to the generally accepted levels of theory in the broader computational chemistry context. As such, we evaluated several different wavefunction sources: HF, LDA (SVWN5) (Dirac, 1930; Bloch, 1929; Vosko et al., 1980), BLYP (Becke, 1988b; Lee et al., 1988; Miehlich et al., 1989), B3LYP (Stephens et al., 1994), ωB97x (Chai & HeadGordon, 2008) and ωB97MV (Mardirossian & HeadGordon, 2016). HF and LDA were chosen because of their ubiquity and simplicity, and because they serve as good reference points for understanding the contributions of relative errors and are well characterized. The other choices represent different `rungs' on the ladder of DFT accuracy: BLYP remains a popular GGA (generalized gradient approximation) density functional approximation, with generally good accuracy for noncovalent interactions when combined with appropriate dispersion corrections (Goerigk et al., 2017). Likewise, dispersioncorrected B3LYP (hybridGGA) has demonstrated accuracy for noncovalent interactions and is the wavefunction source used in the CEB3LYP model from our previous work (Mackenzie et al., 2017). Finally, as previously mentioned, ωB97MV has been widely demonstrated to be a reliable functional for a variety of chemical problems (most importantly noncovalent interactions) (Najibi & Goerigk, 2018), and as such it was chosen for reference calculations for the CE fitting data set. We also investigated ωB97x as a reference point, as it is more widely implemented (rangeseparated hybrid GGA), and given that the default in e.g. ORCA is to use the VV10 (Vydrov & Van Voorhis, 2010) nonlocal correlation functional as a postSCF correction only, there should be only small differences in the MOs between the two methods.
2.4. XDM dispersion implementation
The XDM dispersion model was implemented and tested for correctness against the existing implementation in the postg program (OterodelaRoza & Johnson, 2013; Kannemann & Becke, 2010). The present implementation makes use of the Becke–Roussel (Becke & Roussel, 1989; Proynov et al., 2008) density functional approximation for the exchangehole calculations, with the same sort of grids used in numerical integration for DFT approximations in the code, namely the Becke integration scheme (Becke, 1988a) using Lebedev (1976) angular quadratures and the radial quadrature scheme from Lindh et al. (2001).
For the Hirshfeld (1977) partitioning part of the model, atomic wavefunctions from Koga et al. (1993, 2000) were used for the spherically symmetric gasphase atomic electron densities. The grids used in the numerical integration were generated using the typical Becke partitioning scheme (Becke, 1988a), as the implementation could use the existing framework already present in DFT calculations.
Atomic polarizabilities utilized throughout were the same as those used in the previous CE models (Turner et al., 2014), and the resulting inmolecule polarizabilities are used internally in the XDM methodology as follows:
where α and V are the inmolecule polarizability and Hirshfeld volume, respectively, and α′ and V′ are the corresponding quantities for the free atom. This is a simple approximation that accounts for the (typically contracted) atomic volume in a molecule.
As formulated, the XDM dispersion model has two free parameters associated with the damping at close range, a_{1} and a_{2}. The energy is formulated as follows:
where i and j are atoms, c are the coefficients specific to the atom pair ij and r is the interatomic distance. Subscripts vdW and crit indicate van der Waals and critical radius, respectively. In essence, a_{1} is a scale factor for the critical radius where all three terms (1/r^{6}, 1/r^{8} and 1/r^{10}) are equal, and a_{2} is a constant in distance units.
The above expression is, of course, the total dispersion energy for a system, not the interaction energy between two molecules. The interaction energy E_{int} can itself be calculated either by finding the corresponding energy terms for both monomers and the dimer,
or, equivalently, by simply summing only over i ∈ A and j ∈ B in equation (4).
One consideration here, given that we are not performing an SCF for the dimer AB, is in the choice of XDM coefficients – we can either reuse the coefficients calculated for the monomers A and B (referred to as the monomer dispersion model), or calculate the coefficients for the concatenated dimer wavefunction and use equation (8) (referred to as the dimer dispersion model in this work).
2.5. Other considerations
As briefly mentioned when discussing the fitting procedure, we explored various possible lowcost enhancements and modifications for the accuracy of the energy models in this work. Specifically, we evaluated the possible improvement in performance when fitting the model variants against the S66x8 (Řezáč et al., 2011a,b; Brauer et al., 2016) benchmark set, separately incorporating the geometric counterpoise correction (GCP) term (Kruse & Grimme, 2012) including the accompanying small basisset correction utilized in the HF3c model (Sure & Grimme, 2013), fitting both the XDM a_{1} and a_{2} parameters, using the larger def2TZVP basis set, using a DFT interaction model (see Section S2.1) and fitting a force fieldlike repulsion term similar to that found in the DREIDING model (Mayo et al., 1990) to correct interaction energies at short distances.
Preliminary results demonstrated that none of the above modifications systematically improved the performance of the model energies against their training set when compared with simply refitting scale factors in either the CEB3LYP or CE5p models fitted to the same data. In other words, the only improvements were attributed to the fitting of the parameters to the set being evaluated, rather than the introduction of the different terms. As such, none of these modifications were incorporated into the models proposed in this work.
2.6. Rotation of molecular orbitals using pure spherical basis sets
Since the model shown in this work, and the previous CE models, start with isolated monomer wavefunctions, it is expedient to reuse the same wavefunction (after rotation and translation to the new position) rather than recalculate each monomer multiple times, particularly in the context of molecular crystals where there is typically only one (or a handful of) symmetryunique molecules.
The above rotation was previously implemented only for wavefunctions in a Cartesian spherical harmonic basisset convention, but has now been extended to direct rotation of `pure' spherical harmonic basis sets. This has been implemented using the already known recurrence relations (Ivanic & Ruedenberg, 1998). Facilitating the use of pure spherical basis sets is valuable as it allows for the use of other QM programs (such as ORCA) which may only support the use of pure spherical basis sets, and may help to reduce computation times where many higher angular momentum (d, f, g etc.) functions are included.
2.7. Reference lattice energies for molecular crystals
In previous work we highlighted the surprising success and accuracy of the CEB3LYP model in predicting lattice energies of molecular crystals when using direct summation. To compare the models developed in this work against accurate reference data, along with the previous model, we now examine the performance on the X23 (Reilly & Tkatchenko, 2013) benchmark set, specifically using the revised reference energies in the X23b set (Dolgonos et al., 2019).
We have examined the predicted lattice energies for both the experimental geometry (after normalizing X—H bond lengths to average crystallographic values) and the PBE0+MBD/light optimized geometries given by Dolgonos et al. (2019).
3. Results and discussion
3.1. Fitting results
The previous model fitted four parameters, the scale factors k_{ele}, k_{exchrep}, k_{pol} and k_{disp}. In this work we examine the CE5p model which separately fits five parameters (splitting the exchrep term, resulting in separate k_{exch} and k_{rep} parameters), the CE2p model with two parameters (the scale factors k_{exchrep} and k_{pol}) and the singleparameter model CE1p with only k = k_{rep} = k_{pol}. In part, the last two fits are motivated by the use of a new (and more accurate) dispersion term, and it was thought that fixing the corresponding scale factors for this and the Coulomb term to unity (i.e. unscaled) would be a reasonable starting point and might result in a more transferable fit to other methods and/or basis sets.
To examine transferability, once all energy terms were evaluated for the training set (see Section S2 for full details), rather than simply performing a leastsquares fit for the CE1p and CE2p models to minimize the errors, we examined the error over a range of k in the case of CE1p, and a distribution of k_{exchrep} and k_{pol} in the case of CE2p. The results for this process are shown in Fig. 1 for the case of CE1p and in Fig. S2 for the case of the CE2p model.
In the case of CE1p, there is substantial agreement between different wavefunction sources for the optimal value of the fitted k parameter, with a minimum in the case of ωB97MV at 0.77850434…, i.e. 0.78 (to 2 d.p.). The only method with a significantly shifted minimum value is HF, but even there the difference between the RMSD of 5.0 kJ mol^{−1} for the optimal value (k = 0.86) versus 5.7 kJ mol^{−1} for the HF model is less than 1 kJ mol^{−1}.
For the case of CE2p, all methods yield minima in a similar area and there is a region of the twoparameter space where the error is relatively flat (plots are provided in Fig. S2). Given that, intuitively, ωB97MV should be the most accurate method against the reference (being the same method as used in the reference), we elected to use two parameters from the minima of this approach, located at k_{rep} = 0.485 and k_{pol} = 0.803. However, it must be pointed out that the agreement between the twodimensional minima in CE2p is much worse than in the onedimensional case of CE1p, which indicates that its transferability to different molecular wavefunction sources is likely to be inferior.
The fiveparameter model (CE5p) was also examined, with scale factors fitted using least squares, and the resulting coefficients are given in Table 1. For the fiveparameter fit, values for all scale coefficients were fixed in the range [0, 3], with the exception of k_{rep}, which was bounded in the range [0.6, 3]. These boundaries were introduced because of the potential for overfitting – it was found that freely optimizing the parameters introduced negative scale factors or scale factors close to zero, with very small improvements to the errors for the model. Based on the exploration of the parameter space for the CE1p model, k_{rep} ≥ 0.6 was chosen to allow enough flexibility to alter the k_{rep} parameter while constraining the fit to more (theoretically) reasonable parameters. It can be seen here that the Coulomb scaling coefficient is very close to unity for all wavefunction sources, which is rationalized by the number of charged systems in the training set, where this will be the overwhelmingly dominant term for most of them.

Error distributions for the training set are shown in Fig. 2, where it is notable that, for all sources other than HF, the fiveparameter fit produces only negligible improvements in overall energy errors relative to the transferable one and twoparameter fits. This, coupled with the need to introduce constraints, is indicative that it is likely we are in the realm of overfitting. To examine this further, we must look outside the training set to evaluate the behaviour against external reference data.
3.2. S66x8 set
The primary use case for the proposed energy models, and indeed the major use case for the previous models, is intermolecular interactions in organic molecular crystals, i.e. intermolecular interactions between neutral organic species. To this end, the S66x8 benchmark set (Řezáč et al., 2011a,b; Brauer et al., 2016) constitutes an excellent test set to evaluate the performance of the fitted models over a range of intermolecular interaction distances and understand its behaviour for different interaction types. The overall results are shown in Table 2. Unsurprisingly, the previously fitted model (CEB3LYP) performs very well for neutral systems, just as it did on the training set and in previous work. Likewise, across all fits, the B3LYP, ωB97X and ωB97MV wavefunction sources are significantly better choices than HF, LDA and BLYP. This is not a surprising result, but it is consistent with intuition when it comes to the superior accuracy of (hybrid) DFT functionals.

Perhaps the most important conclusion from the results for this test set is that the fiveparameter fits are only marginally better in performance (or sometimes worse) than the twoparameter fits. The performances for ωB97x and ωB97MV are almost identical, but this is unsurprising when we consider how similar the models are in terms of their wavefunctions (the VV10 nonlocal correlation functional has only been applied as a postSCF correction).
Similar to the results for the training set in the previous section, in Fig. 3 it is apparent that using more parameters does not necessarily lead to wholesale improvement of the reliability of the model, and where it does the improvement is relatively minor when considering potential applications to noncovalent interactions. The best result for CE5p, for example, was actually HF with an RMSD of 2.5 kJ mol^{−1}, versus the best result for CE1p being ωB97X with an RMSD of 3.2 kJ mol^{−1}. This is not totally insignificant, but when put in the context of the training data and the kinds of interactions that are present in the S66x8 data set, and the difference between five parameters (that may be less transferable) and only a single parameter, we believe the trade off of around 0.5 kJ mol^{−1} (which largely corresponds to a shift in mean error) is worthwhile.
With few exceptions (BLYP and HF in some cases), the error distributions are largely symmetric about zero, though there are longer and denser tails in the overbinding direction, which may be explained through the trends at different intermolecular separations in the S66x8 set. Once again, in this context the B3LYP, ωB97MV and ωB97x methods are clearly superior in their reliability, particularly for the more transferable CE1p model.
It is worth comparing our results with the recently published ωB97X3c composite method (Müller et al., 2023), which has also been evaluated against the S66 set (i.e. the subset of the S66x8 set with separation scale factor r = 1.0). Müller et al. (2023) gave ωB97X3c an MAD of roughly 1.2 kJ mol^{−1} over these interactions, compared to, say, the corresponding ωB97X variant of the CE1p model at roughly 1.8 kJ mol^{−1} (the corresponding value for CEB3LYP is 2.5 kJ mol^{−1}). Likewise, natural comparison methods for the approach are the SAPT(DFT) and SAPT0 methods, for which the corresponding MAD values, taken from Xie et al. (2022), are 1.4 kJ mol^{−1} for SAPT(DFT) and 4.1 kJ mol^{−1} for SAPT0. A full list of the MAD values for the S66 set is provided in Table S4. We believe this represents excellent agreement considering the relative cost of the methods, with our model not requiring an SCF for the combined dimer wavefunction.
3.3. Trends in error versus separation distance
In the context of molecular crystals, it is of course important for dimer interaction energies to be accurate over a range of intermolecular separations. Furthermore, characterization of the kinds of errors and their trends with distance and interaction type gives insight into where models are accurate and where their accuracy may be limited, be it systematic or otherwise.
Careful examination of Fig. 4 shows that the distribution of errors surrounding the mean error when using the CE5p fit is marginally narrower, but it appears that the CE5p fit is indeed overfitted, as its systematic overbinding at closerthanequilibrium intermolecular separations (0.9, 0.95 in the separation scale) is more extreme than for both the previous CEB3LYP model and the CE1p model. This trend, coupled with the marginal difference in overall performance for the CE1p fit versus the CE5p fit (3.3 kJ mol^{−1} for CE1p using B3LYP, versus 2.5 kJ mol^{−1} RMSD for the best CE5p model which used HF) is, in our judgement, enough to recommend usage of the CE1p model, with only a single free parameter that is clearly transferable across wavefunction sources.
Errors at closerthanequilibrium separations are of particular relevance to highpressure studies of molecular crystals, and previous work (Eikeland et al., 2017) has shown that the CEB3LYP model demonstrated significant errors compared to counterpoisecorrected B3LYP calculations on dimers in hydroquinone of methanol and acetonitrile. This was at least partly attributed to the k_{exchrep} scaling coefficient affecting repulsion at close separations. To examine this behaviour with the now larger repulsion scaling parameter for the CE1p model relative to the previous CEB3LYP model (k = 0.78 versus k_{exchrep} = 0.6177), we evaluated dimer interaction energies for the same systems as Eikeland et al. (2017) when using ωB97MV/def2svp as the wavefunction source. The results of these calculations, comparing CEB3LYP, CE1p and GFN2xTB, are summarized in Figs. S3 and S4. The results are mixed: some systems show improvements over the previous CEB3LYP model, others do not. The CE1p model cannot be said to be an improvement in this regard over the old model, but neither can it be said to be worse. This is indicative that, despite the previous discussion by Eikeland et al. (2017), the errors are not solely associated with the repulsion scale factor, as there is a significant interplay between the dispersion damping coefficients as well. Ultimately, if we wish to have an improved model at closer intermolecular separations then, at a minimum, the training data set would need to be reexamined, along with changes probably being required in the assumptions of the model (e.g. a departure from global scaling parameters), which is outside the scope of this work. However, it is worth emphasizing that, as interactions get stronger at closer separations, we should expect that the errors get larger for almost any method, and that the present model, where we are assuming relatively little deformation in the molecular electron density, will eventually break down at close enough separations where significant overlap between monomer wavefunctions becomes a reality.
4. Lattice energies for molecular crystals
The previous CE models, particularly CEB3LYP, have already been successfully applied to the prediction of lattice energies for neutral molecular crystals via the direct summation technique. Prediction of such lattice energies is particularly challenging due to the interplay of intermolecular interactions over a variety of distances where systematic errors will tend to manifest themselves largely as wrong absolute values for the lattice energies. At present, the X23 benchmark set (Reilly & Tkatchenko, 2013; Dolgonos et al., 2019), itself an extension of the C21 set (OterodelaRoza & Johnson, 2012a), constitutes the most reliable and robust set of reference values for molecular energies.
Table 3 shows the resulting error statistics in lattice energies for the models explored in this work compared with those for the CEB3LYP model. Optimized geometries that were previously used in the X23 revised benchmark (Dolgonos et al., 2019) were used, as well as experimental geometries with normalized hydrogenbond lengths; the results for the experimental geometries are available in Table S7. When interpreting these data, it is worth noting that the standard error (σ) for enthalpies is typically in the region of 4 kJ mol^{−1} (Chickos, 2003) and the correction procedure involving vibrational contributations to the will further introduce its own errors. It is also worth highlighting that six of the molecules in the X23 set are present in the training set, namely urea (UREAXX), benzene (BENZEN01), imidazole (IMAZOL13), formamide (FORMAM02), succinic acid (SUCACB09) and uracil (URACIL). However, the training set is based only on the nearest dimers (rather than the dozens or hundreds of energies contributing to the lattice energy) and uses a different reference method than in the reference X23 data.

It can be seen in Table 3 that the errors for different models are nowhere near as uniform as they were for the S66x8 benchmark set. This should not be surprising – any systematic over or underestimation for interactions will be amplified when summing across many interactions/dimer pairs. The transferable singleparameter model demonstrates excellent performance and transferability across different wavefunction methods, though the success is significantly diminished when examining the HF results.
There is little doubt that the accuracy of the models presented here, when using direct summation for the lattice energies, is reliant on cancellation of errors (though the same may be said for all other methods). Nevertheless, the results for the CE1p model using B3LYP, with an MAD of only 3.6 kJ mol^{−1}, compares very favourably with the current state of the art (Price et al., 2023), where the best method (a composite of B86bPBE25XDM//B86bPBEXDM approaches including a basisset correction) reported had an MAD of 2 kJ mol^{−1}. For comparison, the same work reported PBE0MBD giving an MAD of 3.5 kJ mol^{−1} when using the `Tight' basis set.
The relative behaviour of the CE1p, CE2p and CE5p models when using the B3LYP functional against CEB3LYP may be seen in Fig. 5, which indicates that the singleparameter model is an improvement over the alternatives across almost every system. The results for cytosine indicate that there may be some systematic error in the models, which merits further investigation into polarization as discussed in the next subsection.
4.1. Crystal polarization effects
There are several sources of possible error in the direct summation method used here for lattice energies, with one of the more significant sources of error being the treatment of polarization. In a crystalline environment, the framework utilized in the CE energy models assumes that the polarization contributions are pairwise additive. However, the true polarization energy is fundamentally manybody in nature because of the quadratic dependence on magnitude of the electric field:
Here i is the ith atom, α_{i} is its polarizability, is the polarization energy for this atom and is the electric field at atom i from its environment, either from the neighbouring monomer (in the case of a pair) or from the crystal environment. While the field is just the pairwise sum over neighbours, the square of this quantity is no longer separable.
The practical difference in the case of crystals may be seen in Fig. 6, where the sum of pairwise polarization terms overestimates the polarization energy of the crystal, especially for cytosine. In contrast, summing the field experienced by a monomer over all interactions in the crystal, before computing the polarization energy, significantly improves most energies for the X23 data set. It should be noted that there is a further effect of using gasphase wavefunctions that may lead to underestimation of the polarization in comparison to the corresponding wavefunction in the solid state. Since gasphase wavefunctions are used throughout the current procedure, we are likewise underestimating the polarization of the initial molecular wavefunction. These two factors will compete and we may (or may not) experience cancellation of errors, but this is an aspect that must be understood if, in future, moleculeincrystal wavefunctions, for example, were to be used with this energy model to estimate lattice energies.
For now, we can examine the influence of evaluating the polarization energy using the theoretically more correct approach of summing over contributions to the field before computing the energy. This manifests in changes in the error statistics: the MAD value shifts to 4.4 kJ mol^{−1} (versus 3.6 kJ mol^{−1}), largely due to the systematic shift in MSD to 2.2 kJ mol^{−1} (versus −1.0 kJ mol^{−1}), while the RMSD only shifts to 5.3 kJ mol^{−1} (versus 4.8 kJ mol^{−1}).
This suggests that if the goal is to predict relative lattice energies, as it so often is in the context of molecular crystals, then using a polarization term based on the electric field of the crystal should be considered. Indeed, there is a minor increase in error (less than 1 kJ mol^{−1}) for the absolute values of the lattice energies compared with the reference values, but the errors are likewise more systematic which tends to improve the relative energies. Furthermore, in ionic systems the effects of the pairwise approximation for polarization are even greater. An illustrative example is cubic NaCl, which has zero net electric field at the ion positions, but when calculated pairwise there is an erroneous net polarization energy of around 400 kJ mol^{−1}.
5. Conclusions
We have implemented, fitted and demonstrated the accuracy of a new singleparameter intermolecular interaction energy model, CE1p, which is transferable across different wavefunction sources. The model constitutes an excellent accuracy/computational cost trade off, with an MAD of 2.0 kJ mol^{−1} across the S66x8 benchmark set, and near stateoftheart performance when predicting lattice energies for the X23 benchmark set, with an MAD of 3.6 kJ mol^{−1}, while using the simplest pairwise sum in real space. While there is little doubt that this model (and others) exhibits error cancellation, we believe it constitutes an efficient and accurate method for quantitative predictions of intermolecular interactions which may also be rationalized and understood through sensible separation of the energy into distinct contributions.
We have also examined the performance of the GFN2xTB tightbinding model as a rapid method for similar purposes, with only relatively minor loss in accuracy for dimer energies in organic neutral systems. In particular, for neutral organic molecules (which constitute a large part of the chemical systems of interest in molecular crystals) using GFN2xTB through visualization software such as CrystalExplorer provides nearinstant feedback and allows for realtime exploratory evaluations, as in most cases it is one to two orders of magnitude faster than our energy model.
Our recommendations for the end user of these methods, who wishes to rationalize the intermolecular interactions in a molecular crystal, are as follows:
(i) First, calculate the interaction energies using GFN2xTB, which is extremely low cost for pairwise interactions. It is particularly worth examining interactions beyond nearestneighbours using such methods, as their longrange behaviour should be reliable. These can be visualized using energy frameworks (Turner et al., 2014).
(ii) Calculate the same energies using the CE1p model, using B3LYP, ωB97x or ωB97MV functionals, depending on which is available in your QM software of choice.
(iii) If these models are not in agreement with respect to the relative strength of interactions, further detailed investigation with higherlevel methods is warranted.
(iv) If lattice energies are of interest, strongly consider the polarity of the molecule(s) in your system, and whether or not the pairwise approximation for polarization is appropriate.
There are several avenues for further improvements and developments, including:
(i) The possibility of fitting a model using DFT exchange correlation rather than HF exchange, given the reduction in computational cost over evaluating the exact exchange matrices for the two monomers (A and B) and the better scaling speedups available for the Coulomb term, particularly if the overall accuracy of the model is not sacrificed.
(ii) Departing from global scaling parameters, particularly at close intermolecular separations, though this would be a significant change to the model.
(iii) Examination of the connection between the orthogonalized electron population (see Section S9) and the repulsion and exchange terms for development of lowcost approximations.
In summary, the new singleparameter CE1p model is generally more accurate, is significantly more transferable (not only to other wavefunction methods but in its coverage of more of the periodic table, which enables application to more chemical problems) and has (qualitatively) a more sound theoretical basis than the existing CEB3LYP model. Hence, we view it as an excellent replacement for the previous method which will no doubt find broad application across a variety of chemical problems, particularly in the domain of molecular crystals.
This work has been implemented and is already available in the open source software occ made available on GitHub (https://github.com/peterspackman/occ) The incorporation of the newly proposed CE1p model, along with access to GFN2xTB, will be available in a forthcoming release of CrystalExplorer and we are optimistic that both will be valuable additions for the community.
6. Related literature
For further literature related to the supporting information, see van Eijck & Kroon (1997) and Thomas et al. (2018).
Supporting information
Additional details, figures and tables. DOI: https://doi.org/10.1107/S2052252523008941/fc5074sup1.pdf
Compressed archive of supporting data files. DOI: https://doi.org/10.1107/S2052252523008941/fc5074sup2.zip
Acknowledgements
PRS and JDG thank the Pawsey Supercomputing Centre and National Computational Infrastructure for provision of computing resources.
Funding information
The following funding is acknowledged: Australian Research Council (grant No. FL180100087).
References
Bannwarth, C., Ehlert, S. & Grimme, S. (2019). J. Chem. Theory Comput. 15, 1652–1671. Web of Science CrossRef CAS PubMed Google Scholar
Becke, A. & Roussel, M. (1989). Phys. Rev. A, 39, 3761–3767. CrossRef CAS Google Scholar
Becke, A. D. (1988a). J. Chem. Phys. 88, 2547–2553. CrossRef CAS Web of Science Google Scholar
Becke, A. D. (1988b). Phys. Rev. A, 38, 3098–3100. CrossRef CAS PubMed Web of Science Google Scholar
Becke, A. D. & Johnson, E. R. (2007). J. Chem. Phys. 127, 154108. Web of Science CrossRef PubMed Google Scholar
Bloch, F. (1929). Z. Phys. 57, 545–555. CrossRef CAS Google Scholar
Boys, S. & Bernardi, F. (1970). Mol. Phys. 19, 553–566. CrossRef CAS Web of Science Google Scholar
Brammer, L. (2017). Faraday Discuss. 203, 485–507. Web of Science CrossRef CAS PubMed Google Scholar
Brauer, B., Kesharwani, M. K., Kozuch, S. & Martin, J. M. L. (2016). Phys. Chem. Chem. Phys. 18, 20905–20925. Web of Science CrossRef CAS PubMed Google Scholar
Chai, J.D. & HeadGordon, M. (2008). J. Chem. Phys. 128, 084106. Web of Science CrossRef PubMed Google Scholar
Chickos, J. S. (2003). Netsu Sokutei, 30(3), 116–124. Google Scholar
Dirac, P. A. M. (1930). Math. Proc. Camb. Phil. Soc. 26, 376–385. CrossRef CAS Google Scholar
Dolgonos, G. A., Hoja, J. & Boese, A. D. (2019). Phys. Chem. Chem. Phys. 21, 24333–24344. CrossRef CAS PubMed Google Scholar
Eijck, B. P. van & Kroon, J. (1997). J. Phys. Chem. B, 101, 1096–1100. Google Scholar
Eikeland, E., Thomsen, M. K., Overgaard, J., Spackman, M. A. & Iversen, B. B. (2017). Cryst. Growth Des. 17, 3834–3846. Web of Science CSD CrossRef CAS Google Scholar
Goerigk, L., Hansen, A., Bauer, C., Ehrlich, S., Najibi, A. & Grimme, S. (2017). Phys. Chem. Chem. Phys. 19, 32184–32215. Web of Science CrossRef CAS PubMed Google Scholar
Grimme, S. (2006). J. Comput. Chem. 27, 1787–1799. Web of Science CrossRef PubMed CAS Google Scholar
Groom, C. R., Bruno, I. J., Lightfoot, M. P. & Ward, S. C. (2016). Acta Cryst. B72, 171–179. Web of Science CrossRef IUCr Journals Google Scholar
Hirshfeld, F. L. (1977). Theor. Chim. Acta, 44, 129–138. CrossRef CAS Web of Science Google Scholar
Ivanic, J. & Ruedenberg, K. (1998). J. Phys. Chem. A, 102, 9099–9100. CrossRef CAS Google Scholar
Johnson, E. R. & Becke, A. D. (2006). J. Chem. Phys. 124, 174104. Web of Science CrossRef PubMed Google Scholar
Kannemann, F. O. & Becke, A. D. (2010). J. Chem. Theory Comput. 6, 1081–1088. CrossRef CAS Google Scholar
Koga, T., Kanayama, K., Watanabe, T., Imai, T. & Thakkar, A. J. (2000). Theor. Chem. Acc. 104, 411–413. Web of Science CrossRef CAS Google Scholar
Koga, T., Tatewaki, H. & Thakkar, A. J. (1993). Phys. Rev. A, 47, 4510–4512. CrossRef CAS PubMed Web of Science Google Scholar
Kruse, H. & Grimme, S. (2012). J. Chem. Phys. 136, 154101. Web of Science CrossRef PubMed Google Scholar
Lebedev, V. I. (1976). USSR Comput. Math. Math. Phys. 16, 10–24. CrossRef Google Scholar
Lee, C., Yang, W. & Parr, R. G. (1988). Phys. Rev. B, 37, 785–789. CrossRef CAS Web of Science Google Scholar
Levin, I. (2020). NIST Inorganic Crystal Structure Database (ICSD). National Institute of Standards and Technology, Maryland, USA. https://data.nist.gov/od/id/mds22147. Google Scholar
Lindh, R., Malmqvist, P.Å & Gagliardi, L. (2001). Theor. Chem. Acc. 106, 178–187. Google Scholar
Mackenzie, C. F., Spackman, P. R., Jayatilaka, D. & Spackman, M. A. (2017). IUCrJ, 4, 575–587. Web of Science CrossRef CAS PubMed IUCr Journals Google Scholar
Mardirossian, N. & HeadGordon, M. (2016). J. Chem. Phys. 144, 214110. CrossRef PubMed Google Scholar
Mayo, S. L., Olafson, B. D. & Goddard, W. A. (1990). J. Phys. Chem. 94, 8897–8909. CrossRef CAS Web of Science Google Scholar
Miehlich, B., Savin, A., Stoll, H. & Preuss, H. (1989). Chem. Phys. Lett. 157, 200–206. CrossRef CAS Web of Science Google Scholar
Müller, M., Hansen, A. & Grimme, S. (2023). J. Chem. Phys. 158, 014103. PubMed Google Scholar
Najibi, A. & Goerigk, L. (2018). J. Chem. Theory Comput. 14, 5725–5738. CrossRef CAS PubMed Google Scholar
Neese, F. (2022). WIREs Comput. Mol. Sci. 12, e1606. Google Scholar
OterodelaRoza, A. & Johnson, E. R. (2012a). J. Chem. Phys. 137, 054103. PubMed Google Scholar
OterodelaRoza, A. & Johnson, E. R. (2012b). J. Chem. Phys. 136, 174109. PubMed Google Scholar
OterodelaRoza, A. & Johnson, E. R. (2013). J. Chem. Phys. 138, 204109. Web of Science PubMed Google Scholar
Peterson, K. A., Figgen, D., Goll, E., Stoll, H. & Dolg, M. (2003). J. Chem. Phys. 119, 11113–11123. Web of Science CrossRef CAS Google Scholar
Price, A. J. A., OterodelaRoza, A. & Johnson, E. R. (2023). Chem. Sci. 14, 1252–1262. CrossRef CAS PubMed Google Scholar
Proynov, E., Gan, Z. & Kong, J. (2008). Chem. Phys. Lett. 455, 103–109. CrossRef PubMed CAS Google Scholar
Reilly, A. M. & Tkatchenko, A. (2013). J. Chem. Phys. 139, 024705. Web of Science CrossRef PubMed Google Scholar
Řezáč, J., Riley, K. E. & Hobza, P. (2011a). J. Chem. Theory Comput. 7, 3466–3470. Google Scholar
Řezáč, J., Riley, K. E. & Hobza, P. (2011b). J. Chem. Theory Comput. 7, 2427–2438. Web of Science PubMed Google Scholar
Spackman, P. R., Turner, M. J., McKinnon, J. J., Wolff, S. K., Grimwood, D. J., Jayatilaka, D. & Spackman, M. A. (2021). J. Appl. Cryst. 54, 1006–1011. Web of Science CrossRef CAS IUCr Journals Google Scholar
Spackman, P. R., Walisinghe, A. J., Anderson, M. W. & Gale, J. D. (2023). Chem. Sci. 14, 7192–7207. CrossRef CAS PubMed Google Scholar
Stephens, P. J., Devlin, F. J., Chabalowski, C. F. & Frisch, M. J. (1994). J. Phys. Chem. 98, 11623–11627. CrossRef CAS Web of Science Google Scholar
Su, P. & Li, H. (2009). J. Chem. Phys. 131, 014102. Web of Science CrossRef PubMed Google Scholar
Sure, R. & Grimme, S. (2013). J. Comput. Chem. 34, 1672–1685. Web of Science CrossRef CAS PubMed Google Scholar
Tan, S. L., Jotani, M. M. & Tiekink, E. R. T. (2019). Acta Cryst. E75, 308–318. Web of Science CrossRef IUCr Journals Google Scholar
Thomas, S. P., Spackman, P. R., Jayatilaka, D. & Spackman, M. A. (2018). J. Chem. Theory Comput. 14, 1614–1623. Web of Science CrossRef CAS PubMed Google Scholar
Turner, M. J., Grabowsky, S., Jayatilaka, D. & Spackman, M. A. (2014). J. Phys. Chem. Lett. 5, 4249–4255. Web of Science CrossRef CAS PubMed Google Scholar
Vosko, S. H., Wilk, L. & Nusair, M. (1980). Can. J. Phys. 58, 1200–1211. Web of Science CrossRef CAS Google Scholar
Vydrov, O. A. & Van Voorhis, T. (2010). J. Chem. Phys. 133, 244103. CrossRef PubMed Google Scholar
Wang, C. & Sun, C. C. (2018). Cryst. Growth Des. 18, 1909–1916. Web of Science CrossRef CAS Google Scholar
Weigend, F. & Ahlrichs, R. (2005). Phys. Chem. Chem. Phys. 7, 3297–3305. Web of Science CrossRef PubMed CAS Google Scholar
Weigend, F., Furche, F. & Ahlrichs, R. (2003). J. Chem. Phys. 119, 12753–12762. Web of Science CrossRef CAS Google Scholar
Xie, Y., Smith, D. G. A. & Sherrill, C. D. (2022). J. Chem. Phys. 157, 024801. CrossRef PubMed Google Scholar
This is an openaccess article distributed under the terms of the Creative Commons Attribution (CCBY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.