Inferring crystal electronic properties from experimental data sets through Semidefinite Programming

TL;DR

This paper introduces a semidefinite programming approach to infer the one-electron reduced density matrix of crystals from experimental data, minimizing assumptions and overcoming previous constraints, with successful application to dry ice.

Contribution

The paper presents a novel semidefinite programming framework for inferring the 1-RDM from experimental data with minimal assumptions, improving upon existing methods.

Findings

Successfully applied to dry ice crystal data

Provides results comparable to ab-initio calculations

Addresses constraints that limited previous models

Abstract

Constructing a quantum description of crystals from scattering experiments is of paramount importance 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. We have 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.

Equations

$$\hat{\Gamma }(\mathbf{r},{\mathbf{r}}^{\prime })=\sum _{i,j}^n{\hat{\mathrm{P}}}_{\mathrm{ij}}{\varphi }_i^{\ast }(\mathbf{r}){\varphi }_j({\mathbf{r}}^{\prime })$$

$$⟨\hat{O}⟩=\int {\left({\hat{\mathcal{O}}}_{{\mathbf{r}}^{\prime }}\hat{\Gamma }(\mathbf{r},{\mathbf{r}}^{\prime })\right)}_{{\mathbf{r}}^{\prime }=\mathbf{r}}d\mathbf{r}$$

$${\hat{\mathrm{O}}}_{\mathrm{ij}}=\int {\varphi }_i^{\ast }(\mathbf{r}){\left({\hat{\mathcal{O}}}_{{\mathbf{r}}^{\prime }}{\varphi }_j({\mathbf{r}}^{\prime })\right)}_{{\mathbf{r}}^{\prime }=\mathbf{r}}d\mathbf{r}$$

$$F(\mathbf{q})=$$

$$=$$

$${\hat{\mathrm{F}}}_{\mathrm{ij}}(\mathbf{q})=$$

$$J^{\mathbf{u}}(q)=$$

$$B(\mathbf{r})=$$

$${\hat{\mathrm{J}}}_{\mathrm{ij}}^{\mathbf{u}}(\mathrm{q})=\int \int \frac{1}{2\pi }{\varphi }_i^{\ast }({\mathbf{r}}^{\prime }){\varphi }_j(t\mathbf{u}+{\mathbf{r}}^{\prime })e^{-itq}dtd{\mathbf{r}}^{\prime }$$

$$\displaystyle\left[\begin{array}[]{cc}\mathrm{I}&{\mathbf{\Delta_{\sigma}O}}\\ ({\mathbf{\Delta_{\sigma}O}})^{T}&t\end{array}\right]$$

$$F_{tot}(\mathbf{q})=$$

$$J_{tot}^{\mathbf{u}}(q)=$$

Full text

I 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); Coleman (1963); McWeeny (1960); Davidson (1976). This quantity provides a quantum description of an average electron and has been proved to be sufficient Lathiotakis and 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, using 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 non-zero temperature.

Several models have been proposed to approximate and refine a 1-RDM from experimental expectation values Deutsch et al. (2012, 2014); Hansen and Coppens (1978); Gillet et al. (2001); Gillet (2007); Gillet and Becker (2004); Gillet et al. (1999); Gueddida et al. (2018); Pillet et al. (2001); Schmider et al. (1992); Clinton and Massa (1972); Tsirelson and 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 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. 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 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 and Boyd (1996), a recent subfield of convex optimization Boyd and Vandenberghe (2004) which is of growing interest in Systems & Control Theory, Geometry and Statistics Wolkowicz et al. (2012).

II Method

II.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 Sec.III), 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 Löwdin orthogonalization procedure Löwdin (1950) for example. Expanding the spin-traced 1-RDM $\widehat{\Gamma}({\mathbf{r}},{\mathbf{r}}^{\prime})$ in such a basis, one reveals its basis set representation: the population matrix $\mathrm{\widehat{P}}$, so that:

[TABLE]

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 $\mathrm{\widehat{P}}$ and state that they must lie in $[0,2]$ (as $N$ is even) and their sum must be equal to $N$.

II.2 Expectation values

Any one-electron expectation value $\langle\widehat{O}\rangle$ can be calculated from its operator $\mathcal{\widehat{O}_{{\mathbf{r}}^{\prime}}}$ applied to the 1-RDM $\widehat{\Gamma}({\mathbf{r}},{\mathbf{r^{\prime}}})$:

[TABLE]

where $\mathcal{\widehat{O}_{{\mathbf{r}}^{\prime}}}$ means that the operator only acts on variable ${\mathbf{r}}^{\prime}$. By defining, the basis set representation of $\mathcal{\widehat{O}}$ as:

[TABLE]

one can conveniently write the expectation value as $\langle\widehat{O}\rangle=\mathrm{tr(\widehat{P}\widehat{O})}$ (using Eq.1), where $\mathrm{tr}$ is the matrix trace operator.

In particular, in position space, the X-ray structure factors $F({\mathbf{q}})$, which are given by:

[TABLE]

have an operator whose basis set representation is:

[TABLE]

In momentum space, the directional Compton profiles $J^{{\mathbf{u}}}(q)$ can be defined through the autocorrelation function $B({\mathbf{r}})$ Pattison and Weyrich (1979); Weyrich et al. (1979); Benesch et al. (1971) as:

[TABLE]

Their operator basis set representation is therefore:

[TABLE]

From Eq.4 and Eq.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.

II.3 Constrained least-squares fitting scheme

In the Bayesian sense, the objective is to infer the most probable population matrix $\mathrm{\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 $\mathrm{\widehat{P}}$, the problem is equivalent to minimizing the so-called $\chi^{2}$ function with respect to the elements of $\mathrm{\widehat{P}}$ Gillet and Becker (2004); Sivia and Skilling (2006). It can be summarized in the following optimization program:

{mini}^P

∑α1σα2— ⟨^Oα⟩- tr(^P^O_α)—^2

\addConstrainttr(^P)= N \addConstraint^P ≽0 \addConstraint2  I-^P ≽0

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

The following passage will cast program (II.3) as a semidefinite optimization program. These steps are quite standard in the field of convex optimization Boyd and Vandenberghe (2004). Introducing a new variable $t$, program (II.3) is equivalent to: {mini} ^P,tt

\addConstrainttr(^P)= N \addConstraint^P ≽0 \addConstraint2  I-^P ≽0 \addConstraintt - ——Δ_σO——^2 ≥0

where ${\mathbf{\Delta_{\sigma}O}}$ is a column vector whose elements are $(\langle\widehat{O}_{\alpha}\rangle-\mathrm{tr(\widehat{P}\widehat{O}_{\alpha})})/\sigma_{\alpha}$ and $||\cdot||$ is the euclidean norm.

Using Schur’s complement Zhang (2006), the last constraint of program (II.3) can be written as a linear matrix inequality:

[TABLE]

where $\mathrm{I}$ is the identity matrix of appropriate dimensions. This inequality is indeed linear with respect to $\widehat{\mathrm{P}}$ as ${\mathbf{\Delta_{\sigma}O}}$ is a linear function of $\widehat{\mathrm{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 and Boyd (1996). Interior-point algorithms can be used to solve this class of problems and no initial guess is required. Treatment of the 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 Yalmip toolbox Löfberg (2004) under Matlab.

III Application to dry ice

Dry ice $\mathrm{CO_{2}}$ is a molecular crystal with four molecules per cubic unit cell (Fig.1).

III.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 et al. (2014a, b). Density Functional Theory and the B3LYP of 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)_{\text{cubic cell}}\in\mathrm{Z}^{3}\,|\,0\leq h\leq 7,-7\leq k\leq 7,-7\leq l\leq 7$, $\sin(\theta_{max})/\lambda\sim 1.08$ Å*-1*) and three directional Compton profiles ($u=(h,k,l)_{\text{cubic cell}}\in\{(0,0,1),(1,1,0),(1,1,1)\}$), with a resolution of $0.15$ a.u. 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^{{\mathbf{u}}}(q)$, it is set to be $\sqrt{J^{{\mathbf{u}}}(q)/\alpha_{{\mathbf{u}}}}$ where $\alpha_{{\mathbf{u}}}$ is such that $\sqrt{J^{{\mathbf{u}}}(0)/\alpha_{{\mathbf{u}}}}=0.03\times J^{{\mathbf{u}}}(0)$. Such distorted DCPs and SFs are illustrated respectively in Fig.2 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.

III.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_{tot}({\mathbf{q}})$ and directional Compton profiles $J^{{\mathbf{u}}}_{tot}(q)$ can be computed from the molecular structure factors and directional Compton profiles $F({\mathbf{q}})$ and $J^{{\mathbf{u}}}(q)$ by:

[TABLE]

where ${\mathbf{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$-$21$G(d)) Binkley et al. (1980); Schuchardt et al. (2007); Feller (1996).

III.3 Results analysis

Program (II.3) 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 (see Ancillary Material for the other two DCPs). 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 of 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) 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$ a.u. - $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 grow, 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 Sec.II.2. Of course, restricting the optimization on the DCPs gives an even worse result.

IV Conclusion

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 modeling of 1-RDM of more general crystalline systems.

While this method is quite general, it still depends on the choice of the atomic basis set. Its result can be refined through optimization of the basis, such as in Gueddida et al. (2018), however the most ideal solution would be an inference process that does not require a basis set altogether. 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.

V Acknowledgments

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, N. Ghermani for invaluable help and suggestions. B.D. wishes to acknowledge G. Valmorbida for his time explaining Semidefinite Programming and X. Adriaens, J-Y Raty for practical suggestions and comments. J-M.G. thanks B. Gillon, M. Souhassou, N. Claiser and C. Lecomte for fruitful discussions and 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.