Support complex values in the QCSchema format

[S-Erik]

Environment

• Qiskit Nature version: 0.7.2
• Python version: 3.11.5
• Operating system: Linux

What is happening?

Summary

The get_overlap_ab_from_qcschema function calculates the molecular orbital (MO) overlap matrix from the atomic orbital (AO) overlap matrix and the MO coefficients. The function returns

return coeff_a.T @ overlap @ coeff_b

where overlap is the AO overlap matrix, and coeff_a and coeff_b are the MO coefficients of the $\alpha$ (up)- and $\beta$ (down)-spin MOs. For real coefficients coeff_a and coeff_b this is correct. For complex coefficients coeff_a and coeff_b the following should be correct (and is also correct for real coefficients):

return coeff_a.conj().T @ overlap @ coeff_b

Mathematical Background

The get_overlap_ab_from_qcschema function calculates the overlap matrix
$$O_{ij}^\mathrm{MO}=c^T_{pi}\cdot O_{pq}^\mathrm{AO}\cdot c_{qj},$$
where we abuse notation by mixing matrix multiplication and index notation. $\cdot$ denotes matrix multiplication, $O_{ij}^\mathrm{MO}$($O_{pq}^\mathrm{AO}$) are the MO(AO) overlap matrix entries, and $c_{pi}$ are the MO coefficients. Further, $i$ and $j$ denote different MOs while $p$ and $q$ denote different AOs. This can be written with sums instead of matrix multiplication as follows
$$O_{ij}^\mathrm{MO}=\sum_{p,q}c_{pi}c_{qj}O_{pq}^\mathrm{AO},$$

But on the other hand $O_{ij}^\mathrm{MO}$ should be defined as
$$O_{ij}^\mathrm{MO}=\braket{i|j}.$$
We can now expand the MOs in terms of AOs:
$$\ket{j}=\sum_q c_{qj}\ket{q}$$
with $c_{qj}=\braket{q|j}$. This implies
$$\bra{j}=\sum_q c^\ast_{qj}\bra{q}.$$
With this we find
$$O_{ij}^\mathrm{MO}=\braket{i|j}=\sum_{p,q}c^\ast_{pi}c_{qj}\braket{p|q}\sum_{p,q}c^\ast_{pi}c_{qj}O_{pq}^\mathrm{AO},$$
where we used $O_{pq}^\mathrm{AO}=\braket{p|q}$. Writing the above equation using matrix multiplication we find
$$O_{ij}^\mathrm{MO}=c^\dagger_{pi}\cdot O_{pq}^\mathrm{AO}\cdot c_{qj},$$
where $c^\dagger_{pi}=(c^T_{pi})^\ast$.
From this it is obvious to me that

return coeff_a.conj().T @ overlap @ coeff_b

should be the correct implementation.

How can we reproduce the issue?

I do not think a minimal working example is necessary here.

What should happen?

Change

return coeff_a.T @ overlap @ coeff_b

to

return coeff_a.conj().T @ overlap @ coeff_b

in get_overlap_ab_from_qcschema.

Any suggestions?

No response

[S-Erik]

Implementing the suggested change seems to cause problems with the AngularMomentum overlap setter function:

    @overlap.setter
    def overlap(self, overlap: np.ndarray | None) -> None:
        self._overlap = overlap

        if overlap is not None:
            norb = self.num_spatial_orbitals
            delta = np.eye(2 * norb)
            delta[:norb, :norb] -= overlap.T @ overlap
            delta[norb:, norb:] -= overlap @ overlap.T
            summed = np.einsum("ij->", np.abs(delta))
            if not np.isclose(summed, 0.0, atol=1e-6):
                LOGGER.warning(
                    "The provided alpha-beta overlap matrix is NOT unitary! This can happen when "
                    "the alpha- and beta-spin orbitals do not span the same space. To provide an "
                    "example of what this means, consider an active space chosen from unrestricted-"
                    "spin orbitals. Computing <S^2> within this active space may not result in the "
                    "same <S^2> value as obtained on the single-reference starting point. More "
                    "importantly, this implies that the inactive subspace will account for the "
                    "difference between these two <S^2> values, possibly resulting in significant "
                    "spin contamination in both subspaces. You should verify whether this is "
                    "intentional/acceptable or whether your choice of active space can be improved."
                    " As a reference, here is the summed-absolute deviation of `S^T @ S` from the "
                    "identity: %s",
                    str(summed),
                )

This function checks if the overlap matrix is unitary, which was added in qiskit_nature verion 0.7.2, see #1292. This is done by using a helper array delta. Since the suggested change turns the overlap in a numpy array of dtype complex while delta is of dtype float, this results in the error
UFuncTypeError: Cannot cast ufunc 'subtract' output from dtype('complex128') to dtype('float64') with casting rule 'same_kind'.

What should happen?

Change

return coeff_a.T @ overlap @ coeff_b

to

return (coeff_a.conj().T @ overlap @ coeff_b).real

in get_overlap_ab_from_qcschema since the overlap matrix is real and casting to real values does not change the matrix.

[mrossinek]

While your detailed analysis is indeed correct, the original reason why this was not implemented via the matrix adjoint, is simply the fact that the QCSchema does not support complex values. At least, this has never been tested or verified.

If a PR were to address this, I would expect the QCSchema data structure to be validated for complex rather than real values, too. Once validated (and unittested), the relevant type hints would need to be updated from float to complex indicating this as officially supported.

Note, that complex values in Python are not natively serializable in json. So the to_json and from_json methods of all QCSchema-related classes will need to account for this, too.

Given that complex values are not advertised as being supported, I am changing this label from a bug to an enhancement.