X-ray constrained spin-coupled technique: theoretical details and further assessment of the method

Genoni, A.; Macetti, G.; Franchini, D.; Pieraccini, S.; Sironi, M.

doi:10.1107/S2053273319011021

research papers

FOUNDATIONS
ADVANCES

ISSN: 2053-2733

Volume 75| Part 6| November 2019| Pages 778-797

https://doi.org/10.1107/S2053273319011021

X-ray constrained spin-coupled technique: theoretical details and further assessment of the method

Alessandro Genoni,^a ^* Giovanni Macetti,^a Davide Franchini,^b Stefano Pieraccini ^b,^c,^d and Maurizio Sironi ^b,^c,^d

^aLaboratoire de Physique et Chimie Théoriques, Université de Lorraine and CNRS, 1 Boulevard Arago, Metz, F-57078, France, ^bDipartimento di Chimica, Università degli Studi di Milano, Via Golgi 19, Milano, I-20133, Italy, ^cIstituto di Scienze e Tecnologie Molecolari (ISTM), CNR, Via Golgi 19, Milano, I-20133, Italy, and ^dConsorzio Interuniversitario Nazionale per la Scienza e Tecnologia dei Materiali (INSTM), UdR Milano, Via Golgi 19, Milano, I-20133, Italy
^*Correspondence e-mail: alessandro.genoni@univ-lorraine.fr

Edited by A. Altomare, Institute of Crystallography - CNR, Bari, Italy (Received 24 May 2019; accepted 7 August 2019; online 24 September 2019)

One of the well-established methods of modern quantum crystallography is undoubtedly the X-ray constrained wavefunction (XCW) approach, a technique that enables the determination of wavefunctions which not only minimize the energy of the system under examination, but also reproduce experimental X-ray diffraction data within the limit of the experimental errors. Initially proposed in the framework of the Hartree–Fock method, the strategy has been gradually extended to other techniques of quantum chemistry, but always remaining limited to a single-determinant ansatz for the wavefunction to extract. This limitation has been recently overcome through the development of the novel X-ray constrained spin-coupled (XCSC) approach [Genoni et al. (2018 ). Chem. Eur. J. 24, 15507–15511] which merges the XCW philosophy with the traditional spin-coupled strategy of valence bond theory. The main advantage of this new technique is the possibility of extracting traditional chemical descriptors (e.g. resonance structure weights) compatible with the experimental diffraction measurements, without the need to introduce information a priori or perform analyses a posteriori. This paper provides a detailed theoretical derivation of the fundamental equations at the basis of the XCSC method and also introduces a further advancement of its original version, mainly consisting in the use of molecular orbitals resulting from XCW calculations at the Hartree–Fock level to describe the inactive electrons in the XCSC computations. Furthermore, extensive test calculations, which have been performed by exploiting high-resolution X-ray diffraction data for salicylic acid and by adopting different basis sets, are presented and discussed. The computational tests have shown that the new technique does not suffer from particular convergence problems. Moreover, all the XCSC calculations provided resonance structure weights, spin-coupled orbitals and global electron densities slightly different from those resulting from the corresponding unconstrained computations. These discrepancies can be ascribed to the capability of the novel strategy to capture the information intrinsically contained in the experimental data used as external constraints.

Keywords: X-ray constrained wavefunction; quantum crystallography; valence bond theory; spin-coupled method.

1. Introduction

As is well known, in theoretical chemistry there exist two main approaches to investigate molecular electronic structure: the valence bond (VB) (Heitler & London, 1927 ; London, 1928 ; Pauling, 1939 ) and the molecular orbital (MO) (Hückel, 1930 , 1931 ; Roothaan, 1951 ; Dewar, 1952 ) theories. The former has been strictly related to the traditional chemical perception since its origin and significantly contributed to the definition of concepts (e.g. Lewis structures, resonance structure, hybridization, local bonds, electronegativity etc.) that, even today, are of customary use among chemists and constitute the basis of the traditional chemical reasoning to interpret bonding and reactivity. In contrast, the latter provides pictures of the electronic structure that are generally delocalized over the whole molecules under examination and that are consequently far from the traditional chemical notions. Despite this fact, the MO-based methods have become more and more common in electronic structure investigations, mainly due to their high predictive power, their intrinsic lower computational cost (at least for the basic strategies) and the ease with which they can be implemented into working computer codes (Shurki, 2006 ; Hiberty & Shaik, 2007 ). Nevertheless, owing to the unquestionable higher chemical interpretability associated with the VB theory, different VB approaches have been continuously proposed over the years (Shurki, 2006; Hiberty & Shaik, 2007; Goddard, 1967 ; Ladner & Goddard, 1969 ; Goddard et al., 1973 ; Gerratt & Lipscomb, 1968 ; Gerratt, 1971 ; Cooper et al., 1986 , 1991 ; Karadakov et al., 1992 ; Cooper et al., 1993 ; Raimondi et al., 1996 ; Thorsteinsson et al., 1996 ; Hirao et al., 1996 ; Voter & Goddard, 1981 ; Hollauer & Nascimento, 1993 ; McDouall, 1992 ; van Lenthe & Balint-Kurti, 1983 ; Hiberty et al., 1994 ; Song et al., 2003 ), each of them with its own features and with its main field of application. In this paper, we will particularly consider the spin-coupled (SC) method (Gerratt & Lipscomb, 1968; Gerratt, 1971; Cooper et al., 1986, 1991, 1993; Karadakov et al., 1992; Raimondi et al., 1996), a technique that uses a very general single-configuration wavefunction by including all the possible spin-coupling modesand without imposing any constraints (e.g. orthogonality constraints) on the orbital expansions. Despite its intrinsically correlated nature, the approach keeps a high degree of chemical interpretability. For example, the coefficients associated with the spin-coupling modes considered in the wavefunction are generally interpreted as a measure of the weights corresponding to the resonance structures of the investigated system, while the spatial extensions of the SC orbitals are generally exploited to get insights into the electron distributions around atoms and, consequently, into the atomic hybridizations.

Therefore, due to the wealth of information that one can generally extract from SC wavefunctions, we have recently combined the SC method of the VB theory with the X-ray constrained wavefunction (XCW) philosophy (Jayatilaka, 1998 ; Jayatilaka & Grimwood, 2001 ; Grimwood & Jayatilaka, 2001 ; Bytheway, Grimwood & Jayatilaka, 2002 ; Bytheway, Grimwood, Figgis et al., 2002 ; Grimwood et al., 2003 ; Hudák et al., 2010 ; Jayatilaka, 2012 ; Bučinský et al., 2016 ), thus giving rise to the novel X-ray constrained spin-coupled (XCSC) technique (Genoni, Franchini et al., 2018 ), with the final goal of proposing another useful method to extract chemical information directly from experimental X-ray diffraction data.

In fact, along with strategies such as the multipole models (Stewart, 1976 ; Hansen & Coppens, 1978 ) and the maximum entropy methods (MEM) (Sakata & Sato, 1990 ; Roversi et al., 1998 ; van Smaalen & Netzel, 2009 ), which are density-based approaches, the XCW technique introduced by Jayatilaka (1998) is currently one of the most important methods in the rising field of quantum crystallography (Genoni, Bučinský et al., 2018; Novara et al., 2018 ; Grabowsky et al., 2017 ; Massa & Matta, 2018 ; Tsirelson, 2018 ) and probably the most popular tool to extract wavefunctions and/or density matrices from experimental X-ray diffraction and/or scattering measurements. Its roots obviously date back to the pioneering methods proposed by Clinton, Massa and their co-workers to obtain N-representable one-electron density matrices from X-ray diffraction data (Clinton & Massa, 1972 ; Clinton et al., 1973 ; Frishberg & Massa, 1981 ; Goldberg & Massa, 1983 ; Massa et al., 1985 ). Inspired by these works, different researchers have devised alternative techniques to determine `experimental' wavefunctions/density matrices (Aleksandrov et al., 1989 ; Howard et al., 1994 ; Snyder & Stevens, 1999 ; Tanaka, 1988 , 2018 ; Hibbs et al., 2005 ; Waller et al., 2006 ; Schmider et al., 1990 , 1992 ; Weyrich, 2006 ; Gillet et al., 2001 ; Gillet & Becker, 2004 ; Gillet, 2007 ; Gueddida et al., 2018 ). Among them, it is worth citing Tanaka's X-ray atomic orbital (XAO) (Tanaka, 1988) and X-ray molecular orbital (XMO) (Tanaka, 2018) methods, the molecular orbital occupation number (MOON) approach (Hibbs et al., 2005; Waller et al., 2006) and, above all, the strategies introduced by Weyrich & Smith (Schmider et al., 1990, 1992; Weyrich, 2006) and by Gillet, Becker and co-workers (Gillet et al., 2001; Gillet & Becker, 2004; Gillet, 2007; Gueddida et al., 2018) to reconstruct both the diagonal and the off-diagonal parts of the one-electron density matrices through joint refinements of X-ray diffraction data, magnetic structure factors and inelastic Compton scattering measurements. However, the above-mentioned Jayatilaka XCW approach (Jayatilaka, 1998; Jayatilaka & Grimwood, 2001; Grimwood & Jayatilaka, 2001) represents the major breakthrough in this research field. In fact, bearing in mind Gilbert's corollary to Coleman's theorem (Gilbert, 1975 ), according to which there are in principle an infinite number of wavefunctions compatible with a given electron density, Jayatilaka proposed a method that can be considered as a practical implementation of Henderson & Zimmermann's idea (Henderson & Zimmerman, 1976 ) and which consists in finding wavefunctions that not only fit a set of diffraction/scattering data within the limit of the experimental uncertainties (as in all the other methods mentioned above), but also minimize the energy of the systems under examination.

Originally devised in the framework of the restricted Hartree–Fock (RHF) formalism (Jayatilaka, 1998; Jayatilaka & Grimwood, 2001; Grimwood & Jayatilaka, 2001; Bytheway, Grimwood & Jayatilaka, 2002; Bytheway, Grimwood, Figgis et al., 2002; Grimwood et al., 2003; Jayatilaka, 2012), the XCW method has been progressively improved by combining it with other strategies of quantum chemistry, such as the unrestricted Hartree–Fock technique (Hudák et al., 2010), relativistic approaches [particularly, the second-order Douglas–Kroll–Hess (Hudák et al., 2010) and the infinite-order two-component (Bučinský et al., 2016) methods] and density functional theory (DFT) (Jayatilaka, 2012). Therefore, the proposed XCW strategies have been mainly developed in combination with MO approaches, thus providing completely delocalized pictures of molecular electronic structures. To recover the traditional chemical interpretability from the obtained `experimental' wavefunctions, suitable a posteriori techniques have been exploited, as done for example by Jayatilaka, Grabowsky and co-workers (Jayatilaka & Grimwood, 2004 ; Grabowsky et al., 2010 , 2011 , 2012 ), who have successfully applied topological strategies [e.g. quantum theory of atoms in molecules (QTAIM) (Bader, 1990 ), the electron localization function (ELF) (Becke & Edgecombe, 1990 ), the electron localizability indicator (ELI) (Kohout, 2004 ) and the localized orbital locator (LOL) (Schmider & Becke, 2000 )] to get back to the traditional bonding patterns for the investigated systems. For the sake of completeness, another possibility suggested by the Grabowsky group is the application of the natural bond orbitals (NBO) method (Weinhold & Landis, 2001 ) and of the natural resonance theory (NRT) (Glendening & Weinhold, 1998 ) to determine the weights of the resonance structures for the molecules under examination (Fugel, Kleemiss et al., 2018 ; Fugel, Beckmann et al., 2018 ).

Before the recent introduction of the XCSC method (Genoni, Franchini et al., 2018), which will be described in detail in this paper, the only techniques that allowed one to directly obtain X-ray constrained wavefunctions already close to the traditional chemical perception without the application of subsequent bonding analyses were the XC-ELMO (Genoni, 2013a ,b ; Dos Santos et al., 2014 ; Genoni & Meyer, 2016 ) and XC-ELMO-VB (Genoni, 2017 ; Casati et al., 2017 ) approaches, even if both of them are characterized by the introduction of preliminary chemical constraints. In fact, the XC-ELMO method really provides X-ray constrained MOs strictly localized on elementary molecular fragments corresponding to atoms, bonds and functional groups, namely the so-called extremely localized molecular orbitals (ELMOs) (Stoll & Wagenblast, 1980 ; Fornili et al., 2003 ; Genoni & Sironi, 2004 ; Genoni, Fornili & Sironi, 2005 ; Genoni, Ghitti et al., 2005 ; Sironi et al., 2007 , 2009 ; Meyer et al., 2016a ,b ; Meyer & Genoni, 2018 ). Nevertheless, it is the user that preliminarily decides the fragmentation of the molecules into subunits according to chemical intuition or the specific computational needs, thus partially biasing the results of the calculations. In a quite analogous way, in the XC-ELMO-VB method the global wavefunction is written as a linear combination of pre-computed ELMO wavefunctions corresponding to the different resonance structures that one wants to take into account, but, during the computations, only the weights of these resonance structures are determined, while the ELMO wavefunctions remain unchanged.

This scenario has thus prompted us to develop a novel XCW technique able to extract chemically meaningful information from X-ray diffraction data without applying a posteriori techniques for the analysis of the obtained `experimental' wavefunctions or without imposing chemical constraints a priori. As mentioned above, this has been accomplished through the development of the new XCSC method (Genoni, Franchini et al., 2018), which is the extension of the XCW approach to the SC strategy of the VB theory. Our very preliminary investigations (Genoni, Franchini et al., 2018) have shown that the determination of XCSC wavefunctions is quite straightforward. The new technique is efficiently able to capture crystal-field effects intrinsically contained in the experimental structure factors used as external constraints in the XCW calculations, particularly revealing slight but non-negligible differences between the electronic structures in the gas phase and in the solid state.

In this paper, after presenting for the first time all the theoretical details of the new method, we will also discuss the results of new test calculations to further evaluate the overall performance of the technique. Finally, in the last section, we will draw some general conclusions and we will indicate possible future perspectives and extensions of the proposed strategy.

2. Theory

In this section we will describe the theory at the basis of the new XCSC method. In particular, we will focus (i) on the basic assumptions of the technique, (ii) on the presentation of the new wavefunction ansatz, (iii) on the expression of the SC electron density and (iv) on the implementation of the new strategy. For the sake of completeness, the detailed mathematical derivation to obtain the analytical expressions of the first and second derivatives of the statistical agreement between theoretical and experimental structure-factor amplitudes is given in Appendix A.

2.1. Basic assumptions of the method

As in any XCW method (Jayatilaka, 1998; Jayatilaka & Grimwood, 2001; Grimwood & Jayatilaka, 2001; Bytheway, Grimwood & Jayatilaka, 2002; Bytheway, Grimwood, Figgis et al., 2002; Grimwood et al., 2003; Hudák et al., 2010; Jayatilaka, 2012; Bučinský et al., 2016; Genoni, 2013a,b; Dos Santos et al., 2014; Genoni & Meyer, 2016; Genoni, 2017; Casati et al., 2017), the new XCSC approach works for molecular crystals and assumes them to be collections of non-interacting molecular units, which are described by formally identical and symmetry-related wavefunctions. On the basis of this working hypothesis and of the additional assumption that each non-interacting unit corresponds to a symmetry-unique portion of the crystal unit cell, the electron density of the crystal unit cell can be expressed as the sum of the N_m crystal unit electron densities $[{\rho }_{j}({\bf r})]$ related to the reference distribution $[{\rho }_{0}({\bf r})]$ through the unit-cell symmetry operations $[\{{{\bf Q}}_{j},{{\bf q}}_{j}\}]$ :

$[{\rho }_{\rm cell}({\bf r}) = \textstyle\sum\limits _{j = 1}^{{N}_{m}}{\rho }_{j}({\bf r}) = \sum\limits _{j = 1}^{{N}_{m}}{\rho }_{0}\left[{{\bf Q}}_{j}^{-1}({\bf r}-{{\bf q}}_{j})\right]. \eqno(1)]$

This equation is exact only provided that $[{\rho }_{0}({\bf r})]$ is not obtained by means of an isolated computation on the reference crystal unit. In the XCW approaches this requirement is satisfied by looking for the wavefunction $[{\Psi }_{0}]$ associated with electron density $[{\rho }_{0}({\bf r})]$ that minimizes the energy of the reference crystal unit and that, at the same time, reproduces a set of experimental structure-factor amplitudes within the limit imposed by the experimental uncertainties. In this way the effects of the surrounding crystal field are intrinsically taken into account in an effective way, although the XCW computations formally consist in single-molecule computations.

In other words, the new XCSC method and all the XCW techniques consist in determining a wavefunction $[{\Psi }_{0}]$ of the reference crystal unit that minimizes the Jayatilaka functional:

$[\eqalignno{J[{\Psi }_{0}]& = {{\langle{\Psi }_{0}|{\hat{\cal H}}_{0}|{\Psi }_{0}\rangle}\over{\langle{\Psi }_{0}|{\Psi }_{0}\rangle}}+\lambda \left({\chi }^{2}[{\Psi }_{0}]-\Delta \right)&\cr &= W[{\Psi }_{0}]+\lambda \left({\chi }^{2}[{\Psi }_{0}]-\Delta \right).& (2)}]$

In equation (2), the first term of the functional is obviously the energy part with $[{\hat{\cal H}}_{0}]$ as the non-relativistic Hamiltonian operator for the reference unit, while the second term of the functional is the part that takes into account the experimental constraints given by X-ray diffraction data, with $[\lambda]$ as an external multiplier that is manually adjusted during the calculations and that gives the strength of the experimental constraints, and $[\Delta]$ as the desired agreement between theoretical and observed structure-factor amplitudes (usually set equal to 1). $[{\chi }^{2}]$ is a measure of the statistical agreement between calculated and experimental values:

$[{\chi }^{2} = {{1}\over{{N}_{r}-{N}_{p}}} \sum _{{\bf h}}{{{\left(\eta |{F}_{{\bf h}}^{\rm calc}|-|{F}_{{\bf h}}^{\rm exp}|\right)}^{2}}\over{{\sigma }_{{\bf h}}^{2}}} \eqno(3)]$

with N_r the number of experimental constraints, N_p the number of adjustable parameters, $[{\bf h}]$ the triad of Miller indices which characterize the reflection, $[{\sigma }_{{\bf h}}]$ the experimental error associated with the generic experimental structure-factor amplitude $[|{F}_{{\bf h}}^{\rm exp}|]$ , and $[\eta]$ an overall $[{\bf h}]$ -independent scale factor which multiplies each computed structure-factor amplitude $[|{F}_{{\bf h}}^{\rm calc}|]$ and which is obtained through minimization of the $[{\chi }^{2}]$ value.

For the sake of completeness, considering the basic assumptions of the XCW philosophy and the fact that structure factors are Fourier transforms of the unit-cell electron density, we remind the reader that the calculated structure factors can be computed as

$[{F}_{{\bf h}}^{\rm calc} = {\rm Tr}\left[{{\bf D}}_{0}{\bf }{{\bf I}}_{{\bf h}}\right] \eqno(4)]$

with $[{{\bf D}}_{0}]$ the one-particle density matrix associated with wavefunction $[{\Psi }_{0}]$ for the reference crystal unit and $[{{\bf I}}_{{\bf h}}]$ the matrix of the Fourier transform integrals of the basis function products summed over the N_m equivalent unit-cell sites, namely

$[\eqalignno{\left[{\bf I}_{\bf h}\right]_{\mu \upsilon } &= \textstyle\sum\limits _{j = 1}^{{N}_{m}}\exp[i2\pi {\bf q}_{j}\cdot ({\bf B}{\bf h})] {T}_{\mu \upsilon }\left[{\bf B}^{-1}{\bf Q}_{j}^{\rm T}{\bf B}{\bf h}\semi {\bf U}^{\mu }, {\bf U}^{\nu }\right]&\cr &\quad\times\textstyle\int \,{\rm d}{\bf r} {\chi }_{\mu }({\bf r}) {\chi }_{\upsilon }({\bf r}) \exp[i2\pi ({\bf Q}_{j}{\bf r})\cdot ({\bf B}{\bf h})]. &(5)}]$

In equation (5), $[{\bf B}]$ is the reciprocal-lattice matrix and $[{T}_{\mu \upsilon }[{\bf B}^{-1}{\bf Q}_{j}^{\rm T}{\bf B}{\bf h}\semi{\bf }{\bf U}^{\mu }, {\bf U}^{\nu }]]$ is a term that accounts for the influence of the thermal motion on the electron density and that, in our new XCSC technique, has been evaluated following the thermal smearing approach proposed by Stewart (1969 ). It is worth noting that this term parametrically depends on the symmetric matrices $[{\bf U}^{\mu }]$ and $[{\bf U}^{\nu }]$ of the anisotropic displacement parameters (ADPs) of the atoms on which basis functions $[{\chi }_{\mu }({\bf r})]$ and $[{\chi }_{\upsilon }({\bf r})]$ are centred. As in any XCW fitting calculation, these ADPs are read as input parameters and never optimized during the computations.

2.2. The SC wavefunction ansatz

In the novel XCSC method, the analytical form for the wavefunction $[{\Psi }_{0}]$ associated with the reference crystal unit is that of the traditional SC wavefunctions for a system of N electrons in the spin state (S, M):

$[\eqalignno{{\Psi }_{0}& = {\cal N}{ \Psi }_{0}^{\rm SC} = {\cal N}\textstyle\sum\limits _{k = 1}^{{f}_{S}^{{N}_{v}}}{c}_{S,k} {\psi }_{S,M\semi k}^{N}&\cr &= {\cal N}\textstyle\sum\limits _{k = 1}^{{f}_{S}^{{N}_{v}}}{c}_{S,k} {\hat{\cal A}}\left({\phi }_{1}^{c}{ \overline{\phi }}_{1}^{c}\ldots {\phi }_{j}^{c}{ \overline{\phi }}_{j}^{c}\ldots {\phi }_{{N}_{1}}^{c}{ \overline{\phi }}_{{N}_{1}}^{c}{\boldvarphi }_{v} {\Theta }_{S,M\semi k}^{{N}_{v}}\right), &(6)}]$

where S and M are the quantum numbers associated with the spin operators $[{\hat S}{}^2]$ and $[{\hat{S}}_{z}]$ , respectively, $[{\cal N}]$ is the normalization constant, f_S^N_v is the number of linearly independent spin eigenfunctions $[\{{\Theta }_{S,M\semi k}^{{N}_{v}}\}]$ and, equivalently, of linearly independent SC structures $[\{{\psi }_{S,M\semi k}^{N}\}]$ with weights given by the spin-coupling coefficients {c_S,k}, and $[{\hat{\cal A}}]$ is the antisymmetrizer defined as

$[{\hat{\cal A}} = {{1}\over{\sqrt{N!}}} \sum _{{\hat{\cal P}}\in {S}_{N}}(-1)^{p} {\hat{\cal P}} \eqno(7)]$

with $[{\hat{\cal P}}]$ as a permutation operator belonging to the symmetric group S_N and with parity (-1)^p. Furthermore, $[\{{\phi }_{i}^{c}\}_{i = 1}^{{N}_{1}}]$ is the set of frozen, doubly occupied `core orbitals' that can be preliminarily obtained by means of unconstrained Hartree–Fock or X-ray constrained Hartree–Fock calculations and that describe the subset of 2N₁ core (or inactive) electrons, while $[{\boldvarphi }_{v}]$ is the product of the N_v `SC orbitals'

$[{\boldvarphi }_{v} \left({{\bf r}}_{1}, {{\bf r}}_{2},\ldots {{\bf r}}_{{N}_{v}}\right) = {\varphi }_{1}\left({{\bf r}}_{1}\right) {\varphi }_{2}\left({{\bf r}}_{2}\right)\ldots { \varphi }_{{N}_{v}}\left({{\bf r}}_{{N}_{v}}\right) \eqno(8)]$

that describe the N_v valence (or active) electrons. As in every SC calculation, these are the only electrons that are really treated at the SC level in order to reduce the computational cost. The `core orbitals' are constructed such that they are orthonormal among each other and orthogonal to the `SC orbitals', on which, in contrast, no orthogonality constraints are imposed:

$[\left\{\matrix{\langle{\phi }_{i}^{c}|{\phi }_{j}^{c}\rangle = {\delta }_{ij}\cr \langle{\phi }_{i}^{c}|{\varphi }_{j}\rangle = 0 \hfill\cr \langle {\varphi }_{i}|{\varphi }_{j}\rangle = {s}_{ij}. }\right.\eqno(9)]$

Here it is worth noting that, since both unconstrained and X-ray constrained Hartree–Fock orbitals are generally delocalized over all the systems under examination, they are usually localized by means of traditional a posteriori techniques (Boys, 1960 ; Foster & Boys, 1960 ; Edmiston & Ruedenberg, 1963 , 1965 ; Pipek & Mezey, 1989 ) to select which of them have to be used to describe the core (inactive) electrons. As a consequence, the `core orbitals' are automatically orthonormal among each other, while, to guarantee the orthogonality between the `core orbitals' and the SC orbitals, the latter are expanded on the set of the remaining localized occupied MOs (namely, those non-selected as `core orbitals') and of all the virtual (Hartree–Fock or X-ray constrained Hartree–Fock) MOs, namely:

$[{\varphi }_{i}({\bf r}) = \textstyle\sum\limits _{\mu = {N}_{1}+1}^{M}{C}_{\mu i} \phi_{\mu }^{\rm (XC)HF}({\bf r}), \eqno(10)]$

where M is the size of the adopted basis set (i.e. the number of atomic basis functions initially used to perform the preliminary Hartree–Fock or X-ray constrained Hartree–Fock calculation).

Because of the non-orthogonality of the SC orbitals, the SC structures in equation (6) are also non-orthogonal. Therefore, the spin-coupling coefficients {c_S,k} do not directly give the weights of the corresponding structures $[\{{\psi }_{S,M\semi k}^{N}\}]$ . One of the possible ways of obtaining the real weights is given by the Chirgwin–Coulson coefficients (Chirgwin & Coulson, 1950 ), which have been used in this work and are defined as

$[{w}_{S,k} = {\left|{c}_{S,k}\right|}^{2}+\textstyle\sum\limits _{j\ne k}{c}_{S,k} {c}_{S,j}{ S}_{kj}, \eqno(11)]$

where S_kj is the overlap integral between SC structures $[{\psi }_{S,M\semi k}^{N}]$ and $[{\psi }_{S,M\semi j}^{N}]$ .

To complete the overview on the SC wavefunction ansatz, it is also important to note that each generic spin eigenfunction can be expressed as a linear combination of spin primitive functions, namely

$[{\Theta }_{S,M\semi k}^{{N}_{v}} = \textstyle\sum\limits _{i = 1}^{{N}_{d}}{d}_{i,k}^{S} {\vartheta }_{i}, \eqno(12)]$

where

$[{N}_{d} = \left(\matrix{{N}_{v}\cr {N}_{\alpha }}\right) = \left(\matrix{{N}_{v}\cr {N}_{\beta }}\right),]$

with $[{N}_{\alpha }]$ and $[{N}_{\beta }]$ as the number of spin functions $[\alpha]$ and $[\beta]$ in the generic spin primitive function, respectively. The values of the coefficients {d_i,k^S} depend on the adopted basis to describe the spin space. In our work we have chosen the one constituted of Rumer spin eigenfunctions (Simonetta et al., 1968 ) and, for this reason, {d_i,k^S} can only be equal to 0 and ±1. Now, if we introduce equation (12) into (6), the SC wavefunction ansatz can also be rewritten like this:

$[{\Psi }_{0} = {\cal N}\textstyle\sum\limits _{i = 1}^{{N}_{d}}{b}_{S,i} {\hat{\cal A}}\left({\phi }_{1}^{c}{ \overline{\phi }}_{1}^{c}\ldots {\phi }_{j}^{c}{ \overline{\phi }}_{j}^{c}\ldots {\phi }_{{N}_{1}}^{c}{ \overline{\phi }}_{{N}_{1}}^{c}{{\boldvarphi }}_{v} {\vartheta }_{i}\right), \eqno(13)]$

where

$[{b}_{S,i} = \textstyle\sum\limits _{k = 1}^{{f}_{S}^{{N}_{v}}}{c}_{S,k} {d}_{i,k}^{S}. \eqno(14)]$

Since the generic spin primitive function $[{\vartheta }_{i}]$ is a product of N_v spin functions $[\alpha]$ and $[\beta]$ , the product $[{{\boldvarphi }}_{v} {\vartheta }_{i}]$ in equation (13) is actually a product of spin orbitals, i.e.

$[{{\boldvarphi }}_{v} {\vartheta }_{i} = {\varphi }_{1} {\varphi }_{2}\ldots { \varphi }_{{N}_{v}}{\vartheta }_{i} = {\varphi }_{1}^{i}{ \varphi }_{2}^{i}\ldots {\varphi }_{{N}_{v} }^{i}, \eqno(15)]$

where $[{ \varphi }_{j}^{i} = {\varphi }_{j} {\omega }_{j}^{i}]$ with $[{\omega }_{j}^{i}]$ as the jth spin function ( $[\alpha]$ or $[\beta]$ ) of the ith spin primitive function $[{\vartheta }_{i}]$ . Therefore, in equation (13), the antisymmetrizer applied to the spin orbitals product generates a Slater determinant and, consequently, the SC ansatz becomes

$[\eqalignno{{\Psi }_{0} &= {\cal N}\sum _{i = 1}^{{N}_{d}}{{{b}_{S,i}}\over{\sqrt{N!}}} \left|{\phi }_{1}^{c}{ \overline{\phi }}_{1}^{c}\ldots {\phi }_{j}^{c}{ \overline{\phi }}_{j}^{c}\ldots {\phi }_{{N}_{1}}^{c}{ \overline{\phi }}_{{N}_{1}}^{c}{\varphi }_{1}^{i}{ \varphi }_{2}^{i}\ldots {\varphi }_{{N}_{v} }^{i}\right|&\cr &= {\cal N}\sum _{i = 1}^{{N}_{d}}{b}_{S,i}{\Omega }_{i}. &(16)}]$

That is, it can be written as a linear combination of Slater determinants $[\{{\Omega }_{i}\}]$ . Furthermore, exploiting equation (16) and bearing in mind the Löwdin rules (Löwdin, 1955 ; McWeeny, 1992 ) for the computation of overlap integrals and matrix elements between Slater determinants constructed with non-orthogonal orbitals, it is easy to show that

$[\specialfonts{\cal N} = {{\frak D}}^{\rm -{{1}/{2}}}, \eqno(17)]$

where $[\specialfonts{\frak D}]$ is the zeroth-order supercofactor given by

$[\specialfonts{\frak D} = \textstyle\sum\limits _{i,j = 1}^{{N}_{d}}{b}_{S,i}{b}_{S,j}\, {\rm det}[{{\bf O}}_{ij}] \eqno(18)]$

with $[{{\bf O}}_{ij}]$ as the matrix of the overlap integrals between the spin orbitals of Slater determinants $[{\Omega }_{i}]$ and $[{\Omega }_{j}]$ appearing in equation (16). Therefore, the wavefunction ansatz for the reference crystal unit can be finally rewritten in this way:

$[\specialfonts\eqalignno{{\Psi }_{0}& = {{\frak D}}^{-{{1}/{2}}}\sum _{i = 1}^{{N}_{d}}{{{b}_{S,i}}\over{\sqrt{N!}}} \left|{\phi }_{1}^{c}{ \overline{\phi }}_{1}^{c}\ldots {\phi }_{j}^{c}{ \overline{\phi }}_{j}^{c}\ldots {\phi }_{{N}_{1}}^{c}{ \overline{\phi }}_{{N}_{1}}^{c}{\varphi }_{1}^{i}{ \varphi }_{2}^{i}\ldots {\varphi }_{{N}_{v} }^{i}\right|&\cr &= {{\frak D}}^{-{{1}/{2}}}\sum _{i = 1}^{{N}_{d}}{b}_{S,i}{\Omega }_{i}. &(19)}]$

2.3. The SC electron density

Exploiting wavefunction $[{\Psi }_{0}]$ , the electron density $[{\rho }_{0}({\bf r})]$ of the reference crystal unit can be simply obtained as

$[\eqalignno{{\rho }_{0}({\bf r})& = \textstyle\int\limits _{ {\bf x}_1 = {\bf x}} \,{\rm d}{\omega }_{1} \,{\rm d}{{\bf x}}_{2} \,{\rm d}{{\bf x}}_{3}\ldots \,{\rm d}{{\bf x}}_{N} {\Psi }_{0}({{\bf x}}_{1},{{\bf x}}_{2},\ldots, {{\bf x}}_{N})&\cr &\quad\times{\Psi }_{0}^{*}({{\bf x}}_{1},{{\bf x}}_{2},\ldots, {{\bf x}}_{N}). &(20)}]$

Substituting (19) in (20) and using again the Löwdin rules, we obtain

$[\specialfonts\eqalignno{{\rho }_{0}({\bf r}) &= {\rho }_{\rm core}({\bf r})+{{{\rho }_{\rm SC}({\bf r})}\over{{\frak D}}} &\cr &= 2 \sum _{t = 1}^{{N}_{1}} \left|{\phi }_{t}^{c}({\bf r})\right|^{2}+{{1}\over{{\frak D}}} \sum _{t,u = 1}^{{N}_{v}}{\varphi }_{t}({\bf r}){ \varphi }_{u}({\bf r}){\frak D}(t|u), &(21)}]$

where $[\specialfonts{\frak D}(t|u)]$ is the first-order supercofactor, which is defined as

$[\specialfonts{\frak D}(t|u) = \textstyle\sum\limits _{i,j = {\rm 1}}^{{N}_{d}}{b}_{S,i}{b}_{S,j} {I}_{tu}^{ij}({\rm -1})^{t+u}{\rm det}\left[{\bf O}_{ij}(t|u)\right] \eqno(22)]$

with I_tu^ij as the spin integral between the spin part of the tth spin orbital of the ith Slater determinant and the spin part of the uth spin orbital of the jth Slater determinant, and $[{{\bf O}}_{ij}(t|u)]$ as the (N − 1)th-order matrix obtained by deleting row t and column u of matrix $[{{\bf O}}_{ij}]$ , which was defined above (the more general definition of the rth-order supercofactors and the main properties of the different kinds of supercofactors are given in Appendix A1).

2.4. Implementation of the XCSC technique

Considering the analytical form of the new wavefunction ansatz [see particularly equation (6)] and of the active SC orbitals [see equation (10)], the new XCSC method consists in determining the coefficients $[\{{C}_{\mu i}\}]$ of the SC orbital expansions and the spin-coupling coefficients {c_S,k} that minimize the Jayatilaka functional shown in equation (2), which can be re-expressed in this way to stress the dependence of the functional on the coefficients that must be determined:

$[\eqalignno{J\left[\{{C}_{\mu i}\}, \{{c}_{S,k}\}\right]& = W\left[\{{C}_{\mu i}\}, \{{c}_{S,k}\}\right]&\cr &\quad +\lambda \left({\chi }^{2}\left[\{{C}_{\mu i}\}, \{{c}_{S,k}\}\right]-\Delta \right).& (23)}]$

To accomplish this task, we started from a working SC program (Cooper et al., 1993) where we have implemented a Newton–Raphson procedure based on the `quadratic hill climbing' algorithm proposed by Goldfeld et al. (1966 ). The algorithm requires the computation of analytical first and second derivatives of functional (23) with respect to SC orbitals and spin-coupling coefficients, whose mathematical derivations and expressions are shown in Appendix A2. Convergence is reached when the norm of the functional gradient is lower than 1×10^-6 and all the Hessian eigenvalues are positive.

3. Test calculations

3.1. Preliminary information

To better assess the performances of the new XCSC method, we have exploited the high-resolution X-ray diffraction data ( $[\sin\theta /\lambda]$ = 1.08 Å⁻¹) for salicylic acid collected at 90 K by Munshi & Guru Row (2006 ). In particular, the molecular geometry obtained from the X-ray diffraction experiment (see Fig. 1) has been exploited to carry out calculations at RHF, X-ray constrained restricted Hartree–Fock (XC-RHF), SC and XCSC levels. Furthermore, since one of the goals of our investigation was also to evaluate how the results obtained through the new XCSC strategy change as a function of the adopted basis sets, we have performed all the above-mentioned computations using the 6-31G, 6-311G, 6-31G(d) and 6-311G(d) sets of basis functions.

Figure 1
X-ray crystal structure of salicylic acid with labels for the carbon atoms of the aromatic ring and SC resonance structures taken into account in the unconstrained SC and XCSC calculations (A, D: Kekulé resonance structures; B, C and E: Dewar resonance structures).

For all the unconstrained and XCSC calculations, only the six π electrons of the aromatic ring have been considered as valence (active) electrons and, consequently, fully treated at the SC level [see the wavefunction ansatz (6)]. The other core (inactive) electrons have been described by frozen doubly occupied `core' MOs and always kept orthogonal to the (active) SC space [see equation (9)]. In particular, two kinds of core orbitals have been considered for our SC and XCSC computations. In one case [and in analogy with the original version of the XCSC approach (Genoni, Franchini et al., 2018)] we have used unconstrained RHF MOs, which allowed us to perform calculations that will be hereafter labelled as SC.0 and XCSC.0. In the other case we have used MOs that resulted from preliminary XC-RHF computations and that enabled us to carry out new types of SC and XCSC calculations. They will be labelled as SC.1 and XCSC.1, respectively. As we will discuss hereafter, this represents an important step forward compared with the previous version of the XCSC strategy (Genoni, Franchini et al., 2018) because the inactive part of the SC wavefunction also intrinsically takes into account the effect of the experimental diffraction data, thus leading to a more balanced description of the electronic structure.

For the sake of completeness, it is also worth noting that, due to their intrinsic delocalized nature, the RHF and XC-RHF occupied MOs have been preliminarily localized exploiting the Pipek–Mezey localization technique (Pipek & Mezey, 1989) to determine which of them to freeze for the description of the inactive electrons and which of them to use in the expansion of the active SC orbitals [see equation (10)].

Furthermore, since for a system of N_v active electrons in the spin state (S, M) the number of possible spin-coupling modes is given by

$[{f}_{S}^{{N}_{v}} = {{(2S+1){N}_{v}!}\over{\left({{1}\over{2}}{N}_{v}+S+1\right)! \left({{1}\over{2}}{N}_{v}-S\right)!}} \eqno(24)]$

and, since in our investigations we treated only singlet states, in all our SC and XCSC calculations we have considered five SC structures [see again equation (6)], which correspond to the five resonance structures depicted in Fig. 1, where structures A and D correspond to the traditional Kekulé resonance structures, while B, C and E are the Dewar ones.

Concerning the XC-RHF and all types of XCSC calculations we also exploited the lattice parameters (namely, unit-cell lengths and angles), ADPs and structure-factor amplitudes deposited with the considered crystallographic structures. Only structure-factor amplitudes characterized by $[|{F}_{{\bf h}}^{\rm exp}|\,\gt\, 3 {\sigma }_{{\bf h}}]$ were considered in the computations, which resulted in the selection of 6157 unique reflections. Furthermore, all the X-ray constrained computations (i.e. XC-RHF, XCSC.0 and XCSC.1) were carried out by gradually increasing the value of the external multiplier λ from 0 (unconstrained computations) and, as we will show below in more detail, to halt the XCSC calculations we have proposed and adopted a novel termination criterion based on the sign of the second derivative of χ² with respect to λ.

While unconstrained RHF and XC-RHF computations were carried out using the quantum chemistry package Gaussian09 (Frisch et al., 2009 ) and the quantum crystallography software TONTO (Jayatilaka & Grimwood, 2003 ), respectively, all kinds of XCSC calculations were performed by exploiting our in-house program mentioned in Section 2.4.

To evaluate the general capabilities of the XCSC technique, in the next subsections we will show the changes in the χ² statistical agreements, in the weights of the resonance structures and in the SC orbitals when the experimental X-ray diffraction data are taken into account in the SC calculations. Furthermore, we will also present how these changes are influenced by the adopted basis set.

For the sake of clarity, we will start discussing the results of the SC.0 and XCSC.0 computations. Afterwards we will focus on the outcomes of the SC.1 and XCSC.1 calculations. Finally, in the last part of this section, we will also show the effects of accounting for the experimental structure factors on the global electron-density distributions, both by visual inspection of density differences and through the evaluation of similarity indices between the obtained electron densities.

3.2. Calculations based on unconstrained RHF MOs

In the first validation tests, both SC and XCSC calculations were performed using the unconstrained RHF MOs to describe the inactive electrons. The χ² and energy values resulting from these computations are reported in Table 1 (`RHF MOs' column), where it is easy to observe that, for all the considered basis sets, the statistical agreement with the experimental data always slightly improves when the XCSC.0 computations are performed. Furthermore, as one would expect, the XCSC.0 description improves (always in terms of agreement with the experimental data) when larger and more flexible basis sets are taken into account, with the χ² value decreasing from 13.51 (6-31G basis set) to 8.07 [6-311G(d) basis set]. Nevertheless, it is also evident that, despite the application of the XCSC approach, for all the basis sets, the values of the statistical agreements remain quite high and close to those associated with the unconstrained RHF and SC.0 computations. The reason is that, in all the XCSC.0 calculations, only six electrons were actually treated at the XCSC level, while the remaining ones (namely, the inactive electrons) were described by doubly occupied MOs that resulted from unconstrained RHF computations and that, consequently, cannot take into account the effect of the experimental X-ray diffraction data. To reach lower values of the χ² statistical agreements by means of XCSC calculations, the only possibility is to exploit XC-RHF MOs to describe the core electrons. This has actually been done through the SC.1 and XCSC.1 computations, which will be discussed in the next subsection. For the sake of completeness, in Table 1 it is also possible to notice that, for all the basis sets, the energies associated with the XCSC.0 wavefunctions are higher than those corresponding to the SC.0 wavefunctions. As already discussed in our preliminary work on the XCSC method (Genoni, Franchini et al., 2018), this arises from the fact that in every XCW technique, the determination of the `experimental' wavefunction is accomplished without introducing new variational parameters. Consequently, the resulting SC orbitals and the spin-coupling coefficients represent a minimum point for the Jayatilaka functional [see equation (23)], but not for the energy of the system.

Table 1
χ² statistical agreement and energy values for all the unconstrained and X-ray constrained calculations performed on salicylic acid

The value of the external parameter λ is also reported for X-ray constrained calculations.

	Calculations with RHF MOs†			Calculations with XC-RHF MOs‡
Method and basis set	$[{\lambda }_{\max}]$	$[{\chi }^{2}]$	Energy (E_h)	$[{\lambda }_{\max}]$	$[{\chi }^{2}]$	Energy (E_h)
6-31G
RHF/XC-RHF		14.05	−492.975931	0.100	5.66	−492.817879
SC.0/SC.1		14.01	−493.037815		5.87	−492.884062
XCSC.0/XCSC.1	0.100	13.51	−493.020372	0.326	5.58	−492.861463

6-311G
RHF/XC-RHF		13.70	−493.088826	0.105	4.93	−492.919497
SC.0/SC.1		13.66	−493.150622		5.14	−492.984446
XCSC.0/XCSC.1	0.105	13.08	−493.129330	0.315	4.86	−492.961320

6-31G(d)
RHF/XC-RHF		9.34	−493.179441	0.150	3.59	−493.053881
SC.0/SC.1		9.14	−493.238706		3.74	−493.120653
XCSC.0/XCSC.1	0.150	8.56	−493.207299	0.610	3.52	−493.096062

6-311G(d)
RHF/XC-RHF		8.85	−493.285839	0.140	3.49	−493.172468
SC.0/SC.1		8.67	−493.344578		3.63	−493.237150
XCSC.0/XCSC.1	0.140	8.07	−493.314518	0.640	3.41	−493.209195

†Results for RHF, SC.0 and XCSC.0 calculations.
‡Results for XC-RHF, SC.1 and XCSC.1 calculations.

As already anticipated, for the termination of the XCSC calculations performed in this study we have adopted a new criterion. This point deserves to be discussed in more detail because, despite many suggestions having been proposed over the years, the determination of the exact value of λ at which to stop the X-ray constrained computations is still an open and debated problem. The new criterion simply consists in halting the X-ray constrained calculations when the curvature of the typical graph of χ² as a function of λ changes, namely when the second derivative $[{\partial }^{2}{\chi }^{2}/\partial {\lambda }^{2}]$ becomes negative (inflection point). In fact, by plotting χ² against λ for all the XCSC.0 calculations, we easily observed clear changes in the curvature of the graphs (see Figs S1–S4 in the supporting information) and, above all, we have also noted that the XCSC computations generally stopped converging a few steps after the inflection point. This observation was further corroborated by the fact that, for λ values only slightly larger than the detected inflection points, even the more traditional XC-RHF computations showed instabilities and difficulties in reaching convergence. In other words, the second derivative $[{\partial }^{2}{\chi }^{2}/\partial {\lambda }^{2}]$ can be considered as a useful descriptor to detect the onset of instability in the XCW procedure. Although at the moment it is only based on empirical observations, the new criterion seems quite promising and robust. In our opinion, it will deserve further investigations in the future by means of more extensive tests and statistical analyses, also applying it to traditional X-ray constrained calculations such as the XC-RHF ones. However, it is also worth pointing out that the use of the adopted termination criterion does not invalidate the new XCSC method and the assessment of its performance discussed in the present paper.

In Section 1 we also mentioned that one of the novelties introduced through the new XCSC techniques is the possibility of extracting resonance structure weights for the systems under examination directly from the collected experimental diffraction data. In Table 2, we have reported the Chirgwin–Coulson weights corresponding to resonance structures of salicylic acid in its singlet state as obtained from the SC.0 and XCSC.0 calculations. The trend is generally analogous for almost all the basis sets. Concerning Kekulé structures A and D, we can see that, while the weight of the former slightly increases in almost all the cases when experimental data are taken into account [the only exception is for basis set 6-31G(d)], the importance of the latter always decreases, with the largest variations observed for basis sets with polarization functions [Δ = −2.5 and Δ = −3.3 for basis sets 6-31G(d) and 6-311G(d), respectively]. Considering the remaining Dewar structures (namely, structures B, C and E), we can observe that structure B loses its weight when X-ray data are included in the calculations and, interestingly, all the X-ray constrained computations provided very similar values. In contrast, structure C is characterized by larger weights in the XCSC.0 calculations and also in this case the resulting values are very similar among the different basis sets. Structure E is the only one for which the trends are different for basis sets with and without polarization functions. In fact, its importance decreases for basis sets 6-31G and 6-311G and, conversely, increases for basis sets 6-31G(d) and 6-311G(d).

Table 2
Chirgwin–Coulson weights (in %) of the salicylic acid resonance structures resulting from the SC.0 and XCSC.0 calculations for all the considered basis sets

	6-31G		6-311G		6-31G(d)		6-311G(d)
Structures	SC.0	XCSC.0	SC.0	XCSC.0	SC.0	XCSC.0	SC.0	XCSC.0
A	16.6	17.2	16.6	18.4	15.7	15.7	15.7	16.9
B	8.8	7.6	8.7	7.6	9.3	7.5	9.1	7.8
C	9.5	12.2	9.4	12.3	9.7	12.9	9.6	12.4
D	54.7	52.9	55.0	52.6	53.8	51.3	54.3	51.0
E	10.4	10.1	10.2	9.1	11.5	12.6	11.3	11.9

To complete the analysis of the results obtained from the XCSC.0 calculations, we have finally decided to focus on the SC orbitals for the six active electrons (see Fig. 2) obtained from the X-ray constrained SC calculations performed on salicylic acid. In particular, we analysed their variations with respect to the corresponding unconstrained SC.0 orbitals, namely we evaluated the effect of the experimental diffraction data on the SC orbitals. From Fig. 3 it is clear that the orbital variations resulting from the calculations with basis sets without polarization functions are very similar to each other. Likewise, the orbital variations obtained from calculations with basis sets including polarization functions present a high degree of similarity.

Figure 2
Three-dimensional plots of the unconstrained SC.0 spin-coupled orbitals (0.06 e bohr⁻³ isosurfaces are plotted). The unconstrained SC.1 spin-coupled orbitals have analogous shapes.

Figure 3
Two-dimensional plots of the differences between the square moduli of the XCSC.0 and SC.0 spin-coupled orbitals for all the considered basis sets. The contours are drawn at $[\pm 1\times {10}^{-3}]$ e bohr⁻³ and at $[\pm 2, 4, 8\times {10}^{n}]$ e bohr⁻³ (with n as an integer ranging from −3 to 0) in a plane parallel to and 0.5 Å above the aromatic ring. Red and blue contours indicate positive and negative values, respectively.

Through a more detailed inspection, it is also possible to see that, for orbitals $[{\varphi }_{2}]$ , $[{\varphi }_{3}]$ and $[{\varphi }_{4}]$ , when basis sets 6-31G and 6-311G are used, the variations mainly consist in a slight delocalization from the atoms on which the unconstrained SC.0 orbitals are mainly localized to the neighbour atoms. Conversely, for the same orbitals, when basis sets 6-31G(d) and 6-311G(d) are exploited, we have an opposite trend, with the SC orbitals tending to further localize on the central atoms. This opposite behaviour can be ascribed to the effect of the experimental diffraction data. In fact, in the absence of polarization functions, the SC.0 orbitals $[{\varphi }_{2}]$ , $[{\varphi }_{3}]$ and $[{\varphi }_{4}]$ are less polarized towards the neighbour atoms and the effect of the X-ray diffraction data in the XCSC.0 computations is a delocalization over the neighbour carbon atoms. In contrast, in the presence of polarization functions, the starting SC.0 orbitals are already polarized towards the other atoms, probably at a higher level than that compatible with the collected experimental data. Therefore, accounting for the X-ray diffraction data in the XCSC.0 calculations, the polarization of the orbitals is tempered mainly through a shift of the electronic clouds from the bonding regions to the central carbon atom (especially for orbitals $[{\varphi }_{3}]$ and $[{\varphi }_{4}]$ ).

For orbitals $[{\varphi }_{1}]$ , $[{\varphi }_{5}]$ and $[{\varphi }_{6}]$ , the variations within the aromatic ring are quite similar for all the basis sets aside from the use of polarization functions in the calculations, even if for basis sets 6-31G(d) and 6-311G(d) these orbitals delocalize towards the carboxylic group when experimental data are taken into account. For orbitals $[{\varphi }_{1}]$ and $[{\varphi }_{5}]$ it is possible to notice stronger localizations on the carbon atoms on which the SC orbitals are mainly localized, with shifts of electron density from the neighbour bonding regions that are more pronounced when polarization functions are used. For orbital $[{\varphi }_{6}]$ the variation is more convoluted, with a shift of the electronic cloud mainly from the carbon atom bearing the hydroxyl group to the neighbour atoms. In this case the effect is more evident when basis sets without polarization functions are used, while it is barely sketched with basis sets 6-31G(d) and 6-311G(d).

3.3. Calculations based on X-ray constrained RHF MOs

To improve the description provided by the XCSC.0 calculations, we have subsequently performed a new series of XCSC computations (XCSC.1 computations) which exploit XC-RHF MOs to describe the inactive electrons.

Energy and χ² values resulting from these calculations are reported in Table 1 (`XC-RHF MOs' column) where it is immediately evident that the statistical agreements associated with the XCSC.1 wavefunctions are significantly better than the corresponding XCSC.0 ones. Furthermore, as for the SC.0 and XCSC.0 computations, the χ² values associated with the XCSC.1 wavefunctions are always slightly lower than the unconstrained SC.1 values and, as expected, the statistical agreement with the experimental data improves as richer basis sets are used, with χ² decreasing from 5.58 (6-31G basis set) to 3.41 [6-311G(d) basis set]. Finally, it is worth noting that, in Table 1, the SC.1 χ² values are always larger than the XC-RHF ones, which is in contrast with the results obtained from the calculations with unconstrained RHF MOs. This different behaviour can be explained considering that, in the SC.1 computations, only the core orbitals fully take into account the effect of the experimental data, while the SC ones (which describe the active electrons) are obtained by simply minimizing the energy of the system and not the Jayatilaka functional. Conversely, in the XC-RHF computation all the MOs result from the minimization of the functional given by equation (2) and, for this reason, they completely take into account the influence of the X-ray diffraction data. Anyway, as seen above, the SC.1 descriptions are improved by carrying out the XCSC.1 calculations, through which also the active SC orbitals are optimized by fully considering the constraints of the experimental data.

Concerning the termination of the XCSC.1 computations, we have used again the criterion based on the sign of second derivative $[{\partial }^{2}{\chi }^{2}/\partial {\lambda }^{2}]$ . In these cases, the curves representing the trend of χ² as a function of λ are much more regular and smoother than those associated with the XCSC.0 calculations, but they are still characterized by inflection points (see Figs S5–S8 in the supporting information), after which convergence rapidly becomes difficult or impossible.

As in the previous subsection, in Table 3 we report the obtained Chirgwin–Coulson weights for the resonance structures of salicylic acid, as obtained from the SC.1 and XCSC.1 calculations. The trends are almost identical for all the basis sets and are generally analogous to those resulting from the calculations based on unconstrained RHF MOs. Moreover, despite unavoidable discrepancies in the absolute values, the results show quite consistently how the resonance structure weights change when the experimental X-ray diffraction data are used as external constraints in the computations. In particular, Kekulé resonance structure A gains weight when X-ray data are taken into account, with the only exception of basis set 6-31G, for which a 0.3 drop has been observed. Conversely, the other Kekulé structure (namely, structure D) consistently loses its importance for all the considered basis sets, also with significant reductions compared with the SC.1 values [e.g. Δ = −4.3 for basis set 6-311G(d)]. Concerning the Dewar resonance structures, when experimental constraints are introduced in the calculations, the weight of B slightly decreases, while both C and E acquire importance, with the only exception of the 6-311G(d) calculations, for which the Chirgwin–Coulson coefficient for structure E decreases from 11.3 to 10.9.

Table 3
Chirgwin–Coulson weights (in %) of the salicylic acid resonance structures resulting from the SC.1 and XCSC.1 calculations for all the considered basis sets

	6-31G		6-311G		6-31G(d)		6-311G(d)
Structures	SC.1	XCSC.1	SC.1	XCSC.1	SC.1	XCSC.1	SC.1	XCSC.1
A	15.1	14.8	15.8	17.7	14.8	17.8	15.2	20.0
B	9.3	7.8	9.4	7.7	9.9	8.6	9.7	8.4
C	9.5	15.2	9.3	13.0	9.4	11.1	9.1	10.3
D	55.6	50.3	55.5	50.4	54.5	50.4	54.7	50.4
E	10.5	11.9	10.1	11.1	11.4	12.0	11.3	10.9

Also in this case we have completed the analysis of the SC.1 and XCSC.1 computations by analysing the obtained SC orbitals and, in particular, by focusing on the orbital variations when the X-ray data are directly taken into account in the calculations (see Fig. 4). In most cases the trends are analogous to those observed in Fig. 3 for the SC.0 and XCSC.0 calculations, even if with different intensities for the concentrations and depletions of electronic charge. Furthermore, as for the calculations based on unconstrained RHF MOs, it is possible to notice high similarities between the orbital variations resulting from the computations with the 6-31G and 6-311G basis sets and, at the same time, between the orbital variations obtained from the calculations exploiting basis sets that include polarization functions.

Figure 4
Two-dimensional plots of the differences between the square moduli of the XCSC.1 and SC.1 spin-coupled orbitals for all the considered basis sets. The contours are drawn at $[\pm 1\times {10}^{-3}]$ e bohr⁻³ and at $[\pm 2, 4, 8\times {10}^{n}]$ e bohr⁻³ (with n as an integer ranging from −3 to 0) in a plane parallel to and 0.5 Å above the aromatic ring. Red and blue contours indicate positive and negative values, respectively.

However, the most important and evident difference between the XCSC.1–SC.1 and XCSC.0–SC.0 orbital variations can be observed in the region of the carboxylic group, especially for orbitals $[{\varphi }_{1}]$ , $[{\varphi }_{5}]$ and $[{\varphi }_{6}]$ . We have previously pointed out that, using basis sets 6-31G(d) and 6-311G(d), the XCSC.0 calculations provided orbitals $[{\varphi }_{1}]$ , $[{\varphi }_{5}]$ and $[{\varphi }_{6}]$ characterized by significant charge redistributions on the carbon atom of the carboxylic group compared with the unconstrained computation (see Fig. 3). Conversely, in Fig. 4, we can observe that these charge redistributions disappear or significantly reduce when we consider the XCSC.1–SC.1 variations of the same orbitals. For orbital $[{\varphi }_{1}]$ , if we consider the results of the computations with basis sets 6-31G and 6-311G, we can even notice a slight depletion of electronic charge close to the carboxylic group, which is not detected when the XCSC.0–SC.0 orbital variations are considered. All these differences in the region of the carboxylic group actually derive from the different treatment of the inactive electrons in the XCSC calculations. In fact, in our computations, only the six π electrons of the aromatic ring were fully treated at the SC level, while the remaining ones (including those associated with the carboxylic group) were described by frozen doubly occupied MOs. In the XCSC.0 calculations, the MOs describing the electrons of the carboxylic group are obtained from unconstrained RHF computations and do not take into account the influence of the experimental X-ray diffraction data, which contain information on the full electron density of the system (therefore, also on the electron-density `basins' corresponding to the carboxylic group). Hence, the active SC orbitals (and the corresponding electron-density distributions, see next subsection) resulting from the XCSC.0 calculations delocalize towards the carboxylic group to compensate for the fact that the information contained in the X-ray diffraction data for the carboxylic group is not already captured by the frozen MOs used to describe the inactive electrons. In contrast, when XCSC.1 calculations are carried out, the inactive electrons are described by orbitals resulting from XC-RHF computations and, therefore, the information contained in the X-ray data for the inactive electrons is already taken into account. This is probably the reason why in the active XCSC.1 orbitals (and in the related electron densities) the degree of the delocalization towards the carboxylic group is significantly lower.

The previous observation clearly reveals that, unlike the original version of the XCSC method (Genoni, Franchini et al., 2018), the new XCSC technique based on XC-RHF MOs enables one to correctly deconvolute the X-ray constrained treatment of the inactive electrons from that of the active electrons, thus leading to more balanced descriptions of the electronic structures and to better electron-density distributions for the systems under examination.

3.4. Analysis of the obtained electron densities

Since X-ray structure factors used as constraints in XCW methods are strictly related to electron density, in our investigation we also analysed the charge-density distributions resulting from the performed calculations. In particular, we determined how the inclusion of the experimental structure factors in the XCSC computations influences the electron distributions. Therefore, we have evaluated and plotted the difference maps between corresponding XCSC and SC electron densities, both for calculations using unconstrained RHF MOs and for calculations using XC-RHF MOs. These maps are depicted in Fig. 5, where it can be easily noticed that the obtained differences are strictly related to the orbital variations shown in Figs. 3 and 4. For this reason, as for the orbital variations, in this case we also observe high similarities between the difference maps associated with basis sets without polarization functions and between those corresponding to basis sets including polarization functions. This is true both for the $[{\rho }_{\rm XCSC.0}-{\rho }_{\rm SC.0}]$ maps (see Fig. 5, left panel) and for the $[{\rho }_{\rm XCSC.1}-{\rho }_{\rm SC.1}]$ maps (see Fig. 5, right panel).

Figure 5
Two-dimensional plots of the differences between X-ray constrained and unconstrained SC electron densities for all the considered basis sets. The contours are drawn at $[\pm 1\times {10}^{-3}]$ e bohr⁻³ and at $[\pm 2, 4, 8\times {10}^{n}]$ e bohr⁻³ (with n as an integer ranging from −3 to 0) in a plane parallel to and 0.5 Å above the aromatic ring. Red and blue contours indicate positive and negative values, respectively.

Analysing Fig. 5 in more detail, the difference maps obtained with basis sets 6-31G and 6-311G are generally characterized by a charge-density accumulation on atoms C1 and C5 (see Fig. 1 for the labels of the atoms belonging to the aromatic ring) and in the C3–C4 bonding region. At the same time, they also present charge depletions on atoms C2, C3, C4 and C6. All these charge accumulations and depletions decrease in intensity when XC-RHF MOs are used in the SC and XCSC computations. Furthermore, we can observe that, while in the $[{\rho }_{\rm XCSC.0}-{\rho }_{\rm SC.0}]$ maps there is a slight accumulation of charge on the carbon atom belonging to the carboxylic group, in the $[{\rho }_{\rm XCSC.1}-{\rho }_{\rm SC.1}]$ maps it is actually possible to observe a slight charge depletion for the same atom. An analogous trend can be observed for atom C5 when basis set 6-311G is used. These trends and the observation that the intensities of charge accumulations and depletions decrease in the $[{\rho }_{\rm XCSC.1}-{\rho }_{\rm SC.1}]$ maps can be ascribed again to the fact that, in the SC.1 and XCSC.1 calculations, XC-RHF orbitals are used to describe the inactive electrons. As explained in the previous subsection, these orbitals already take into account the effects of the experimental data and enable one to deconvolute the treatment of the active electrons from that of the inactive electrons, thus leading to better descriptions of the electronic structures and to better electron densities.

However, unlike the difference densities obtained for basis sets without polarization functions, those corresponding to basis sets 6-31G(d) and 6-311G(d) are more similar, regardless of the use of XC-RHF orbitals in the computations. In particular, charge accumulations are mainly observed for almost all the atoms belonging to the aromatic ring and for the carbon atom of the carboxylic group, while charge depletions are seen in the bonding regions and for atom C6. However, also in this case, in the $[{\rho }_{\rm XCSC.1}-{\rho }_{\rm SC.1}]$ maps we have noted a slight decrease of the intensities (especially for the carboxylic group region), which is again due to the use of X-ray constrained orbitals to treat the inactive electrons.

Finally, to perform a more global comparison, we have also computed similarity indices between the obtained electron densities. In particular, we have adopted the Walker–Mezey similarity measure $[L({\rho }_{x},{\rho }_{y},a,{a}')]$ (Walker & Mezey, 1994 ), which allows the comparison of two charge distributions $[{\rho }_{x}]$ and $[{\rho }_{y}]$ in a three-dimensional region bound by two isosurfaces characterized by the values a and a′ (hereafter indicated in e bohr⁻³; see the supporting information for more details about this similarity index). By changing a and a′, the similarities between the electron densities can be evaluated in different three-dimensional basins. Therefore, for all our comparisons, two Walker–Mezey similarity indices have been taken into account: (i) the $[L({\rho }_{x},{\rho }_{y}, 0.1, 10)]$ index to compare the electron densities in regions close to the nuclei (core regions); (ii) the $[L({\rho }_{x},{\rho }_{y}, 0.001, 0.1)]$ index to compare the charge distributions in the valence regions.

At first, we focused on the similarities between all the XCSC.1 electron densities. Both in the core and valence regions (see Tables 4 and 5, respectively), the Walker–Mezey similarity index confirms what was already observed in Fig. 5, namely the fact that charge distributions obtained with the same kind of basis set (i.e. with or without polarization functions) are generally more similar to each other. Completely analogous results have been obtained for the XCSC.0 electron distributions (see Tables S1 and S2 in the supporting information).

Table 4
Values of the Walker–Mezey similarity index $[L({\rho }_{x},{\rho }_{y}, 0.1, 10)]$ corresponding to the comparisons of the XCSC.1 electron densities obtained with the different basis sets

$[L({\rho }_{x},{\rho }_{y}, 0.1, 10)]$	6-31G	6-311G	6-31G(d)	6-311G(d)
6-31G	100.00
6-311G	99.24	100.00
6-31G(d)	97.31	97.52	100.00
6-311G(d)	97.44	97.71	99.47	100.00

Table 5
Values of the Walker–Mezey similarity index $[L({\rho }_{x},{\rho }_{y}, 0.001, 0.1)]$ corresponding to the comparisons of the XCSC.1 electron densities obtained with the different basis sets

$[L({\rho }_{x},{\rho }_{y}, 0.001, 0.1)]$	6-31G	6-311G	6-31G(d)	6-311G(d)
6-31G	100.00
6-311G	95.38	100.00
6-31G(d)	93.36	92.09	100.00
6-311G(d)	92.47	93.57	95.89	100.00

Then we compared the complete set of electron densities obtained for each basis set. For the sake of clarity, we here report and discuss only the results obtained for basis set 6-31G(d) (see Tables 6 and 7), but completely analogous results have also been obtained for the other sets of basis functions (see Tables S3–S8 in the supporting information). In Table 6, which shows the values for similarity index $[L({\rho }_{x},{\rho }_{y}, 0.1, 10)]$ , we can see that all the electron distributions resulting from calculations based on unconstrained RHF MOs (RHF, SC.0 and XCSC.0) are very similar (values greater than 99%). This also holds true for the electron densities associated with computations based on XC-RHF molecular orbitals (XC-RHF, SC.1 and XCSC.1). This observation is obviously in line with the fact that, for the two different families of computations, the real core electrons (which are a subset of the inactive electrons) have been treated in exactly the same way.

Table 6
Values of the Walker–Mezey similarity index $[L({\rho }_{x},{\rho }_{y}, 0.1, 10)]$ corresponding to the comparisons of all the electron densities obtained with basis set 6-31G(d)

$[L({\rho }_{x},{\rho }_{y}, 0.1, 10)]$	RHF	SC.0	XCSC.0	XC-RHF	SC.1	XCSC.1
RHF	100.00
SC.0	99.61	100.00
XCSC.0	99.15	99.34	100.00
XC-RHF	97.04	97.25	97.56	100.00
SC.1	97.16	97.39	97.57	99.65	100.00
XCSC.1	96.88	97.09	97.45	99.75	99.47	100.00

Table 7
Values of the Walker–Mezey similarity index $[L({\rho }_{x},{\rho }_{y}, 0.001, 0.1)]$ corresponding to the comparisons of all the electron densities obtained with basis set 6-31G(d)

$[L({\rho }_{x},{\rho }_{y}, 0.001, 0.1)]$	RHF	SC.0	XCSC.0	XC-RHF	SC.1	XCSC.1
RHF	100.00
SC.0	98.59	100.00
XCSC.0	96.52	97.57	100.00
XC-RHF	93.50	94.23	93.84	100.00
SC.1	94.19	94.81	93.24	98.63	100.00
XCSC.1	92.79	93.53	93.91	98.77	97.63	100.00

Finally, in Table 7, we show the values for index $[L({\rho }_{x},{\rho }_{y}, 0.001, 0.1)]$ , which measures the degree of similarity in the valence region. Also, in this case, despite lower values for the similarity indices, it is easy to distinguish the two different groups of electron densities (in one group: RHF, SC.0 and XCSC.0; in the other one: XC-RHF, SC.1 and XCSC.1). It is interesting to notice that, for the first group, the similarity with the RHF electron density clearly decreases from SC.0 to XCSC.0. This is a clear indication of the different way in which we determined the orbitals that describe the active electrons for the system under examination (i.e. the six π electrons of the aromatic ring). In fact, in the RHF and SC.0 calculations, these orbitals have been obtained by simply minimizing the energy of the molecule, while in the XCSC.0 case they resulted from the minimization of the Jayatilaka functional. Conversely, in the second group of electron distributions, the similarity with the XC-RHF charge slightly increases passing from SC.1 to XCSC.1. This is due to the fact that, in the XC-RHF and XCSC.1 wavefunctions, the active electrons are described by orbitals that include the effect of the experimental X-ray diffraction data, while in the SC.1 computations the SC orbitals for the active electrons are obtained by only minimizing the energy functional, without introducing the constraint of the experimental structure factors. The effect of this different treatment is also obviously reflected in the lower value of the similarity index between the SC.1 and XCSC.1 electron-density distributions.

4. Conclusions and perspectives

In this paper, we have presented the novel XCSC method, a technique that introduces the typical high chemical interpretability of the SC strategy of quantum chemistry into the XCW approach of quantum crystallography. In particular, we have shown the fundamental equations for the new method and we have discussed the results of further test calculations that were performed to better evaluate the capabilities of the proposed strategy.

The tests have revealed that the XCSC computations are quite straightforward and that the new strategy does not present particular convergence problems compared with the corresponding unconstrained SC technique. In analogy with results obtained by means of other XCW methods, our calculations have also shown that the novel strategy provides better results and better statistical agreements with the experimental X-ray diffraction data as larger and more flexible basis sets are used. Furthermore, as we initially expected, the new XCSC technique has always provided better (i.e. lower) χ² values compared with all the other strategies taken into account for the comparison.

One of the advantages of the novel XCSC method is that it extracts the weights of the resonance structures for the system under investigation which are compatible with the X-ray structure factors observed experimentally. Notwithstanding some unavoidable variations in the absolute values, the XCSC calculations generally provided resonance weights that changed consistently and with very similar trends when the X-ray diffraction data have been considered, practically regardless of the basis set used.

Furthermore, to evaluate the effects of including the experimental structure factors in the calculations, we have also compared the obtained SC orbitals and electron densities, particularly focusing on the direct comparison between corresponding unconstrained and X-ray constrained quantities. The results indicate that the XCSC procedure entails clear redistributions of the electron density, which are generally similar for all the performed calculations. However, they sometimes also show some discrepancies depending on the adopted basis set and, above all, on the treatment of the inactive electrons in the computations. In fact, compared with the original version of the XCSC method (Genoni, Franchini et al., 2018), in this paper we have introduced a further advancement of the technique consisting in using MOs obtained from a preliminary XC-RHF calculation to describe the inactive electrons in the XCSC procedure. The analysis of the obtained orbitals and electron densities has highlighted that the use of this new strategy enables one to better deconvolute the description of the inactive electrons from that of the active electrons, which consequently gives better and more balanced descriptions of the electronic structures for the systems under examination. Therefore, this will definitely be the procedure to follow for future XCSC computations and, for this reason, the new XCSC technique can also be fully classified as the first post-XC-RHF method.

Finally, in the present study, we have proposed a new possible criterion to stop the XCW procedure, which consists in halting the computations when a clear inflection point is detected in the graph representing the trend of the χ² statistical agreement as a function of the external multiplier λ. So far, the new criterion has been tested only for XCSC calculations and, at least for them, it seems robust and consistent. Further computational tests will obviously be necessary on other types of XCW methods to test its general applicability.

In conclusion, in this paper we have introduced a new XCW strategy able to improve the description provided by the traditional X-ray constrained Hartree–Fock approach. Moreover, one of the main advantages of the novel strategy is that traditional chemical information (i.e. resonance structure weights, but also spatial distributions of the electronic clouds around atoms) can be directly obtained without imposing any preliminary constraint a priori or applying techniques a posteriori to the obtained wavefunction. However, notwithstanding the promising results obtained with the new method, here we do not claim the superiority of our novel approach over the existing ones. On the contrary, we want to point out that the proposed XCSC strategy will be only one of the different tools currently available in quantum crystallography to investigate chemical/physical problems/phenomena in the solid state from different perspectives. In fact, each technique has its own features and provides specific answers. Therefore, as already highlighted by Fugel, Beckmann et al. (2018), we believe that only the application of different and complementary methods can lead to a more global and complete description of the physical reality.

Concerning the future perspectives of the method, two aspects should be considered. First of all, it is worth pointing out that the proposed XCSC technique is a real many-determinant XCW strategy. For this reason, it might be fruitfully used to obtain two-particle density matrices from experimental scattering data and also to shed further light on the capability of the XCW approach to capture electron correlation effects, a problem that, so far, has been investigated only by inspecting how the wavefunction fitting reflects on the electron density (Genoni et al., 2017 ; Grabowsky, 2017 ), but without taking into account quantities more strictly related to electron correlation (e.g. two-particle quantities such as the intracule densities). In the context of these investigations, we are also planning to evaluate the capabilities of the XCSC method when the external constraints are only the low-angle reflections, which are those reflections that actually contain information on the electron correlation effects (Genoni et al., 2017). Finally, considering that SC wavefunctions are exact spin eigenfunctions also in the case of open-shell systems, another tantalizing future perspective for the XCSC approach is its extension to the joint refinement of X-ray diffraction and polarized neutron diffraction data, with the final goal of obtaining experimental spin densities, which could be analysed through advanced quantum chemical topology techniques (Gatti et al., 2015 , 2017 ; Macetti et al., 2018 ). This will be in line with techniques already developed in this context, such as those proposed in the framework of the multipole models (Deutsch et al., 2012 , 2014 ; Voufack et al., 2017 ) and those already based on a wavefunction ansatz (Gueddida et al., 2018).

APPENDIX A

Further theoretical details

A1. Higher-order supercofactors

As we will show in Section A2, supercofactors play a fundamental role in the XCSC method for the computation of the first and second derivatives of the Jayatilaka functional and can actually be defined in a very general way for the generic order r:

$[\specialfonts\eqalignno{&{\frak D}({t}_{1} {t}_{2}\ldots { t}_{r}| {u}_{1} {u}_{2}\ldots { u}_{r}) &\cr &= \textstyle\sum\limits _{i,j = 1}^{{N}_{d}}{b}_{S,i}{b}_{S,j} {I}_{{t}_{1}{u}_{1}}^{ij } {I}_{{t}_{2}{u}_{2}}^{ij }\ldots {I}_{{t}_{r}{u}_{r}}^{ij } {\cal{E}}_{{t}_{1}{t}_{2}\ldots {t}_{r}}{\cal{E}}_{{u}_{1}{u}_{2}\ldots {u}_{r}} (-1)^{\textstyle\sum _{k = 1}^{r}{t}_{k}+{u}_{k}}&\cr &\quad\times {\rm det}\left[{{\bf O}}_{ij}({t}_{1}{t}_{2}\ldots {t}_{r}|{u}_{1}{u}_{2}\ldots {u}_{r})\right], &(25)}]$

where the symbols {I_{t_ku_k}^ij} are analogous to those in equation (22), $[{{\bf O}}_{ij}({t}_{1}{t}_{2}\ldots {t}_{r}|{u}_{1}{u}_{2}\ldots {u}_{r})]$ is the (N − r)th-order matrix obtained by deleting rows t₁, t₂, …, t_r and columns u₁, u₂, …, u_r of matrix $[{{\bf O}}_{ij}]$ and $[{\cal{E}}_{{v}_{1}{v}_{2}\ldots {v}_{r}}]$ is a kind of Levi-Civita symbol related to the number of permutations that lead from the collection of integer numbers ( $[{v}_{1}{v}_{2}\ldots {v}_{r})]$ to the ordered collection ( $[{v}_{1}'{v}_{2}'\ldots {v}_{r}']$ ) with $[{v}_{1}'\,\lt\, v_{2}'\,\lt\, \ldots \,\lt\, {v}_{r}']$ . In particular, the symbol is defined as

$[{\cal{E}}_{{v}_{1}{v}_{2}\ldots {v}_{r}} = \left\{\matrix{+1 &{\rm for\ an\ even\ number\ of\ permutations }\cr -1 & {\rm for\ an\ odd\ number\ of\ permutations}.\hfill }\right. \eqno(26)]$

By applying the Laplace theorem for the computation of the determinants, it is possible to show that every supercofactor of order r can be expressed in terms of the supercofactors of order r + 1 through the following recurrence relation:

$[\specialfonts\eqalignno{&{\frak D}({t}_{1} {t}_{2}\ldots { t}_{r}| {u}_{1} {u}_{2}\ldots { u}_{r})&\cr &= \textstyle\sum\limits _{w\ne {t}_{1},{t}_{2},\ldots, {t}_{r}}^{N}\langle{\phi }_{w} | {\phi }_{v}\rangle {\frak D}({t}_{1} {t}_{2}\ldots { t}_{r}\ w | {u}_{1} {u}_{2}\ldots { u}_{r} v),& (27)}]$

where $[{\phi }_{v}]$ and $[{\phi }_{w}]$ can be `core' or SC orbitals in wavefunction (19), while v is a dummy index that can be arbitrarily chosen provided that it is different from the indices appearing in the `ket part' of the rth-order supercofactor. For the development of the traditional and the new XCSC methods, it is also important to introduce two other classes of supercofactors: (i) the indexed supercofactors and (ii) the symmetrized indexed supercofactors. For the generic order r, the former are defined as

$[\specialfonts\eqalignno{&{\frak D}_k (t_1 t_2 \ldots t_r | u_1 u_2 \ldots u_r ) \cr &= \textstyle\sum\limits _{i,j = 1}^{N_d} b_{S,i} d_{j,k} I_{t_1 u_1}^{ij} I_{t_2 u_2}^{ij} \ldots I_{t_r u_r}^{ij} {\cal E}_{t_1 t_2 \ldots t_r} {\cal E}_{u_1 u_2 \ldots u_r} (-1)^{\textstyle\sum _{k=1}^r t_k + u_k} \cr &\quad \times {\rm det} \left [ {\bf O}_{ij} (t_1 t_2 \ldots t_r | u_1 u_2 \ldots u_r )\right ], &(28)}]$

where d_j,k represents the weight of the jth Slater determinant in the kth SC structure. Consequently, always for the generic rth order, the latter are

$[\specialfonts\eqalignno{&{\frak D}_{k}'' ({t}_{1} {t}_{2}\ldots { t}_{r}| {u}_{1} {u}_{2}\ldots { u}_{r})&\cr &= {{\frak D}}_{k}({t}_{1} {t}_{2}\ldots { t}_{r}| {u}_{1} {u}_{2}\ldots { u}_{r})+{{\frak D}}_{k}({u}_{1} {u}_{2}\ldots { u}_{r}| {t}_{1} {t}_{2}\ldots { t}_{r}).&\cr &&(29)}]$

It is easy to show that this relation exists between ordinary supercofactors and indexed supercofactors:

$[\specialfonts{\frak D}{}({t}_{\rm 1} {t}_{\rm 2}\ldots { t}_{r}| {u}_{\rm 1} {u}_{\rm 2}\ldots { u}_{r}) = \textstyle\sum\limits _{k = {\rm 1}}^{{f}_{S}^{{N}_{v}}}{c}_{S,k}{{\frak D}}_{k}({t}_{\rm 1} {t}_{\rm 2}\ldots { t}_{r}| {u}_{\rm 1} {u}_{\rm 2}\ldots { u}_{r}). \eqno(30)]$

Furthermore, the following two other properties are also valid and exploited in the derivation of the working equations for the XCSC technique:

$[\specialfonts{\frak D}\left({\hat Q}[{t}_{\rm 1} {t}_{\rm 2}\ldots { t}_{r}] | {\hat Q} [{u}_{\rm 1} {u}_{\rm 2}\ldots { u}_{r}]\right) = {\frak D}({t}_{\rm 1} {t}_{\rm 2}\ldots { t}_{r}| {u}_{\rm 1} {u}_{\rm 2}\ldots { u}_{r}), \eqno(31)]$

where $[{\hat Q}]$ is a generic permutation of the symmetric group of order r and

$[\specialfonts{\frak D}({u}_{\rm 1} {u}_{\rm 2}\ldots { u}_{r}| {t}_{\rm 1} {t}_{\rm 2}\ldots { t}_{r}) = {\frak D}({t}_{\rm 1} {t}_{\rm 2}\ldots { t}_{r}| {u}_{\rm 1} {u}_{\rm 2}\ldots { u}_{r}). \eqno(32)]$

However, while equation (31) is valid for every class of supercofactors, relation (32) can be used only for the ordinary and the symmetrized indexed ones.

A2. Derivatives of the $[{{\boldchi }}^{2}]$ statistical agreement

To obtain the derivatives of functional (23), we clearly need both the derivatives of its energy part W and the derivatives of the statistical agreement between experimental and theoretical structure-factor amplitudes, where the latter must obviously be multiplied by λ and added to the former. The energy derivatives are identical to those already exploited in the traditional SC method (Cooper et al., 1993) and, therefore, in this subsection, we will only focus on the formal derivation of the first and second derivatives of χ² with respect to the coefficients of the SC orbitals and with respect to the spin-coupling coefficients.

Exploiting the fact that $[|{F}_{{\bf h}}^{\rm calc}| = [{F}_{\bf h}^{\rm calc} ({F}_{\bf h}^{\rm calc})^{*}]^{1/2}]$ , the first derivative of χ² with respect to a generic variable x can be written as

$[{{\partial {\chi }^{2}}\over{\partial x}} = \sum _{{\bf h}}{{{K}_{{\bf h}}}\over{2}} \left[{{\partial {F}_{{\bf h}}^{\rm calc}}\over{\partial x}} {\left({F}_{{\bf h}}^{\rm calc}\right)}^{*}+{F}_{{\bf h}}^{\rm calc}{{\partial {\left({F}_{{\bf h}}^{\rm calc}\right)}^{*}}\over{\partial x}}\right], \eqno(33)]$

where x can be a generic SC orbital coefficient $[{C}_{\mu k}]$ or a generic spin-coupling coefficient c_s,k and where

$[{K}_{{\bf h}} = {{2\eta }\over{{N}_{r}-{N}_{p}}} {{\eta \left|{F}_{{\bf h}}^{\rm calc}\right|-\left|{F}_{{\bf h}}^{\rm exp}\right|}\over{{\sigma }_{{\bf h}}^{2} \left|{F}_{{\bf h}}^{\rm calc}\right|}}. \eqno(34)]$

Furthermore, since we assume that we are working with real orbitals, both for the first derivatives with respect to $[{C}_{\mu k}]$ and for the first derivatives with respect to c_s,k, it is possible to show that

$[{{\partial {\left({F}_{{\bf h}}^{\rm calc}\right)}^{*}}\over{\partial x}} = {\left({{\partial {F}_{{\bf h}}^{\rm calc}}\over{\partial x}}\right)}^{*} \eqno(35)]$

and, consequently, equation (33) can be simply rewritten as

$[\eqalignno{{{\partial {\chi }^{2}}\over{\partial x}} &= \sum _{{\bf h}}{K}_{{\bf h}}\Biggl[{\rm Re}\left\{{F}_{{\bf h}}^{\rm calc}\right\} {\rm Re}\left\{{{\partial {F}_{{\bf h}}^{\rm calc}}\over{\partial x}}\right\} &\cr &\quad+ {\rm Im}\left\{{F}_{{\bf h}}^{\rm calc}\right\} {\rm Im}\left\{{{\partial {F}_{{\bf h}}^{\rm calc}}\over{\partial x}}\right\}\Biggr]. &(36)}]$

Now, introducing the definition of the structure-factor operator

$[{\hat I}_{{\bf h}} = \textstyle\sum\limits _{j = 1}^{{N}_{m}}\exp[i2\pi ({{\bf Q}}_{j}{\bf r}+ {{\bf q}}_{j}) \cdot {\bf Bh}] = {\hat I}_{{\bf h}, {\rm R}}+i {\hat I}_{{\bf h},{\rm C}}, \eqno(37)]$

where $[{\hat I}_{{\bf h}, {\rm R}}]$ and $[{\hat I}_{{\bf h},{\rm C}}]$ (real and imaginary parts of $[{\hat I}_{{\bf h}}]$ , respectively) are Hermitian operators, and exploiting equations (1) and (21), we obtain

$[\specialfonts\eqalignno{{F}_{{\bf h}}^{\rm calc}& = {F}_{{\bf h}}^{\rm core}+{{{F}_{{\bf h}}^{\rm SC}}\over{{\frak D}}} &\cr &= 2\sum _{t = 1}^{{N}_{1}}\langle {\phi }_{t}^{c} |{ \hat I}_{{\bf h}} |{ \phi }_{t}^{c}\rangle +{{1}\over{{\frak D}}} \sum _{t,u = 1}^{{N}_{v}}\langle {\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle {\frak D}(t|u). &(38)}]$

Therefore, considering that $[{F}_{{\bf h}}^{\rm core}]$ depends neither on the SC orbitals nor on the spin-coupling coefficients, the first derivatives of the calculated structure factors can be expressed in this way:

$[\specialfonts{{\partial {F}_{{\bf h}}^{\rm calc}}\over{\partial x}} = -{{{F}_{{\bf h}}^{\rm SC}}\over{{{\frak D}}^{\rm 2}}} {{\partial {\frak D}}\over{\partial x}}+ {{\rm 1}\over{{\frak D}}} {{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial x}} \eqno(39)]$

and the expression for the first derivatives of the statistical agreement becomes

$[\specialfonts\eqalignno{{{\partial {\chi }^{2}}\over{\partial x}} &= \sum _{{\bf h}}{K}_{{\bf h}}\Biggl[{\rm Re}\{{F}_{{\bf h}}^{\rm calc}\}\left(-{{{\rm Re}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial x}}+ {{1}\over{{\frak D}}} {\rm Re}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial x}}\right\}\right)&\cr &\quad +{\rm Im}\{{F}_{{\bf h}}^{\rm calc}\}\left(-{{{\rm Im}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial x}}+ {{1}\over{{\frak D}}} {\rm Im}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial x}}\right\}\right)\Biggr], &\cr &&(40)}]$

where we have exploited the fact that $[\specialfonts{\frak D}]$ is a real quantity.

From equation (40) it clearly follows that, in order to compute the first derivatives of χ², other than the calculated structure factors given by equation (38), it is also necessary to know the first derivatives $[\specialfonts\partial {\frak D}/\partial x]$ and $[\specialfonts\partial {F}_{{\bf h}}^{\rm SC}/\partial x]$ .

First of all, let us consider the first derivatives with respect to the coefficients of the SC orbitals and, for the sake of example, let us consider in detail the mathematical derivation of the analytical expression for $[\specialfonts\partial {F}_{{\bf h}}^{\rm SC}/\partial {C}_{\mu k}]$ . To this purpose, let us expand $[{F}_{{\bf h}}^{\rm SC}]$ in this way:

$[\specialfonts\eqalignno{{F}_{{\bf h}}^{\rm SC} &= \textstyle\sum\limits _{t\ne k}^{{N}_{v}}\sum\limits _{u\ne k}^{{N}_{v}}\langle {\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle{\frak D}(t|u)+\sum\limits _{u\ne k}^{{N}_{v}}\langle {\varphi }_{k}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle{\frak D}(k|u)&\cr &\quad +\textstyle\sum\limits _{t\ne k}^{{N}_{v}}\langle {\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{k}\rangle{\frak D}(t|k)+\langle {\varphi }_{k}|{ \hat I}_{{\bf h}} |{ \varphi }_{k}\rangle{\frak D}(k|k). &(41)}]$

By using the properties of supercofactors [particularly, equation (32)] and bearing in mind that we assume we are working with real orbitals, we easily obtain

$[\specialfonts\eqalignno{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial {C}_{\mu k}}}& = \sum _{t\ne k}^{{N}_{v}}\sum _{u\ne k}^{{N}_{v}}\langle {\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle {{\partial {\frak D}(t|u)}\over{\partial {C}_{\mu k}}} +\sum _{u\ne k}^{{N}_{v}}\langle {\varphi }_{k}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle {{\partial {\frak D}(k|u)}\over{\partial {C}_{\mu k}}}&\cr &\quad +\sum _{t\ne k}^{{N}_{v}}\langle {\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{k}\rangle {{\partial {\frak D}(t|k)}\over{\partial {C}_{\mu k}}}+2\sum _{u = 1}^{{N}_{v}}\langle {\chi }_{\mu }|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle{\frak D}(k|u). &\cr &&(42)}]$

Now, exploiting recurrence relation (27), $[\specialfonts{\frak D}(t|u)]$ can be simply expressed as

$[\specialfonts{\frak D}(t|u) = \textstyle\sum\limits _{r\ne u}^{{N}_{v}}\langle{\varphi }_{k}|{\varphi }_{r}\rangle{\frak D}(tk|ur) \eqno(43)]$

and $[\specialfonts\partial {\frak D}(t|u)/\partial {C}_{\mu k}]$ becomes

$[\specialfonts\eqalignno{{{\partial {\frak D}(t|u)}\over{\partial {C}_{\mu k}}}& =\langle {\chi }_{\mu }|{ \varphi }_{k}\rangle {\frak D}(tk|uk)+\langle{\varphi }_{k}|{\chi }_{\mu }\rangle{\frak D}(tk|uk)&\cr &\quad +\sum _{r\ne u,k}^{{N}_{v}}\langle{\chi }_{\mu }|{ \varphi }_{r}\rangle{\frak D}(tk|ur)+\sum _{r\ne u,k}^{{N}_{v}}\langle{\varphi }_{k}|{\varphi }_{r}\rangle {{\partial {\frak D}(tk|ur)}\over{\partial {C}_{\mu k}}}.&\cr &&(44)}]$

Using again recurrence relation (27) for $[\specialfonts{\frak D}(tk|ur)]$ , we have

$[\specialfonts{\frak D}(tk|ur) = \textstyle\sum\limits _{s\ne t,k}^{{N}_{v}}\langle{\varphi }_{s}|{\varphi }_{k}\rangle {\frak D}(tks|urk)\eqno (45)]$

and, therefore, equation (44) can be rewritten like this:

$[\specialfonts\eqalignno{{{\partial {\frak D}(t|u)}\over{\partial {C}_{\mu k}}}& =\langle {\chi }_{\mu }|{ \varphi }_{k}\rangle{\frak D}(tk|uk)+\langle{\varphi }_{k}|{\chi }_{\mu }\rangle{\frak D}(tk|uk)&\cr &\quad +\sum _{r\ne u,k}^{{N}_{v}}\langle{\chi }_{\mu }|{ \varphi }_{r}\rangle{\frak D}(tk|ur)&\cr &\quad +\sum _{s\ne t,k}^{{N}_{v}}\langle{\varphi }_{s}|{\chi }_{\mu }\rangle\sum _{r\ne u,k}^{{N}_{v}}\langle{\varphi }_{k}|{\varphi }_{r}\rangle{\frak D}(tks|urk). &(46)}]$

Now, since $[\specialfonts{\frak D}(tks|urk)]$ is equivalent to $[\specialfonts{\frak D}(tsk|ukr)]$ [see equation (31)], by exploiting the recurrence relation for supercofactors in the reverse way, we obtain

$[\specialfonts\textstyle\sum\limits _{r\ne u,k}^{{N}_{v}}\langle{\varphi }_{k}|{\varphi }_{r}\rangle{\frak D}(tsk|ukr) = {\frak D}(ts|uk). \eqno(47)]$

Therefore, introducing (47) into (46), we can write

$[\specialfonts{{\partial {\frak D}(t|u)}\over{\partial {C}_{\mu k}}} = \sum _{r\ne u}^{{N}_{v}}\langle{\chi }_{\mu }|{ \varphi }_{r}\rangle {\frak D}(tk|ur)+\sum _{s\ne t}^{{N}_{v}}\langle{\chi }_{\mu }|{\varphi }_{s}\rangle {\frak D}(ts|uk). \eqno(48)]$

Furthermore, always exploiting recurrence relation (27), the other derivatives of the first-order supercofactor appearing in (42) can be expressed as

$[\specialfonts{{\partial {\frak D}(k|u)}\over{\partial {C}_{\mu k}}} = \sum _{r\ne k}^{{N}_{v}}\langle{\chi }_{\mu }|{ \varphi }_{r}\rangle{\frak D}(kr|uk) \eqno(49)]$

and

$[\specialfonts{{\partial {\frak D}(t|k)}\over{\partial {C}_{\mu k}}} = \sum _{r\ne u}^{{N}_{v}}\langle{\chi }_{\mu }|{ \varphi }_{r}\rangle{\frak D}(tk|kr). \eqno(50)]$

Now, using relations (48), (49) and (50) in (42), we obtain

$[\specialfonts\eqalignno{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial {C}_{\mu k}}} &= \sum _{t\ne k}^{{N}_{v}}\sum _{u = 1}^{{N}_{v}}\sum _{r\ne u}^{{N}_{v}}\langle {\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle \langle{\chi }_{\mu }|{ \varphi }_{r}\rangle{\frak D}(tk|ur)&\cr &\quad +\sum _{t = 1}^{{N}_{v}}\sum _{u\ne k}^{{N}_{v}}\sum _{r\ne t}^{{N}_{v}}\langle {\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle \langle{\chi }_{\mu }|{ \varphi }_{r}\rangle{\frak D}(tr|uk)&\cr &\quad +2\sum _{t = 1}^{{N}_{v}}\langle {\chi }_{\mu }|{ \hat I}_{{\bf h}} |{ \varphi }_{t}\rangle {\frak D}(k|t). &(51)}]$

Finally, considering that $[\specialfonts{\frak D}(tr|uk)]$ is equivalent to $[\specialfonts{\frak D}(uk|tr)]$ [see relation (32)] and exchanging the labels t and u in the second term of the right-hand side of the previous equation, the expression of $[\partial {F}_{{\bf h}}^{\rm SC}/\partial {C}_{\mu k}]$ can be simply written in this way:

$[\specialfonts\eqalignno{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial {C}_{\mu k}}}& = 2\Biggl[\sum _{t\ne k}^{{N}_{v}}\sum _{u = 1}^{{N}_{v}}\sum _{r\ne u}^{{N}_{v}}\langle {\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle \langle{\chi }_{\mu }|{ \varphi }_{r}\rangle{\frak D}(tk|ur)&\cr &\quad +\sum _{t = 1}^{{N}_{v}}\langle {\chi }_{\mu }|{ \hat I}_{{\bf h}} |{ \varphi }_{t}\rangle{\frak D}(k|t)\Biggr]. &(52)}]$

By means of a similar mathematical procedure, namely by expanding supercofactors in terms of the higher-order ones and by exploiting the supercofactor properties [in particular properties (31) and (32)], it is also possible to obtain the analytical expression for the first derivatives of $[\specialfonts{\frak D}]$ with respect to the coefficients of the SC orbitals:

$[\specialfonts{{\partial {\frak D}}\over{\partial {C}_{\mu k}}} = {\rm 2}\sum _{r = {\rm 1}}^{{N}_{v}}\langle{\chi }_{\mu }|{ \varphi }_{r}\rangle{\frak D}(k|r). \eqno(53)]$

Now, let us take into account the first derivatives with respect to the spin-coupling coefficients {c_S,k}. In this regard, it is worth observing that the derivative of any rth-order ordinary supercofactor gives a symmetrized indexed supercofactor of the same order:

$[\specialfonts{{\partial {\frak D}({t}_{\rm 1} {t}_{\rm 2}\ldots { t}_{r}| {u}_{\rm 1} {u}_{\rm 2}\ldots { u}_{r})}\over{\partial {c}_{S,k}}} = {{\frak D}}_{k}'' ({t}_{\rm 1} {t}_{\rm 2}\ldots { t}_{r}| {u}_{\rm 1} {u}_{\rm 2}\ldots { u}_{r}) .\eqno(54)]$

Therefore, exploiting the previous property and using the analytical expressions of $[{F}_{{\bf h}}^{\rm SC}]$ and $[\specialfonts{\frak D}]$ in equations (38) and (18), respectively, it is easy to see that

$[\specialfonts{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial {c}_{S,k}}} = \sum _{t,u = {\rm 1}}^{{N}_{v}}\langle {\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle {{\frak D}}_{k}'' (t| u) \eqno(55)]$

and

$[\specialfonts{{\partial {\frak D}}\over{\partial {c}_{S,k}}} = {{\frak D}}_{k}''. \eqno(56)]$

As mentioned above, to determine the coefficients of the SC orbitals and the spin-coupling coefficients that minimize functional (23), the second derivatives of χ² are also necessary. To obtain their analytical form, let us start from equation (33) and let us derive it with respect to the generic variable y. In compact form, we can write

$[ {{{\partial }^{2}{\chi }^{2}}\over{\partial y \partial x}} = \sum _{{\bf h}}{K}_{{\bf h}}{{\partial {A}_{{\bf h}}}\over{\partial y}}+{{\partial {K}_{{\bf h}}}\over{\partial y}} {A}_{{\bf h}}, \eqno(57)]$

where $[{A}_{{\bf h}}]$ is given by

$[\specialfonts\eqalignno{{A}_{{\bf h}} &= {{1}\over{2}} \left[{{\partial {F}_{{\bf h}}^{\rm calc}}\over{\partial x}} {({F}_{{\bf h}}^{\rm calc})}^{*}+{F}_{{\bf h}}^{\rm calc}{{\partial {({F}_{{\bf h}}^{\rm calc})}^{*}}\over{\partial x}}\right]&\cr & = {\rm Re}\{{F}_{{\bf h}}^{\rm calc}\}\left(-{{{\rm Re}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial x}}+ {{1}\over{{\frak D}}} {\rm Re}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial x}}\right\}\right)&\cr &\quad +{\rm Im}\{{F}_{{\bf h}}^{\rm calc}\}\left(-{{{\rm Im}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial x}}+ {{1}\over{{\frak D}}} {\rm Im}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial x}}\right\}\right) &\cr &&(58)}]$

and where $[\partial {K}_{{\bf h}}/\partial y]$ is

$[\specialfonts\eqalignno{{{\partial {K}_{{\bf h}}}\over{\partial y}} &= {\Gamma }_{{\bf h}}\Biggl[{{{\rm Re}\{{F}_{{\bf h}}^{\rm calc}\}}\over{{|{F}_{{\bf h}}^{\rm calc}|}^{3}}}\left(-{{{\rm Re}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial y}}+ {{1}\over{{\frak D}}} {\rm Re}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial y}}\right\}\right)&\cr &\quad +{{{\rm Im}\{{F}_{{\bf h}}^{\rm calc}\}}\over{{|{F}_{{\bf h}}^{\rm calc}|}^{3}}}\left(-{{{\rm Im}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial y}}+ {{1}\over{{\frak D}}} {\rm Im}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial y}}\right\}\right)\Biggr] &\cr &&(59)}]$

with

$[{\Gamma }_{{\bf h}} = {{2\eta |{F}_{{\bf h}}^{\rm exp}|}\over{({N}_{r}-{N}_{p}) {\sigma }_{{\bf h}}^{2}}}. \eqno(60)]$

Therefore, equation (58) and equation (59) depend only on $[{F}_{{\bf h}}^{\rm calc}]$ and on the first derivatives of $[{F}_{{\bf h}}^{\rm SC}]$ and $[\specialfonts{\frak D}]$ , which are already known. In contrast, if we consider the derivative $[\partial {A}_{{\bf h}}/\partial y]$ , we have

$[\eqalignno{{{\partial {A}_{{\bf h}}}\over{\partial y}}& = {{1}\over{2}} \Biggl[{{{\partial }^{2}{F}_{{\bf h}}^{\rm calc}}\over{\partial y \partial x}} {({F}_{{\bf h}}^{\rm calc})}^{*}+{{\partial {F}_{{\bf h}}^{\rm calc}}\over{\partial x}} {{\partial {({F}_{{\bf h}}^{\rm calc})}^{*}}\over{\partial y}}&\cr &\quad +{{\partial {F}_{{\bf h}}^{\rm calc}}\over{\partial y}} {{\partial {({F}_{{\bf h}}^{\rm calc})}^{*}}\over{\partial x}}+F_{{\bf h}}^{\rm calc} {{{\partial }^{2}{({F}_{{\bf h}}^{\rm calc})}^{*}}\over{\partial y\partial x}}\Biggr]. &\cr &&(61)}]$

Also, in this case, since we assume we are working with real orbitals, both for the second derivatives with respect to the coefficients of the SC orbitals and for the second derivatives with respect to the spin-coupling coefficients, it is possible to show that

$[{{{\partial }^{2}{({F}_{{\bf h}}^{\rm calc})}^{*}}\over{\partial y\partial x}} = {\left({{{\partial }^{2}{F}_{{\bf h}}^{\rm calc}}\over{\partial y \partial x}}\right)}^{*}. \eqno(62)]$

Therefore, taking into account relations (35) and (62), we obtain

$[\eqalignno{{{\partial {A}_{{\bf h}}}\over{\partial y}}& = { \rm Re}\{{F}_{{\bf h}}^{\rm calc}\} {\rm Re}\left\{{{{\partial }^{2}{F}_{{\bf h}}^{\rm calc}}\over{\partial y \partial x}}\right\}+{\rm Re}\left\{{{\partial {F}_{{\bf h}}^{\rm calc}}\over{\partial x}}\right\} {\rm Re}\left\{{{\partial {F}_{{\bf h}}^{\rm calc}}\over{\partial y}}\right\}&\cr &\quad +{\rm Im}\{{F}_{{\bf h}}^{\rm calc}\} {\rm Im}\left\{{{{\partial }^{2}{F}_{{\bf h}}^{\rm calc}}\over{\partial y \partial x}}\right\}+{\rm Im}\left\{{{\partial {F}_{{\bf h}}^{\rm calc}}\over{\partial x}}\right\} {\rm Im}\left\{{{\partial {F}_{{\bf h}}^{\rm calc}}\over{\partial y}}\right\} .&\cr &&(63)}]$

Now, starting from equation (39), we have

$[\specialfonts\eqalignno{{{{\partial }^{2}{F}_{{\bf h}}^{\rm calc}}\over{\partial y\partial x}} &= 2{{{F}_{{\bf h}}^{\rm SC}}\over{{{\frak D}}^{3}}} {{\partial {\frak D}}\over{\partial y}} {{\partial {\frak D}}\over{\partial x}} -{{1}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial x}} {{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial y}}-{{{F}_{{\bf h}}^{\rm SC}}\over{{{\frak D}}^{2}}} {{{\partial }^{2}{\frak D}}\over{\partial y\partial x}}&\cr &\quad -{{1}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial y}} {{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial x}}+{{1}\over{{\frak D}}} {{{\partial }^{2}{F}_{{\bf h}}^{\rm SC}}\over{\partial y\partial x}} &(64)}]$

and substituting (39) and (64) into (63), $[\partial {A}_{{\bf h}}/\partial y]$ becomes

$[\specialfonts\eqalignno{&{{\partial {A}_{{\bf h}}}\over{\partial y}} = {\rm Re}\{{F}_{{\bf h}}^{\rm calc}\}\Biggl({{2 {\rm Re}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{3}}} {{\partial {\frak D}}\over{\partial y}} {{\partial {\frak D}}\over{\partial x}} -{{1}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial x}} {\rm Re}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial y}}\right\}&\cr &\quad -{{{\rm Re}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{2}}} {{{\partial }^{2}{\frak D}}\over{\partial y\partial x}} -{{1}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial y}} {\rm Re}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial x}}\right\}+{{1}\over{{\frak D}}} {\rm Re}\left\{{{{\partial }^{2}{F}_{{\bf h}}^{\rm SC}}\over{\partial y\partial x}}\right\}\Biggr)&\cr &\quad + {\rm Im}\{{F}_{{\bf h}}^{\rm calc}\}\Biggl({{2 {\rm Im}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{3}}} {{\partial {\frak D}}\over{\partial y}} {{\partial {\frak D}}\over{\partial x}} -{{1}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial x}} {\rm Im}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial y}}\right\}&\cr &\quad -{{{\rm Im}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{2}}} {{{\partial }^{2}{\frak D}}\over{\partial y\partial x}}-{{1}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial y}} {\rm Im}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial x}}\right\}+{{1}\over{{\frak D}}} {\rm Im}\left\{{{{\partial }^{2}{F}_{{\bf h}}^{\rm SC}}\over{\partial y\partial x}}\right\}\Biggr)&\cr &\quad +\left(-{{{\rm Re}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial x}}+ {{1}\over{{\frak D}}} {\rm Re}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial x}}\right\}\right)&\cr &\quad\times\left(-{{{\rm Re}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial y}}+ {{1}\over{{\frak D}}} {\rm Re}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial y}}\right\}\right)&\cr &\quad +\left(-{{{\rm Im}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial x}}+ {{1}\over{{\frak D}}} {\rm Im}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial x}}\right\}\right)&\cr &\quad\times\left(-{{{\rm Im}\{{F}_{{\bf h}}^{\rm SC}\}}\over{{{\frak D}}^{2}}} {{\partial {\frak D}}\over{\partial y}}+ {{1}\over{{\frak D}}} {\rm Im}\left\{{{\partial {F}_{{\bf h}}^{\rm SC}}\over{\partial y}}\right\}\right). &(65)}]$

Therefore, it is clear that $[\partial {A}_{{\bf h}}/\partial y]$ does not depend only on $[{F}_{{\bf h}}^{\rm calc}]$ and on the first derivatives of $[{F}_{{\bf h}}^{\rm SC}]$ and $[\specialfonts{\frak D}]$ , but also on the second derivatives $[{\partial }^{2}{F}_{{\bf h}}^{\rm SC}/\partial y\partial x]$ and $[\specialfonts{\partial }^{2}{\frak D}/\partial y\partial x]$ , which will be discussed below.

Let us consider the orbital–orbital second derivatives of $[{F}_{{\bf h}}^{\rm SC}]$ and $[\specialfonts{\frak D}]$ , namely the second derivatives of $[{F}_{{\bf h}}^{\rm SC}]$ and $[\specialfonts{\frak D}]$ only with respect to the coefficients of the SC orbitals. These second derivatives can be obtained by following the mathematical procedure already adopted to derive the expression of $[\partial {F}_{{\bf h}}^{\rm SC}/\partial {C}_{\mu k}]$ , namely by exploiting the properties of the supercofactors and, in particular, by using recurrence relation (27). Nevertheless, since these derivations require a large number of steps, here, for the sake of simplicity, we will report only the final analytical expressions of the second derivatives. It is necessary to distinguish between two different cases: (i) the simplest one, in which the derivatives are with respect to two coefficients of the same SC orbital (case h = k) and (ii) the one in which the derivatives are with respect to two coefficients of different SC orbitals (case $[h\ne k]$ ).

Let us consider the case h = k. It is possible to show that

$[\specialfonts\eqalignno{&{{{\partial }^{2}{F}_{{\bf h}}^{\rm SC}}\over{\partial {C}_{\nu k}\partial {C}_{\mu k}}}&\cr & = 2\Biggl[\sum _{t\ne k}^{{N}_{v}}\sum _{u\ne k}^{{N}_{v}}\sum _{r\ne u,k}^{{N}_{v}}\sum _{s\ne t,k}^{{N}_{v}}\langle {\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle \langle{\chi }_{\mu }|{ \varphi }_{r}\rangle \langle{\chi }_{\nu }|{ \varphi }_{s}\rangle{\frak D}(tks|urk)&\cr &\quad +\sum _{t\ne k}^{{N}_{v}}\sum _{u\ne k}^{{N}_{v}}\langle {\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle \langle{\chi }_{\mu }|{\chi }_{\nu }\rangle {\frak D}(tk|uk) &\cr &\quad +\sum _{t\ne k}^{{N}_{v}}\sum _{r\ne k}^{{N}_{v}}\left(\langle {\chi }_{\nu }|{ \hat I}_{{\bf h}} |{ \varphi }_{t}\rangle \langle{\chi }_{\mu }|{ \varphi }_{r}\rangle +\langle {\chi }_{\mu }|{ \hat I}_{{\bf h}} |{ \varphi }_{t}\rangle \langle{\chi }_{\nu }|{ \varphi }_{r}\rangle\right){\frak D}(tk|kr)&\cr &\quad +\langle {\chi }_{\mu }|{ \hat I}_{{\bf h}} |{\chi }_{\nu }\rangle{\frak D}(k|k)\Biggr] &(66)}]$

and

$[\specialfonts\eqalignno{{{{\partial }^{2}{\frak D}}\over{\partial {C}_{\nu k}\partial {C}_{\mu k}}} &= {\rm 2}\Biggl[\sum _{r\ne k}^{{N}_{v}}\sum _{s\ne k}^{{N}_{v}}\langle{\chi }_{\mu }|{ \varphi }_{r}\rangle \langle{\chi }_{\nu }|{ \varphi }_{s}\rangle{\frak D}(ks|rk)&\cr &\quad+\langle{\chi }_{\mu }|{\chi }_{\nu }\rangle{\frak D}(k|k)\Biggr]. &(67)}]$

For the case $[h\ne k]$ the derivation is slightly more cumbersome and we obtain

$[\specialfonts\eqalignno{&{{{\partial }^{2}{F}_{{\bf h}}^{\rm SC}}\over{\partial {C}_{\nu h}\partial {C}_{\mu k}}}&\cr & = 2\Biggl[\sum _{t\ne k,h}^{{N}_{v}}\sum _{u = 1}^{{N}_{v}}\sum _{r\ne u}^{{N}_{v}}\sum _{s\ne u,r}^{{N}_{v}}\langle {\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle \langle{\chi }_{\mu }|{ \varphi }_{r}\rangle \langle{\chi }_{\nu }|{ \varphi }_{s}\rangle{\frak D}(tkh|urs)&\cr &\quad +\sum _{t\ne k}^{{N}_{v}}\sum _{u\ne h}^{{N}_{v}}\sum _{r\ne u,h}^{{N}_{v}}\sum _{s\ne t,k}^{{N}_{v}}\langle {\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle \langle{\chi }_{\mu }|{ \varphi }_{r}\rangle \langle{\chi }_{\nu }|{ \varphi }_{s}\rangle{\frak D}(tks|urh)&\cr &\quad +\sum _{t\ne k}^{{N}_{v}}\sum _{u\ne h}^{{N}_{v}}\langle{\varphi }_{t}|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle \langle{\chi }_{\mu }|{\chi }_{\nu }\rangle{\frak D}(tk|uh)&\cr &\quad +\sum _{t\ne k}^{{N}_{v}}\sum _{r\ne h}^{{N}_{v}}\langle {\chi }_{\nu }|{ \hat I}_{{\bf h}} |{ \varphi }_{t}\rangle \langle{\chi }_{\mu }|{ \varphi }_{r}\rangle{\frak D}(tk|hr) &\cr &\quad +\sum _{u = 1}^{{N}_{v}}\sum _{r\ne u}^{{N}_{v}}\langle {\chi }_{\nu }|{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle \langle{\chi }_{\mu }|{ \varphi }_{r}\rangle{\frak D}(hk|ur)&\cr &\quad +\sum _{t = 1}^{{N}_{v}}\sum _{s\ne t}^{{N}_{v}}\langle {\chi }_{\mu }|{ \hat I}_{{\bf h}} |{ \varphi }_{t}\rangle \langle{\chi }_{\nu }|{ \varphi }_{s}\rangle{\frak D}(kh|ts)&\cr &\quad +\sum _{t\ne h}^{{N}_{v}}\sum _{s\ne k}^{{N}_{v}}\langle {\chi }_{\mu }|{ \hat I}_{{\bf h}} |{ \varphi }_{t}\rangle \langle{\chi }_{\nu }|{ \varphi }_{s}\rangle{\frak D}(ks|th)+\langle {\chi }_{\mu }|{ \hat I}_{{\bf h}} |{\chi }_{\nu }\rangle{\frak D}(k|h)\Biggr] &\cr &&(68)}]$

and

$[\specialfonts\eqalignno{&{{\partial ^{2}{\frak D}}\over{\partial {C}_{\nu h}\partial {C}_{\mu k}}} =2\Biggl[\sum _{r=1}^{{N}_{v}}\sum _{s \ne r}^{{N}_{v}}\langle {\chi }_{\mu}|{ \varphi }_{r}\rangle \langle{\chi }_{\nu }|{ \varphi }_{s}\rangle {\frak D}(kh|rs)&\cr &\quad +\sum_{r\ne h}^{N_v}\sum_{s\ne k}^{N_v}\langle \chi_\mu | \varphi_r \rangle \langle \chi_\nu |\varphi_s\rangle {\frak D}(ks|rh)+\langle\chi_\mu |\chi_\nu\rangle {\frak D}(k|h)\Biggr].&\cr &&(69)}]$

Now, let us take into account the spin–orbital second derivatives of $[F_{\bf h}^{\rm SC}]$ and $[\specialfonts{\frak D}]$ . Unlike the orbital–orbital second derivatives, the derivation is really straightforward. In fact, exploiting property (54) and applying it to equations (52) and (53), we easily obtain

$[\specialfonts\eqalignno{{{\partial ^{2}F_{\bf h}^{\rm SC}}\over{\partial {c}_{S,h}\partial {C}_{\mu k}}}& =2\Biggl[\sum _{t\ne k}^{{N}_{v}}\sum _{u=1}^{{N}_{v}}\sum_{r\ne u}^{N_v}\langle \varphi_t |{ \hat I}_{{\bf h}} |{ \varphi }_{u}\rangle \langle{\chi }_{\mu }|{ \varphi }_{r}\rangle {\frak D}''_h(tk|ur)&\cr &\quad +\sum_{t=1}^{N_v}\langle \chi_\mu | {\hat I}_{\bf h}|\varphi_t \rangle {\frak D}''_{h}(k|t)\Biggr]&(70)}]$

and

$[\specialfonts{{\partial ^{2}{\frak D}}\over{\partial {c}_{S,h}\partial {C}_{\mu k}}} ={\rm 2}\sum_{r=1}^{N_v}\langle \chi_\mu |\varphi_r\rangle {\frak D}''_{h}(k|r).\eqno(71)]$

Finally, to obtain the analytical expressions of the spin–spin second derivatives of $[F_{\bf h}^{\rm SC}]$ and $[\specialfonts{\frak D}]$ , namely the analytical expressions of the second derivatives only with respect to spin-coupling coefficients, we have used a different approach. To this purpose, let us start considering the spin–spin second derivatives of $[\specialfonts{\frak D}]$ . It is easy to see that, using equations (6) and (17), $[\specialfonts{\frak D}]$ can also be expressed in this way:

$[\specialfonts{\frak D} =\langle \Psi_0^{\rm SC}|\Psi_0^{\rm SC}\rangle = \textstyle\sum\limits_{p=1}^{f_S^{N_v}} \sum\limits_{q=1}^{f_S^{N_v}} c_{S,p}c_{S,q}\langle \psi^N_{S,M\semi p}|\psi^N_{S,M\semi q}\rangle .\eqno(72)]$

Therefore, it is straightforward to show that

$[\specialfonts{{\partial ^{2}{\frak D}}\over{\partial {c}_{S,h}\partial {c}_{S,k}}}={\rm 2} \langle \psi^N_{S,M\semi h}|\psi^N_{S,M\semi k}\rangle={\rm 2}\textstyle\sum\limits_{i,j=1}^{N_d}d_{i,h}d_{j,k}\langle \Omega_i |\Omega_j\rangle, \eqno(73)]$

where we have exploited the fact that each generic SC structure $[\psi _{S,M\semi k}^N]$ can be written as a linear combination of Slater determinants $[\{ \Omega _i\}]$ . The coefficients of the linear combination { d_i,k } are identical to the coefficients in the expansion of the corresponding spin eigenfunction in terms of spin primitive functions [see equation (12)].

In order to compute the spin–spin second derivatives of $[F_{\bf h}^{\rm SC}]$ , first of all it is worth nothing that, after introducing the N-electron structure-factor operator

$[{\hat{\Im}}_{\bf h}=\textstyle\sum\limits_{a=1}^N{\hat I}_{\bf h}({\bf x}_a),\eqno(74)]$

$[F_{\bf h}^{\rm calc}]$ can be rewritten like this:

$[\specialfonts\eqalignno{F_{\bf h}^{\rm calc}&= \langle \Psi_0|{\hat{\Im}}_{\bf h}|\Psi_0\rangle =F_{\bf h}^{\rm core}+{{F_{\bf h}^{\rm SC}}\over {\frak D}}&\cr &= {1\over {\frak D}} \sum_{p=1}^{f_S^{N_v}}\sum_{q=1}^{f_S^{N_v}} c_{S,p}c_{S,q}\langle \psi^N_{S,M\semi p}|{\hat{\Im}}_{\bf h}| \psi^N_{S,M\semi q}\rangle .&(75)}]$

Therefore,

$[\specialfonts F_{\bf h}^{\rm SC}=\textstyle\sum\limits_{p={\rm 1}}^{f_S^{N_v}}\sum\limits_{q={\rm 1}}^{f_S^{N_v}} c_{S,p}c_{S,q}\langle \psi^N_{S,M\semi p}|{\hat{\Im}}_{\bf h}| \psi^N_{S,M\semi q}\rangle -F_{\bf h}^{\rm core}{\frak D}\eqno(76)]$

and it is easy to obtain

$[\specialfonts\eqalignno{{{\partial ^{2}F_{\bf h}^{\rm SC}}\over{\partial {c}_{S,h}\partial {c}_{S,k}}}&=2\langle \psi_{S,M\semi h}^N |{\hat{\Im}}_{\bf h}|\psi_{S,M\semi k}^N\rangle -F_{\bf h}^{\rm core} {{\partial^2{\frak D}}\over{\partial c_{S,h}\partial c_{S,k}}}&\cr &= 2\sum_{i,j=1}^{N_d} d_{i,h}d_{j,k}\langle \Omega_i |{\hat{\Im}}_{\bf h}|\Omega_j \rangle -F_{\bf h}^{\rm core}{{\partial^2 {\frak D}}\over{\partial c_{S,h}\partial c_{S,k}}}&(77)}]$

with $[\specialfonts{\partial ^2}{\frak D}/\partial {c_{S,h}}\partial {c_{S,k}}]$ given by equation (73).

Supporting information

Supporting figures and tables, and theoretical details of the Walker and Mezey similarity index. DOI: https://doi.org/10.1107/S2053273319011021/ae5070sup1.pdf

Funding information

AG and GM acknowledge the French Research Agency (ANR) for financial support in the form of the Young Researcher Project QuMacroRef (grant No. ANR-17-CE29-0005-01 to AG).

References

Aleksandrov, Y. V., Tsirelson, V. G., Reznik, I. M. & Ozerov, R. P. (1989). Phys. Status Solidi B, 155, 201–207. CrossRef CAS Web of Science Google Scholar
Bader, R. F. W. (1990). Atoms in Molecules: a Quantum Theory. Oxford University Press. Google Scholar
Becke, A. D. & Edgecombe, K. E. (1990). J. Chem. Phys. 92, 5397–5403. CrossRef CAS Web of Science Google Scholar
Boys, S. F. (1960). Rev. Mod. Phys. 32, 296–299. CrossRef CAS Web of Science Google Scholar
Bučinský, L., Jayatilaka, D. & Grabowsky, S. (2016). J. Phys. Chem. A, 120, 6650–6669. Web of Science PubMed Google Scholar
Bytheway, I., Grimwood, D. J., Figgis, B. N., Chandler, G. S. & Jayatilaka, D. (2002). Acta Cryst. A58, 244–251. Web of Science CrossRef CAS IUCr Journals Google Scholar
Bytheway, I., Grimwood, D. J. & Jayatilaka, D. (2002). Acta Cryst. A58, 232–243. Web of Science CrossRef CAS IUCr Journals Google Scholar
Casati, N., Genoni, A., Meyer, B., Krawczuk, A. & Macchi, P. (2017). Acta Cryst. B73, 584–597. Web of Science CSD CrossRef IUCr Journals Google Scholar
Chirgwin, B. H. & Coulson, C. A. (1950). Proc. R. Soc. London A, 201, 196–209. CAS Google Scholar
Clinton, W. L., Frishberg, C. A., Massa, L. J. & Oldfield, P. L. (1973). Int. J. Quantum Chem. 7, 505–514. CrossRef Google Scholar
Clinton, W. L. & Massa, L. J. (1972). Phys. Rev. Lett. 29, 1363–1366. CrossRef CAS Web of Science Google Scholar
Cooper, D. L., Gerratt, J. & Raimondi, M. (1986). Nature, 323, 699–701. CrossRef CAS Web of Science Google Scholar
Cooper, D. L., Gerratt, J. & Raimondi, M. (1991). Chem. Rev. 91, 929–964. CrossRef CAS Web of Science Google Scholar
Cooper, D. L., Gerratt, J., Raimondi, M., Sironi, M. & Thorsteinsson, T. (1993). Theor. Chim. Acta, 85, 261–270. CrossRef CAS Web of Science Google Scholar
Deutsch, M., Claiser, N., Pillet, S., Chumakov, Yu., Becker, P., Gillet, J.-M., Gillon, B., Lecomte, C. & Souhassou, M. (2012). Acta Cryst. A68, 675–686. Web of Science CSD CrossRef IUCr Journals Google Scholar
Deutsch, M., Gillon, B., Claiser, N., Gillet, J.-M., Lecomte, C. & Souhassou, M. (2014). IUCrJ, 1, 194–199. Web of Science CrossRef CAS PubMed IUCr Journals Google Scholar
Dewar, M. J. S. (1952). J. Am. Chem. Soc. 74, 3341–3345. CrossRef CAS Web of Science Google Scholar
Dos Santos, L. H. R., Genoni, A. & Macchi, P. (2014). Acta Cryst. A70, 532–551. Web of Science CSD CrossRef IUCr Journals Google Scholar
Edmiston, C. & Ruedenberg, K. (1963). Rev. Mod. Phys. 35, 457–464. CrossRef CAS Web of Science Google Scholar
Edmiston, C. & Ruedenberg, K. (1965). J. Chem. Phys. 43, S97–S116. CrossRef CAS Web of Science Google Scholar
Fornili, A., Sironi, M. & Raimondi, M. (2003). J. Mol. Struct. Theochem. 632, 157–172. Web of Science CrossRef CAS Google Scholar
Foster, J. M. & Boys, S. F. (1960). Rev. Mod. Phys. 32, 300–302. CrossRef CAS Web of Science Google Scholar
Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., Scalmani, G., Barone, V., Mennucci, B., Petersson, G. A., Nakatsuji, H., Caricato, M., Li, X., Hratchian, H. P., Izmaylov, A. F., Bloino, J., Zheng, G., Sonnenberg, J. L., Hada, M., Ehara, M., Toyota, K., Fukuda, R., Hasegawa, J., Ishida, M., Nakajima, T., Honda, Y., Kitao, O., Nakai, H., Vreven, T., Montgomery, J. A. Jr, Peralta, J. E., Ogliaro, F., Bearpark, M., Heyd, J. J., Brothers, E., Kudin, K. N., Staroverov, V. N., Kobayashi, R., Normand, J., Raghavachari, K., Rendell, A., Burant, J. C., Iyengar, S. S., Tomasi, J., Cossi, M., Rega, N., Millam, J. M., Klene, M., Knox, J. E., Cross, J. B., Bakken, V., Adamo, C., Jaramillo, J., Gomperts, R., Stratmann, R. E., Yazyev, O., Austin, A. J., Cammi, R., Pomelli, C., Ochterski, J. W., Martin, R. L., Morokuma, K., Zakrzewski, V. G., Voth, G. A., Salvador, P., Dannenberg, J. J., Dapprich, S., Daniels, A. D., Farkas, Ö., Foresman, J. B., Ortiz, J. V., Cioslowski, J. & Fox, D. J. (2009). Gaussian09, Revision D.01. Gaussian, Inc., Wallingford, Connecticut, USA. Google Scholar
Frishberg, C. & Massa, L. J. (1981). Phys. Rev. B, 24, 7018–7024. CrossRef CAS Web of Science Google Scholar
Fugel, M., Beckmann, J., Jayatilaka, D., Gibbs, G. V. & Grabowsky, S. (2018). Chem. Eur. J. 24, 6248–6261. Web of Science CrossRef CAS PubMed Google Scholar
Fugel, M., Kleemiss, F., Malaspina, L. A., Pal, R., Spackman, P. R., Jayatilaka, D. & Grabowsky, S. (2018). Aust. J. Chem. 71, 227–237. Web of Science CrossRef ICSD CAS Google Scholar
Gatti, C., Macetti, G. & Lo Presti, L. (2017). Acta Cryst. B73, 565–583. Web of Science CrossRef IUCr Journals Google Scholar
Gatti, C., Orlando, A. M. & Lo Presti, L. (2015). Chem. Sci. 6, 3845–3852. Web of Science CrossRef CAS PubMed Google Scholar
Genoni, A. (2013a). J. Phys. Chem. Lett. 4, 1093–1099. Web of Science CrossRef CAS PubMed Google Scholar
Genoni, A. (2013b). J. Chem. Theory Comput. 9, 3004–3019. Web of Science CrossRef CAS PubMed Google Scholar
Genoni, A. (2017). Acta Cryst. A73, 312–316. Web of Science CrossRef IUCr Journals Google Scholar
Genoni, A., Bučinský, L., Claiser, N., Contreras-García, J., Dittrich, B., Dominiak, P. M., Espinosa, E., Gatti, C., Giannozzi, P., Gillet, J.-M., Jayatilaka, D., Macchi, P., Madsen, A. Ø., Massa, L. J., Matta, C. F., Merz, K. M. Jr, Nakashima, P. N. H., Ott, H., Ryde, U., Schwarz, K., Sierka, M. & Grabowsky, S. (2018). Chem. Eur. J. 24, 10881–10905. Web of Science CrossRef CAS PubMed Google Scholar
Genoni, A., Dos Santos, L. H. R., Meyer, B. & Macchi, P. (2017). IUCrJ, 4, 136–146. Web of Science CrossRef CAS PubMed IUCr Journals Google Scholar
Genoni, A., Fornili, A. & Sironi, M. (2005). J. Comput. Chem. 26, 827–835. Web of Science CrossRef PubMed CAS Google Scholar
Genoni, A., Franchini, D., Pieraccini, S. & Sironi, M. (2018). Chem. Eur. J. 24, 15507–15511. Web of Science CrossRef CAS PubMed Google Scholar
Genoni, A., Ghitti, M., Pieraccini, S. & Sironi, M. (2005). Chem. Phys. Lett. 415, 256–260. Web of Science CrossRef CAS Google Scholar
Genoni, A. & Meyer, B. (2016). Adv. Quantum Chem. 73, 333–362. Web of Science CrossRef CAS Google Scholar
Genoni, A. & Sironi, M. (2004). Theor. Chem. Acc. 112, 254–262. Web of Science CrossRef CAS Google Scholar
Gerratt, J. (1971). Adv. Atom Mol. Phys. 7, 141–221. CrossRef Google Scholar
Gerratt, J. & Lipscomb, W. N. (1968). Proc. Natl Acad. Sci. USA, 59, 332–335. CrossRef PubMed CAS Web of Science Google Scholar
Gilbert, T. L. (1975). Phys. Rev. B, 12, 2111–2120. CrossRef Web of Science Google Scholar
Gillet, J.-M. (2007). Acta Cryst. A63, 234–238. Web of Science CrossRef CAS IUCr Journals Google Scholar
Gillet, J.-M. & Becker, P. J. (2004). J. Phys. Chem. Solids, 65, 2017–2023. Web of Science CrossRef CAS Google Scholar
Gillet, J.-M., Becker, P. J. & Cortona, P. (2001). Phys. Rev. B, 63, 235115. Web of Science CrossRef Google Scholar
Glendening, E. D. & Weinhold, F. (1998). J. Comput. Chem. 19, 593–609. CrossRef CAS Google Scholar
Goddard, W. A. III (1967). Phys. Rev. 157, 81–93. CrossRef CAS Web of Science Google Scholar
Goddard, W. A. III, Dunning, T. H. Jr, Hunt, W. J. & Hay, P. J. (1973). Acc. Chem. Res. 6, 368–376. CrossRef CAS Web of Science Google Scholar
Goldberg, M. J. & Massa, L. J. (1983). Int. J. Quantum Chem. 24, 113–126. CrossRef CAS Web of Science Google Scholar
Goldfeld, S. M., Quandt, R. & Trotter, F. (1966). Econometrica, 34, 541–551. CrossRef Web of Science Google Scholar
Grabowsky, S. (2017). Acta Cryst. A73, C568. Web of Science CrossRef IUCr Journals Google Scholar
Grabowsky, S., Genoni, A. & Bürgi, H.-B. (2017). Chem. Sci. 8, 4159–4176. Web of Science CrossRef CAS PubMed Google Scholar
Grabowsky, S., Jayatilaka, D., Mebs, S. & Luger, P. (2010). Chem. Eur. J. 16, 12818–12821. Web of Science CrossRef CAS PubMed Google Scholar
Grabowsky, S., Luger, P., Buschmann, J., Schneider, T., Schirmeister, T., Sobolev, A. N. & Jayatilaka, D. (2012). Angew. Chem. Int. Ed. 51, 6776–6779. Web of Science CSD CrossRef ICSD CAS Google Scholar
Grabowsky, S., Weber, M., Jayatilaka, D., Chen, Y.-S., Grabowski, M. T., Brehme, R., Hesse, M., Schirmeister, T. & Luger, P. (2011). J. Phys. Chem. A, 115, 12715–12732. Web of Science CSD CrossRef CAS PubMed Google Scholar
Grimwood, D. J., Bytheway, I. & Jayatilaka, D. (2003). J. Comput. Chem. 24, 470–483. Web of Science CrossRef PubMed CAS Google Scholar
Grimwood, D. J. & Jayatilaka, D. (2001). Acta Cryst. A57, 87–100. Web of Science CrossRef CAS IUCr Journals Google Scholar
Gueddida, S., Yan, Z., Kibalin, I., Voufack, A. B., Claiser, N., Souhassou, M., Lecomte, C., Gillon, B. & Gillet, J.-M. (2018). J. Chem. Phys. 148, 164106. Web of Science CrossRef PubMed Google Scholar
Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909–921. CrossRef CAS IUCr Journals Web of Science Google Scholar
Heitler, W. & London, F. (1927). Z. Phys. 44, 455–472. CrossRef CAS Google Scholar
Henderson, G. A. & Zimmerman, R. K. Jr (1976). J. Chem. Phys. 65, 619–622. CrossRef CAS Web of Science Google Scholar
Hibbs, D. E., Howard, S. T., Huke, J. P. & Waller, M. P. (2005). Phys. Chem. Chem. Phys. 7, 1772–1778. Web of Science CrossRef PubMed CAS Google Scholar
Hiberty, P. C., Humbel, S., Byrman, C. P. & van Lenthe, J. H. (1994). J. Chem. Phys. 101, 5969–5976. CrossRef CAS Web of Science Google Scholar
Hiberty, P. C. & Shaik, S. (2007). J. Comput. Chem. 28, 137–151. Web of Science CrossRef PubMed CAS Google Scholar
Hirao, K., Nakano, H., Nakayama, K. & Dupuis, M. (1996). J. Chem. Phys. 105, 9227–9239. CrossRef CAS Web of Science Google Scholar
Hollauer, E. & Nascimento, M. A. C. (1993). J. Chem. Phys. 99, 1207–1214. CrossRef CAS Web of Science Google Scholar
Howard, S. T., Huke, J. P., Mallinson, P. R. & Frampton, C. S. (1994). Phys. Rev. B, 49, 7124–7136. CrossRef CAS Web of Science Google Scholar
Hückel, E. (1930). Z. Phys. 60, 423–456. Google Scholar
Hückel, E. (1931). Z. Phys. 72, 310–337. Google Scholar
Hudák, M., Jayatilaka, D., Perašínová, L., Biskupič, S., Kožíšek, J. & Bučinský, L. (2010). Acta Cryst. A66, 78–92. Web of Science CrossRef IUCr Journals Google Scholar
Jayatilaka, D. (1998). Phys. Rev. Lett. 80, 798–801. Web of Science CrossRef CAS Google Scholar
Jayatilaka, D. (2012). Modern Charge-Density Analysis, edited by C. Gatti & P. Macchi, pp. 213–257. Dordrecht: Springer Netherlands. Google Scholar
Jayatilaka, D. & Grimwood, D. (2004). Acta Cryst. A60, 111–119. Web of Science CrossRef CAS IUCr Journals Google Scholar
Jayatilaka, D. & Grimwood, D. J. (2001). Acta Cryst. A57, 76–86. Web of Science CrossRef CAS IUCr Journals Google Scholar
Jayatilaka, D. & Grimwood, D. J. (2003). Computational Science – ICCS 2003, edited by P. M. A. Sloot, D. Abramson, A. V. Bogdanov, J. J. Dongarra, A. Y. Zomaya & Y. E. Gorbachev, pp. 142–151. Berlin, Heidelberg: Springer-Verlag. Google Scholar
Karadakov, P. B., Gerratt, J., Cooper, D. L. & Raimondi, M. (1992). J. Chem. Phys. 97, 7637–7655. CrossRef CAS Web of Science Google Scholar
Kohout, M. (2004). Int. J. Quantum Chem. 97, 651–658. Web of Science CrossRef CAS Google Scholar
Ladner, R. C. & Goddard, W. A. III (1969). J. Chem. Phys. 51, 1073–1087. CrossRef CAS Web of Science Google Scholar
London, F. (1928). Z. Phys. 46, 455–477. CrossRef CAS Google Scholar
Löwdin, P.-O. (1955). Phys. Rev. 97, 1474–1489. Google Scholar
Macetti, G., Lo Presti, L. & Gatti, C. (2018). J. Comput. Chem. 39, 587–603. Web of Science CrossRef CAS PubMed Google Scholar
Massa, L., Goldberg, M., Frishberg, C., Boehme, R. F. & La Placa, S. J. (1985). Phys. Rev. Lett. 55, 622–625. CrossRef PubMed CAS Web of Science Google Scholar
Massa, L. & Matta, C. F. (2018). J. Comput. Chem. 39, 1021–1028. Web of Science CrossRef CAS PubMed Google Scholar
McDouall, J. J. W. (1992). Theor. Chim. Acta, 83, 339–350. CrossRef CAS Web of Science Google Scholar
McWeeny, R. (1992). Methods of Molecular Quantum Mechanics, 2nd ed. London: Academic Press. Google Scholar
Meyer, B. & Genoni, A. (2018). J. Phys. Chem. A, 122, 8965–8981. Web of Science CrossRef CAS PubMed Google Scholar
Meyer, B., Guillot, B., Ruiz-Lopez, M. F. & Genoni, A. (2016a). J. Chem. Theory Comput. 12, 1052–1067. Web of Science CrossRef CAS PubMed Google Scholar
Meyer, B., Guillot, B., Ruiz-Lopez, M. F., Jelsch, C. & Genoni, A. (2016b). J. Chem. Theory Comput. 12, 1068–1081. Web of Science CrossRef CAS PubMed Google Scholar
Munshi, P. & Guru Row, T. N. (2006). Acta Cryst. B62, 612–626. Web of Science CSD CrossRef CAS IUCr Journals Google Scholar
Novara, R. F., Genoni, A. & Grabowsky, S. (2018). ChemViews, doi:10.1002/chemv.201800066. Google Scholar
Pauling, L. (1939). The Nature of the Chemical Bond. An Introduction to Modern Structural Chemistry. Ithaca: Cornell University Press. Google Scholar
Pipek, J. & Mezey, P. G. (1989). J. Chem. Phys. 90, 4916–4926. CrossRef CAS Web of Science Google Scholar
Raimondi, M., Sironi, M., Gerratt, J. & Cooper, D. L. (1996). Int. J. Quantum Chem. 60, 225–233. CrossRef CAS Web of Science Google Scholar
Roothaan, C. C. J. (1951). Rev. Mod. Phys. 23, 69–89. CrossRef CAS Web of Science Google Scholar
Roversi, P., Irwin, J. J. & Bricogne, G. (1998). Acta Cryst. A54, 971–996. Web of Science CrossRef CAS IUCr Journals Google Scholar
Sakata, M. & Sato, M. (1990). Acta Cryst. A46, 263–270. CrossRef CAS Web of Science IUCr Journals Google Scholar
Schmider, H., Smith, V. H. Jr & Weyrich, W. (1990). Trans. Am. Crystallogr. Assoc. 26, 125–140. CAS Google Scholar
Schmider, H., Smith, V. H. Jr & Weyrich, W. (1992). J. Chem. Phys. 96, 8986–8994. CrossRef CAS Web of Science Google Scholar
Schmider, H. L. & Becke, A. D. (2000). J. Mol. Struct. Theochem, 527, 51–61. Web of Science CrossRef CAS Google Scholar
Shurki, A. (2006). Theor. Chem. Acc. 116, 253–261. Web of Science CrossRef CAS Google Scholar
Simonetta, M., Gianinetti, E. & Vandoni, I. (1968). J. Chem. Phys. 48, 1579–1594. CrossRef CAS Web of Science Google Scholar
Sironi, M., Genoni, A., Civera, M., Pieraccini, S. & Ghitti, M. (2007). Theor. Chem. Acc. 117, 685–698. Web of Science CrossRef CAS Google Scholar
Sironi, M., Ghitti, M., Genoni, A., Saladino, G. & Pieraccini, S. (2009). J. Mol. Struct. Theochem. 898, 8–16. Web of Science CrossRef CAS Google Scholar
Smaalen, S. van & Netzel, J. (2009). Phys. Scr. 79, 048304. Web of Science CrossRef Google Scholar
Snyder, J. A. & Stevens, E. D. (1999). Chem. Phys. Lett. 313, 293–298. Web of Science CrossRef CAS Google Scholar
Song, L., Wu, W., Hiberty, P. C., Danovich, D. & Shaik, S. (2003). Chem. Eur. J. 9, 4540–4547. Web of Science CrossRef PubMed CAS Google Scholar
Stewart, R. F. (1969). J. Chem. Phys. 51, 4569–4577. CrossRef CAS Web of Science Google Scholar
Stewart, R. F. (1976). Acta Cryst. A32, 565–574. CrossRef CAS IUCr Journals Web of Science Google Scholar
Stoll, H., Wagenblast, G. & Preuss, H. (1980). Theor. Chim. Acta, 57, 169–178. CrossRef CAS Web of Science Google Scholar
Tanaka, K. (1988). Acta Cryst. A44, 1002–1008. CrossRef CAS Web of Science IUCr Journals Google Scholar
Tanaka, K. (2018). Acta Cryst. A74, 345–356. Web of Science CrossRef IUCr Journals Google Scholar
Thorsteinsson, T., Cooper, D. L., Gerratt, J., Karadakov, P. B. & Raimondi, M. (1996). Theor. Chim. Acta, 93, 343–366. CrossRef CAS Google Scholar
Tsirelson, V. (2018). J. Comput. Chem. 39, 1029–1037. Web of Science CrossRef CAS PubMed Google Scholar
Van Lenthe, J. H. & Balint-Kurti, G. G. (1983). J. Chem. Phys. 78, 5699–5713. CrossRef CAS Web of Science Google Scholar
Voter, A. F. & Goddard, W. A. III (1981). J. Chem. Phys. 75, 3638–3639. CrossRef CAS Web of Science Google Scholar
Voufack, A. B., Claiser, N., Lecomte, C., Pillet, S., Pontillon, Y., Gillon, B., Yan, Z., Gillet, J.-M., Marazzi, M., Genoni, A. & Souhassou, M. (2017). Acta Cryst. B73, 544–549. Web of Science CSD CrossRef IUCr Journals Google Scholar
Walker, P. D. & Mezey, P. G. (1994). J. Am. Chem. Soc. 116, 12022–12032. CrossRef CAS Web of Science Google Scholar
Waller, M. P., Howard, S. T., Platts, J. A., Piltz, R. O., Willock, D. J. & Hibbs, D. E. (2006). Chem. Eur. J. 12, 7603–7614. Web of Science CrossRef PubMed CAS Google Scholar
Weinhold, F. & Landis, C. R. (2001). Chem. Educ. Res. Pract. 2, 91–104. CrossRef CAS Google Scholar
Weyrich, W. (2006). Lect. Ser. Comput. Comput. Sci. 5, 1–3. Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.

FOUNDATIONS
ADVANCES

ISSN: 2053-2733

Volume 75| Part 6| November 2019| Pages 778-797

https://doi.org/10.1107/S2053273319011021