research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Journal logoFOUNDATIONS
ADVANCES
ISSN: 2053-2733

Inferring the one-electron reduced density matrix of molecular crystals from experimental data sets through semidefinite programming

CROSSMARK_Color_square_no_text.svg

aCentraleSupélec School, Paris-Saclay University, Gif-sur-Yvette, 91190, France, and bStructures, Properties and Modeling of Solids Laboratory, CentraleSupélec School, CNRS-UMR8580, Paris-Saclay University, Gif-sur-Yvette, 91190, France
*Correspondence e-mail: benjamin.debruyne@centraliens.net

Edited by P. M. Dominiak, University of Warsaw, Poland (Received 18 October 2019; accepted 25 November 2019)

Constructing a quantum description of crystals from scattering experiments is of great significance to explain their macroscopic properties and to evaluate the pertinence of theoretical ab initio models. While reconstruction methods of the one-electron reduced density matrix have already been proposed, they are usually tied to strong assumptions that limit and may introduce bias in the model. The goal of this paper is to infer a one-electron reduced density matrix (1-RDM) with minimal assumptions. It has been found that the mathematical framework of semidefinite programming can achieve this goal. Additionally, it conveniently addresses the nontrivial constraints on the 1-RDM which were major hindrances for the existing models. The framework established in this work can be used as a reference to interpret experimental results. This method has been applied to the crystal of dry ice and provides very satisfactory results when compared with periodic ab initio calculations.

1. Introduction

The computation of one-electron expectation values such as the mean position, the mean momentum or the mean kinetic energy of electrons in a crystal does not require more than the mere knowledge of the one-electron reduced density matrix (1-RDM) (Löwdin, 1955[Löwdin, P.-O. (1955). Phys. Rev. 97, 1474-1489.]; Coleman, 1963[Coleman, A. J. (1963). Rev. Mod. Phys. 35, 668-686.]; McWeeny, 1960[McWeeny, R. (1960). Rev. Mod. Phys. 32, 335-369.]; Davidson, 1976[Davidson, E. (1976). Quantum Chemistry. New York: Academic Press.]). This quantity provides a quantum description of an average electron and has been proved to be sufficient (Lathiotakis & Marques, 2008[Lathiotakis, N. N. & Marques, M. A. L. (2008). J. Chem. Phys. 128, 184103.]; Gilbert, 1975[Gilbert, T. L. (1975). Phys. Rev. B, 12, 2111-2120.]). Furthermore, the electron density in position and momentum spaces can easily be derived from such a quantity. It is therefore a useful tool for describing electronic properties at a quantum level. Additionally, use of the 1-RDM is well suited to represent mixed states systems using statistical ensembles of pure states. This is generally the case for crystals at nonzero temperature.

Several models have been proposed to approximate and refine a 1-RDM from experimental expectation values (Deutsch et al., 2012[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.], 2014[Deutsch, M., Gillon, B., Claiser, N., Gillet, J.-M., Lecomte, C. & Souhassou, M. (2014). IUCrJ, 1, 194-199.]; Hansen & Coppens, 1978[Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909-921.]; Gillet et al., 1999[Gillet, J.-M., Fluteaux, C. & Becker, P. J. (1999). Phys. Rev. B, 60, 2345-2349.], 2001[Gillet, J.-M., Becker, P. J. & Cortona, P. (2001). Phys. Rev. B, 63, 235115.]; Gillet, 2007[Gillet, J.-M. (2007). Acta Cryst. A63, 234-238.]; Gillet & Becker, 2004[Gillet, J.-M. & Becker, P. J. (2004). J. Phys. Chem. Solids, 65, 2017-2023.]; Gueddida et al., 2018[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.]; Pillet et al., 2001[Pillet, S., Souhassou, M., Pontillon, Y., Caneschi, A., Gatteschi, D. & Lecomte, C. (2001). New J. Chem. 25, 131-143.]; Schmider et al., 1992[Schmider, H., Smith, V. H. & Weyrich, W. (1992). J. Chem. Phys. 96, 8986-8994.]; Clinton & Massa, 1972[Clinton, W. L. & Massa, L. J. (1972). Phys. Rev. Lett. 29, 1363-1366.]; Tsirelson & Ozerov, 1996[Tsirelson, V. G. & Ozerov, R. P. (1996). Electron Density and Bonding in Crystals: Principles, Theory and X-ray Diffraction Experiments in Solid State Physics and Chemistry. New York: Taylor and Francis.]). The complementarity between position and momentum space expectation values in the description of the 1-RDM is now well accepted (Cooper et al., 2004[Cooper, M. J., Mijnarends, P. E., Mijnarends, P., Shiotani, N., Sakai, N. & Bansil, A. (2004). X-ray Compton Scattering. Oxford Series on Synchrotron Radiation, Vol. 5. Oxford University Press.]; Pisani, 2012[Pisani, C. (2012). Quantum-mechanical ab-initio Calculation of the Properties of Crystalline Materials. Lecture Notes in Chemistry, Vol. 67. Berlin: Springer Science & Business Media.]). For this reason, deep inelastic X-ray scattering data known as `directional Compton scattering profiles' (DCPs) have been taken into account in addition to X-ray or polarized neutron diffraction structure factors (SFs) to refine a variety of models. The former are related to 2D projections of electron density in momentum space, while the latter are linked to the Fourier coefficients of the electron density in position space. However, almost all of these models require an initial guess or assumption on the electronic configuration. When these are inappropriate or too simple, there is a risk that the model, hence the results, will be affected by a severe bias. Notable exceptions of model-free approaches are the maximum entropy (MEM) and charge flipping methods (Oszlányi & Sütő, 2004[Oszlányi, G. & Sütő, A. (2004). Acta Cryst. A60, 134-141.]; Dobrzynski & Zukowski, 1999[Dobrzynski, L. & Zukowski, E. (1999). J. Phys. Condens. Matter, 11, 8049-8060.]; Binkley et al., 1980[Binkley, J. S., Pople, J. A. & Hehre, W. J. (1980). J. Am. Chem. Soc. 102, 939-947.]). The purpose of this work is to investigate and assess a new method to obtain a 1-RDM from expectation values with minimal bias.

In order to serve as a reference, an initial periodic ab initio calculation [at the density functional theory (DFT) level] has been conducted from which the reference 1-RDM was extracted. From the same calculation, a limited number of structure factors and directional Compton profiles were generated. Once a random noise was added, these deteriorated data constituted our pseudo-experimental data.

The method explicitly takes into account the so-called N-representability conditions (Coleman, 1963[Coleman, A. J. (1963). Rev. Mod. Phys. 35, 668-686.]), which ensure that the inferred 1-RDM is quantum mechanically acceptable, i.e. that there exists a many-electron wavefunction from which the 1-RDM can be derived. Addressing these nontrivial conditions is made possible by the use of semidefinite programming (Vandenberghe & Boyd, 1996[Vandenberghe, L. & Boyd, S. (1996). SIAM Rev. 38, 49-95. ]), a recent subfield of convex optimization (Boyd & Vandenberghe, 2004[Boyd, S. & Vandenberghe, L. (2004). Convex Optimization. New York: Cambridge University Press.]) which is of growing interest in systems and control theory, geometry and statistics (Wolkowicz et al., 2012[Wolkowicz, H., Saigal, R. & Vandenberghe, L. (2012). Handbook of Semidefinite Programming: Theory, Algorithms, and Applications. International Series in Operations Research and Management Science, Vol. 27. New York: Springer Science & Business Media.]).

2. Method

2.1. Molecular spin-traced 1-RDM

In the following section, for simplicity, we will restrict our treatment to a crystal with a single molecule per cell that has N paired electrons. The method can be generalized to several molecules by either assigning a 1-RDM to each molecule provided that they can be considered electronically isolated from each other (as in Section 3[link]), or defining one 1-RDM for a group of interacting molecules. Additionally, spin orbitals can be employed to construct two spin-resolved 1-RDMs when the system bears unpaired electrons.

Let [\{\chi_i\}_{i\in\{1,\ldots,n\}}] be a set of atomic orbitals describing the electrons of each atom taken as an independent system. From [\{\chi_i\}_{i\in\{1,\ldots,n\}}], one can deduce an orthogonal basis set [\{\phi_i\}_{i\in\{1,\ldots,n\}}] for the molecule, using the Löwdin orthogonalization procedure (Löwdin, 1950[Löwdin, P. (1950). J. Chem. Phys. 18, 365-375.]) for example. Expanding the spin-traced 1-RDM [\widehat{\Gamma} ({\bf r},{\bf r}')] in such a basis, one reveals its basis-set representation: the population matrix [{\widehat{P}}], so that

[\widehat{\Gamma} ({\bf r},{\bf r}') = \textstyle\sum\limits_{i,j}^{n} \widehat{P}_{ij} \phi^*_i({\bf r}) \phi_j({\bf r'}).\eqno(1)]

Although it is not necessary to use an orthogonal basis, it is done here because the N-representability conditions are conveniently expressed in such a basis. In general, these conditions are expressed on the eigenvalues of the spin-traced 1-RDM. In this case, they are translated into conditions on the eigenvalues of [{\widehat{P}}] and state that they must lie in [0, 2] (as N is even) and their sum must be equal to N.

2.2. Expectation values

Any one-electron expectation value [\langle \widehat{O} \rangle] can be calculated from its operator [{ \widehat{\cal O}}_{{\bf r}'}] applied to the 1-RDM [\widehat{\Gamma}({\bf r},{\bf r'})]:

[\langle \widehat{O} \rangle = \textstyle\int \left[{ \widehat{\cal O}}_{{\bf r}'} \widehat{\Gamma}({\bf r},{\bf r'})\right]_{{\bf r'} = {\bf r}} \,{\rm d}{\bf r}\eqno(2)]

where [{ \widehat{\cal O}}_{{\bf r}'}] means that the operator only acts on variable [{\bf r}']. By defining the basis-set representation of [{ \widehat{\cal O}}] as

[{ \widehat{O}}_{ij} = \textstyle\int\phi^*_i({\bf r}) \left[{ \widehat{\cal O}}_{{\bf r}'} \phi_j({\bf r}')\right]_{{\bf r}' = {\bf r}} \,{\rm d}{\bf r}\eqno(3)]

one can conveniently write the expectation value as [\langle \widehat{O} \rangle] = [{\rm tr}\,(\widehat{P}\widehat{O})] [using equation (1[link])], where tr is the matrix trace operator.

In particular, in position space, the X-ray structure factors [F({\bf q})], which are given by

[F({\bf q}) = \textstyle \int \widehat{\Gamma}({\bf r},{\bf r}) \exp(i{\bf q}\cdot{\bf r}) \,{\rm d}{\bf r} \eqno(4)]

[= \textstyle\sum\limits_{i,j}^n {\widehat{P}}_{ij} \int \phi_i^*({\bf r}) \phi_j({\bf r}) \exp(i{\bf q}\cdot{\bf r}) \,{\rm d}{\bf r},\eqno(5)]

have an operator whose basis-set representation is

[{\widehat{F}}_{ij}({\bf q}) = \textstyle \int \phi^*_i({\bf r}) \phi_j({\bf r}) \exp(i{\bf q}\cdot{\bf r})\,{\rm d}{\bf r}.\eqno(6)]

In momentum space, the directional Compton profiles [J^{\bf u}(q)] can be defined through the autocorrelation function [B({\bf r})] (Pattison & Weyrich, 1979[Pattison, P. & Weyrich, W. (1979). J. Phys. Chem. Solids, 40, 213-222.]; Weyrich et al., 1979[Weyrich, W., Pattison, P. & Williams, B. (1979). Chem. Phys. 41, 271-284.]; Benesch et al., 1971[Benesch, R., Singh, S. & Smith, V. Jr (1971). Chem. Phys. Lett. 10, 151-153.]) as

[J^{\bf u}(q) = \int {{1}\over{2\pi}} B(t{\bf u}) \exp(-i t q) \,{\rm d}t \eqno(7)]

[B({\bf r}) = \textstyle \int \widehat{\Gamma}({\bf r}',{\bf r}+{\bf r}') \,{\rm d}{\bf r}'. \eqno(8)]

Their operator basis-set representation is therefore

[{\widehat{J}}_{ij}^{\ \bf u}(q) = \int \int {{1}\over{2\pi}} \phi^*_i({\bf r}') \phi_j(t{\bf u}+{\bf r}') \exp(-i t q)\, {\rm d}t\, {\rm d}{\bf r}' .\eqno(9)]

From equations (4[link]) and (7[link])–(8[link]), one can appreciate the complementarity of both expectation values as they, respectively, shed light upon the diagonal and the off-diagonal directions of the 1-RDM

2.3. Constrained least-squares fitting scheme

In the Bayesian sense, the objective is to infer the most probable population matrix [{\widehat{P}}] so that it fits given independent expectation values [\langle \widehat{O}_\alpha \rangle]. In the following, the expectation values [\langle \widehat{O}_\alpha \rangle] are SFs and DCPs data. Supposing the latter follow Gaussian error distributions with standard deviations [\sigma_\alpha] and no a priori knowledge is given on [{\widehat{P}}], the problem is equivalent to minimizing the so-called [\chi^2] function with respect to the elements of [{\widehat{P}}] (Gillet & Becker, 2004[Gillet, J.-M. & Becker, P. J. (2004). J. Phys. Chem. Solids, 65, 2017-2023.]; Sivia & Skilling, 2006[Sivia, D. & Skilling, J. (2006). Data Analysis: A Bayesian Tutorial. Oxford University Press.]). It can be summarized in the following optimization program:

[\matrix{{\rm minimize}\atop{\widehat{P}}& \sum_\alpha {{1}\over{\sigma_\alpha^2}}\left| \langle \widehat{O}_\alpha \rangle - {\rm tr}(\widehat{P}\widehat{O}_\alpha)\right|^2 \cr {\rm subject\ to} & {\rm tr}(\widehat{P}) = N,\hfill\cr \cr &\widehat{P}\succeq 0,\hfill\cr \cr &2I-\widehat{P}\succeq 0\hfill}\eqno(10)]

where I is the identity matrix and the notation [A \succeq 0] means that A is a symmetric positive semidefinite matrix, i.e. its eigenvalues are non-negative. The last two constraints are mathematically equivalent to the condition that the eigen­values of [{\widehat{P}}] must lie in [0, 2].

The following passage will cast program (10[link]) as a semidefinite optimization program. These steps are quite standard in the field of convex optimization (Boyd & Vandenberghe, 2004[Boyd, S. & Vandenberghe, L. (2004). Convex Optimization. New York: Cambridge University Press.]). Introducing a new variable t, program (10[link]) is equivalent to

[\matrix{{\rm minimize}\atop{\widehat{P},t}& t \hfill\cr {\rm subject\ to} & {\rm tr}(\widehat{P}) = N,\hfill\cr\cr &\widehat{P}\succeq 0,\hfill\cr \cr &2I-\widehat{P}\succeq 0,\hfill\cr\cr &t-||\boldDelta_{\sigma}{\bf O}||^2 \ge 0\hfill}\eqno(11)]

where [\boldDelta_{\sigma}{\bf O}] is a column vector whose elements are [[\langle \widehat{O}_\alpha \rangle] [- {\rm tr}\,(\widehat{P}\widehat{O}_\alpha)]]/[\sigma_\alpha] and [||\cdot||] is the Euclidean norm.

Using Schur's complement (Zhang, 2006[Zhang, F. (2006). The Schur Complement and its Applications. Numerical Algorithms, Vol. 4. New York: Springer Science & Business Media.]), the last constraint of program (11[link]) can be written as a linear matrix inequality:

[\left[\matrix{ I & \boldDelta_\sigma {\bf O} \cr (\boldDelta_\sigma {\bf O})^{\rm T} & t }\right] \succeq 0 \eqno(12)]

where I is the identity matrix of appropriate dimensions. This inequality is indeed linear with respect to [\widehat{P}] as [\boldDelta_\sigma {\bf O}] is a linear function of [\widehat{P}].

This type of program, where the objective function is linear and the constraints are linear combinations of symmetric matrices that must be positive semidefinite, has been extensively studied and is referred to as the class of semidefinite programming (Vandenberghe & Boyd, 1996[Vandenberghe, L. & Boyd, S. (1996). SIAM Rev. 38, 49-95. ]). Interior-point algorithms can be used to solve this class of problems and no initial guess is required. Treatment of the two-electron reduced density matrix (2-RDM) by semidefinite programming has already been reported in the context of variational computation of molecules (Mazziotti, 2007[Mazziotti, D. A. (2007). ESAIM: M2AN, 41, 249-259.]).

In the present work, this program has been addressed by using the optimization software Mosek (ApS, 2017[ApS, M. (2017). The MOSEK optimization toolbox for MATLAB manual. Version 8.1. https://docs.mosek.com/8.1/toolbox/index.html.]) interfaced by the Yalmip toolbox (Löfberg, 2004[Löfberg, J. (2004). In Proceedings of the CACSD Conference. Taipei, Taiwan.]) under Matlab.

3. Application to dry ice

Dry ice CO2 is a molecular crystal with four molecules per cubic unit cell (Fig. 1[link]).

[Figure 1]
Figure 1
Unit cell of dry ice: space group [Pa\bar{3}], a = 5.63 Å, C—O bond length 1.05 Å (De Smedt & Keesom, 1924[De Smedt, J. & Keesom, W. (1924). Proc. K. Akad. Wetenschappen Amsterdam, 27, 839-846.]).

3.1. Expectation values generation

For the following example, structure factors and directional Compton profiles have been generated using the Crystal14 periodic ab initio software (Dovesi, Orlando et al., 2014[Dovesi, R., Orlando, R., Erba, A., Zicovich-Wilson, C. M., Civalleri, B., Casassa, S., Maschio, L., Ferrabone, M., Pierre, M. D. L., D'Arco, P., Noel, Y., Causa, M., Rerat, M. & Kirtman, B. (2014). Int. J. Quantum Chem. 114, 1287.]; Dovesi, Saunders et al., 2014[Dovesi, R., Saunders, V. R., Roetti, C., Orlando, R., Zicovich-Wilson, C. M., Pascale, F., Civalleri, B., Doll, K., Harrison, N. M., Bush, I. J., D'Arco, P., Llunell, M., Causa, M. & Noel, Y. (2014). Crystal14 User's Manual. University of Torino, Italy.]). DFT and the B3LYP of the hybrid exchange and correlation functional have been chosen as a theoretical framework. Large polarized and diffuse atomic basis sets (triple-zeta valence with polarization quality) (Peintinger et al., 2013[Peintinger, M. F., Oliveira, D. V. & Bredow, T. (2013). J. Comput. Chem. 34, 451-459.]; Civalleri et al., 2012[Civalleri, B., Presti, D., Dovesi, R. & Savin, A. (2012). Chem. Modell. 9, 168-185.]) for both types of atoms have been used.

In the following, 1800 structure factors [(h, k, l)cubic cell [ \in Z^3] | 0 ≤ h ≤ 7, −7 ≤ k ≤ 7, −7 ≤ l ≤ 7, [\sin(\theta_{\max})/\lambda] ∼ 1.08 Å−1] and three directional Compton profiles [u = (h, k, l)cubic cell [\in] {(0, 0, 1), (1, 1, 0), (1, 1, 1)}], with a resolution of 0.15 a.u. (a.u. = atomic unit) and limited to 6 a.u., were computed.

To prove the robustness of the method, Gaussian errors have been added to the data. For each structure factor, the standard deviation is 3% of its modulus and for each directional Compton profile [J^{\bf u}(q)], it is set to be [[ J^{\bf u}(q)/\alpha_{\bf u}]^{1/2}] where [\alpha_{\bf u}] is such that [[ J^{\bf u}(0)/\alpha_{\bf u}]^{1/2} = 0.03 \times J^{\bf u}(0)]. Such distorted DCPs and SFs are illustrated, respectively, in Fig. 2[link] and in Fig. 3[link] by means of a Fourier density map. In the following, the resulting distorted DCPs and SFs will be qualified as `pseudo-experimental' data.

[Figure 2]
Figure 2
Directional Compton profile [J^{\bf u}(q)] (in red) and [{\rm tr}[\widehat{P}\widehat{J}^{\ \bf u}(q)]] (in blue) for dry ice in the crystallographic direction [{\bf u} = (1,1,1)] in the conventional cell. The spectrum is in atomic units and the profile is normalized to one electron.
[Figure 3]
Figure 3
Density map reconstructed from truncated Fourier series with coefficients [{\rm tr}\,[\widehat{P}\widehat{F}({\bf q})]] (left) and [F({\bf q})] (right) in a plane including the O—C—O bonding. Contours at intervals of ±0.01 × 2n a.u.−3 (n = 0–20): positive and negative contours are blue solid lines and red dashed lines, respectively. (Although density is a non-negative valued function, negative regions appear, as the Fourier series representing density is truncated. This is not an issue as the goal is to compare the structure factors.)

3.2. Independent molecule model

As the four CO2 molecules in the unit cell are identical and sufficiently distant from each other, each molecule can be described by the same molecular spin-traced 1-RDM in a different orientation set of local axes. Consequently, the total structure factors [F_{\rm tot}({\bf q})] and directional Compton profiles [J^{\bf u}_{\rm tot}(q)] can be computed from the molecular structure factors and directional Compton profiles [F({{\bf q}})] and [J^{\bf u}(q)] by

[F_{\rm tot}({\bf q}) = F({\bf q}) + \textstyle\sum\limits_{m = 2}^4 \exp[-i (\widehat{\Omega}_m{\bf q})\cdot{\bf r}_m] F(\widehat{\Omega}_m{\bf q})\eqno(13)]

[J^{\bf u}_{\rm tot}(q) = J^{\bf u}(q) + \textstyle\sum\limits_{m = 2}^4 J^{(\widehat{\Omega}_m{\bf u})}(q)\eqno(14)]

where [{\bf r}_{m}] and [\widehat{\Omega}_{m}] are, respectively, the translation vector and the inverse of the rotation matrix, bringing the first molecule to molecule [m\,(m\in{2,3,4})].

To assess the robustness of the method, the basis set [\{\chi _{i}\} _{{i\in\{ 1,\ldots,n\}}}] used to represent the spin-traced 1-RDM has been chosen to have fewer degrees of freedom and diffuseness than the one used to generate the expectation values [3-21G(d)] (Binkley et al., 1980[Binkley, J. S., Pople, J. A. & Hehre, W. J. (1980). J. Am. Chem. Soc. 102, 939-947.]; Schuchardt et al., 2007[Schuchardt, K. L., Didier, B. T., Elsethagen, T., Sun, L., Gurumoorthi, V., Chase, J., Li, J. & Windus, T. L. (2007). J. Chem. Inf. Model. 47, 1045-1052.]; Feller, 1996[Feller, D. (1996). J. Comput. Chem. 17, 1571-1586.]).

3.3. Results analysis

Program (11[link]) has been successfully solved for the case of dry ice. The DCPs and SFs computed with the optimized population matrix are near identical to their reference. In Fig. 2[link], one DCP derived from the 1-RDM model is plotted together with its pseudo-experimental reference for comparison (the Compton profiles anisotropies are not compared as they are negligible compared with the added random noise). The same comparison is made for the SFs in a Fourier density map in Fig. 3[link].

The inferred and the periodic ab initio spin-traced 1-RDM are in close agreement along the O—C—O bond (Fig. 4[link]). Although slight differences are observed in the off-diagonal regions, corresponding to the subtle interactions between both bonds, the general features have been accurately reproduced.

[Figure 4]
Figure 4
Spin-traced 1-RDM [\widehat{\Gamma}({\bf r},{\bf r}^{\prime})] contour maps for two different segments [as the 1-RDM is a six-variable function, no convenient graphical representation exists apart from restricting the variation of the two position vectors of [\Gamma({\bf r}, {\bf r}')] along a path]. For each segment, the position vectors [{\bf r}] (horizontal axis) and [{\bf r}^{\prime}] (vertical axis) are restricted to vary along the segment. Upper panel: along the O—C—O bonding. Lower panel: along a segment parallel and 1 a.u. away from the O—C—O bonding. Left column: inferred from position and momentum space expectation values. Right column: periodic ab initio computation. Contours at intervals of ±0.01 × 2n a.u.−3 (n = 0–20): positive and negative contours are blue solid lines and red dashed lines, respectively.

In a plane comprising the atoms of the molecule, the overall expected picture of the deformation density map, i.e. the difference between the total density and the non-interacting atom density, is recovered with minor discrepancies on the oxygen atoms and around the carbon atom (Fig. 5[link]). The fact that the axial symmetry is not obtained originates from the lack of symmetry constraints and the limited amount of experimental information (Fig. 3[link]). It could possibly be recovered by providing additional knowledge (symmetry constraints on the basis set or the population matrix for example) to the model or using more expectation values.

[Figure 5]
Figure 5
Deformation density contour map in a plane including the O—C—O bonding. Contours at intervals of ±0.01 × 2n a.u.−3(n = 0–20): positive and negative contours are blue solid lines and red dashed lines, respectively. Left: inferred from position and momentum expectation values. Right: periodic ab initio computation.

The off-diagonal regions in Fig. 4[link] are highly sensitive to the amount of noise added to the DCPs and the sharp contrast around the O—O interaction (region 5–1 a.u.) is quickly lost as the standard deviation is increased. This sensitivity might be particularly high for the case of dry ice as limited information can be deduced from DCPs because of their relatively low anistropies. Additionally, as the noise added to the SFs grows, further discrepancies appear quite naturally on the deformation density map.

Furthermore, restricting the optimization on the SFs only severely impacts the results (Fig. 6[link]) and therefore clearly illustrates the complementarity of both momentum and position expectation values as mentioned in Section 2.2[link]. Of course, restricting the optimization on the DCPs gives an even worse result.

[Figure 6]
Figure 6
SF-only inferred deformation density map in a plane including the O—C—O bonding (left) and spin-traced 1-RDM [\widehat{\Gamma}({\bf r}, {\bf r}')] along the O—C—O bonding (right). Contours at intervals of ±0.01 × 2n a.u.−3 (n = 0–20): positive and negative contours are blue solid lines and red dashed lines, respectively.

4. Conclusions

With the aim of inferring a 1-RDM from structure factors and directional Compton profiles with minimal prior knowledge, a method based on semidefinite programming was proposed. The effectiveness of this method has been evaluated on the crystal of dry ice taking periodic ab initio calculations as reference. In this example, the method was in very good agreement with the reference, showing that the use of both structure factors and directional Compton profiles provides sufficient information to infer the 1-RDM in a given atomic basis set.

Such a method could be used as a reference to interpret experimental results. For now, it is only applicable to mol­ecular crystals but it could possibly, in the future, be extended to the modelling of the 1-RDM of more general crystalline systems by using the periodic 1-RDM of the crystal as opposed to its molecular block-diagonal version as presented above.

While this method is quite general, it still depends on the choice of the atomic basis set and the atomic orbitals for each atom. This constitutes a prejudice for the electronic configuration as it truncates the Hilbert space available for the 1-RDM to be refined with respect to the experimental data. Its result can be refined through optimization of the basis, such as in the work of Gueddida et al. (2018[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.]); however, the most ideal solution would be an inference process that does not require a basis set at all. At this stage, further work is required to achieve this as, to the best of our knowledge, the N-representability conditions are more conveniently expressed in a given basis.

Ultimately, one would like to infer a 2-RDM from experimental data but technical and conceptual obstacles make the N-representability difficult to account for in this case.

Acknowledgements

The authors gratefully acknowledge Z. Yan whose program was used to compute the periodic ab initio 1-RDM along segments. Special thanks are also addressed to S. Gueddida, P. Becker and N. Ghermani for invaluable help and suggestions. BDB wishes to acknowledge G. Valmorbida for his time explaining semidefinite programming, and X. Adriaens and J.-Y. Raty for practical suggestions and comments. J-MG thanks B. Gillon, M. Souhassou, N. Claiser and C. Lecomte for fruitful discussions and help with experimental issues. We express our warm thanks to Julie McDonald and Michelle Seeto for carefully proofreading the English in this paper. The computing cluster of CentraleSupélec has been used for this work.

References

First citationApS, M. (2017). The MOSEK optimization toolbox for MATLAB manual. Version 8.1. https://docs.mosek.com/8.1/toolbox/index.htmlGoogle Scholar
First citationBenesch, R., Singh, S. & Smith, V. Jr (1971). Chem. Phys. Lett. 10, 151–153.  CrossRef CAS Google Scholar
First citationBinkley, J. S., Pople, J. A. & Hehre, W. J. (1980). J. Am. Chem. Soc. 102, 939–947.  CrossRef CAS Web of Science Google Scholar
First citationBoyd, S. & Vandenberghe, L. (2004). Convex Optimization. New York: Cambridge University Press.  Google Scholar
First citationCivalleri, B., Presti, D., Dovesi, R. & Savin, A. (2012). Chem. Modell. 9, 168–185.  CrossRef CAS Google Scholar
First citationClinton, W. L. & Massa, L. J. (1972). Phys. Rev. Lett. 29, 1363–1366.  CrossRef CAS Web of Science Google Scholar
First citationColeman, A. J. (1963). Rev. Mod. Phys. 35, 668–686.  CrossRef Web of Science Google Scholar
First citationCooper, M. J., Mijnarends, P. E., Mijnarends, P., Shiotani, N., Sakai, N. & Bansil, A. (2004). X-ray Compton Scattering. Oxford Series on Synchrotron Radiation, Vol. 5. Oxford University Press.  Google Scholar
First citationDavidson, E. (1976). Quantum Chemistry. New York: Academic Press.  Google Scholar
First citationDe Smedt, J. & Keesom, W. (1924). Proc. K. Akad. Wetenschappen Amsterdam, 27, 839–846.  CAS Google Scholar
First citationDeutsch, 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
First citationDeutsch, 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
First citationDobrzynski, L. & Zukowski, E. (1999). J. Phys. Condens. Matter, 11, 8049–8060.  CrossRef CAS Google Scholar
First citationDovesi, R., Orlando, R., Erba, A., Zicovich-Wilson, C. M., Civalleri, B., Casassa, S., Maschio, L., Ferrabone, M., Pierre, M. D. L., D'Arco, P., Noel, Y., Causa, M., Rerat, M. & Kirtman, B. (2014). Int. J. Quantum Chem. 114, 1287.  CrossRef Google Scholar
First citationDovesi, R., Saunders, V. R., Roetti, C., Orlando, R., Zicovich-Wilson, C. M., Pascale, F., Civalleri, B., Doll, K., Harrison, N. M., Bush, I. J., D'Arco, P., Llunell, M., Causa, M. & Noel, Y. (2014). Crystal14 User's Manual. University of Torino, Italy.  Google Scholar
First citationFeller, D. (1996). J. Comput. Chem. 17, 1571–1586.  CrossRef CAS Google Scholar
First citationGilbert, T. L. (1975). Phys. Rev. B, 12, 2111–2120.  CrossRef Web of Science Google Scholar
First citationGillet, J.-M. (2007). Acta Cryst. A63, 234–238.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationGillet, J.-M. & Becker, P. J. (2004). J. Phys. Chem. Solids, 65, 2017–2023.  Web of Science CrossRef CAS Google Scholar
First citationGillet, J.-M., Becker, P. J. & Cortona, P. (2001). Phys. Rev. B, 63, 235115.  Web of Science CrossRef Google Scholar
First citationGillet, J.-M., Fluteaux, C. & Becker, P. J. (1999). Phys. Rev. B, 60, 2345–2349.  CrossRef CAS Google Scholar
First citationGueddida, 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
First citationHansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909–921.  CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationLathiotakis, N. N. & Marques, M. A. L. (2008). J. Chem. Phys. 128, 184103.  Web of Science CrossRef PubMed Google Scholar
First citationLöfberg, J. (2004). In Proceedings of the CACSD Conference. Taipei, Taiwan.  Google Scholar
First citationLöwdin, P. (1950). J. Chem. Phys. 18, 365–375.  Google Scholar
First citationLöwdin, P.-O. (1955). Phys. Rev. 97, 1474–1489.  Google Scholar
First citationMazziotti, D. A. (2007). ESAIM: M2AN, 41, 249–259.  Google Scholar
First citationMcWeeny, R. (1960). Rev. Mod. Phys. 32, 335–369.  CrossRef Web of Science Google Scholar
First citationOszlányi, G. & Sütő, A. (2004). Acta Cryst. A60, 134–141.  Web of Science CrossRef IUCr Journals Google Scholar
First citationPattison, P. & Weyrich, W. (1979). J. Phys. Chem. Solids, 40, 213–222.  CrossRef CAS Web of Science Google Scholar
First citationPeintinger, M. F., Oliveira, D. V. & Bredow, T. (2013). J. Comput. Chem. 34, 451–459.  Web of Science CrossRef CAS PubMed Google Scholar
First citationPillet, S., Souhassou, M., Pontillon, Y., Caneschi, A., Gatteschi, D. & Lecomte, C. (2001). New J. Chem. 25, 131–143.  Web of Science CSD CrossRef CAS Google Scholar
First citationPisani, C. (2012). Quantum-mechanical ab-initio Calculation of the Properties of Crystalline Materials. Lecture Notes in Chemistry, Vol. 67. Berlin: Springer Science & Business Media.  Google Scholar
First citationSchmider, H., Smith, V. H. & Weyrich, W. (1992). J. Chem. Phys. 96, 8986–8994.  CrossRef CAS Web of Science Google Scholar
First citationSchuchardt, K. L., Didier, B. T., Elsethagen, T., Sun, L., Gurumoorthi, V., Chase, J., Li, J. & Windus, T. L. (2007). J. Chem. Inf. Model. 47, 1045–1052.  Web of Science CrossRef PubMed CAS Google Scholar
First citationSivia, D. & Skilling, J. (2006). Data Analysis: A Bayesian Tutorial. Oxford University Press.  Google Scholar
First citationTsirelson, V. G. & Ozerov, R. P. (1996). Electron Density and Bonding in Crystals: Principles, Theory and X-ray Diffraction Experiments in Solid State Physics and Chemistry. New York: Taylor and Francis.  Google Scholar
First citationVandenberghe, L. & Boyd, S. (1996). SIAM Rev. 38, 49–95.   CrossRef Google Scholar
First citationWeyrich, W., Pattison, P. & Williams, B. (1979). Chem. Phys. 41, 271–284.  CrossRef CAS Google Scholar
First citationWolkowicz, H., Saigal, R. & Vandenberghe, L. (2012). Handbook of Semidefinite Programming: Theory, Algorithms, and Applications. International Series in Operations Research and Management Science, Vol. 27. New York: Springer Science & Business Media.  Google Scholar
First citationZhang, F. (2006). The Schur Complement and its Applications. Numerical Algorithms, Vol. 4. New York: Springer Science & Business Media.  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.

Journal logoFOUNDATIONS
ADVANCES
ISSN: 2053-2733
Follow Acta Cryst. A
Sign up for e-alerts
Follow Acta Cryst. on Twitter
Follow us on facebook
Sign up for RSS feeds