research papers
Inferring the one-electron reduced density matrix of molecular crystals from experimental data sets through semidefinite programming
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
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.
Keywords: one-electron reduced density matrix; X-ray diffraction; Compton diffusion; semidefinite programming.
1. Introduction
The computation of one-electron expectation values such as the mean position, the mean momentum or the mean ; Coleman, 1963; McWeeny, 1960; Davidson, 1976). This quantity provides a quantum description of an average electron and has been proved to be sufficient (Lathiotakis & Marques, 2008; Gilbert, 1975). 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.
of electrons in a crystal does not require more than the mere knowledge of the one-electron reduced density matrix (1-RDM) (Löwdin, 1955Several models have been proposed to approximate and refine a 1-RDM from experimental expectation values (Deutsch et al., 2012, 2014; Hansen & Coppens, 1978; Gillet et al., 1999, 2001; Gillet, 2007; Gillet & Becker, 2004; Gueddida et al., 2018; Pillet et al., 2001; Schmider et al., 1992; Clinton & Massa, 1972; Tsirelson & Ozerov, 1996). The complementarity between position and momentum space expectation values in the description of the 1-RDM is now well accepted (Cooper et al., 2004; Pisani, 2012). 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 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 (MEM) and methods (Oszlányi & Sütő, 2004; Dobrzynski & Zukowski, 1999; Binkley et al., 1980). 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), 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), a recent subfield of convex optimization (Boyd & Vandenberghe, 2004) which is of growing interest in systems and control theory, geometry and statistics (Wolkowicz et al., 2012).
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), 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 be a set of atomic orbitals describing the electrons of each atom taken as an independent system. From , one can deduce an orthogonal basis set for the molecule, using the Löwdin orthogonalization procedure (Löwdin, 1950) for example. Expanding the spin-traced 1-RDM in such a basis, one reveals its basis-set representation: the population matrix , so that
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 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
can be calculated from its operator applied to the 1-RDM :where means that the operator only acts on variable . By defining the basis-set representation of as
one can conveniently write the )], where tr is the matrix trace operator.
as = [using equation (1In particular, in position space, the X-ray structure factors , which are given by
have an operator whose basis-set representation is
In momentum space, the directional Compton profiles can be defined through the autocorrelation function (Pattison & Weyrich, 1979; Weyrich et al., 1979; Benesch et al., 1971) as
Their operator basis-set representation is therefore
From equations (4) and (7)–(8), 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 so that it fits given independent expectation values . In the following, the expectation values are SFs and DCPs data. Supposing the latter follow Gaussian error distributions with standard deviations and no a priori knowledge is given on , the problem is equivalent to minimizing the so-called function with respect to the elements of (Gillet & Becker, 2004; Sivia & Skilling, 2006). It can be summarized in the following optimization program:
where I is the identity matrix and the notation 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 eigenvalues of must lie in [0, 2].
The following passage will cast program (10) as a semidefinite optimization program. These steps are quite standard in the field of convex optimization (Boyd & Vandenberghe, 2004). Introducing a new variable t, program (10) is equivalent to
where is a column vector whose elements are [ ]/ and is the Euclidean norm.
Using Schur's complement (Zhang, 2006), the last constraint of program (11) can be written as a linear matrix inequality:
where I is the identity matrix of appropriate dimensions. This inequality is indeed linear with respect to as is a linear function of .
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). 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).
In the present work, this program has been addressed by using the optimization software Mosek (ApS, 2017) interfaced by the Yalmip toolbox (Löfberg, 2004) under Matlab.
3. Application to dry ice
Dry ice CO2 is a molecular crystal with four molecules per cubic (Fig. 1).
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, Saunders et al., 2014). 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; Civalleri et al., 2012) for both types of atoms have been used.
In the following, 1800 structure factors [(h, k, l)cubic cell | 0 ≤ h ≤ 7, −7 ≤ k ≤ 7, −7 ≤ l ≤ 7, ∼ 1.08 Å−1] and three directional Compton profiles [u = (h, k, l)cubic cell {(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 and in Fig. 3 by means of a Fourier density map. In the following, the resulting distorted DCPs and SFs will be qualified as `pseudo-experimental' data.
the standard deviation is 3% of its modulus and for each directional Compton profile , it is set to be where is such that . Such distorted DCPs and SFs are illustrated, respectively, in Fig. 2 | Figure 2 The spectrum is in and the profile is normalized to one electron. |
3.2. Independent molecule model
As the four CO2 molecules in the 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 and directional Compton profiles can be computed from the molecular structure factors and directional Compton profiles and by
where and are, respectively, the translation vector and the inverse of the rotation matrix, bringing the first molecule to molecule .
To assess the robustness of the method, the basis set used to represent the spin-traced 1-RDM has been chosen to have fewer et al., 1980; Schuchardt et al., 2007; Feller, 1996).
and diffuseness than the one used to generate the expectation values [3-21G(d)] (Binkley3.3. Results analysis
Program (11) 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, 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.
The inferred and the periodic ab initio spin-traced 1-RDM are in close agreement along the O—C—O bond (Fig. 4). 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.
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). 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). 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.
The off-diagonal regions in Fig. 4 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) and therefore clearly illustrates the complementarity of both momentum and position expectation values as mentioned in Section 2.2. Of course, restricting the optimization on the DCPs gives an even worse result.
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 molecular 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 et al. (2018); 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.
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 GueddidaUltimately, 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
ApS, M. (2017). The MOSEK optimization toolbox for MATLAB manual. Version 8.1. https://docs.mosek.com/8.1/toolbox/index.html. Google Scholar
Benesch, R., Singh, S. & Smith, V. Jr (1971). Chem. Phys. Lett. 10, 151–153. CrossRef CAS Google Scholar
Binkley, J. S., Pople, J. A. & Hehre, W. J. (1980). J. Am. Chem. Soc. 102, 939–947. CrossRef CAS Web of Science Google Scholar
Boyd, S. & Vandenberghe, L. (2004). Convex Optimization. New York: Cambridge University Press. Google Scholar
Civalleri, B., Presti, D., Dovesi, R. & Savin, A. (2012). Chem. Modell. 9, 168–185. CrossRef CAS Google Scholar
Clinton, W. L. & Massa, L. J. (1972). Phys. Rev. Lett. 29, 1363–1366. CrossRef CAS Web of Science Google Scholar
Coleman, A. J. (1963). Rev. Mod. Phys. 35, 668–686. CrossRef Web of Science Google Scholar
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. Google Scholar
Davidson, E. (1976). Quantum Chemistry. New York: Academic Press. Google Scholar
De Smedt, J. & Keesom, W. (1924). Proc. K. Akad. Wetenschappen Amsterdam, 27, 839–846. CAS 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
Dobrzynski, L. & Zukowski, E. (1999). J. Phys. Condens. Matter, 11, 8049–8060. CrossRef CAS Google Scholar
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. CrossRef Google Scholar
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. Google Scholar
Feller, D. (1996). J. Comput. Chem. 17, 1571–1586. CrossRef CAS 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
Gillet, J.-M., Fluteaux, C. & Becker, P. J. (1999). Phys. Rev. B, 60, 2345–2349. CrossRef CAS 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
Lathiotakis, N. N. & Marques, M. A. L. (2008). J. Chem. Phys. 128, 184103. Web of Science CrossRef PubMed Google Scholar
Löfberg, J. (2004). In Proceedings of the CACSD Conference. Taipei, Taiwan. Google Scholar
Löwdin, P. (1950). J. Chem. Phys. 18, 365–375. Google Scholar
Löwdin, P.-O. (1955). Phys. Rev. 97, 1474–1489. Google Scholar
Mazziotti, D. A. (2007). ESAIM: M2AN, 41, 249–259. Google Scholar
McWeeny, R. (1960). Rev. Mod. Phys. 32, 335–369. CrossRef Web of Science Google Scholar
Oszlányi, G. & Sütő, A. (2004). Acta Cryst. A60, 134–141. Web of Science CrossRef IUCr Journals Google Scholar
Pattison, P. & Weyrich, W. (1979). J. Phys. Chem. Solids, 40, 213–222. CrossRef CAS Web of Science Google Scholar
Peintinger, M. F., Oliveira, D. V. & Bredow, T. (2013). J. Comput. Chem. 34, 451–459. Web of Science CrossRef CAS PubMed Google Scholar
Pillet, 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
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. Google Scholar
Schmider, H., Smith, V. H. & Weyrich, W. (1992). J. Chem. Phys. 96, 8986–8994. CrossRef CAS Web of Science Google Scholar
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. Web of Science CrossRef PubMed CAS Google Scholar
Sivia, D. & Skilling, J. (2006). Data Analysis: A Bayesian Tutorial. Oxford University Press. Google Scholar
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. Google Scholar
Vandenberghe, L. & Boyd, S. (1996). SIAM Rev. 38, 49–95. CrossRef Google Scholar
Weyrich, W., Pattison, P. & Williams, B. (1979). Chem. Phys. 41, 271–284. CrossRef CAS Google Scholar
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. Google Scholar
Zhang, 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.