Add converter function for PySCF's CISD to PennyLane statevector

doc/releases/changelog-dev.md

@@ -37,6 +37,9 @@ def circuit():
 >>> circuit(shots=1)
 array([False, False])
 
+* Functions added to convert wavefunctions obtained from `PySCF` to a state vector.
+  [(#4427)](https://github.com/PennyLaneAI/pennylane/pull/4427)
+
 <h3>Improvements 🛠</h3>
 
 * Transform Programs, `qml.transforms.core.TransformProgram`, can now be called on a batch of circuits
@@ -278,6 +281,7 @@ array([False, False])
 This release contains contributions from (in alphabetical order):
 
 Isaac De Vlugt,
+Stepan Fomichev,
 Lillian M. A. Frederiksen,
 Soran Jahangiri,
 Edward Jiang,

pennylane/qchem/__init__.py

@@ -25,7 +25,7 @@
 )
 from .basis_data import load_basisset
 from .basis_set import BasisFunction, atom_basis_data, mol_basis_data
-from .convert import import_operator
+from .convert import import_operator, import_state
 from .dipole import dipole_integrals, fermionic_dipole, dipole_moment
 from .factorization import basis_rotation, factorize
 from .hamiltonian import electron_integrals, fermionic_hamiltonian, diff_hamiltonian

pennylane/qchem/convert.py

@@ -15,13 +15,14 @@
 This module contains the functions for converting an external operator to a Pennylane operator.
 """
 import warnings
+from itertools import product
 
 # pylint: disable=import-outside-toplevel
 import pennylane as qml
 from pennylane import numpy as np
-from pennylane.wires import Wires
 from pennylane.operation import Tensor, active_new_opmath
 from pennylane.pauli import pauli_sentence
+from pennylane.wires import Wires
 
 
 def _process_wires(wires, n_wires=None):
@@ -354,3 +355,277 @@ def import_operator(qubit_observable, format="openfermion", wires=None, tol=1e01
         return qml.dot(*_openfermion_to_pennylane(qubit_observable, wires=wires))
 
     return qml.Hamiltonian(*_openfermion_to_pennylane(qubit_observable, wires=wires))
+
+
+def _excitations(electrons, orbitals):
+    r"""Generate all possible single and double excitations from a Hartree-Fock reference state.
+
+    This function is a more performant version of `qchem.excitations`, where the order of the
+    generated excitations is consistent with PySCF.
+
+    Single and double excitations can be generated by acting with the operators
+    :math:`\hat T_1` and :math:`\hat T_2` on the Hartree-Fock reference state:
+
+    .. math::
+
+        && \hat{T}_1 = \sum_{p \in \mathrm{occ} \\ q \in \mathrm{unocc}}
+        \hat{c}_q^\dagger \hat{c}_p \\
+        && \hat{T}_2 = \sum_{p>q \in \mathrm{occ} \\ r>s \in
+        \mathrm{unocc}} \hat{c}_r^\dagger \hat{c}_s^\dagger \hat{c}_p \hat{c}_q.
+
+
+    In the equations above the indices :math:`p, q, r, s` run over the
+    occupied (occ) and unoccupied (unocc) *spin* orbitals and :math:`\hat c` and
+    :math:`\hat c^\dagger` are the electron annihilation and creation operators,
+    respectively.
+
+    Args:
+        electrons (int): number of electrons
+        orbitals (int): number of spin orbitals
+
+    Returns:
+        tuple(list, list): lists with the indices of the spin orbitals involved in the excitations
+
+    **Example**
+
+    >>> electrons = 2
+    >>> orbitals = 4
+    >>> _excitations(electrons, orbitals)
+    ([[0, 2], [0, 3], [1, 2], [1, 3]], [[0, 1, 2, 3]])
+    """
+    singles_p, singles_q = [], []
+    doubles_pq, doubles_rs = [], []
+
+    for i in range(electrons):
+        singles_p += [i]
+        doubles_pq += [[k, i] for k in range(i)]
+    for j in range(electrons, orbitals):
+        singles_q += [j]
+        doubles_rs += [[k, j] for k in range(electrons, j)]
+
+    singles = [[p] + [q] for p in singles_p for q in singles_q]
+    doubles = [pq + rs for pq in doubles_pq for rs in doubles_rs]
+
+    return singles, doubles
+
+
+def _excited_configurations(electrons, orbitals, excitation):
+    r"""Generate excited configurations from a Hartree-Fock reference state.
+
+    This function generates excited configurations in the form of integers representing a binary
+    string, e.g., :math:`|1 1 0 1 0 0 \rangle is represented by :math:`int('110100', 2) = 52`.
+
+    The excited configurations are generated from a Hartree-Fock (HF) reference state. The HF state
+    is assumed to have the form :math:`|1 1 ...1 0 ... 0 0 \rangle` where the number of :math:`1`
+    and :math:`0` elements are the number of occupied and unoccupied spin orbitals, respectively.
+    The string representation of the state is obtained by converting the occupation-number vector to
+    a string, e.g., ``111000`` to represent :math:`|1 1 1 0 0 0 \rangle.
+
+    Each excited configuration has a sign, :math:`+1` or :math:`-1`, that is obtained by reordering
+    the creation operators.
+
+    Args:
+        electrons (int): number of electrons
+        orbitals (int): number of spin orbitals
+        excitation (int): number of excited electrons
+
+    Returns:
+        tuple(array, array): arrays of excited configurations and signs obtained by reordering the
+         creation operators
+
+    **Example**
+
+    >>> electrons = 3
+    >>> orbitals = 5
+    >>> excitation = 2
+    >>> _excited_configurations(electrons, orbitals, excitation)
+    ([28, 26, 25], [1, -1, 1])
+    """
+    if excitation not in [1, 2]:
+        raise ValueError(
+            "Only single (excitation = 1) and double (excitation = 2) excitations are supported."
+        )
+
+    hf_state = qml.qchem.hf_state(electrons, orbitals)
+    singles, doubles = _excitations(electrons, orbitals)
+    states, signs = [], []
+
+    if excitation == 1:
+        for s in singles:
+            state = hf_state.copy()
+            state[s] = state[[s[1], s[0]]]  # apply single excitation
+            states += [state]
+            signs.append((-1) ** len(range(s[0], electrons - 1)))
+
+    if excitation == 2:
+        for d in doubles:
+            state = hf_state.copy()
+            state[d] = state[[d[2], d[3], d[0], d[1]]]  # apply double excitation
+            states += [state]
+            order_pq = len(range(d[0], electrons - 1))
+            order_rs = len(range(d[1], electrons - 1))
+            signs.append((-1) ** (order_pq + order_rs + 1))
+
+    states_str = ["".join([str(i) for i in state]) for state in states]
+    states_int = [int(state[::-1], 2) for state in states_str]
+
+    return states_int, signs
+
+
+def _ucisd_state(cisd_solver, tol=1e-15):
+    r"""Construct a wavefunction from PySCF's `UCISD` solver object.
+
+    The generated wavefunction is a dictionary where the keys represent a configuration, which
+    corresponds to a Slater determinant, and the values are the CI coefficients of the that Slater
+    determinant. Each dictionary key is a tuple of two integers. The binary representation of these
+    integers correspond to a specific configuration such that the first number represents the
+    configuration of the alpha electrons and the second number represents the configuration of the
+    beta electrons. For instance, the Hartree-Fock state :math:`|1 1 0 0 \rangle` will be
+    represented by the flipped binary string `0011`. This string can be splitted to `01` and `01` for
+    the alpha and beta electrons. The integer corresponding to `01` is `1`. Then the dictionary
+    representation of a state with `0.99` contribution of the Hartree-Fock state and `0.01`
+    contribution from the doubly-excited state will be `{(1, 1): 0.99, (2, 2): 0.01}`.
+
+    Args:
+        cisd_solver (PySCF UCISD Class instance): the class object representing the UCISD calculation in PySCF
+        tol (float):  the tolerance to which the wavefunction is being built -- Slater determinants
+         with coefficients smaller than this are discarded. Default is 1e-15 (all coefficients are included).
+
+    Returns:
+        dict: Dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b`
+        having binary representation corresponding to the Fock occupation vector in alpha and beta
+        spin sectors, respectively, and coeff being the CI coefficients of those configurations.
+
+    **Example**
+
+    >>> from pyscf import gto, scf, ci
+    >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g')
+    >>> myhf = scf.UHF(mol).run()
+    >>> myci = ci.UCISD(myhf).run()
+    >>> wf_cisd = _ucisd_state(myci, tol=1e-1)
+    >>> print(wf_cisd)
+    {(1, 1): -0.9942969785398778, (2, 2): 0.10664669927602159}
+    """
+    mol = cisd_solver.mol
+    cisdvec = cisd_solver.ci
+
+    norb = mol.nao
+    nelec = mol.nelectron
+    nelec_a = int((nelec + mol.spin) / 2)
+    nelec_b = int((nelec - mol.spin) / 2)
+
+    nvir_a, nvir_b = norb - nelec_a, norb - nelec_b
+
+    size_a, size_b = nelec_a * nvir_a, nelec_b * nvir_b
+    size_aa = int(nelec_a * (nelec_a - 1) / 2) * int(nvir_a * (nvir_a - 1) / 2)
+    size_bb = int(nelec_b * (nelec_b - 1) / 2) * int(nvir_b * (nvir_b - 1) / 2)
+    size_ab = nelec_a * nelec_b * nvir_a * nvir_b
+
+    sizes = [1, size_a, size_b, size_aa, size_ab, size_bb]
+    cumul = np.cumsum(sizes)
+    idxs = [0] + [slice(cumul[ii], cumul[ii + 1]) for ii in range(len(cumul) - 1)]
+    c0, c1a, c1b, c2aa, c2ab, c2bb = [cisdvec[idx] for idx in idxs]
+
+    # numbers representing the Hartree-Fock vector, e.g., bin(ref_a)[::-1] = 1111...10...0
+    ref_a = int(2**nelec_a - 1)
+    ref_b = int(2**nelec_b - 1)
+
+    dict_fcimatr = dict(zip(list(zip([ref_a], [ref_b])), [c0]))
+
+    # alpha -> alpha excitations
+    c1a_configs, c1a_signs = _excited_configurations(nelec_a, norb, 1)
+    dict_fcimatr.update(dict(zip(list(zip(c1a_configs, [ref_b] * size_a)), c1a * c1a_signs)))
+
+    # beta -> beta excitations
+    c1b_configs, c1b_signs = _excited_configurations(nelec_b, norb, 1)
+    dict_fcimatr.update(dict(zip(list(zip([ref_a] * size_b, c1b_configs)), c1b * c1b_signs)))
+
+    # alpha, alpha -> alpha, alpha excitations
+    c2aa_configs, c2aa_signs = _excited_configurations(nelec_a, norb, 2)
+    dict_fcimatr.update(dict(zip(list(zip(c2aa_configs, [ref_b] * size_aa)), c2aa * c2aa_signs)))
+
+    # alpha, beta -> alpha, beta excitations
+    rowvals, colvals = np.array(list(product(c1a_configs, c1b_configs)), dtype=int).T.numpy()
+    dict_fcimatr.update(
+        dict(zip(list(zip(rowvals, colvals)), c2ab * np.kron(c1a_signs, c1b_signs).numpy()))
+    )
+
+    # beta, beta -> beta, beta excitations
+    c2bb_configs, c2bb_signs = _excited_configurations(nelec_b, norb, 2)
+    dict_fcimatr.update(dict(zip(list(zip([ref_a] * size_bb, c2bb_configs)), c2bb * c2bb_signs)))
+
+    # filter based on tolerance cutoff
+    dict_fcimatr = {key: value for key, value in dict_fcimatr.items() if abs(value) > tol}
+
+    return dict_fcimatr
+
+
+def import_state(solver, tol=1e-15):
+    r"""Convert an external wavefunction to a PennyLane state vector.
+
+    Args:
+        solver: external wavefunction object that will be converted
+        tol (float): the tolerance for discarding Slater determinants with small coefficients
+
+    Raises:
+        ValueError: if ``method`` is not supported
+
+    Returns:
+        statevector (array): normalized state vector of length 2**len(number_of_spinorbitals)
+
+    **Example**
+
+    >>> from pyscf import gto, scf, ci
+    >>> mol = gto.M(atom=[['H', (0, 0, 0)], ['H', (0,0,0.71)]], basis='sto6g')
+    >>> myhf = scf.UHF(mol).run()
+    >>> myci = ci.UCISD(myhf).run()
+    >>> wf_cisd = qml.qchem.convert.import_state(myci, tol=1e-1)
+    >>> print(wf_cisd)
+    [ 0.        +0.j  0.        +0.j  0.        +0.j  0.1066467 +0.j
+      0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
+      0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
+     -0.99429698+0.j  0.        +0.j  0.        +0.j  0.        +0.j]
+    """
+    try:
+        import pyscf
+    except ImportError as Error:
+        raise ImportError(
+            "This feature requires pyscf. It can be installed with: pip install pyscf"
+        ) from Error
+
+    if isinstance(solver, pyscf.ci.ucisd.UCISD):
+        wf_dict = _ucisd_state(solver, tol=tol)
+    else:
+        raise ValueError("The supported option is 'ucisd' for unrestricted CISD calculations.")
+    wf = _wfdict_to_statevector(wf_dict, solver.mol.nao)
+
+    return wf
+
+
+def _wfdict_to_statevector(wf_dict, norbs):
+    r"""Convert a wavefunction in sparse dictionary format to a PennyLane statevector.
+
+    In the sparse dictionary format, the keys (int_a, int_b) are integers whose binary
+    representation shows the Fock occupation vector for alpha and beta electrons and values are the
+    CI coefficients.
+
+    Args:
+        wf_dict (dict): the sparse dictionary format of a wavefunction
+        norbs (int): number of molecular orbitals
+
+    Returns:
+        statevector (array): normalized state vector of length 2^(number_of_spinorbitals)
+    """
+    statevector = np.zeros(2 ** (2 * norbs), dtype=complex)
+
+    for (int_a, int_b), coeff in wf_dict.items():
+        bin_a = bin(int_a)[2:][::-1]
+        bin_b = bin(int_b)[2:][::-1]
+        bin_a += "0" * (norbs - len(bin_a))
+        bin_b += "0" * (norbs - len(bin_b))
+        bin_ab = "".join(i + j for i, j in zip(bin_a, bin_b))
+        statevector[int(bin_ab, 2)] += coeff
+
+    statevector = statevector / np.sqrt(np.sum(statevector**2))
+
+    return statevector