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 crystal structure, molecule/ion pairs were obtained by generating a cluster of nearest-neighbour interactions around each symmetry-unique molecule/ion in the crystal. A complete list of the crystal structures, including Cambridge Structural Database (CSD) refcodes (Groom et al., 2016) and Inorganic Crystal Structure Database (ICSD) identifiers (Levin, 2020), is available in the supporting information.

Reference ωB97M-V/def2-QZVP interaction energies were then calculated via evaluation of the full dimer SCF energy, subtracting the monomer energies (in the monomer basis). The ωB97M-V method (Mardirossian & Head-Gordon, 2016) was selected because of its general accuracy and transferability across many problems, including non-covalent interactions (Najibi & Goerigk, 2018). The large def2-QZVP (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 gas-phase dimer does not necessarily correspond well to that present when the dimer is embedded in a crystalline environment with competing charge-transfer 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 closed-shell 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_exch-rep = E_rep + E_exch was used in place of expressing the two components separately. In that work, the E_exch-rep 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_exch-rep, 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₁ = 0.65 and a₂ = 1.70 Å, which are relatively close to those utilized for the B86bPBE-25 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) CE-1p, a global single-parameter fit for the repulsion term and polarization term (i.e. k_rep = k_pol, using the same parameter regardless of wavefunction source).

(ii) CE-2p, a two-parameter fit for the exchange-repulsion term and polarization term, fitting k_exch-rep and k_pol separately, and using the same parameters regardless of wavefunction source.

(iii) CE-5p, a five-parameter 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 long-range 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 def2-SVP (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 & Head-Gordon, 2008) and ωB97M-V (Mardirossian & Head-Gordon, 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 non-covalent interactions when combined with appropriate dispersion corrections (Goerigk et al., 2017). Likewise, dispersion-corrected B3LYP (hybrid-GGA) has demonstrated accuracy for non-covalent interactions and is the wavefunction source used in the CE-B3LYP model from our previous work (Mackenzie et al., 2017). Finally, as previously mentioned, ωB97M-V has been widely demonstrated to be a reliable functional for a variety of chemical problems (most importantly non-covalent inter­actions) (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 (range-separated hybrid GGA), and given that the default in e.g. ORCA is to use the VV10 (Vydrov & Van Voorhis, 2010) non-local correlation functional as a post-SCF 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 (Otero-de-la-Roza & 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 exchange-hole 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 gas-phase 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 in-molecule polarizabilities are used internally in the XDM methodology as follows:

where α and V are the in-molecule 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₁ and a₂. 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₁ is a scale factor for the critical radius where all three terms (1/r⁶, 1/r⁸ and 1/r¹⁰) are equal, and a₂ 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 re-use 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 low-cost 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 basis-set correction utilized in the HF-3c model (Sure & Grimme, 2013), fitting both the XDM a₁ and a₂ parameters, using the larger def2-TZVP basis set, using a DFT interaction model (see Section S2.1) and fitting a force field-like 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 re-fitting scale factors in either the CE-B3LYP or CE-5p 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 re-use the same wavefunction (after rotation and translation to the new position) rather than recalculate each monomer multiple times, particularly in the context of mol­ecular crystals where there is typically only one (or a handful of) symmetry-unique molecules.

The above rotation was previously implemented only for wavefunctions in a Cartesian spherical harmonic basis-set 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 CE-B3LYP 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.1. Fitting results

The previous model fitted four parameters, the scale factors k_ele, k_exch-rep, k_pol and k_disp. In this work we examine the CE-5p model which separately fits five parameters (splitting the exch-rep term, resulting in separate k_exch and k_rep parameters), the CE-2p model with two parameters (the scale factors k_exch-rep and k_pol) and the single-parameter model CE-1p 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 least-squares fit for the CE-1p and CE-2p models to minimize the errors, we examined the error over a range of k in the case of CE-1p, and a distribution of k_exch-rep and k_pol in the case of CE-2p. The results for this process are shown in Fig. 1 for the case of CE-1p and in Fig. S2 for the case of the CE-2p model.

Figure 1 Error distributions for HF, LDA, BLYP, B3LYP, ωB97x and ωB97M-V methods over a range of k = krep = kpol values. While there are clearly some differences between the models, it should be seen that the minimum for ωB97M-V at k = 0.78 (highlighted by the grey dot-dashed line) still produces good results for the other models.

In the case of CE-1p, there is substantial agreement between different wavefunction sources for the optimal value of the fitted k parameter, with a minimum in the case of ωB97M-V 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⁻¹ for the optimal value (k = 0.86) versus 5.7 kJ mol⁻¹ for the HF model is less than 1 kJ mol⁻¹.

For the case of CE-2p, all methods yield minima in a similar area and there is a region of the two-parameter space where the error is relatively flat (plots are provided in Fig. S2). Given that, intuitively, ωB97M-V 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 two-dimensional minima in CE-2p is much worse than in the one-dimensional case of CE-1p, which indicates that its transferability to different molecular wavefunction sources is likely to be inferior.

The five-parameter model (CE-5p) was also examined, with scale factors fitted using least squares, and the resulting coefficients are given in Table 1. For the five-parameter 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 over-fitting – 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 CE-1p 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 five-parameter fit produces only negligible improvements in overall energy errors relative to the transferable one and two-parameter 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.

Figure 2 Kernel density estimate plots, showing the errors for the predicted energies versus the reference (ωB97M-V/def2-QZVP) dimer interaction energies in the training set when using the previous CE-B3LYP model, and GFN2-xTB, along with the fitted CE-1p, CE-2p and CE-5p using the six wavefunction sources explored in this work (top to bottom): HF, LDA (SVWN5), BLYP, B3LYP, ωB97X, ωB97M, with RMSD values provided to the right-hand side of each plot. Each wavefunction source has been shifted on the y axis to facilitate visual inspection of the distributions.

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 (CE-B3LYP) 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 ωB97M-V 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 five-parameter fits are only marginally better in performance (or sometimes worse) than the two-parameter fits. The performances for ωB97x and ωB97M-V are almost identical, but this is unsurprising when we consider how similar the models are in terms of their wavefunctions (the VV10 non-local correlation functional has only been applied as a post-SCF 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 non-covalent interactions. The best result for CE-5p, for example, was actually HF with an RMSD of 2.5 kJ mol⁻¹, versus the best result for CE-1p being ωB97X with an RMSD of 3.2 kJ mol⁻¹. 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⁻¹ (which largely corresponds to a shift in mean error) is worthwhile.

Figure 3 Kernel density estimate plots for the S66x8 set using the CE-1p, CE-2p and CE-5p fits, the previous CE-B3LYP model and the six wavefunction sources used in this work (top to bottom): HF, LDA (SVWN5), BLYP, B3LYP, ωB97X, ωB97M-V, with RMSD values provided on the right-hand side of each plot. Each wavefunction source has been shifted on the y axis to facilitate visual inspection of the distributions.

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 over-binding direction, which may be explained through the trends at different intermolecular separations in the S66x8 set. Once again, in this context the B3LYP, ωB97M-V and ωB97x methods are clearly superior in their reliability, particularly for the more transferable CE-1p model.

It is worth comparing our results with the recently published ωB97X-3c 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 ωB97X-3c an MAD of roughly 1.2 kJ mol⁻¹ over these interactions, compared to, say, the corresponding ωB97X variant of the CE-1p model at roughly 1.8 kJ mol⁻¹ (the corresponding value for CE-B3LYP is 2.5 kJ mol⁻¹). 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⁻¹ for SAPT(DFT) and 4.1 kJ mol⁻¹ 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 CE-5p fit is marginally narrower, but it appears that the CE-5p fit is indeed over-fitted, as its systematic over-binding at closer-than-equilibrium intermolecular separations (0.9, 0.95 in the separation scale) is more extreme than for both the previous CE-B3LYP model and the CE-1p model. This trend, coupled with the marginal difference in overall performance for the CE-1p fit versus the CE-5p fit (3.3 kJ mol⁻¹ for CE-1p using B3LYP, versus 2.5 kJ mol⁻¹ RMSD for the best CE-5p model which used HF) is, in our judgement, enough to recommend usage of the CE-1p model, with only a single free parameter that is clearly transferable across wavefunction sources.

Figure 4 Trends in errors across different scaled separations for the S66x8 benchmark set when using different models in this work: CE-B3LYP (top), and CE-1p (upper middle), CE-2p (lower middle) and CE-5p (bottom) for the B3LYP wavefunction source. The solid lines indicate the mean error value, while the shaded regions represent the 95% confidence interval. The separation scale is relative to the equilibrium separation, i.e. a value of 1.0 is at equilibrium, while a value of 0.9 is 90% of the equlibrium separation. Values between explicitly calculated separations (0.9, 0.95, 1.0, 1.05, 1.1, 1.25, 1.5, 2.0) have been linearly interpolated in order to show the behavioural trends.

Errors at closer-than-equilibrium separations are of particular relevance to high-pressure studies of molecular crystals, and previous work (Eikeland et al., 2017) has shown that the CE-B3LYP model demonstrated significant errors compared to counterpoise-corrected B3LYP calculations on dimers in hydroquinone clathrates of methanol and acetonitrile. This was at least partly attributed to the k_exch-rep scaling coefficient affecting repulsion at close separations. To examine this behaviour with the now larger repulsion scaling parameter for the CE-1p model relative to the previous CE-B3LYP model (k = 0.78 versus k_exch-rep = 0.6177), we evaluated dimer inter­action energies for the same systems as Eikeland et al. (2017) when using ωB97M-V/def2-svp as the wavefunction source. The results of these calculations, comparing CE-B3LYP, CE-1p and GFN2-xTB, are summarized in Figs. S3 and S4. The results are mixed: some systems show improvements over the previous CE-B3LYP model, others do not. The CE-1p 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 re-examined, 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.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 many-body 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 gas-phase wavefunctions that may lead to underestimation of the polarization in comparison to the corresponding wavefunction in the solid state. Since gas-phase 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, molecule-in-crystal wavefunctions, for example, were to be used with this energy model to estimate lattice energies.

Figure 6 Errors in predicted X23 energies (difference from reference values) for the CE-1p method using B3LYP/def2-SVP as the wavefunction source when using two different polarization models. Here `Pairwise field' corresponds to the usual method of calculating the polarization term for the lattice energy, i.e. each pair is calculated individually, whereas `Crystal field' evaluates the electric field from all neighbours at the atomic positions of the symmetry-unique molecule(s) in the crystal. See equation (9).

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⁻¹ (versus 3.6 kJ mol⁻¹), largely due to the systematic shift in MSD to 2.2 kJ mol⁻¹ (versus −1.0 kJ mol⁻¹), while the RMSD only shifts to 5.3 kJ mol⁻¹ (versus 4.8 kJ mol⁻¹).

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⁻¹) 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⁻¹.