Source code for pennylane.qchem.convert

# Copyright 2018-2022 Xanadu Quantum Technologies Inc.

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     http://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
This module contains the functions for converting an external operator to a Pennylane operator.
"""
import warnings
from itertools import product

import numpy as np

# pylint: disable= import-outside-toplevel,no-member,too-many-function-args
import pennylane as qml
from pennylane.wires import Wires


def _process_wires(wires, n_wires=None):
    r"""Checks and consolidates custom user wire mapping into a consistent, direction-free, ``Wires``
    format. Used for converting between OpenFermion qubit numbering and PennyLane wire labels.

    Since OpenFermion's qubit numbering is always consecutive int, simple iterable types such as
    list, tuple, or Wires can be used to specify the 2-way `qubit index` <-> `wire label` mapping
    with indices representing qubits. Dict can also be used as a mapping, but does not provide any
    advantage over lists other than the ability to do partial mapping/permutation in the
    `qubit index` -> `wire label` direction.

    It is recommended to pass Wires/list/tuple `wires` since it's direction-free, i.e. the same
    `wires` argument can be used to convert both ways between OpenFermion and PennyLane. Only use
    dict for partial or unordered mapping.

    Args:
        wires (Wires, list, tuple, dict): User wire labels.
            For types Wires, list, or tuple, each item in the iterable represents a wire label
            corresponding to the qubit number equal to its index.
            For type dict, only int-keyed dict (for qubit-to-wire conversion) or
            consecutive-int-valued dict (for wire-to-qubit conversion) is accepted.
            If None, will be set to consecutive int based on ``n_wires``.
        n_wires (int): Number of wires used if known. If None, will be inferred from ``wires``; if
            ``wires`` is not available, will be set to 1.

    Returns:
        Wires: Cleaned wire mapping with indices corresponding to qubits and values
            corresponding to wire labels.

    **Example**

    >>> # consec int wires if no wires mapping provided, ie. identity map: 0<->0, 1<->1, 2<->2
    >>> _process_wires(None, 3)
    Wires([0, 1, 2])

    >>> # List as mapping, qubit indices with wire label values: 0<->w0, 1<->w1, 2<->w2
    >>> _process_wires(['w0','w1','w2'])
    Wires(['w0', 'w1', 'w2'])

    >>> # Wires as mapping, qubit indices with wire label values: 0<->w0, 1<->w1, 2<->w2
    >>> _process_wires(Wires(['w0', 'w1', 'w2']))
    Wires(['w0', 'w1', 'w2'])

    >>> # Dict as partial mapping, int qubits keys to wire label values: 0->w0, 1 unchanged, 2->w2
    >>> _process_wires({0:'w0',2:'w2'})
    Wires(['w0', 1, 'w2'])

    >>> # Dict as mapping, wires label keys to consec int qubit values: w2->2, w0->0, w1->1
    >>> _process_wires({'w2':2, 'w0':0, 'w1':1})
    Wires(['w0', 'w1', 'w2'])
    """

    # infer from wires, or assume 1 if wires is not of accepted types.
    if n_wires is None:
        n_wires = len(wires) if isinstance(wires, (Wires, list, tuple, dict)) else 1

    # defaults to no mapping.
    if wires is None:
        return Wires(range(n_wires))

    if isinstance(wires, (Wires, list, tuple)):
        # does not care about the tail if more wires are provided than n_wires.
        wires = Wires(wires[:n_wires])

    elif isinstance(wires, dict):
        if all(isinstance(w, int) for w in wires.keys()):
            # Assuming keys are taken from consecutive int wires. Allows for partial mapping.
            n_wires = max(wires) + 1
            labels = list(range(n_wires))  # used for completing potential partial mapping.
            for k, v in wires.items():
                if k < n_wires:
                    labels[k] = v
            wires = Wires(labels)
        elif set(range(n_wires)).issubset(set(wires.values())):
            # Assuming values are consecutive int wires (up to n_wires, ignores the rest).
            # Does NOT allow for partial mapping.
            wires = {v: k for k, v in wires.items()}  # flip for easy indexing
            wires = Wires([wires[i] for i in range(n_wires)])
        else:
            raise ValueError("Expected only int-keyed or consecutive int-valued dict for `wires`")

    else:
        raise ValueError(
            f"Expected type Wires, list, tuple, or dict for `wires`, got {type(wires)}"
        )

    if len(wires) != n_wires:
        # check length consistency when all checking and cleaning are done.
        raise ValueError(f"Length of `wires` ({len(wires)}) does not match `n_wires` ({n_wires})")

    return wires


def _openfermion_to_pennylane(qubit_operator, wires=None, tol=1.0e-16):
    r"""Convert OpenFermion ``QubitOperator`` to a 2-tuple of coefficients and
    PennyLane Pauli observables.

    Args:
        qubit_operator (QubitOperator): fermionic-to-qubit transformed operator in terms of
            Pauli matrices
        wires (Wires, list, tuple, dict): Custom wire mapping used to convert the qubit operator
            to an observable terms measurable in a PennyLane ansatz.
            For types Wires/list/tuple, each item in the iterable represents a wire label
            corresponding to the qubit number equal to its index.
            For type dict, only int-keyed dict (for qubit-to-wire conversion) is accepted.
            If None, will use identity map (e.g. 0->0, 1->1, ...).
        tol (float): whether to keep the imaginary part of the coefficients if they are smaller
            than the provided tolerance.

    Returns:
        tuple[array[float], Iterable[pennylane.operation.Operator]]: coefficients and their
        corresponding PennyLane observables in the Pauli basis

    **Example**

    >>> q_op = 0.1 * QubitOperator('X0') + 0.2 * QubitOperator('Y0 Z2')
    >>> q_op
    0.1 [X0] +
    0.2 [Y0 Z2]
    >>> _openfermion_to_pennylane(q_op, wires=['w0','w1','w2','extra_wire'])
    (tensor([0.1, 0.2], requires_grad=False), [X('w0'), Y('w0') @ Z('w2')])
    """
    n_wires = (
        1 + max(max(i for i, _ in t) if t else 1 for t in qubit_operator.terms)
        if qubit_operator.terms
        else 1
    )
    wires = _process_wires(wires, n_wires=n_wires)

    if not qubit_operator.terms:  # added since can't unpack empty zip to (coeffs, ops) below
        return np.array([0.0]), [qml.Identity(wires[0])]

    xyz2pauli = {"X": qml.X, "Y": qml.Y, "Z": qml.Z}

    def _get_op(term, wires):
        """A function to compute the PL operator associated with the term string."""
        if len(term) > 1:
            return qml.prod(*[xyz2pauli[op[1]](wires=wires[op[0]]) for op in term])

        if len(term) == 1:
            return xyz2pauli[term[0][1]](wires=wires[term[0][0]])

        return qml.Identity(wires[0])

    coeffs, ops = zip(
        *[(coef, _get_op(term, wires)) for term, coef in qubit_operator.terms.items()]
        # example term: ((0,'X'), (2,'Z'), (3,'Y'))
    )
    coeffs = np.array(coeffs)

    if (np.abs(coeffs.imag) < tol).all():
        coeffs = coeffs.real

    return coeffs, list(ops)


def _ps_to_coeff_term(ps, wire_order):
    """Convert a non-empty pauli sentence to a list of coeffs and terms."""
    ops_str = []
    pws, coeffs = zip(*ps.items())

    for pw in pws:
        if len(pw) == 0:
            ops_str.append("")
            continue
        wires, ops = zip(*pw.items())
        ops_str.append(" ".join([f"{op}{wire_order.index(wire)}" for op, wire in zip(ops, wires)]))

    return coeffs, ops_str


# pylint:disable=too-many-branches
def _pennylane_to_openfermion(coeffs, ops, wires=None, tol=1.0e-16):
    r"""Convert a 2-tuple of complex coefficients and PennyLane operations to
    OpenFermion ``QubitOperator``.

    Args:
        coeffs (array[complex]):
            coefficients for each observable, same length as ops
        ops (Iterable[pennylane.operation.Operations]): list of PennyLane operations that
            have a valid PauliSentence representation.
        wires (Wires, list, tuple, dict): Custom wire mapping used to convert to qubit operator
            from an observable terms measurable in a PennyLane ansatz.
            For types Wires/list/tuple, each item in the iterable represents a wire label
            corresponding to the qubit number equal to its index.
            For type dict, only consecutive-int-valued dict (for wire-to-qubit conversion) is
            accepted. If None, will map sorted wires from all `ops` to consecutive int.
        tol (float): whether to keep the imaginary part of the coefficients if they are smaller
            than the provided tolerance.

    Returns:
        QubitOperator: an instance of OpenFermion's ``QubitOperator``.

    **Example**

    >>> coeffs = np.array([0.1, 0.2, 0.3, 0.4])
    >>> ops = [
    ...     qml.prod(qml.X('w0')),
    ...     qml.prod(qml.Y('w0'), qml.Z('w2')),
    ...     qml.sum(qml.Z('w0'), qml.s_prod(-0.5, qml.X('w0'))),
    ...     qml.prod(qml.X('w0'), qml.Z('w1')),
    ... ]
    >>> _pennylane_to_openfermion(coeffs, ops, wires=Wires(['w0', 'w1', 'w2']))
    (-0.05+0j) [X0] +
    (0.4+0j) [X0 Z1] +
    (0.2+0j) [Y0 Z2] +
    (0.3+0j) [Z0]
    """
    try:
        import openfermion
    except ImportError as Error:
        raise ImportError(
            "This feature requires openfermion. "
            "It can be installed with: pip install openfermion"
        ) from Error

    all_wires = Wires.all_wires([op.wires for op in ops], sort=True)

    if wires is not None:
        qubit_indexed_wires = _process_wires(
            wires,
        )
        if not set(all_wires).issubset(set(qubit_indexed_wires)):
            raise ValueError("Supplied `wires` does not cover all wires defined in `ops`.")
    else:
        qubit_indexed_wires = all_wires

    coeffs = np.array(coeffs)
    if (np.abs(coeffs.imag) < tol).all():
        coeffs = coeffs.real

    q_op = openfermion.QubitOperator()
    for coeff, op in zip(coeffs, ops):

        if (ps := op.pauli_rep) is None:
            raise ValueError(
                f"Expected a Pennylane operator with a valid Pauli word representation, but got {op}."
            )

        if len(ps) > 0:
            sub_coeffs, op_strs = _ps_to_coeff_term(ps, wire_order=qubit_indexed_wires)
            for c, op_str in zip(sub_coeffs, op_strs):
                # This is how one makes QubitOperator in OpenFermion
                q_op += complex(coeff * c) * openfermion.QubitOperator(op_str)

    return q_op


def _openfermion_pennylane_equivalent(
    openfermion_qubit_operator, pennylane_qubit_operator, wires=None
):
    r"""Check equivalence between OpenFermion :class:`~.QubitOperator` and Pennylane VQE
    ``Hamiltonian`` (Tensor product of Pauli matrices).

    Equality is based on OpenFermion :class:`~.QubitOperator`'s equality.

    Args:
        openfermion_qubit_operator (QubitOperator): OpenFermion qubit operator represented as
            a Pauli summation
        pennylane_qubit_operator (pennylane.Hamiltonian): PennyLane
            Hamiltonian object
        wires (Wires, list, tuple, dict): Custom wire mapping used to convert to qubit operator
            from an observable terms measurable in a PennyLane ansatz.
            For types Wires/list/tuple, each item in the iterable represents a wire label
            corresponding to the qubit number equal to its index.
            For type dict, only int-keyed dict (for qubit-to-wire conversion) is accepted.
            If None, will map sorted wires from all Pennylane terms to consecutive int.

    Returns:
        (bool): True if equivalent
    """
    coeffs, ops = pennylane_qubit_operator.terms()
    return openfermion_qubit_operator == _pennylane_to_openfermion(
        np.array(coeffs), ops, wires=wires
    )


def import_operator(qubit_observable, format="openfermion", wires=None, tol=1e010):
    r"""Convert an external operator to a PennyLane operator.

    We currently support `OpenFermion <https://quantumai.google/openfermion>`__ operators: the function accepts most types of
    OpenFermion qubit operators, such as those corresponding to Pauli words and sums of Pauli words.

    Args:
        qubit_observable: external qubit operator that will be converted
        format (str): the format of the operator object to convert from
        wires (.Wires, list, tuple, dict): Custom wire mapping used to convert the external qubit
            operator to a PennyLane operator.
            For types ``Wires``/list/tuple, each item in the iterable represents a wire label
            for the corresponding qubit index.
            For type dict, only int-keyed dictionaries (for qubit-to-wire conversion) are accepted.
            If ``None``, the identity map (e.g., ``0->0, 1->1, ...``) will be used.
        tol (float): Tolerance in `machine epsilon <https://numpy.org/doc/stable/reference/generated/numpy.real_if_close.html>`_
            for the imaginary part of the coefficients in ``qubit_observable``.
            Coefficients with imaginary part less than :math:`(2.22 \cdot 10^{-16}) \cdot \text{tol}` are considered to be real.

    Returns:
        (.Operator): PennyLane operator representing any operator expressed as linear combinations of
        Pauli words, e.g.,
        :math:`\sum_{k=0}^{N-1} c_k O_k`

    **Example**

    >>> h_pl = import_operator(h_of, format='openfermion')
    >>> print(h_pl)
    (-0.0548 * X(0 @ X(1) @ Y(2) @ Y(3))) + (0.14297 * Z(0 @ Z(1)))
    """
    if format not in ["openfermion"]:
        raise TypeError(f"Converter does not exist for {format} format.")

    coeffs = np.array([np.real_if_close(coef, tol=tol) for coef in qubit_observable.terms.values()])

    if any(np.iscomplex(coeffs)):
        warnings.warn(
            f"The coefficients entering the QubitOperator must be real;"
            f" got complex coefficients in the operator"
            f" {list(coeffs[np.iscomplex(coeffs)])}"
        )

    return qml.dot(*_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(list, list): lists 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 import_state(solver, tol=1e-15):
    r"""Convert an external wavefunction to a state vector.

    The sources of wavefunctions that are currently accepted are listed below.

        * The PySCF library (the restricted and unrestricted CISD/CCSD
          methods are supported). The `solver` argument is then the associated PySCF CISD/CCSD Solver object.
        * The library Dice implementing the SHCI method. The `solver` argument is then the tuple(list[str], array[float]) of Slater determinants and their coefficients.
        * The library Block2 implementing the DMRG method. The `solver` argument is then the tuple(list[int], array[float]) of Slater determinants and their coefficients.

    Args:
        solver: external wavefunction object
        tol (float): the tolerance for discarding Slater determinants based on their coefficients

    Raises:
        ValueError: if external object type is not supported

    Returns:
        array: normalized state vector of length :math:`2^M`, where :math:`M` is the number of spin orbitals

    **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.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]
    """

    method = str(solver.__str__)

    if "CISD" in method and "UCISD" not in method:
        wf_dict = _rcisd_state(solver, tol=tol)
    elif "UCISD" in method:
        wf_dict = _ucisd_state(solver, tol=tol)
    elif "CCSD" in method and "UCCSD" not in method:
        wf_dict = _rccsd_state(solver, tol=tol)
    elif "UCCSD" in method:
        wf_dict = _uccsd_state(solver, tol=tol)
    elif "tuple" in method and len(solver) == 2:
        if isinstance(solver[0][0], str):
            wf_dict = _shci_state(solver, tol=tol)
        elif isinstance(solver[0][0][0], int):
            wf_dict = _dmrg_state(solver, tol=tol)
        else:
            raise ValueError(
                "For tuple input, the supported objects are"
                " tuple(list[str], array[float]) for SHCI calculations with Dice library and "
                "tuple(list[int], array[float]) for DMRG calculations with the Block2 library."
            )
    else:
        raise ValueError(
            "The supported objects are RCISD, UCISD, RCCSD, and UCCSD for restricted and"
            " unrestricted configuration interaction and coupled cluster calculations, and"
            " tuple(list[str], array[float]) for SHCI calculations with Dice library and "
            "tuple(list[int], array[float]) for DMRG calculations with the Block2 library."
        )
    if "tuple" in method:
        wf = _wfdict_to_statevector(wf_dict, len(solver[0][0]))
    else:
        wf = _wfdict_to_statevector(wf_dict, solver.mol.nao)

    return wf


def _wfdict_to_statevector(fcimatr_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:
        fcimatr_dict (dict[tuple(int,int),float]): the sparse dictionary format of a wavefunction
        norbs (int): number of molecular orbitals

    Returns:
        array: normalized state vector of length :math:`2^M`, where :math:`M` is the number of spin orbitals
    """
    statevector = np.zeros(2 ** (2 * norbs), dtype=complex)

    for (int_a, int_b), coeff in fcimatr_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


def _rcisd_state(cisd_solver, tol=1e-15):
    r"""Construct a wavefunction from PySCF's ``RCISD`` 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 Slater
    determinant. Each dictionary key is a tuple of two integers. The binary representation of these
    integers correspond to a specific configuration: 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`` which is split to ``01`` and ``01`` for
    the alpha and beta electrons. The integer corresponding to ``01`` is ``1`` and the dictionary
    representation of the Hartree-Fock state will be ``{(1, 1): 1.0}``. The dictionary
    representation of a state with ``0.99`` contribution from the Hartree-Fock state and ``0.01``
    contribution from the doubly-excited state, i.e., :math:`|0 0 1 1 \rangle`, will be
    ``{(1, 1): 0.99, (2, 2): 0.01}``.

    Args:
        cisd_solver (PySCF CISD Class instance): the class object representing the CISD calculation in PySCF
        tol (float): the tolerance for discarding Slater determinants based on their coefficients

    Returns:
        fcimatr_dict (dict[tuple(int,int),float]): 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', symmetry='d2h')
    >>> myhf = scf.RHF(mol).run()
    >>> myci = ci.CISD(myhf).run()
    >>> wf_cisd = _rcisd_state(myci, tol=1e-1)
    >>> print(wf_cisd)
    {(1, 1): -0.9942969785398775, (2, 2): 0.10664669927602162}
    """
    mol = cisd_solver.mol
    cisdvec = cisd_solver.ci

    norb = mol.nao
    nelec = mol.nelectron
    nocc, nvir = nelec // 2, norb - nelec // 2

    c0, c1, c2 = (
        cisdvec[0],
        cisdvec[1 : nocc * nvir + 1],
        cisdvec[nocc * nvir + 1 :].reshape(nocc, nocc, nvir, nvir),
    )
    c2ab = c2.transpose(0, 2, 1, 3).reshape(nocc * nvir, -1)

    # numbers representing the Hartree-Fock vector, e.g., bin(ref_a)[::-1] = 1111...10...0
    ref_a = int(2**nocc - 1)
    ref_b = ref_a

    fcimatr_dict = dict(zip(list(zip([ref_a], [ref_b])), [c0]))

    # alpha -> alpha excitations
    c1a_configs, c1a_signs = _excited_configurations(nocc, norb, 1)
    fcimatr_dict.update(
        dict(zip(list(zip(c1a_configs, [ref_b] * len(c1a_configs))), c1 * c1a_signs))
    )
    # beta -> beta excitations
    fcimatr_dict.update(
        dict(zip(list(zip([ref_a] * len(c1a_configs), c1a_configs)), c1 * c1a_signs))
    )

    # check if double excitations within one spin sector (aa->aa and bb->bb) are possible
    if nocc > 1 and nvir > 1:
        # get rid of excitations from identical orbitals, double-count the allowed ones
        c2_tr = c2 - c2.transpose(1, 0, 2, 3)
        # select only unique excitations, via lower triangle of matrix (already double-counted)
        ooidx, vvidx = np.tril_indices(nocc, -1), np.tril_indices(nvir, -1)
        c2aa = c2_tr[ooidx][:, vvidx[0], vvidx[1]].ravel()

        # alpha, alpha -> alpha, alpha excitations
        c2aa_configs, c2aa_signs = _excited_configurations(nocc, norb, 2)
        fcimatr_dict.update(
            dict(zip(list(zip(c2aa_configs, [ref_b] * len(c2aa_configs))), c2aa * c2aa_signs))
        )
        # beta, beta -> beta, beta excitations
        fcimatr_dict.update(
            dict(zip(list(zip([ref_a] * len(c2aa_configs), c2aa_configs)), c2aa * c2aa_signs))
        )

    # alpha, beta -> alpha, beta excitations
    # generate all possible pairwise combinations of _single_ excitations of alpha and beta sectors
    fcimatr_dict.update(
        dict(
            zip(
                list(product(c1a_configs, c1a_configs)),
                np.einsum("i,j,ij->ij", c1a_signs, c1a_signs, c2ab, optimize=True).ravel(),
            )
        )
    )

    # filter based on tolerance cutoff
    fcimatr_dict = {key: value for key, value in fcimatr_dict.items() if abs(value) > tol}

    # convert sign parity from chemist to physicist convention (interleaving spin operators
    # rather than commuting all spin-up operators to the left)
    fcimatr_dict = _sign_chem_to_phys(fcimatr_dict, norb)

    return fcimatr_dict


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 Slater
    determinant. Each dictionary key is a tuple of two integers. The binary representation of these
    integers correspond to a specific configuration: 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`` which is split to ``01`` and ``01`` for
    the alpha and beta electrons. The integer corresponding to ``01`` is ``1`` and the dictionary
    representation of the Hartree-Fock state will be ``{(1, 1): 1.0}``. The dictionary
    representation of a state with `0.99` contribution from the Hartree-Fock state and ``0.01``
    contribution from the doubly-excited state, i.e., :math:`|0 0 1 1 \rangle`, 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 for discarding Slater determinants based on their coefficients

    Returns:
        fcimatr_dict (dict[tuple(int,int),float]): 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', symmetry='d2h')
    >>> 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 = size_a * size_b

    cumul = np.cumsum([0, 1, size_a, size_b, size_ab, size_aa, size_bb])
    c0, c1a, c1b, c2ab, c2aa, c2bb = [
        cisdvec[cumul[idx] : cumul[idx + 1]] for idx in range(len(cumul) - 1)
    ]
    c2ab = (
        c2ab.reshape(nelec_a, nelec_b, nvir_a, nvir_b)
        .transpose(0, 2, 1, 3)
        .reshape(nelec_a * nvir_a, -1)
    )

    # 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)

    fcimatr_dict = dict(zip(list(zip([ref_a], [ref_b])), c0))

    # alpha -> alpha excitations
    c1a_configs, c1a_signs = _excited_configurations(nelec_a, norb, 1)
    fcimatr_dict.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)
    fcimatr_dict.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)
    fcimatr_dict.update(dict(zip(list(zip(c2aa_configs, [ref_b] * size_aa)), c2aa * c2aa_signs)))

    # alpha, beta -> alpha, beta excitations
    fcimatr_dict.update(
        dict(
            zip(
                list(product(c1a_configs, c1b_configs)),
                np.einsum("i,j,ij->ij", c1a_signs, c1b_signs, c2ab, optimize=True).ravel(),
            )
        )
    )

    # beta, beta -> beta, beta excitations
    c2bb_configs, c2bb_signs = _excited_configurations(nelec_b, norb, 2)
    fcimatr_dict.update(dict(zip(list(zip([ref_a] * size_bb, c2bb_configs)), c2bb * c2bb_signs)))

    # filter based on tolerance cutoff
    fcimatr_dict = {key: value for key, value in fcimatr_dict.items() if abs(value) > tol}

    # convert sign parity from chemist to physicist convention (interleaving spin operators
    # rather than commuting all spin-up operators to the left)
    fcimatr_dict = _sign_chem_to_phys(fcimatr_dict, norb)

    return fcimatr_dict


def _rccsd_state(ccsd_solver, tol=1e-15):
    r"""Construct a wavefunction from PySCF's ``RCCSD`` 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 Slater
    determinant. Each dictionary key is a tuple of two integers. The binary representation of these
    integers correspond to a specific configuration: 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` which is split to ``01`` and ``01`` for
    the alpha and beta electrons. The integer corresponding to ``01`` is ``1`` and the dictionary
    representation of the Hartree-Fock state will be ``{(1, 1): 1.0}``. The dictionary
    representation of a state with ``0.99`` contribution from the Hartree-Fock state and ``0.01``
    contribution from the doubly-excited state, i.e., :math:`|0 0 1 1 \rangle`, will be
    ``{(1, 1): 0.99, (2, 2): 0.01}``.

    In the current version, the exponential ansatz :math:`\exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}}`
    is expanded to second order, with only single and double excitation terms collected and kept.
    In the future this will be amended to also collect terms from higher order. The expansion gives

    .. math::
        \exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}} = \left[ 1 + \hat{T}_1 +
        \left( \hat{T}_2 + 0.5 * \hat{T}_1^2 \right) \right] \ket{\text{HF}}

    The coefficients in this expansion are the CI coefficients used to build the wavefunction
    representation.

    Args:
        ccsd_solver (PySCF RCCSD Class instance): the class object representing the RCCSD calculation in PySCF
        tol (float): the tolerance for discarding Slater determinants with small coefficients

    Returns:
        fcimatr_dict (dict[tuple(int,int),float]): dictionary of the form ``{(int_a, int_b) :coeff}``, with integers ``int_a, int_b``
        having binary represention 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, cc
    >>> mol = gto.M(atom=[['Li', (0, 0, 0)], ['Li', (0,0,0.71)]], basis='sto6g', symmetry="d2h")
    >>> myhf = scf.RHF(mol).run()
    >>> mycc = cc.CCSD(myhf).run()
    >>> wf_ccsd = _rccsd_state(mycc, tol=1e-1)
    >>> print(wf_ccsd)
    {(7, 7): -0.8886969878256522, (11, 11): 0.30584590248164206,
     (19, 19): 0.30584590248164145, (35, 35): 0.14507552651170982}
    """

    mol = ccsd_solver.mol

    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

    # build the full, unrestricted representation of the coupled cluster amplitudes
    t1a = ccsd_solver.t1
    t1b = t1a
    t2aa = ccsd_solver.t2 - ccsd_solver.t2.transpose(1, 0, 2, 3)
    t2ab = ccsd_solver.t2.transpose(0, 2, 1, 3)
    t2bb = t2aa

    # add in the disconnected part ( + 0.5 T_1^2) of double excitations
    t2aa = t2aa - 0.5 * np.kron(t1a, t1a).reshape(nelec_a, nvir_a, nelec_a, nvir_a).transpose(
        0, 2, 1, 3
    )
    t2bb = t2bb - 0.5 * np.kron(t1b, t1b).reshape(nelec_b, nvir_b, nelec_b, nvir_b).transpose(
        0, 2, 1, 3
    )
    # align the entries with how the excitations are ordered when generated by _excited_configurations()
    t2ab = t2ab - 0.5 * np.kron(t1a, t1b).reshape(nelec_a, nvir_a, nelec_b, nvir_b)

    # 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)

    fcimatr_dict = dict(zip(list(zip([ref_a], [ref_b])), [1.0]))

    # alpha -> alpha excitations
    t1a_configs, t1a_signs = _excited_configurations(nelec_a, norb, 1)
    fcimatr_dict.update(
        dict(zip(list(zip(t1a_configs, [ref_b] * len(t1a_configs))), t1a.ravel() * t1a_signs))
    )

    # beta -> beta excitations
    t1b_configs, t1b_signs = _excited_configurations(nelec_b, norb, 1)
    fcimatr_dict.update(
        dict(zip(list(zip([ref_a] * len(t1b_configs), t1b_configs)), t1b.ravel() * t1b_signs))
    )

    # alpha, alpha -> alpha, alpha excitations
    if nelec_a > 1 and nvir_a > 1:
        t2aa_configs, t2aa_signs = _excited_configurations(nelec_a, norb, 2)
        # select only unique excitations, via lower triangle of matrix
        ooidx = np.tril_indices(nelec_a, -1)
        vvidx = np.tril_indices(nvir_a, -1)
        t2aa = t2aa[ooidx][:, vvidx[0], vvidx[1]]
        fcimatr_dict.update(
            dict(
                zip(list(zip(t2aa_configs, [ref_b] * len(t2aa_configs))), t2aa.ravel() * t2aa_signs)
            )
        )

    if nelec_b > 1 and nvir_b > 1:
        t2bb_configs, t2bb_signs = _excited_configurations(nelec_b, norb, 2)
        # select only unique excitations, via lower triangle of matrix
        ooidx = np.tril_indices(nelec_b, -1)
        vvidx = np.tril_indices(nvir_b, -1)
        t2bb = t2bb[ooidx][:, vvidx[0], vvidx[1]]
        fcimatr_dict.update(
            dict(
                zip(list(zip([ref_a] * len(t2bb_configs), t2bb_configs)), t2bb.ravel() * t2bb_signs)
            )
        )

    # alpha, beta -> alpha, beta excitations
    fcimatr_dict.update(
        dict(
            zip(
                list(product(t1a_configs, t1b_configs)),
                np.einsum(
                    "i,j,ij->ij",
                    t1a_signs,
                    t1b_signs,
                    t2ab.reshape(nelec_a * nvir_a, -1),
                    optimize=True,
                ).ravel(),
            )
        )
    )

    # renormalize, to get the HF coefficient (CC wavefunction not normalized)
    norm = np.sqrt(np.sum(np.array(list(fcimatr_dict.values())) ** 2))
    fcimatr_dict = {key: value / norm for (key, value) in fcimatr_dict.items()}

    # filter based on tolerance cutoff
    fcimatr_dict = {key: value for key, value in fcimatr_dict.items() if abs(value) > tol}

    # convert sign parity from chemist to physicist convention (interleaving spin operators
    # rather than commuting all spin-up operators to the left)
    fcimatr_dict = _sign_chem_to_phys(fcimatr_dict, norb)

    return fcimatr_dict


def _uccsd_state(ccsd_solver, tol=1e-15):
    r"""Construct a wavefunction from PySCF's ``UCCSD`` 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 Slater
    determinant. Each dictionary key is a tuple of two integers. The binary representation of these
    integers correspond to a specific configuration: 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`` which is split to ``01`` and ``01`` for
    the alpha and beta electrons. The integer corresponding to ``01`` is ``1`` and the dictionary
    representation of the Hartree-Fock state will be ``{(1, 1): 1.0}``. The dictionary
    representation of a state with ``0.99`` contribution from the Hartree-Fock state and ``0.01``
    contribution from the doubly-excited state, i.e., :math:`|0 0 1 1 \rangle`, will be
    ``{(1, 1): 0.99, (2, 2): 0.01}``.

    In the current version, the exponential ansatz :math:`\exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}}`
    is expanded to second order, with only single and double excitation terms collected and kept.
    In the future this will be amended to also collect terms from higher order. The expansion gives

    .. math::
        \exp(\hat{T}_1 + \hat{T}_2) \ket{\text{HF}} = \left[ 1 + \hat{T}_1 +
                                    \left( \hat{T}_2 + 0.5 * \hat{T}_1^2 \right) \right] \ket{\text{HF}}

    The coefficients in this expansion are the CI coefficients used to build the wavefunction
    representation.

    Args:
        ccsd_solver (PySCF UCCSD Class instance): the class object representing the UCCSD calculation in PySCF
        tol (float): the tolerance for discarding Slater determinants with small coefficients

    Returns:
        fcimatr_dict (dict[tuple(int,int),float]): dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b`
        having binary represention 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, cc
    >>> mol = gto.M(atom=[['Li', (0, 0, 0)], ['Li', (0,0,0.71)]], basis='sto6g', symmetry="d2h")
    >>> myhf = scf.UHF(mol).run()
    >>> mycc = cc.UCCSD(myhf).run()
    >>> wf_ccsd = _uccsd_state(mycc, tol=1e-1)
    >>> print(wf_ccsd)
    {(7, 7): -0.8886970081919591, (11, 11): 0.3058459002168582,
     (19, 19): 0.30584590021685887, (35, 35): 0.14507552387854625}
    """

    mol = ccsd_solver.mol

    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

    t1a, t1b = ccsd_solver.t1
    t2aa, t2ab, t2bb = ccsd_solver.t2
    # add in the disconnected part ( + 0.5 T_1^2) of double excitations
    t2aa = t2aa - 0.5 * np.kron(t1a, t1a).reshape(nelec_a, nvir_a, nelec_a, nvir_a).transpose(
        0, 2, 1, 3
    )
    t2bb = t2bb - 0.5 * np.kron(t1b, t1b).reshape(nelec_b, nvir_b, nelec_b, nvir_b).transpose(
        0, 2, 1, 3
    )
    # align the entries with how the excitations are ordered when generated by _excited_configurations()
    t2ab = t2ab.transpose(0, 2, 1, 3) - 0.5 * np.kron(t1a, t1b).reshape(
        nelec_a, nvir_a, nelec_b, nvir_b
    )

    # 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)

    fcimatr_dict = dict(zip(list(zip([ref_a], [ref_b])), [1.0]))

    # alpha -> alpha excitations
    t1a_configs, t1a_signs = _excited_configurations(nelec_a, norb, 1)
    fcimatr_dict.update(
        dict(zip(list(zip(t1a_configs, [ref_b] * len(t1a_configs))), t1a.ravel() * t1a_signs))
    )

    # beta -> beta excitations
    t1b_configs, t1b_signs = _excited_configurations(nelec_b, norb, 1)
    fcimatr_dict.update(
        dict(zip(list(zip([ref_a] * len(t1b_configs), t1b_configs)), t1b.ravel() * t1b_signs))
    )

    # alpha, alpha -> alpha, alpha excitations
    if nelec_a > 1 and nvir_a > 1:
        t2aa_configs, t2aa_signs = _excited_configurations(nelec_a, norb, 2)
        # select only unique excitations, via lower triangle of matrix
        ooidx = np.tril_indices(nelec_a, -1)
        vvidx = np.tril_indices(nvir_a, -1)
        t2aa = t2aa[ooidx][:, vvidx[0], vvidx[1]]
        fcimatr_dict.update(
            dict(
                zip(list(zip(t2aa_configs, [ref_b] * len(t2aa_configs))), t2aa.ravel() * t2aa_signs)
            )
        )

    # beta, beta -> beta, beta excitations
    if nelec_b > 1 and nvir_b > 1:
        t2bb_configs, t2bb_signs = _excited_configurations(nelec_b, norb, 2)
        # select only unique excitations, via lower triangle of matrix
        ooidx = np.tril_indices(nelec_b, -1)
        vvidx = np.tril_indices(nvir_b, -1)
        t2bb = t2bb[ooidx][:, vvidx[0], vvidx[1]]
        fcimatr_dict.update(
            dict(
                zip(list(zip([ref_a] * len(t2bb_configs), t2bb_configs)), t2bb.ravel() * t2bb_signs)
            )
        )

    # alpha, beta -> alpha, beta excitations
    fcimatr_dict.update(
        dict(
            zip(
                list(product(t1a_configs, t1b_configs)),
                np.einsum(
                    "i,j,ij->ij",
                    t1a_signs,
                    t1b_signs,
                    t2ab.reshape(nelec_a * nvir_a, -1),
                    optimize=True,
                ).ravel(),
            )
        )
    )

    # renormalize, to get the HF coefficient (CC wavefunction not normalized)
    norm = np.sqrt(np.sum(np.array(list(fcimatr_dict.values())) ** 2))
    fcimatr_dict = {key: value / norm for (key, value) in fcimatr_dict.items()}

    # filter based on tolerance cutoff
    fcimatr_dict = {key: value for key, value in fcimatr_dict.items() if abs(value) > tol}

    # convert sign parity from chemist to physicist convention (interleaving spin operators
    # rather than commuting all spin-up operators to the left)
    fcimatr_dict = _sign_chem_to_phys(fcimatr_dict, norb)

    return fcimatr_dict


def _dmrg_state(wavefunction, tol=1e-15):
    r"""Construct a wavefunction from the DMRG wavefunction obtained from the Block2 library.

    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 Slater
    determinant. Each dictionary key is a tuple of two integers. The binary representation of these
    integers correspond to a specific configuration: 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`` which is split to ``01`` and ``01`` for
    the alpha and beta electrons. The integer corresponding to ``01`` is ``1`` and the dictionary
    representation of the Hartree-Fock state will be ``{(1, 1): 1.0}``. The dictionary
    representation of a state with ``0.99`` contribution from the Hartree-Fock state and ``0.01``
    contribution from the doubly-excited state, i.e., :math:`|0 0 1 1 \rangle`, will be
    ``{(1, 1): 0.99, (2, 2): 0.01}``.

    The determinants and coefficients should be supplied externally. They should be calculated by
    using Block2 DMRGDriver's `get_csf_coefficients()` method.

    Args:
        wavefunction tuple(list[int], array[float]): determinants and coefficients in physicist notation
        tol (float): the tolerance for discarding Slater determinants with small coefficients

    Returns:
        fcimatr_dict (dict[tuple(int,int),float]): 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**

    >>> import numpy as np
    >>> wavefunction = ( [[0, 3], [3, 0]], np.array([-0.10660077,  0.9943019 ]))
    >>> wf_dmrg = _dmrg_state(wavefunction, tol=1e-1)
    >>> print(wf_dmrg)
    {(2, 2): -0.10660077, (1, 1): 0.9943019}
    """
    dets, coeffs = wavefunction

    row, col = [], []

    for det in dets:
        stra, strb = _sitevec_to_fock(det, format="dmrg")
        row.append(stra)
        col.append(strb)

        # slater determinant sign convention of block2 (physicist,
        # interleave spin-up/down operators) is consistent with pennylane

    ## create the FCI matrix as a dict
    fcimatr_dict = dict(zip(list(zip(row, col)), coeffs))

    # filter based on tolerance cutoff
    fcimatr_dict = {key: value for key, value in fcimatr_dict.items() if abs(value) > tol}

    return fcimatr_dict


def _shci_state(wavefunction, tol=1e-15):
    r"""Construct a wavefunction from the SHCI wavefunction obtained from the Dice library.

    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 Slater
    determinant. Each dictionary key is a tuple of two integers. The binary representation of these
    integers correspond to a specific configuration: 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`` which is split to ``01`` and ``01`` for
    the alpha and beta electrons. The integer corresponding to ``01`` is ``1`` and the dictionary
    representation of the Hartree-Fock state will be ``{(1, 1): 1.0}``. The dictionary
    representation of a state with ``0.99`` contribution from the Hartree-Fock state and ``0.01``
    contribution from the doubly-excited state, i.e., :math:`|0 0 1 1 \rangle`, will be
    ``{(1, 1): 0.99, (2, 2): 0.01}``.

    The determinants and coefficients should be supplied externally. They are typically stored under
    SHCI.outputfile.

    Args:
        wavefunction tuple(list[str], array[float]): determinants and coefficients in chemist notation
        tol (float): the tolerance for discarding Slater determinants with small coefficients
    Returns:
        fcimatr_dict (dict[tuple(int,int),float]): 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**

    >>> import numpy as np
    >>> wavefunction = ( ['20', '02'], np.array([-0.9943019036, 0.1066007711]))
    >>> wf_shci = _shci_state(wavefunction, tol=1e-1)
    >>> print(wf_shci)
    {(1, 1): -0.9943019036, (2, 2): 0.1066007711}
    """

    dets, coeffs = wavefunction
    norb = len(dets[0])

    xa = []
    xb = []

    for det in dets:
        bin_a, bin_b = _sitevec_to_fock(list(det), "shci")

        xa.append(bin_a)
        xb.append(bin_b)

    ## create the FCI matrix as a dict
    fcimatr_dict = dict(zip(list(zip(xa, xb)), coeffs))

    # filter based on tolerance cutoff
    fcimatr_dict = {key: value for key, value in fcimatr_dict.items() if abs(value) > tol}

    # convert sign parity from chemist to physicist convention (interleaving spin operators
    # rather than commuting all spin-up operators to the left)
    fcimatr_dict = _sign_chem_to_phys(fcimatr_dict, norb)

    return fcimatr_dict


def _sitevec_to_fock(det, format):
    r"""Convert a Slater determinant from site vector to occupation number vector representation.

    Args:
        det (list(int) or list(str)): determinant in site vector representation
        format (str): the format of the determinant

    Returns:
        tuple: tuple of integers representing binaries that correspond to occupation vectors in
            alpha and beta spin sectors

    **Example**

    >>> det = [1, 2, 1, 0, 0, 2]
    >>> _sitevec_to_fock(det, format = 'dmrg')
    (5, 34)

    >>> det = ["a", "b", "a", "0", "0", "b"]
    >>> _sitevec_to_fock(det, format = 'shci')
    (5, 34)
    """

    if format == "dmrg":
        format_map = {0: "00", 1: "10", 2: "01", 3: "11"}
    elif format == "shci":
        format_map = {"0": "00", "a": "10", "b": "01", "2": "11"}

    strab = [format_map[key] for key in det]

    stra = "".join(i[0] for i in strab)
    strb = "".join(i[1] for i in strab)

    inta = int(stra[::-1], 2)
    intb = int(strb[::-1], 2)

    return inta, intb


def _sign_chem_to_phys(fcimatr_dict, norb):
    r"""Convert the dictionary-form wavefunction from chemist sign convention
    for ordering the creation operators by spin (i.e. all spin-up operators
    on the left) to the physicist convention native to PennyLane, which
    storing spin operators as interleaved for the same spatial orbital index.

    Note that convention change in the opposite direction -- starting from physicist
    and going to chemist -- can be accomplished with the same function
    (the sign transformation is reversible).

    Args:
        fcimatr_dict (dict[tuple(int, int), float]): dictionary of the form `{(int_a, int_b) :coeff}`, with integers `int_a, int_b`
        having binary represention corresponding to the Fock occupation vector in alpha and beta
        spin sectors, respectively, and coeff being the CI coefficients of those configurations
        norb (int): total number of spatial orbitals of the underlying system

    Returns:
        signed_dict (dict): the same dictionary-type wavefunction with appropriate signs converted

    **Example**

    >>> fcimatr_dict = {(3, 1): 0.96, (6, 1): 0.1, \
                        (3, 4): 0.1, (6, 4): 0.14, (5, 2): 0.19}
    >>> _sign_chem_to_phys(fcimatr_dict, 3)
    {(3, 1): -0.96, (6, 1): 0.1, (3, 4): 0.1, (6, 4): 0.14, (5, 2): -0.19}
    """

    signed_dict = {}
    for key, elem in fcimatr_dict.items():
        lsta, lstb = bin(key[0])[2:][::-1], bin(key[1])[2:][::-1]
        # highest energy state is on the right -- pad to the right
        lsta = np.array([int(elem) for elem in lsta] + [0] * (norb - len(lsta)))
        lstb = np.array([int(elem) for elem in lstb] + [0] * (norb - len(lstb)))
        which_occ = np.where(lsta == 1)[0]
        parity = (-1) ** np.sum([np.sum(lstb[: int(ind)]) for ind in which_occ])
        signed_dict[key] = parity * elem
    return signed_dict

_modules/pennylane/qchem/convert

 

Download Python script

 

Download Notebook

 

View on GitHub