A critical investigation of proposed electrostatic corrections to quantum mechanical volumes: the importance of variation and the irrelevance of imbalance

ABSTRACT

The crystal density of neutral and ionic molecular crystals is remarkably well approximated by the enclosed volume of molecular surfaces, where these surfaces are defined as regions of constant and small electron density. Several workers have proposed that estimates may be improved if one includes quantities extracted from the electrostatic potential on the surface of the molecule. The variation of the potential and the imbalance of positive and negative values have been considered to be of importance. In this study we demonstrate that whereas variation is important for improving crystal density predictions, imbalance is not. We present a density functional theory study on a set of 44 neutral molecular crystals. Ten-fold cross-validations were performed on models that incorporate variation, imbalance and combinations of both. Geometries were optimised using B3LYP and basis sets of type 6-31G(d). Electron densities and electrostatic potentials were computed with B3LYP and M05. Regardless of functional, models that correct for variation yield a relative decrease of 15%–18% in root-mean-square error of prediction. This correction appears to sharpen the error distribution about zero. Models based on imbalance yield no improvement, and we argue that it plays an insignificant role.

KEYWORDS:

• Crystal density
• quantum mechanical volume
• electrostatic potential corrections
• molecular surface
• density functional theory

1. Introduction

The density is an important parameter for the detonation performance of a material. To a large extent, it determines the detonation velocity and detonation pressure [1], quantities of great interest in the design and development of energetic materials [2]. Theoretical prediction of the density is an attractive goal, not only because of its importance, but also because it avoids the high costs associated with synthesising and testing novel compounds.

Predicting the density of a crystal has turned out to be very difficult. One possible strategy is to make an attempt to determine the crystal structure, and subsequently deduce the density. Unfortunately, crystal structure prediction remains, in general, an unsolved problem [3,4]. A less ambitious but promising strategy consists of approximating the effective volume a molecule occupies in the crystal. This approach has been remarkably successful in the quantum mechanical setting, where estimates of the effective volume are given by enclosed volumes of surfaces of constant and small electron density [5–8]. Denoting the crystal density by ρ, we can define the molecular volume V by (1) where M is the mass of a single molecule. Here V is the effective volume of a molecule in the crystal, implying that it may differ from the intuitive measure of volume as space. In fact, V will be smaller or larger depending on factors affecting the crystal structure, for instance the molecular shape and intermolecular interactions. To a first approximation, V is simply the size of the molecule, and this can most definitely be estimated by quantum mechanics.

Bader et al. [9] proposed that the surface of molecules may be represented by regions with a small and constant electron density , defined as the probability density to find an electron at the position . They suggested to apply surfaces on which or 0.0010 a.u. The latter, (2) has mostly been settled on in the literature [5–8,10–15], although it has been remarked that other surfaces would serve equally well as the molecular surface [16–18]. Rice and coworkers [5] found that V is exceptionally well approximated by the volume V_encl enclosed by . They achieved a root-mean-square error (RMSE) within 4% of experimentally obtained densities for a set of 180 neutral molecular crystals comprised of carbon, hydrogen, nitrogen and oxygen. The molecular surface was computed by density functional theory (DFT), applying the B3LYP functional and basis sets of type 6-31G(d,p).

Although no systematic errors were found in the Rice et al. [5] study, it is clear that the arbitrary definition of as well as the particular choice of functional and basis set might lead to such errors. Politzer and coworkers [6] observed that V_encl was systematically underestimated with the B3PW91 functional. Consequently, in order to correct for potential systematic errors, we will add a constant γ to V_encl: (3) In the following, we will denote approximations of V by . It is implicitly assumed in Equation (3(3) ) that if a systematic error is present, it is independent of the volume. Another option is to scale V_encl by some constant α, i.e. , making the assumption that the error is proportional to the volume. One may also combine the two in , but we have chosen to avoid this in order to keep the number of parameters to a minimum. The constant γ is absent in previous studies [6–8], making the comparison to more complicated models unfair. Moreover, these complicated models have, in addition to proposed electrostatic corrections, both a scaling constant α and an overall constant γ. These facts naturally raise the question of whether their effectiveness is due to α and γ, or as reported, electrostatics.

Although the molecular volume V is well approximated by V_encl, Politzer and coworkers [6] pointed out that the role of intermolecular electrostatic interactions are completely ignored. Building on previous work, they proposed that information about such interactions are present in the electrostatic potential , the electric potential felt by a unit charge positioned at : (4) Restricting the potential to on , a number of studies have been conducted in which statistics from assumed to be relevant for predicting some physical or chemical property have been extracted. Such properties include solubilities [15,19], impact sensitivities [13,14], partition coefficients [20], crystal densities [6–8], chemical inhibition [16,18], critical points and boiling points [21], as well as sublimation enthalpies and solvation free energies [17]. The procedure has been conceptually generalised as (5) where Φ⁽ⁱ⁾ are parameters computed from restricted to . Applied parameters Φ⁽ⁱ⁾ in Equation (5(5) ) include the extrema of , average values, variances and other measures of spread, and areas of negative and positive regions [16]. The functional relationship f is unknown and has been termed a general interaction properties function (GIPF) [11]. While it is our view that modelling based on is a sensible approach that should be further explored, we consider much of the research based on Equation (5(5) ) to be inconclusive, because the performance of a model has often been equated with the extent to which it fits the data-set, with no consideration of the model’s complexity [11–14,16,20–22]. The complexity may be taken into account by considering information criteria, such the Akaike information criterion (AIC) or Bayesian information criterion (BIC) [23], providing more suitable quantities in the comparison of unvalidated models. Models should nevertheless be validated by testing them on unseen data-sets, for instance by test-set/training-set division, cross-validation or boot-strapping [24].

In this paper we will investigate proposed electrostatic corrections to V_encl and identify the cause of the observed improvement. As such, the work presented provides a firm foundation on which to formulate more accurate models and to interpret the reported results. We begin by describing the quantities that have been applied as corrections. Brinck and coworkers [25] introduced the ‘charge separation’ quantity Π, defined by (6) A non-zero value of Π has been interpreted as a prerequisite for intermolecular interactions [6]. In this equation, N is the number of points on the discretised surface , and is the average value of on . The quantity Π can also be identified as the mean absolute deviation from , a measure of the spread about the mean. A related and more commonly used measure of variation is the standard deviation estimate [26], (7) a quantity that we will consider as an alternative to Π.¹Politzer et al. [19] proposed another measure of spread, , in which negative and positive variances were treated separately: (8) In this expression, the mean positive/negative value on is denoted .

While the spread is a parameter that quantifies the variability of on , other aspects may be important for intermolecular interactions. It has been suggested that one such aspect is the imbalance of negative and positive values, the argument being that if significant imbalance is present, the probability for favourable intermolecular interactions diminishes [7,11,17]. Murray and coworkers [21] proposed to use the quantity (9) for this purpose, a reasonable choice if σ²₋ and σ²₊ correspond to ‘amounts’ of negative and positive values on . A perhaps simpler way to construct a parameter signifying imbalance is to make use of well-known measures of the asymmetry of distributions. To define such a quantity, we thus apply the statistical notion of the coefficient of skewness [27], (10) It is our understanding of the term ‘imbalance’ that it should correspond to the ‘absolute skewness of about the origin’. Consequently, we define an imbalance parameter Θ by first replacing by zero in s, and then computing its absolute value. We will consider Θ as a possible correction to quantum mechanical volumes. To incorporate both variation and imbalance, Murray et al. [21] combined these two aspects of in the parameter .

Politzer and coworkers [6] reasoned that the likelihood of attractive intermolecular interactions is maximised when the imbalance is close to zero (ν is large) and the variation is considerable ( is large). With this interpretation, the molecules for which is large are the ones that one would expect to be more closely packed in the solid than suggested by their volume. They therefore investigated two models for predicting the crystal density, one correcting for variation, and one correcting for both variation and imbalance: (11) (12) Fitting these equations to a set of 36 CHNO compounds, the training-set RMSE obtained was 0.053 g/cm³ for Equation (11(11) ), and 0.047 g/cm³ for Equation (12(12) ). The uncorrected V_encl yielded an RMSE of 0.068 g/cm³. In an evaluation [8], Equation (12(12) ) was validated on a training-set and test-set consisting of 180 and 38 CHNO compounds, respectively. They found a mean absolute error (MAE) of 0.035 g/cm³, compared to 0.050 g/cm³ for the uncorrected V_encl. To our knowledge, Π has never been validated. But it does deserve such a test, given that the training-set RMSE obtained by Politzer and coworkers [6] was very similar for Π and .

While it has been noted [16–18] that the particular choice of surface should be of minor importance, this claim has important mathematical consequences severely restricting the space of permissible parameters Φ used in linear models. This is the subject of a recent paper [28], in which we present a proof of the following. Consider a linear model of the form (13) For such a model to make the same predictions on a range of surfaces of electron density ε, the parameter values of Φ computed on and must be linearly related to each other. This requirement is independent of the modelled property, and hence of considerable importance for all attempts to apply molecular surface parameters in modelling. We tested the mathematical requirement in the particular case of imbalance and variation, showing that whereas variation (Φ = σ, Π) satisfies it in the electron density range [0.0001 a.u., 0.0040 a.u.], imbalance (Φ = ν, Θ) does not. This implies, in particular, that the model in Equation (16b(16b) ) will fail to make the same predictions on these surfaces. The surface is illustrated in Figure 1.

Figure 1. The electron density isosurface for compound 1 with the B3LYP functional. The illustration was produced with the GaussView 4.1 software [29].

2. Procedure

We study a set of 44 neutral molecular crystals. Among these, 23 have empirical formula of the form C_xH_yN_zO_w, six of the form C_xH_yO_z and four can be written C_xN_yO_z. The remaining 11 compounds also contain carbon, hydrogen, nitrogen and oxygen, but some of them contain either chlorine, phosphor or sulphur. These latter compounds were included to shed some light on the composition-dependence of the studied models, and, in particular, to investigate the suggestion by Rice et al. that the density of halogen-containing compounds might be much less accurately estimated by their volume (10%–40%) [5]. An overview of the compounds is given in the Supporting Information Table 1. The experimental densities were obtained from the SciFinder database [30].

Table 1. Model prediction errors (RMSEP). All volumes are given in units of cm³/mol.

	Isosurface	Π	σ		Θ
B3LYP	7.40	6.27	6.04	6.59	7.63
M05	7.49	6.39	6.16	6.75	7.72

Table 2. Listed are the B3LYP/6-31G(d) enclosed volumes V_encl in our study, and the B3LYP/6-31G(d,p) enclosed volumes obtained by Rice et al. [5,8] using the geometry in the crystal as the initial guess. All values are given in cm³/mol.

Structure	Formula	Compound	Vencl
	C3H5N7O5	1	126.0	125.1
	C12H11NO3	26	151.6	152.3
	C6H4N2O5	33	108.1	108.2
	CH2N6O2	40	74.2	73.7

Table 3. Fitted parameters for the enclosed volume estimates and the models including Π and σ corrections. The AICc values are listed with respect to the model the lowest value. The unit of β is cm³/(mol Hartree) and the unit of γ is cm³/mol.

		β	γ	AICc − AICcmin
B3LYP/6-31G(d)	Enclosed (14(14) )	–	2.332	15.8
	Π (15a(15a) )	−703.735	18.380	1.9
	σ (15b(15b) )	−720.314	22.296	0.0

Molecular geometry optimisations were performed with the GAUSSIAN03 Revision E.01 software [31], applying DFT with the B3LYP functional and basis sets of type 6-31G(d). Adopting the approach of Qiu et al. [32], we constructed the initial geometries using no information of the molecular geometry in the crystal. This procedure was different from that of Rice et al. [5,8] who selected the initial geometry as that of the crystal, thereby minimising the error caused by an inaccurate guess of the initial geometry, leading to a good volume estimate on which to add further corrections. It is of course also of interest to be able to predict the density for novel or notional compounds for which there exists no crystalline data. The viability of this approach rests on the assumption that geometries are, for the most part, of sufficiently small importance for the resulting volume estimate.

We performed single-point calculations with the B3LYP and M05 functionals in the search for possible functional dependencies of the volume estimates and the electrostatic parameters. The electron density and the electrostatic potential were evaluated on a 100 by 100 by 100 grid. Points were identified as being on the surface with a tolerance of 10⁻⁴ a.u., yielding 6000–10, 000 points for each molecule. The statistical quantities Π, σ²₋, σ₊², σ and Θ were subsequently evaluated on this surface, and UCSF Chimera [33] was applied to find the enclosed volume V_encl.

Five models for the molecular volume V are presented. While earlier studies have proposed models where corrections are linear in ρ [6–8], we will continue to assume linearity in our models of V. First, we study the uncorrected model (14) Equation (14(14) ) will provide a good reference to ascertain the efficiency of the more complicated models. These have been selected to identify the factor responsible for the reported improvement from electrostatic considerations [6–8]. We consider two models that correct for the variation of : (15a) (15b) We also investigate whether the imbalance of positive and negative values is a relevant factor, both in combination with variation () and as a potential correction in its own right (Θ): (16a) (16b) The models which include the quantities Π and suggested by Politzer et al. [6] are Equations (15a(15a) ) and (16a(16a) ), respectively. Although the predictions of the model based on Θ will be altered with changes in the surface [28], we have included it for the purpose of ascertaining its usefulness on the particular surface.

Due to the relatively small set of 44 compounds, we performed a 10-fold cross-validation for each model. Dividing the set into 10 separate folds, the models were trained using 9 out of 10 folds in 10 consecutive runs. In each run i, we computed the root-mean-square error of prediction (RMSEP), defined by [24] (17) In this equation, N is the number of predicted compounds in the given run, and V^exp the experimental molecular volume, as determined by Equation (1(1) ) from the experimental density ρ^exp. The model RMSEP was then taken to be the average error of prediction of the 10 consecutive runs: (18)

3. Results and discussion

The main results of our work can be divided into two parts. We begin by studying the predictive ability of the models that include electrostatic corrections based on variation, imbalance and the combination of the two. Subsequently we take a closer look at the relation of to variation measures, such as , Π and σ, arguing that the success of is due to variation and not imbalance.

3.1. Predictive ability

In this section, we review the results of the 10-fold cross-validation. The predictive ability of the models, given in terms of their RMSEP values, is shown in Table 1.

The best performing models are those based on the variation parameters Π and σ, i.e. Equations (15a(15a) ) and (15b(15b) ). Compared to the uncorrected enclosed volume estimates, the Π and σ parameters show a 15% and 18% relative decrease in RMSEP, respectively, regardless of functional. Moreover, it is evident that the imbalance-variation parameter leads to no improvement over the simpler variation parameters Π and σ. Since applying imbalance as a linear term leads to inconsistent predictions on different surfaces [28], it is not surprising that Θ performs poorly, yielding no improvement over uncorrected volume estimates. Employing the surface invariant skewness s′ about zero (see [28]) in leads to no improvement over Θ = |s′|; however, we believe this is due to its irrelevance to the problem at hand. This is in contrast to Θ, where the poor performance is likely attributable to its non-linear variability with changes in the molecular surface. In conclusion, our study provides evidence in favour of the view that the variation of on the molecular surface is related to the error of naive volume estimates. Furthermore, the findings clearly indicate that the imbalance of negative and positive values is of no relevance, as this factor does not improve upon variation parameters and appears to be of no use on its own.

In order to investigate the observed changes in RMSEP further, we inspected the error distributions of the models, see Figure 2. Observe that the effect of Π and σ sharpens the peak of the error distribution, causing more compounds to be estimated within 5 cm³/mol of experiment. The tails are left mostly unchanged. Considering B3LYP, the σ correction increases the number of compounds within 5 cm³/mol from 22 to 27 for Π and to 29 for σ. The figures are similar for the M05 functional. The similar accuracy of the two functionals is due to the fact that they agree to a high extent on the value of the parameters (R² > 0.99) and that their volume estimates V_encl are very similar on the B3LYP geometries (standard deviation equal to 0.25 cm³/mol). It should be noted, however, that the functionals might disagree to a significant extent if the optimisations are done independently. Thus, if one observes a considerable functional dependence, one may expect that this is mainly due to their influence on the optimised geometry.

Figure 2. Error distributions of the models. The M05 distributions are illustrated in dark grey, while those of B3LYP are light grey. The errors are given in cm³/mol.

A concern that must be addressed is that optimised geometries of isolated molecules will differ from the molecular geometries in the crystal. To address this issue, we compared our values with those obtained by Rice et al. [5,8], as they ran the optimisation using the geometry of the compound in the crystal as the initial guess. Their values for the four compounds in common with our study are in good agreement with our data, the deviations being less than 1 cm³/mol, see Table 2. We also achieve comparable accuracy, with our estimates yielding for Equations (14(14) ) and (15b(15b) ) a mean average error of 0.070 and 0.054 g/cm³ from the experimental densities, compared to the respective values 0.050 g/cm³ (V_encl) and 0.035 g/cm³ () obtained by Rice and coworkers [8]. Although a more extensive study is required to draw definite conclusions, we note that the considerable agreement and comparable accuracy suggest that crystal densities are in many cases effectively predicted without any knowledge of the molecular geometry in the solid. It may be that our results underestimate the importance of the initial geometry in general, however. We note that many of the compounds studied have limited conformational freedom to change as a consequence of the altered environment in the solid.

Since the parameters considered in previous studies [6–8] have only been applied to CHNO compounds, we decided to look closer at the σ values of non-CHNO compounds, see Figure 3. The effectiveness of variation (here given by σ) is evident from the figure, as it is seen that the higher the variation of the potential , the more likely V_encl is an overestimate. This is consistent with the interpretation given by Politzer et al. [6], who suggested that variation of makes attractive intermolecular interactions more likely, causing the molecules to be closer in the solid than one would expect from their volume. While the non-CHNO compounds appear to be accurately described by both and σ, as compared to the CHNO systems, we note that the three compounds of highest σ all contain sulphur. A more extensive study is thus needed to ascertain the domain of validity of the variation parameter σ.

Figure 3. The error of the enclosed volume estimates V_encl plotted against σ. Non-CHNO compounds are shown in dark grey, and CHNO compounds in light grey. The values of σ and V_encl were computed with the M05 functional and are given in units of Hartree and cm^3/mol, respectively.

For future reference and use, we have fitted the uncorrected models as well as those corrected with Π and σ for B3LYP, see Table 3. It should be kept in mind that these models, being based on all compounds, have not been validated. For this reason, we have also listed their AIC values, corrected for sample size (AICc). A rule-of-thumb lends ‘substantial support’ (<2) to the Π and σ corrected models, but ‘essentially no support’ (>10) to the uncorrected model [23].

3.2. The possible irrelevance of imbalance

We have demonstrated that the imbalance-variation parameter leads to no noticeable improvement over the simpler variation parameters Π and σ. Our work suggests that the findings of Rice and Byrd [8] is due to the close relation of to , a measure of variation. A comparison of the two reveals R² to be about 0.75, see Figure 4.

Figure 4. The values of plotted against the variation quantities σ² (top) and (bottom). All values are given in Hartree².

Furthermore, since is a measure of spread, it is also related to Π and σ. Considering the relation of to σ², we found a correlation coefficient of R² = 0.40–0.45. Though not particularly high, it might be sufficiently high for both of them to correct for the same aspect of , namely variation. Politzer et al. [6] reported R² = 0.68 between Π and , but concluded that this value was sufficiently low for them to signify different aspects of . While it is true that ν might be a small correction, there is as yet no evidence for this, neither in our results nor in the work of Politzer and coworkers [6]. The slight decrease in the training set RMSE is not significant (0.053–0.047 g/cm³). We therefore conclude that ν is adding noise, not correcting.

This discussion has presupposed that Θ or ν are suitable measures of imbalance. We note that they are at the very least consistent with each other, as they have been shown to be in fair agreement [28]. By this we mean that, if the absolute skewness Θ is zero, it is likely that the imbalance parameter ν is equal to 0.25 (indicating no imbalance). Similarly, as Θ increases, the value of ν is seen to decrease.

4. Conclusions and summary

In agreement with previous studies, the enclosed volume of surfaces of small electron density is seen to efficiently predict the crystal density of neutral molecular crystals. Moreover, we have found that, even when the geometry is not taken to be the one experimentally found in the solid, the predictive ability is largely maintained. Our work shows that the presence of other atoms than carbon, hydrogen, nitrogen and oxygen does not result in a breakdown of the method, as the eleven such compounds studied are equally well approximated as the CHNO systems. Our study presents evidence in favour of the view that variation of the electrostatic potential (measured by the standard deviation and related quantities) on the surface of the molecule is responsible for small but possibly significant corrections to quantum mechanical volumes. Applying two measures of variation, we obtained a relative decrease in RMSEP of 15%–18%. The corrections sharpen the error distribution peak, yielding more predictions within 5 cm³/mol of experiment.

The imbalance of the electrostatic potential on the surface of molecules is demonstrated to be irrelevant for the prediction of crystal densities. First, when used as a linear term, it leads to no improvement over the uncorrected estimates. In view of the fact that imbalance yields inconsistent predictions on different surfaces [28], this finding is not surprising. Second, when combined with the variation of the electrostatic potential, it does not improve predictions compared to the simpler corrections based on variation alone. The previously reported success of imbalance combined with variation is therefore likely due to variation alone, i.e. the role of imbalance is to add noise rather than to correct.

Supplemental data for this article can be accessed here: http://dx.doi.org/10.1080/00268976.2016.1158420.

Disclosure statement

No potential conflict of interest was reported by the authors.

Notes

1. We have used the unbiased estimator of the standard deviation, dividing by N − 1 instead of N. This choice does not affect the results since N is greater than 6000 in our calculations.