qiskit-community / qiskit-nature

Qiskit Nature is an open-source, quantum computing, framework for solving quantum mechanical natural science problems.
https://qiskit-community.github.io/qiskit-nature/
Apache License 2.0
306 stars 207 forks source link

Support complex values in the QCSchema format #1351

Open S-Erik opened 8 months ago

S-Erik commented 8 months ago

Environment

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}=\sumq c{qj}\ket{q}$$ with $c_{qj}=\braket{q|j}$. This implies $$\bra{j}=\sumq 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 commented 8 months ago

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 commented 8 months ago

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.