Remarks on X-ray constrained/restrained wavefunction fitting

X-ray constrained/restrained wavefunctions (XCWs/XRWs) result from a combination of theory and experiment and are therefore affected by experimental errors and model uncertainties. The present XCW/XRW procedure does not take this into account, thus limiting the meaning and significance of the obtained wavefunctions.

1. Introduction

In 1998, Jayatilaka (1998) presented the X-ray constrained wavefunction (XCW) fitting approach, a method that combines quantum chemistry and X-ray diffraction to produce what the author called an experimental wavefunction. The intention is to find a wavefunction Ψ_opt and the corresponding electron density ρ_opt, (i) that is as close as possible to the crystal density ρ_cryst, i.e. explains the observed structure factor amplitudes to within a weighted mean square deviation Δ between experiment and model and (ii) that simultaneously minimizes the quantum mechanical energy of the investigated system.

At present, XCW fitting (Jayatilaka, 1998; Jayatilaka & Grimwood, 2001; Grimwood & Jayatilaka, 2001; Bytheway, Grimwood & Jayatilaka, 2002; Bytheway, Grimwood, Figgis et al., 2002; Grimwood et al., 2003) is the last in a series of steps to model ρ_cryst. The procedure usually begins with Independent Atom Modelling (IAM, with tabulated spherical atomic densities and form factors). It may be followed by interpreting the diffraction data with more refined descriptions of the atomic and bonding electron densities, such as multipole models (MMs) (Stewart, 1976; Hansen & Coppens, 1978) or quantum mechanical Hirshfeld atom refinements (HARs) (Jayatilaka & Dittrich, 2008; Capelli et al., 2014; Kleemiss et al., 2021). They are usually based on Hartree–Fock (HF) or density functional theory (DFT) calculations, with or without consideration of the environment of the structural fragment that may also be treated quantum mechanically (Jayatilaka & Dittrich, 2008; Wieduwilt et al., 2021; Ruth et al., 2022). HAR at correlated levels has also been attempted (Wieduwilt et al., 2020). At the last stage, the XCW procedure tries to extract information not included in the preceding steps, typically the polarization due to the crystal environment (Ernst et al., 2020) as well as correlation (Genoni et al., 2017) and relativistic effects (Bučinský et al., 2016).

The XCW procedure thus attempts to extract the last little bit of information from a diffraction experiment. In other words, it tries to improve an already quite detailed model of the electron density, whose R factor is typically as low as 1–3%. The likelihood that this information is of similar magnitude as the systematic and random errors in the diffraction data and inadequacies of the quantum chemical model is thus not negligible and creates some unresolved problems concerning the meaning of the results obtained from XCW fitting. Some of these problems, not all of which have been addressed explicitly, are discussed below.

2. A summary of the X-ray constrained/restrained wavefunction approach

As mentioned in the introduction, the XCW fitting approach aims to determine a wavefunction that minimizes the quantum mechanical energy of the considered system and simultaneously maximizes the agreement between experimental and calculated structure factor amplitudes (Jayatilaka, 1998; Jayatilaka & Grimwood, 2001; Grimwood & Jayatilaka, 2001; Bytheway et al., 2002; Bytheway, Grimwood, Figgis et al., 2002; Grimwood et al., 2003; Jayatilaka, 2012). The problem is mathematically equivalent to finding the wavefunction that minimizes the energy under the constraint

where GoF² is the squared goodness-of-fit given by the following expression:

is the r-th structure factor amplitude calculated from a model electron density ρ(Ψ); is the r-th observed structure factor amplitude with as its estimated standard uncertainty; η is a scale factor that puts the calculated structure factor amplitudes on the same scale as the experimental ones; N_refl is the number of experimental observations, and N_par the number of parameters. Furthermore, as indicated by Jayatilaka and collaborators in the seminal papers on the technique, in principle Δ should be set equal to 1.0 (Jayatilaka, 1998; Jayatilaka & Grimwood, 2001; Grimwood & Jayatilaka, 2001; Bytheway, Grimwood & Jayatilaka, 2002; Bytheway, Grimwood, Figgis et al., 2002; Grimwood et al., 2003). This implies that the objective of XCW fitting is the determination of an `experimental wavefunction' that reproduces – on average – the experimental structure factor amplitudes to within one standard deviation of the experimental measurements.

To minimize the energy of the system under the constraint expressed by relation (1), the following functional has been introduced:

The first term on the right-hand side is the quantum mechanical energy of the system [i.e. = , with as the Hamiltonian operator], and λ is a Lagrange multiplier that has the dimension of an energy. To achieve the aim of XCW fitting, functional J must be made stationary with respect to both the wavefunction Ψ and the Lagrange multiplier λ. However, in the current implementation of the XCW technique, the second stationary condition (i.e. the one with respect to λ) is never considered and, for this reason, the functional that is really minimized is the following one:

This implies that the present XCW method simply attempts to determine the wavefunction that minimizes the energy of the system under examination and that reproduces as much as possible the measured X-ray diffraction data. λ thus loses its original meaning of Lagrange multiplier; instead, it can be considered as the weight of the experimental data in the computations and simply becomes an external parameter that is gradually changed during the XCW procedure (see below). The term λGoF²(Ψ) in equation (4) is a restraint in the crystallographic sense and the method should better be called X-ray restrained wavefunction (XRW) approach (Jayatilaka, 2012; Grabowsky et al., 2017; Ernst et al., 2020; Macetti et al., 2021) or X-ray regularization procedure (Davidson et al., 2022a). For the sake of completeness, it is mentioned here that one of the present authors has very recently proposed to reformulate the Jayatilaka technique by explicitly considering the stationary condition of functional (3) with respect to λ (Genoni, 2022). Only in that case does λ regain its meaning as a Lagrange multiplier, and the approach can then be called X-ray constrained wavefunction (XCW) method.

In this work, we focus on the original and current version of the Jayatilaka strategy, but refer to it as X-ray restrained wavefunction (XRW) technique. Furthermore, although multi-determinant versions of this method have already been developed (Genoni, 2017; Casati et al., 2017; Genoni et al., 2018, 2019), for simplicity and without losing generality, we limit our discussion to the original and most used version of the Jayatilaka strategy, namely the case in which the wavefunction to be determined has the analytical form of a single Slater determinant.

For the single Slater determinant wavefunction ansatz, it is possible to show (Jayatilaka, 1998; Jayatilaka & Grimwood, 2001; Grimwood et al., 2003) that finding the wavefunction which minimizes functional (4) corresponds to determining the molecular orbitals {φ_i} that satisfy the following modified Hartree–Fock equations (for a 2N-electron closed-shell system):

is the usual Fock operator used in quantum chemistry, while the XRW operator is

with

and

and in equation (7) are the real and imaginary parts, respectively, of the one-electron scattering operator

where N_m is the number of symmetry-unique positions in the unit-cell, {R_k, r_k} are the roto-translations associated with the unit-cell symmetry operations (R_k and r_k being a rotation matrix and a translation vector, respectively), B is the reciprocal lattice matrix, and h_r is the vector of Miller indices associated with the r-th reflection.

Note that the term can be seen as a perturbation with weight λ to the traditional Fock operator of quantum chemistry. However, depends on the experimentally determined structure factor amplitudes. Given that experimental measurements are subject to random measurement errors, is subject to random error too and will differ somewhat for different X-ray experiments performed on the same system under the same experimental conditions. In other words, not only accounts for shortcomings of the model of the electron density but is also characterized by a random component.

The interpretation of as a perturbation operator can be traced back to a pioneering idea by Weiss from the 1960s (Weiss, 1966). He envisaged the possibility of correcting the deficiencies of the Hartree–Fock model by recovering effects of electron correlation from the experimental X-ray diffraction data. More recently, the relation between Weiss' original ideas and the XRW approach has been highlighted by Macchi (2022) in a review on origins, state of the art and perspectives of quantum crystallography.

Expressing the molecular orbitals {φ_i} in equations (5) and (7) in terms of a set of basis functions (as it is usually done in quantum chemistry), the modified Hartree–Fock equations expressed by relation (5) straightforwardly transform to modified Roothaan–Hall matrix equations:

with F as the typical Fock matrix of quantum chemistry calculations, C as the matrix whose columns contain the coefficients that expand the molecular orbitals {φ_i} in the chosen basis, E as the diagonal matrix of the orbital energies, and v as the matrix associated with the operator , which obviously has the same dimensions as the Fock matrix and is weighted by the external multiplier λ.

In the current version of the XRW approach, equations (10) are solved self-consistently for different values of λ by following a `trial and error' procedure. This means that the molecular orbitals obtained through the computation with λ⁽ⁱ⁾ at the i-th step are used as guess for the calculation with λ⁽ⁱ⁺¹⁾ = λ⁽ⁱ⁾ + Δλ. The procedure is usually iterated, but without a clear and definitive criterion to stop. Several different halting criteria have been proposed over the years. For a review see Davidson et al. (2022a); for two recent proposals see Davidson et al. (2022b). For reasons that will be discussed in §3, none of these seems convincing. The consequence is that although the XRW process should be halted at GoF² ∼ 1.0, the final value of GoF² often remains greater than the desired value or, under some circumstances, the calculations may go on beyond the desired limit to a final GoF² value that is lower than 1.0.

The present version of the XRW procedure assumes that the nuclear positions and the atomic displacement parameters are known (to within their uncertainties). As also mentioned in the Introduction, these parameters may be obtained with the help of HAR (Capelli et al., 2014) since structural parameters resulting from Hirshfeld atom refinements of X-ray diffraction data are almost as accurate as those obtained through neutron diffraction experiments (Woińska et al., 2016). HAR also provides a starting point for the XRW procedure in the form of the final HAR wavefunction Ψ and, consequently, of the final electron density ρ, which hereafter will be called Ψ_ref and ρ_ref, respectively, with a corresponding squared goodness-of-fit [GoF²]_ref. Note, however, that, given accurate atomic positions and ADPs, the quantities Ψ_ref, ρ_ref and [GoF²]_ref might also be obtained from a level of theory and a basis set different from those used for the Hirshfeld atom refinement. No problem so far.

The subsequent XRW procedure modifies Ψ_ref (and ρ_ref) to Ψ_opt (and ρ_opt) until one obtains the best fitting of the experimental crystal density ρ_exp or, equivalently, the best fitting of its Fourier transforms . In other words, the XRW method tries to extract information from the measured structure factor amplitudes that is not already contained in , such as polarization, electron correlation and relativistic effects. An XRW fitting adjusts – at least in principle – all molecular orbital coefficients. However, for a 2N-electron closed shell system described by a single Slater determinant, there are N × (N_bf − N) independent parameters to be determined (with N_bf as the number of basis functions used in the calculation). This number is usually about an order of magnitude larger than the number of observations N_refl. It is therefore important to choose the λ value or, alternatively, the target value for GoF² in such a way that the procedure produces the best approximation Ψ_opt to Ψ_cryst, with Ψ_opt as little contaminated as possible by imponderables of the diffraction experiment or the quantum mechanical method used, i.e. systematic and random error.¹ The problem of choosing λ has been termed the halting problem. For reasons given in §3, it is our opinion that no tried-and-tested, convincing solution of this problem has been presented so far.

3. Some difficulties

To gain more insight into the XRW fitting process, a few cases with specific assumptions about , and are analyzed. Some of these assumptions are unrealistic but provide a reference for gaining insights into the fitting procedure applied to more realistic cases.

(A) First of all, let us assume the unrealistic case with:

(i) = , where is the random error in ;

(ii) = ;

(iii) are trustworthy and normally distributed estimates of the random errors of ;

In this case, GoF² (Ψ_ref) = 1 and it is trustworthy. Therefore, no further fitting is indicated.

(B) Making the same assumptions as in case (A) except that the uncertainties unknowingly underestimate the random errors of .

In this situation, GOF² (Ψ_ref) > 1 and further fitting is indicated. However, any such fitting would lead to wavefunctions and thus to structure factor amplitudes which include contributions from both and . In other words: for any non-zero value of λ, the XRW procedure merely fits experimental noise! For different experiments of the same substance under the same conditions, the random errors will be distributed differently and, consequently, the fit will look different. Underestimating is not uncommon in actual experiments. Note that this conclusion would not hold if, under the given assumptions, = 0. However, this seems unlikely and has certainly not been shown to be the case.

(C) Another artificial case assumes that

(i) (i.e. no random errors are included);

(ii) for all reflections;

(iii) Ψ_ref is based on a quantum mechanical model with a basis set that is insufficiently flexible to reproduce Ψ_cryst.

This situation corresponds to the case of theoretically generated structure factor amplitudes and the fitting leads to = and = = , where is the part of that cannot be modelled with . Given the shortcomings of the basis set, ≠ 0; GoF² ≠ 0 as well, but its expectation value is unknown because it depends on the unknown limitations of the model used. GoF² is expected to decrease with increasing the number of basis functions (and thus the number of parameters that can be optimized).

For this kind of problem, a way to find an optimal λ has been proposed, albeit based on electron densities rather than on their Fourier transforms, the structure factors (Tozer et al., 1996). Equations analogous to equations (10) are solved by combining Ψ_ref with three different values of λ and extracting an optimal λ from the three results. One of the procedures proposed in Davidson et al. (2022b), although inspired by Tozer et al. (1996), is different. It updates Ψ_ref at each fitting step until a halting point is reached. It has not been shown that the two procedures are equivalent, i.e. lead to the same result. In any case, the important point is that the choice of Ψ_ref affects the GoF² independently of and in an unpredictable way.

(D) For a more realistic case of XRW fitting, it is assumed that:

(i) ;

(ii) ;

(iii) Ψ_ref is based on a quantum mechanical model with a basis set that is likely to be insufficiently flexible to reproduce Ψ_cryst;

(iv) the are trustworthy and normally distributed estimates of the random errors of .

In this case, the fitting leads to and , in analogy to case (C). However, in contrast to that situation, the modifications of in correct not only for inadequacies of Ψ_ref, but also absorb into the model parts of , i.e. some of the experimental noise, similarly to case (B). The two contributions are of unknown magnitudes and thus has an undefined uncertainty on two accounts. Conversely accounts for those parts of and that could not be absorbed into . It is thus unclear again to know the expectation value of GoF² and λ.

(E) The most realistic case in XRW fitting makes the same assumptions as in (D) except that:

(i) , where represents a systematic error of due to inadequacies of the processing of the primary data, e.g. the frames from the diffraction experiment;

(ii) the values are not trustworthy estimates of the random errors of because they are too small, too large or show an unusual (non-normal) distribution.

In this situation, the fitting leads to and , in analogy to cases (C) and (D). However, unlike (C) and (D), the calculated correct not only for inadequacies of Ψ_ref, but also absorb into the model parts of , i.e. some of the systematic and random errors. The three contributions are again of unknown magnitudes and thus has undefined meaning and uncertainty. Conversely accounts for those parts of , and that could not be absorbed into . If, in addition, the quality of the values is unknown, estimating the expectation values of GoF² and λ is further aggravated.

In a recent paper, Davidson et al. (2022a) consider the XRW fitting approach as a regularization problem, i.e. the process of adding information in order to solve an ill-posed problem or to prevent overfitting [here we use the term `regularization' as defined in https://en.wikipedia.org/wiki/Regularization_(mathematics)]. In the present case there are two ways of looking at this problem: either as adding experimental information to a quantum chemical model that is insufficient to determine Ψ_cryst on its own, or as adding theoretical quantum mechanical information to a model of the experiment that is usually strongly underdetermined because the number of coefficients in the wavefunction exceeds the number of observed structure factor amplitudes by far. In both cases the question arises as to how much extra information is required to get an optimal result, i.e. what is the optimal magnitude of λ and how should an optimal mix of calculated and experimental information be defined anyway? This is an alternative formulation of what has been called the halting problem (Davidson et al., 2022a,b). If one can reasonably assume that the model of equation (4) is sufficiently flexible to determine Ψ_cryst and that the values are truthful estimates of the experimental uncertainties, one might choose a value of λ that results in GoF² = 1. However, these conditions are rarely fulfilled implying that the choice of an appropriate final value for GoF² remains an open problem in the XRW scheme.

In both ways of looking at the problem, systematic and random errors also introduce uncertainty in the wavefunction, Ψ_opt, not only in the choice of λ or of the final GoF² value. Since the goal of the exercise is to find an `experimental' X-ray restrained wavefunction and thus an experimental electron density function ρ_opt, an additional question arises: what are the uncertainties in Ψ_op or ρ_opt? Which parts of ρ_opt are reliable and which are not? The present XRW fitting algorithm provides no information on this question. In contrast, least-squares modelling of the crystal electron density ρ_opt through more traditional methods (e.g. multipole models) determines a set of parameters, their standard uncertainties, and mutual correlations. There are statistical criteria for judging the significance of introducing additional parameters into the model, such as Hamilton's R factor ratio test (Hamilton, 1965). Such criteria help preventing overfitting by the model. The uncertainty of ρ_opt can be estimated through propagation of error or analogous procedures. None of these options are available in the present XRW algorithm. Note that the problem raised in this paragraph would also affect the reformulation of the Jayatilaka approach recently envisaged by one of the present authors (Genoni, 2022), although, in that case, if convergence is achieved, one can obtain the desired value Δ of statistical agreement between observed and calculated structure factor amplitudes.

In summary, there are several difficulties with the outcome of an XRW fitting and the meaning of an X-ray restrained wavefunction:

(i) an X-ray restrained wavefunction models random error to an unknown degree, while the more traditional least-squares optimization determines expectation values of parameters together with their uncertainties. XRW fitting with its quantum mechanical minimization condition determines parameters too, but not necessarily their expectation values and without uncertainty.

(ii) an X-ray restrained wavefunction models systematic errors in the structure factors or inadequacies of the atomic basis set to an unknown degree (as does a least-squares procedure), thereby obscuring the statistical meaning of GoF².

(iii) no procedure is available at present to quantify, at least approximately, the relative and absolute magnitudes of random and systematic errors, and of basis-set inadequacies in the optimized density ρ_opt.

(iv) another non-negligible drawback associated with both the XRW and XCW approaches has been alluded to before: it is the problem of correctly defining GoF² and its expectation value. Even among developers of the XRW/XCW approach, it is still unclear today which value to assign to N_par in equation (2). In some circumstances, it is simply set equal to the number of adjustable parameters used in the `wavefunction refinement' (such as λ when it is manually adjusted in the XRW calculations); in other cases, for example when the XRW computation follows a Hirshfeld atom refinement, N_par includes the number of atomic positions and ADPs determined through the preliminary structural refinement; finally, in one of the seminal papers of the X-ray restrained wavefunction approach (Grimwood et al., 2003), it is pointed out that N_par should vary as a function of λ, approaching the number of molecular orbital coefficients when λ becomes large [which in most practical applications would lead to a negative value in the denominator of the first factor on the right-hand side of equation (2)]. This uncertainty in the definition of N_par further aggravates establishing the final and desired value of GoF² and, consequently, determining the correct value of λ.

Therefore, due to above-mentioned problems, an X-ray restrained wavefunction loses – to an unknown degree – its intended meaning as a representation of ρ_opt, because it cannot be decided whether the XRW density ρ_λ = ρ_opt, i.e. it cannot be decided whether or not the final electron density distribution resulting from the XRW process corresponds to the best possible representation of the real electron density ρ_cryst. For the same reasons the halting problem is ill defined.

4. Possible remedies

What can be done to alleviate these uncertainties?

(a) The data: utmost care is needed to get minimally contaminated by systematic errors. It should be clearly separated from diffuse scattering and very carefully corrected for the properties of the sample (absorption, extinction, etc.), the detector (non-linear response, oblique incidence, etc.) and other imponderables of the experiment. Such corrections are ideally based on physical rather than empirical models of the factors to be accounted for. Good data, ideally coming from more than one sample, tend to minimize effects of some systematic errors, . Data should be tested for consistency between and within samples.

(b) The values should be trustworthy and their distribution tested for consistency between and within samples. Experiments that produce small, trustworthy values minimize effects of . Good values help to make GoF² quantities more meaningful, of course provided that the number of extra (effective) parameters coming from the XRW fitting can be estimated. So far this is not done. should be tested for outliers with defined criteria before and after XRW fitting, e.g. with normal probability plots.

(c) Separate analysis of experimental data should be done for different choices of data: splitting data into parts (cross validation; Krause et al., 2017), looking at different samples separately, and combining all data. It is expected that such a procedure leads to different because the influence of will be different. Comparisons and analysis of different for the same compound could be a first step towards characterizing experimental uncertainty.

(d) To gain a fuller picture of the influence of experimental uncertainty, including the part coming from the atomic positional and displacement parameters, the following iterative procedure might be used in validation experiments: the initial XRW electron density can serve to define Hirshfeld atoms which are used in turn for refining atomic positions and displacement parameters. This might be followed by another XRW fit, etc., until convergence is reached, following the spirit of the XWR (X-ray wavefunction refinement) approach already envisaged and proposed by Grabowsky and collaborators (Grabowsky et al., 2012; Woińska et al., 2017).

(e) In addition, the results should be tested by analysing the same data starting from different models Ψ_ref (e.g. HF, DFT and different basis sets). If the different XRW fittings give comparable results for , the analysis can be considered robust.

(f) One could also think of comparing results from alternative expressions for GoF², e.g. one that includes the phases of the structure factors:

or some function of their Fourier transform, i.e. some function of the difference density. Since the XRW fitting comes into play at the very end of the ρ_cryst modelling, it can reasonably be assumed that the phases of and are the same (in the centrosymmetric case) or very close to each other (in the non-centrosymmetric case). This definition would be more akin to a restraint based on the difference electron density, provided one includes a complete set of structure factors up to the resolution limit of the experimental data and GoF² is divided by the volume V of the unit cell.

(g) For many of the open problems, a better understanding of the XRW procedure could be obtained from a well crafted study with known, predefined , , , and combined in various ways, e.g. as mentioned for cases (A)–(E) above.

Footnotes

¹According to the IUCr definition, a systematic error is a contribution to the difference between an estimate and the true value of a quantity due to deficiencies of the model (https://www.iucr.org/resources/commissions/crystallographic nomenclature/statdes/terms.html). In the case of structure factors, such contributions may arise from two sources. On the one hand, from an inadequate model of the diffraction experiment needed to relate the raw intensity observations to what is commonly called the observed structure factor amplitude, . The symbol is somewhat misleading as it does not represent a direct experimental estimate, but rather a derived quantity that depends on assumptions and on estimates or measurements of factors affecting the measurement process (background, X-ray absorption, X-ray beam polarization, extinction, etc.). If such factors are insufficiently accounted for, the derived quantity will be unreliable to an unknown degree. On the other hand, an inadequate model of the electron density distribution, e.g. an insufficiently flexible basis used in the quantum mechanical calculation, can also cause systematic differences between and .