code: 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:
          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', 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),
      )
  
      # 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
  
      dict_fcimatr = dict(zip(list(zip([ref_a], [ref_b])), [c0]))
  
      # alpha -> alpha excitations
      c1a_configs, c1a_signs = _excited_configurations(nocc, norb, 1)
      dict_fcimatr.update(
          dict(
              zip(
                  list(zip(c1a_configs, [ref_b] * len(c1a_configs))),
                  c1 * c1a_signs,
              )
          )
      )
      # beta -> beta excitations
      dict_fcimatr.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
          )
          dict_fcimatr.update(
              dict(
                  zip(
                      list(
                          zip(c2aa_configs, [ref_b] * len(c2aa_configs))
                      ),
                      c2aa * c2aa_signs,
                  )
              )
          )
          # beta, beta -> beta, beta excitations
          dict_fcimatr.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
      rowvals, colvals = np.array(
          list(product(c1a_configs, c1a_configs)), dtype=int
      ).T
      c2ab = (c2.transpose(0, 2, 1, 3).reshape(nocc * nvir, -1)).ravel()
      dict_fcimatr.update(
          dict(
              zip(
                  list(zip(rowvals, colvals)),
                  c2ab * np.kron(c1a_signs, c1a_signs),
              )
          )
      )
  
      # filter based on tolerance cutoff
      dict_fcimatr = {
          key: value
          for key, value in dict_fcimatr.items()
          if abs(value) > tol
      }
  
      return dict_fcimatr
  
imports: import warnings
  from itertools import product
  import numpy
  import pennylane
  from pennylane.operation import Tensor, active_new_opmath
  from pennylane.pauli import pauli_sentence
  from pennylane.wires import Wires
label_desc: Use image node_id 12 for calling a global function with example usage: _rcisd_state(cisd_solver, tol) and returns: dict_fcimatr