PennyLaneAI / pennylane

PennyLane is a cross-platform Python library for quantum computing, quantum machine learning, and quantum chemistry. Train a quantum computer the same way as a neural network.
https://pennylane.ai
Apache License 2.0
2.18k stars 569 forks source link

Measuring fermionic reduced density matrices #4074

Open Emieeel opened 1 year ago

Emieeel commented 1 year ago

Feature details

We're working on algorithms for quantum chemistry and I am using PennyLane to calculate 1- and 2-particle reduced density matrices that are auto-differentiable w.r.t. PQC parameters:

$$ \gamma{pq}(\boldsymbol{\theta}) = \langle \psi(\boldsymbol{\theta})|E{pq} | \psi(\boldsymbol{\theta}) \rangle \quad \quad \Gamma{pqrs}(\boldsymbol{\theta}) = \langle \psi(\boldsymbol{\theta})|e{pqrs} | \psi(\boldsymbol{\theta}) \rangle $$

where

$$ E{pq} = \sum{\sigma} a^\dagger{p\sigma} a{q\sigma} \quad \text{and} \quad e{pqrs} = \sum{\sigma \tau} a^\dagger{p \sigma} a^\dagger{r \tau} a{s \tau} a{q \sigma} = E{pq}E{rs} - \delta{qr}E{ps} $$

I put them here in the restricted format with the sum over spin (and "chemist" ordering of indices), but one could also estimate them in unrestricted form. It would be great to be able to have them directly from a QNode output and integrate it in PennyLane. Something like:

dev = qml.device("default.qubit", wires=N)
@qml.qnode(dev)
def circuit_ansatz(theta, *args):
    qml.UCCSD(theta, *args)
    return qml.fermion_rdm()

I anticipate this will be useful in a lot of quantum computing for quantum chemistry/physics projects. Any standard quantum chemistry code allows you to extract them from a CI-solver, so a quantum analogue would be great. Apart from being able to estimate any fermionic two-body observable, such as electronic structure Hamiltonians but also dipole moments, etc. with them, they are fundamental for a lot of different classical (post-)processing techniques like orbital optimization. PennyLane would be a nice platform to enable extracting gradients or Hessians of the RDMs with respect to circuit parameters.

Implementation

The simple implementation that I have is to estimate the RDMs as a function of a state, coming from some circuit QNode using qml.state(). I do this by extracting $E{pq}$ and $e{pqrs}$ as sparse matrices using OpenFermion and computing all the expectation values w.r.t. the state.

One way to go could be to implement a custom measurement process in PennyLane, such that one can just call a function like qml.fermion_rdm() and have a qnode that output the RDMs directly (containing the measurements of the observables $E{pq}$ and $e{pqrs}$). This could use the process_state that computes expectation values of the sparse matrices as I did above. Maybe there is a a smarter/more efficient way to do it already using existing machinery (using qml.expval?) and/or integrating it into qml.qchem somehow.

Feedback for how to proceed would be greatly appreciated. I would be happy to work/contribute on this if there is interest.

How important would you say this feature is?

3: Very important! Blocking work.

Additional information

No response

isaacdevlugt commented 1 year ago

Hey @Emieeel! Interesting... We can look into this as a possible pennylane feature! In the meantime, have you checked out our newer custom measurements features? We released this in v0.28:

https://docs.pennylane.ai/en/stable/development/release_notes.html#release-0-28-0

Might be worth seeing if you can implement what you need with what's there!

josh146 commented 1 year ago

@Emieeel in addition, we are working to ensure that devices will be able to utilize the custom processing that is included in the new custom measurement processes :) This is definitely where it would be super valuable to get your feedback, since it will help guide our development of measurement processes!

A couple of questions regarding fermion_rdm:

Emieeel commented 1 year ago

Hi @isaacdevlugt and @josh146, thanks for your response. Good to know that there is interest!

I indeed was referring to the new custom measurements features, I think it should be quite straightforward to implement a state measurement with the code I already have. A thing to note is that it needs OpenFermion right now to extract sparse matrices to act on the state. I used sparse matrices as that is a huge save in memory and time. I extract them like so:

def e_pq(p, q, restricted=True):
    r"""
    Can generate either spin-unrestricted single excitation operator:

    .. math::
        E_{pq} = a_{p}^\dagger a_{q}
    where :math:`p` and :math:`q` are composite spatial/spin indices,
    or spin-restricted single excitation operator:

    .. math::
            E_{pq} = \sum_\sigma a_{p \sigma}^\dagger a_{q \sigma}
    where :math:`p` and :math:`q` are spatial indices. For the spin-
    restricted case, One can either select up-then-down convention,
    or up-down-up-down.
    """
    if restricted:
        operator = (openfermion.FermionOperator(f'{2*p}^ {2*q}') +
                    openfermion.FermionOperator(f'{2*p+1}^ {2*q+1}'))
    else:
        operator = openfermion.FermionOperator(f'{p}^ {q}')

    return operator

def e_pqrs(p, q, r, s, n_modes, restricted=True, up_then_down=False):
    r"""
    Can generate either spin-unrestricted double excitation operator:

    .. math::
        e_{pqrs} = a_{p}^\dagger a_{q}^\dagger a_{r} a_{s}
    where :math:`p` and :math:`q` are composite spatial/spin indices,
    or spin-restricted double excitation operator:

    .. math::
        e_{pqrs} = \sum_{\sigma \tau} a_{p \sigma}^\dagger a_{
            r \tau}^\dagger a_{s \tau} a_{q \sigma}
                 = E_{pq}E_{rs} -\delta_{qr}E_{ps}
    where the indices are spatial indices.

    Indices in the restricted case are in chemist order, meant to be
    contracted with two-electron integrals in chemist order to obtain the
    Hamiltonian or to obtain the two-electron RDM.
    """
    if restricted:
        operator = e_pq(p, q, restricted) * e_pq(
            r, s, restricted)
        if q == r:
            operator += - e_pq(p, s, restricted)
    else:
        operator = openfermion.FermionOperator(f'{p}^ {q}^ {r} {s}')
    return operator

def initialize_e_pq(ncas, restricted=True):
    """Initialize full e_pq operator"""
    if restricted:
        num_ind = ncas
    else:
        num_ind = 2 * ncas
    return [[openfermion.get_sparse_operator(
        e_pq(p, q, restricted), n_qubits=2*ncas)
        for q in range(num_ind)] for p in range(num_ind)]

def initialize_e_pqrs(ncas, restricted=True):
    """Initialize full e_pqrs operator"""
    if restricted:
        num_ind = ncas
    else:
        num_ind = 2 * ncas
    return [[[[openfermion.get_sparse_operator(
        e_pqrs(p, q, r, s, restricted), n_qubits=2*ncas)
        for s in range(num_ind)] for r in range(num_ind)]
        for q in range(num_ind)] for p in range(num_ind)]

Is there a native PennyLane solution to this? Another solution is not converting the operators to sparse matrices but use PL operators such that we can calculate the elements with qml.expval, could that be more efficient? The code above already becomes a bit slow when you have an active space of >4 (spatial) orbitals.

To answer @josh146: